update: allow -DFLAT=1 to set a flat spacetime without black hole during compilation

This is because when I pass *x then it's just a copy.

src/common.h

@@ -1,12 +1,34 @@
+#ifndef __COMMON_H
+#define __COMMON_H
 #include <stdio.h>
+#include <stdlib.h>
 #include <gsl/gsl_errno.h>
 #include <math.h>
 
+#define POW2(x) ((x)*(x))
+
+#ifndef FLAT
+#define FLAT 0
+#endif
 #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 f(r) ((1.0-(2*M)/((double)(r))))
-#define bmin (1.5*sqrt(3.0))
+#define bmin (FLAT ? 1e-10 : 1.5*sqrt(3.0)+5e-10)
+#define cotpsi_max (sqrt(R0*R0/(bmin*bmin*f(R0)) - 1))
+#define tanpsi_min (1.0/cotpsi_max)
 #define THETAERROR 100000
+#define color_index(i,j,c) (j*W*3+i*3+c)
 
+#define cutperc (FLAT ? 1.0 : 0.99)
+#define SCALE 1
+
+typedef struct {
+  int W;
+  int H;
+  double w; //width of the sensor. focal length=1
+  int (*angle_to_pixel)(double *, double *);
+} System;
+
+#endif // !__COMMON_H

src/init.c

@@ -1,7 +1,11 @@
 #include <gsl/gsl_matrix_double.h>
 #include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_matrix.h>
+#include <stdio.h>
 #include "common.h"
+#include "init.h"
+
+#define MAX_DEPTH 15
 
 typedef struct {
   int capacity;
@@ -24,11 +28,11 @@ static int jac(double phi, const double y[], double *dfdy, double dfdphi[], void
   gsl_matrix *m = &dfdy_mat.matrix;
   gsl_matrix_set(m, 0, 0, 0.0);
   gsl_matrix_set(m, 0, 1, 1.0);
-  gsl_matrix_set(m, 1, 0, -3*M*y[0]-1.0);
+  gsl_matrix_set(m, 1, 0, 6*M*y[0]-1.0);
   gsl_matrix_set(m, 1, 1, 0.0);
 
   dfdphi[0] = 0.0;
-  dfdphi[0] = 0.0;
+  dfdphi[1] = 0.0;
 
   return GSL_SUCCESS;
 }
@@ -59,16 +63,16 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2,
        pow(u2,2))*v1*v2);
 }
 
-// \chi(b) is the totol deflection angle with parameter b
+// \chi(tanpsi) is the totol deflection angle with parameter tanpsi=-(r √f dphi)/dr
-double chi(double b){
+double chi(double tanpsi){
-  double x;
+  gsl_odeiv2_system sys = {func, jac, 2, NULL};
-  gsl_odeiv2_system sys = {func, jac, 2, &x};
   gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10);
   
   double u0 = 1.0/R0;
   double y[2] = {
     u0,
-    sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0)
+    //sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0)
+    u0*sqrt(f(R0))/tanpsi
   };
   double phi = 0;
   double phi_max = 8*PI;
@@ -106,39 +110,59 @@ void sample_push(SampleData *sample, double x, double y){
     sample->x = realloc(sample->x, sizeof(double) * newcap);
     sample->y = realloc(sample->y, sizeof(double) * newcap);
     sample->capacity = newcap;
+    //printf("enlarge sample capacity: %d\n", newcap);
   }
   sample->x[sample->size] = x;
   sample->y[sample->size] = y;
   sample->size += 1;
 }
 
-static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample){
+static void refine_interval(double xl, double xr, double yl, double yr, double rela_err_limit, SampleData *sample, int depth){
+  //printf("calculating interval (%.16g, %.16g)\n", xl, xr);
   double xm = (xl + xr)/2;
   double ym = chi(xm);
+  int need_refine = (fabs(ym-(yl+yr)/2) >= fabs(rela_err_limit * ym)) && depth <= MAX_DEPTH;
+  if (need_refine) {
+    //printf(" ");
+    refine_interval(xl, xm, yl, ym, rela_err_limit, sample, depth + 1);
+  }
   sample_push(sample, xm, ym);
-  if (fabs((ym-(yl+yr)/2)/ym) >= rela_err_limit) {
+  if (need_refine) {
-    refine_interval(xl, xm, yl, ym, rela_err_limit, sample);
+    //printf(" ");
-    refine_interval(xm, xr, ym, yr, rela_err_limit, sample);
+    refine_interval(xm, xr, ym, yr, rela_err_limit, sample, depth + 1);
   }
   return ;
 }
 
-static int init(double bmax, double rela_err_limit){
+int init(double rmax, double rela_err_limit, int *size, double **x, double **y){
+  if (*x) free(*x);
+  if (*y) free(*y);
   SampleData sample;
-  sample.capacity = 2;
+  sample.capacity = 1;
-  sample.size=2;
+  sample.size=1;
-  sample.x = malloc(sizeof(double) * 2);
+  sample.x = malloc(sizeof(double) * 1);
-  sample.y = malloc(sizeof(double) * 2);
+  sample.y = malloc(sizeof(double) * 1);
-  sample.x[0] = bmin;
+  sample.x[0] = tanpsi_min;
-  sample.x[1] = bmax;
+  sample.y[0] = chi(tanpsi_min);
-  sample.y[0] = chi(bmin);
+  double chi_tanpsimax = chi(rmax);
-  sample.y[1] = chi(bmax);
 
-  refine_interval(bmin, bmax, sample.y[0], sample.y[1], rela_err_limit, &sample);
+  refine_interval(tanpsi_min, rmax, sample.y[0], chi_tanpsimax, rela_err_limit, &sample, 1);
+
+  sample_push(&sample, rmax, chi_tanpsimax);
 
   //sort
 
   //output
+  //printf("Total sample points capacity: %d\n", sample.capacity);
+  //printf("Total sample points number: %d\n", sample.size);
+  *size = sample.size;
+  *x = malloc(sizeof(double) * sample.size);
+  *y = malloc(sizeof(double) * sample.size);
+  for (int i = 0; i < sample.size; i++) {
+    //printf("%.16g %.16g\n", sample.x[i], sample.y[i]);
+    (*x)[i] = sample.x[i];
+    (*y)[i] = sample.y[i];
+  }
 
   free(sample.x);
   free(sample.y);

src/init.h

@@ -1 +1,9 @@
-double chi(double b);
+#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/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*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);
+  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

@@ -0,0 +1,124 @@
+#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 <math.h>
+#include <stdio.h>
+#include "render.h"
+#include "common.h"
+
+typedef struct {
+  int c;
+  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;
+  double b = R0/sqrt(f(R0)*(1+cotpsi2));
+  return b;
+}
+
+double integrand(double *xy, size_t dim, void *integrand_params_void) {
+  double tanpsi = hypot(xy[0], xy[1]);
+  double cotpsi = 1.0/tanpsi;
+  if (tanpsi <= tanpsi_min) return 0;
+  Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
+  double theta = gsl_spline_eval(integrand_params->spline_data->spline, tanpsi, 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);
+  //printf("(x,y)=(%g,%g), tanpsi=%g, theta=%g, %g\n", xy[0], xy[1], tanpsi, theta, POW2(sin(theta))*(2+POW2(cotpsi)+POW2(tanpsi)));
+  return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, tanpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*(2+POW2(tanpsi)+POW2(cotpsi));
+}
+
+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};
+  double err_dw = pixel_render_err*dw*dw;
+
+  //printf("integrating for (i,j)=(%d, %d)\n", i, j);
+
+  //TODO: fill the MC intergral
+  // linear rgb = ∫ rgb(χ, ϕ) sin^2(χ)/sin^2(ψ) |χ'(ψ)| dx dy
+  // χ'(ψ) = χ'(tan(ψ)) * tan'(ψ) = - χ'(tan(ψ))/cos^2(ψ) = - χ'(tan(ψ))(1+tan^2(ψ))
+  gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2);
+
+  double color;
+  double err;
+  Integrand_params integrand_params = {0, system, &spline_data};
+  gsl_monte_function F;
+  F.f = integrand;
+  F.params = &integrand_params;
+  F.dim = 2;
+  //printf("i=%d, j=%d, rgb=( ", i, j);
+  for (int c = 0; c < 3; c++) {
+    integrand_params.c = c;
+    if (i==0 && j==0) {
+      //printf("(x,y) = (%g, %g), intergrand[%d]=%g\n", xl[0], xl[1], c, integrand(xl, 2, &integrand_params));
+    }
+    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 > err_dw) {
+      //int ncalls = 100*POW2(err/err_dw);
+      //printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
+      //if (ncalls < 0 || ncalls > pixel_render_max) ncalls = pixel_render_max;
+      int ncalls = pixel_render_max;
+      //printf("(i,j)=(%d,%d), err=%g, err_dw=%g, ncalls=%d\n", i, j, err, err_dw, ncalls);
+      gsl_monte_miser_integrate(&F, xl, xu, 2, ncalls, r, miser_state, &color, &err);
+    }
+    rgb[c] = color/(dw*dw);
+  }
+  //printf("finished for (i,j)=(%d, %d)\n", i, j);
+
+  return 0;
+}
+
+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_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) {
+  double h = (system->H)*(system->w)/(system->W);
+  //double tanpsi2 = (system->w)*(system->w)+h*h;
+  //double cotpsi2 = 1/tanpsi2;
+  //double bmax = R0/sqrt(f(R0)*(1+cotpsi2));
+  double rmax = hypot(system->w, h);
+ 
+  double *x=NULL;
+  double *y = NULL;
+  int size;
+  int W = system->W;
+  int H = system->H;
+
+  init(rmax, chi_rela_err, &size, &x, &y);
+
+  gsl_spline *spline = gsl_spline_alloc(gsl_interp_steffen, size);
+  gsl_spline_init(spline, x, y, size);
+
+  #pragma omp parallel
+  {
+    gsl_interp_accel *acc = gsl_interp_accel_alloc();
+    Spline_data spline_data = {spline, acc};
+    gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
+ 
+    #pragma omp for
+    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);
+      }
+    }
+
+    gsl_interp_accel_free(acc);
+    gsl_rng_free(r);
+  }
+  printf("rendering finished.\n");
+  return 0;
+}

src/render.h

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

src/test.c

@@ -0,0 +1,22 @@
+#include "common.h"
+#include "init.h"
+#include <stdio.h>
+
+int main(){
+  int size;
+  double *x = NULL;
+  double *y = NULL;
+  init(15.0, 1e-1, &size, &x, &y);
+  printf("Done!\n");
+  printf("Sample size: %d\n", size);
+  printf("x: %p\n", x);
+  printf("y: %p\n", y);
+  for (int i = 0; i < size; i++){
+    printf("%g %g\n", x[i], y[i]);
+  }
+  //chi(bmin);
+
+  free(x);
+  free(y);
+  return 0;
+}

src/test_png.c

@@ -0,0 +1,73 @@
+#include "render.h"
+#include "write_png.h"
+
+int angle_to_pixel_white(double *angle, double *rgb) {
+  rgb[0] = 1.0;
+  rgb[1] = 1.0;
+  rgb[2] = 1.0;
+  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 = 1600;
+  int H = 900;
+  double w = 2;
+  System system = {W, H, w, test};
+
+  double *buffer = malloc(sizeof(double)*W*H*3);
+
+  render(&system, buffer, 10000, 1.0/256, 1e-1);
+  write_png("test_1.png", buffer, W, H);
+  free(buffer);
+  return 0;
+}

src/write_png.c

@@ -0,0 +1,81 @@
+#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;
+  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) {
+  #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++) {
+        buffer[color_index(i, j, c)] = linear_to_srgb(SCALE*buffer[color_index(i, j, 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;
+  //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[color_index(i,j,c)]*255);
+        if (v < 0) v = 0;
+        if (v > 255) v = 255;
+        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);
+      }
+    }
+  }
+
+  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,3 @@
+#include "common.h"
+
+int write_png(char *filename, double *buffer, int W, int H);