update: rerun MC integral

src/quick_select.c

@@ -46,7 +46,7 @@ double percentile(const double *arr, size_t n, double perc) {
   }
 
   double *arr_copy = malloc(sizeof(double)*n);
-  memcpy(arr_copy, arr, n);
+  memcpy(arr_copy, arr, n*sizeof(double));
   // nearest-rank: ceil(0.9*n) - 1
   size_t k = (size_t)(perc*n) - 1;
   double ans = select_k(arr_copy, n, k);

src/render.c

@@ -89,6 +89,11 @@ 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);
+    if (err > pixel_render_err) {
+      int ncalls = 100*(err/pixel_render_err)*(err/pixel_render_err);
+      if (ncalls > pixel_render_max) ncalls = pixel_render_max;
+      gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
+    }
     rgb[c] = color;
   }
   //printf(")\n");

src/write_png.c

@@ -25,6 +25,7 @@ int buffer_normalize_srgb(double *buffer, int W, int H) {
   max = percentile(buffer, W*H*3, cutperc);
   printf("max in the buffer: %g\n", max);
   if (max == 0) return 1;
+  #pragma omp parallel for
   for (int j = 0; j < H; j++) {
     for (int i = 0; i < W; i++) {
       for (int c = 0; c < 3; c++) {
@@ -39,7 +40,7 @@ int buffer_normalize_srgb(double *buffer, 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;
-  printf("bufsize: %d\n", bufsize);
+  //printf("bufsize: %d\n", bufsize);
   uint8_t *img = (uint8_t *)malloc(bufsize*sizeof(uint8_t));
   if (!img) {return 1;}