new: add write_png

src/init.h

@@ -1,2 +1,9 @@
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 double chi(double cotpsi);
 int init(double rmax, double rela_err_limit, int *size, double **x, double **y);
+
+typedef struct {
+  gsl_spline *spline;
+  gsl_interp_accel *acc;
+} Spline_data ;

src/render.c

@@ -1,8 +1,16 @@
 #include <gsl/gsl_interp.h>
+#include <gsl/gsl_monte.h>
+#include <gsl/gsl_rng.h>
 #include <gsl/gsl_spline.h>
+#include <gsl/gsl_monte_miser.h>
 #include "common.h"
 #include "init.h"
 
+typedef struct {
+  const System *system;
+  const Spline_data *spline_data;
+} Integrand_params;
+
 inline double xy_to_b(double x, double y) {
   double tanpsi2 = x*x+y*y;
   double cotpsi2 = 1/tanpsi2;
@@ -10,19 +18,79 @@ inline double xy_to_b(double x, double y) {
   return b;
 }
 
-int MC_pixel_render(const System *system, int i, int j, const gsl_spline *spline, double *rgb, int pixel_render_max, double pixel_render_err){
+double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  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);
+}
+
+double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  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);
+}
+
+double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
+  double cotpsi = 1.0/hypot(xy[0], xy[1]);
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc);
+  double phi = atan2(xy[1], xy[0]);
+  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);
+}
+
+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){
   double dw = (system->w)/(system->W); // width and hight of one pixel
   double x_lu = i*dw - (system->w)/2;
   double y_lu = -j*dw + (system->H)*dw/2;
+  double xl[2] = {x_lu, y_lu - dw};
+  double xu[2] = {x_lu + dw, y_lu};
 
   //TODO: fill the MC intergral
-  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) χ'(ψ) dx dy
+  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ) |χ'(ψ)| dx dy
+  // χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ)
+  gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
+
+  double color;
+  double err;
+  Integrand_params integrand_params = {system, &spline_data};
+  gsl_monte_function F;
+  F.params = &integrand_params;
+  F.dim = 2;
+  for (int c = 0; c < 3; c++) {
+    switch (c) {
+      case 0:
+        F.f = integrand_r;
+        break;
+      case 1:
+        F.f = integrand_g;
+        break;
+      case 2:
+        F.f = integrand_b;
+        break;
+    }
+    gsl_monte_miser_init(miser_state);
+    gsl_monte_miser_integrate(&F, xl, xu, 2, 100, r, miser_state, &color, &err);
+    rgb[c] = color;
+  }
 
   return 0;
 }
 
-int pixel_render(const System *system, int i, int j, const gsl_spline *spline, double *rgb, int pixel_render_max, double pixel_render_err){
+int 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){
-  return MC_pixel_render(system, i, j, spline, rgb, pixel_render_max, pixel_render_err);
+  return MC_pixel_render(system, i, j, spline_data, rgb, pixel_render_max, pixel_render_err, r);
 }
 
 int render(System *system, double **buffer, int pixel_render_max, double pixel_render_err, double chi_rela_err) {
@@ -31,6 +99,7 @@ int render(System *system, double **buffer, int pixel_render_max, double pixel_r
   double cotpsi2 = 1/tanpsi2;
   //double bmax = R0/sqrt(f(R0)*(1+cotpsi2));
   double rmax = sqrt(tanpsi2);
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  
   double *x=NULL;
   double *y = NULL;
@@ -40,11 +109,12 @@ int render(System *system, double **buffer, int pixel_render_max, double pixel_r
 
   gsl_interp_accel *acc = gsl_interp_accel_alloc();
   gsl_spline *spline = gsl_spline_alloc(gsl_interp_steffen, size);
+  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, buffer[j*(system->W) + i], pixel_render_max, pixel_render_err);
+      pixel_render(system, i, j, spline_data, buffer[j*(system->W) + i], pixel_render_max, pixel_render_err, r);
     }
   }
 

src/write_png.c

@@ -0,0 +1,82 @@
+#include "write_png.h"
+
+static inline double linear_to_srgb(double x){
+  if (x <= 0.0) return 0.0;
+  if (x >= 1.0) return 1.0;
+  if (x <= 0.0031308) return 12.92 * x;
+  return 1.055 * pow(x, 1.0/2.4) - 0.055;
+}
+
+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];
+      }
+    }
+  }
+  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]);
+      }
+    }
+  }
+  return 0;
+}
+
+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);
+  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]);
+        if (v < 0) v = 0;
+        if (v > 255) v = 255;
+        img[j*W*3+i*3+c] = (uint8_t) v;
+      }
+    }
+  }
+
+  int flag = 0;
+  spng_ctx *ctx = spng_ctx_new(SPNG_CTX_ENCODER);
+  if(!ctx){ fprintf(stderr, "spng_ctx_new failed\n"); free(img); return 1; }
+
+  FILE *f = fopen(filename, "wb");
+  if(!f){ fprintf(stderr, "fopen failed\n"); spng_ctx_free(ctx); free(img); return 1; }
+  spng_set_png_file(ctx, f);
+
+  struct spng_ihdr ihdr = {0};
+  ihdr.width = W;
+  ihdr.height = H;
+  ihdr.bit_depth = 8;
+  ihdr.color_type = SPNG_COLOR_TYPE_TRUECOLOR; // RGB
+  ihdr.interlace_method = SPNG_INTERLACE_NONE;
+  ihdr.compression_method = 0;
+  ihdr.filter_method = 0;
+  if( (flag= spng_set_ihdr(ctx, &ihdr)) ){
+    fprintf(stderr, "spng_set_ihdr: %s\n", spng_strerror(flag));
+    fclose(f); spng_ctx_free(ctx); free(img); return 1;
+  }
+
+  #ifndef SPNG_SRGB_INTENT_PERCEPTUAL
+  #define SPNG_SRGB_INTENT_PERCEPTUAL 0
+  #endif
+  spng_set_srgb(ctx, SPNG_SRGB_INTENT_PERCEPTUAL);
+
+  flag = spng_encode_image(ctx, img, bufsize, SPNG_FMT_PNG, SPNG_ENCODE_FINALIZE);
+  if(flag){
+    fprintf(stderr, "spng_encode_image: %s\n", spng_strerror(flag));
+  }
+  fclose(f);
+  spng_ctx_free(ctx);
+  free(img);
+  return flag? 1 : 0;
+}

src/write_png.h

@@ -0,0 +1,6 @@
+#include <spng.h>
+#include <math.h>
+#include <stdint.h>
+
+
+int write_png(char *filename, double **buffer, int W, int H);