update: more tests

src/common.h

@@ -1,16 +1,22 @@
+#ifndef __COMMON_H
+#define __COMMON_H
 #include <stdio.h>
 #include <stdlib.h>
 #include <gsl/gsl_errno.h>
 #include <math.h>
 
+#define FLAT 0
 #define PI 3.1415926535897932384626433832795028841971693993751058
 #define Rs 1
-#define M (0.5*Rs)
+#define M (FLAT ? 0 : 0.5*Rs)
 #define R0 (15*Rs)
-#define f(r) ((1.0-Rs/((double)r)))
-#define bmin (1.5*sqrt(3.0)+5e-10)
+#define f(r) ((1.0-(2*M)/((double)r)))
+#define bmin (FLAT ? 1e-10 : 1.5*sqrt(3.0)+5e-10)
 #define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1))
 #define THETAERROR 100000
+#define color_index(i,j,c) (j*W*3+i*3+c)
+
+#define cutperc (FLAT ? 1.0 : 0.99)
 
 typedef struct {
   int W;
@@ -19,3 +25,4 @@ typedef struct {
   int (*angle_to_pixel)(double *, double *);
 } System;
 
+#endif // !__COMMON_H

src/quick_select.c

@@ -0,0 +1,55 @@
+#include <string.h>
+#include <stdlib.h>
+#include "quick_select.h"
+
+static inline void swap_double(double *a, double *b) {
+  double t = *a;
+  *a = *b;
+  *b = t;
+}
+
+// Lomuto partition
+static size_t partition(double *arr, size_t lo, size_t hi) {
+  double pivot = arr[hi];
+  size_t i = lo;
+  for (size_t j = lo; j < hi; ++j) {
+    if (arr[j] < pivot) {
+      swap_double(&arr[i], &arr[j]);
+      ++i;
+    }
+  }
+  swap_double(&arr[i], &arr[hi]);
+  return i;
+}
+
+// 选择第 k 小的元素 (0-based)，平均 O(n)，最坏 O(n^2)
+static double select_k(double *arr, size_t n, size_t k) {
+    size_t lo = 0;
+    size_t hi = n - 1;
+
+  for (;;) {
+    size_t p = partition(arr, lo, hi);
+    if (p == k) {
+      return arr[p];
+    } else if (k < p) {
+      hi = p - 1;
+    } else {
+      lo = p + 1;
+    }
+  }
+  return arr[k];
+}
+
+double percentile(const double *arr, size_t n, double perc) {
+  if (n == 0) {
+    return 0.0;
+  }
+
+  double *arr_copy = malloc(sizeof(double)*n);
+  memcpy(arr_copy, arr, n);
+  // nearest-rank: ceil(0.9*n) - 1
+  size_t k = (size_t)(perc*n) - 1;
+  double ans = select_k(arr_copy, n, k);
+  free(arr_copy);
+  return ans;
+}

src/quick_select.h

@@ -0,0 +1,2 @@
+#include <stddef.h>
+double percentile(const double *arr, size_t n, double perc);

src/render.c

@@ -3,6 +3,7 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_spline.h>
 #include <gsl/gsl_monte_miser.h>
+#include <stdio.h>
 #include "render.h"
 
 typedef struct {
@@ -26,7 +27,7 @@ double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
   double rgb[3] = {0};
   double theta_phi[2] = {theta, phi};
   integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[0]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  return rgb[0]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
 }
 
 double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
@@ -38,7 +39,7 @@ double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
   double rgb[3] = {0};
   double theta_phi[2] = {theta, phi};
   integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[1]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  return rgb[1]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
 }
 
 double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
@@ -50,7 +51,7 @@ double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
   double rgb[3] = {0};
   double theta_phi[2] = {theta, phi};
   integrand_params->system->angle_to_pixel(theta_phi, rgb);
-  return rgb[2]*gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  return rgb[2]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc));
 }
 
 int MC_pixel_render(const System *system, int i, int j, const Spline_data spline_data, double *rgb, int pixel_render_max, double pixel_render_err, gsl_rng *r){
@@ -71,6 +72,7 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
   gsl_monte_function F;
   F.params = &integrand_params;
   F.dim = 2;
+  //printf("i=%d, j=%d, rgb=( ", i, j);
   for (int c = 0; c < 3; c++) {
     switch (c) {
       case 0:
@@ -85,8 +87,10 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
     }
     gsl_monte_miser_init(miser_state);
     gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err);
+    //printf("%g +- %g ", color, err);
     rgb[c] = color;
   }
+  //printf(")\n");
 
   return 0;
 }
@@ -95,7 +99,7 @@ int pixel_render(const System *system, int i, int j, const Spline_data spline_da
   return MC_pixel_render(system, i, j, spline_data, rgb, pixel_render_max, pixel_render_err, r);
 }
 
-int render(System *system, double (*buffer)[3], int pixel_render_max, double pixel_render_err, double chi_rela_err) {
+int render(System *system, double *buffer, int pixel_render_max, double pixel_render_err, double chi_rela_err) {
   double h = (system->H)*(system->w)/(system->W);
   double tanpsi2 = (system->w)*(system->w)+h*h;
   double cotpsi2 = 1/tanpsi2;
@@ -106,6 +110,8 @@ int render(System *system, double (*buffer)[3], int pixel_render_max, double pix
   double *x=NULL;
   double *y = NULL;
   int size;
+  int W = system->W;
+  int H = system->H;
 
   init(rmax, chi_rela_err, &size, &x, &y);
 
@@ -114,9 +120,9 @@ int render(System *system, double (*buffer)[3], int pixel_render_max, double pix
   Spline_data spline_data = {spline, acc};
   gsl_spline_init(spline, x, y, size);
 
-  for(int j = 0; j < system->H; j++) {
-    for (int i = 0; i < system->W; i++) {
-      pixel_render(system, i, j, spline_data, buffer[j*(system->W) + i], pixel_render_max, pixel_render_err, r);
+  for(int j = 0; j < H; j++) {
+    for (int i = 0; i < W; i++) {
+      pixel_render(system, i, j, spline_data, buffer+color_index(i,j,0), pixel_render_max, pixel_render_err, r);
     }
   }
 

src/render.h

@@ -1,4 +1,4 @@
 #include "common.h"
 #include "init.h"
 
-int render(System *system, double (*buffer)[3], int pixel_render_max, double pixel_render_err, double chi_rela_err);
+int render(System *system, double *buffer, int pixel_render_max, double pixel_render_err, double chi_rela_err);

src/test_png.c

@@ -1,4 +1,3 @@
-#include <stdio.h>
 #include "render.h"
 #include "write_png.h"
 
@@ -9,15 +8,66 @@ int angle_to_pixel_white(double *angle, double *rgb) {
   return 0;
 }
 
+int test(double *thetaphi, double *rgb){
+  int i = floor(8.0*thetaphi[0]/PI);
+  int j = floor(8.0*(thetaphi[1]+PI)/PI);
+  int k = (i+j)&1;
+  int l = floor(2*thetaphi[1]/PI);
+  rgb[0] = k;
+  rgb[1] = k;
+  rgb[2] = k;
+  switch (l%4) {
+    case 1:
+      rgb[1] = rgb[2] = 0;
+      break;
+    case 2:
+      rgb[0] = rgb[2] = 0;
+      break;
+    case 3:
+      rgb[0] = rgb[1] = 0;
+      break;
+    default:
+      break;
+  }
+  return 0;
+}
+
+int test2(double *thetaphi, double *rgb){
+  double tantheta = tan(PI-thetaphi[0]);
+  int i = (int)floor(4*tantheta*cos(thetaphi[1]));
+  int j = (int)floor(4*tantheta*sin(thetaphi[1]));
+  int k = (i+j)&1;
+  rgb[0] = 1-k;
+  rgb[1] = 1-k;
+  rgb[2] = 1-k;
+  return 0;
+}
+
+int test3(double *thetaphi, double *rgb){
+  double x = sin(thetaphi[0])*sin(thetaphi[1]);
+  double y = sin(thetaphi[0])*cos(thetaphi[1]);
+  double z = cos(thetaphi[0]);
+  double theta1 = atan2(hypot(z, x), y);
+  double phi1 = atan2(-z, x);
+  int i = (int)floor(16*theta1/PI);
+  int j = (int)floor(16*phi1/PI);
+  int k = (i+j)&1;
+  rgb[0] = 1-k;
+  rgb[1] = 1-k;
+  rgb[2] = 1-k;
+  return 0;
+}
+
 int main(){
-  int W = 1;
-  int H = 1;
-  double w = 1;
-  System system = {W, H, w, angle_to_pixel_white};
+  int W = 1600;
+  int H = 900;
+  double w = 2;
+  System system = {W, H, w, test};
 
-  double buffer[W*H][3];
+  double *buffer = malloc(sizeof(double)*W*H*3);
 
   render(&system, buffer, 10000, 1.0/256, 1e-1);
-  write_png("test_white.png", buffer, W, H);
+  write_png("test_1.png", buffer, W, H);
+  free(buffer);
   return 0;
 }

src/write_png.c

@@ -1,7 +1,10 @@
 #include <spng.h>
 #include <math.h>
 #include <stdint.h>
+#include <stdio.h>
 #include "write_png.h"
+#include "common.h"
+#include "quick_select.h"
 
 static inline double linear_to_srgb(double x){
   if (x <= 0.0) return 0.0;
@@ -10,40 +13,44 @@ static inline double linear_to_srgb(double x){
   return 1.055 * pow(x, 1.0/2.4) - 0.055;
 }
 
-int buffer_normalize_srgb(double (*buffer)[3], int W, int H) {
+int buffer_normalize_srgb(double *buffer, int W, int H) {
   double max = 0;
-  for (int j = 0; j < H; j++) {
-    for (int i = 0; i < W; i++) {
-      for (int c = 0; c < 3; c++) {
-        if (max < buffer[j*W + i][c]) max = buffer[j*W + i][c];
-      }
-    }
-  }
+  //for (int j = 0; j < H; j++) {
+  //  for (int i = 0; i < W; i++) {
+  //    for (int c = 0; c < 3; c++) {
+  //      if (max < buffer[color_index(i, j, c)]) max = buffer[color_index(i, j, c)];
+  //    }
+  //  }
+  //}
+  max = percentile(buffer, W*H*3, cutperc);
+  printf("max in the buffer: %g\n", max);
   if (max == 0) return 1;
   for (int j = 0; j < H; j++) {
     for (int i = 0; i < W; i++) {
       for (int c = 0; c < 3; c++) {
-        buffer[j*W + i][c] /= max;
-        buffer[j*W+i][c] = linear_to_srgb(buffer[j*W+i][c]);
+        buffer[color_index(i, j, c)] /= max;
+        buffer[color_index(i, j, c)] = linear_to_srgb(buffer[color_index(i, j, c)]);
       }
     }
   }
   return 0;
 }
 
-int write_png(char *filename, double (*buffer)[3], int W, int H) {
+int write_png(char *filename, double *buffer, int W, int H) {
   buffer_normalize_srgb(buffer, W, H);
   size_t bufsize = W*H*3;
-  uint8_t *img = (uint8_t *)malloc(bufsize);
+  printf("bufsize: %d\n", bufsize);
+  uint8_t *img = (uint8_t *)malloc(bufsize*sizeof(uint8_t));
   if (!img) {return 1;}
 
   for (int j = 0; j < H; j++) {
     for (int i = 0; i < W; i++) {
       for (int c = 0; c < 3; c++) {
-        int v = (int)lrint(buffer[j*W*i][c]);
+        int v = (int)lrint(buffer[color_index(i,j,c)]*255);
         if (v < 0) v = 0;
         if (v > 255) v = 255;
-        img[j*W*3+i*3+c] = (uint8_t) v;
+        img[color_index(i, j, c)] = (uint8_t) v;
+        //printf("i: %d, j: %d, c:%d, buffer=%g, v=%d\n", i,j,c,buffer[color_index(i, j, c)],v);
       }
     }
   }

src/write_png.h

@@ -1,3 +1,3 @@
+#include "common.h"
 
-
-int write_png(char *filename, double (*buffer)[3], int W, int H);
+int write_png(char *filename, double *buffer, int W, int H);