wyj/blackhole-ray-trace
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update: add MC integral render

Browse Source
This commit is contained in:
Yingjie Wang
2025-12-07 04:50:42 -05:00
parent 0c491d64ce
commit a17bff625c
2 changed files with 82 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/init.h
View File

@@ -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 ;

80
src/render.c
View File

@@ -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){
return MC_pixel_render(system, i, j, spline, rgb, pixel_render_max, 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_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);
}
}