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

update: use a single integrand function

Browse Source
This commit is contained in:
Yingjie Wang
2025-12-07 22:47:53 -05:00
parent 480c49dcf0
commit de6d3bc04a
1 changed files with 7 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
src/render.c
View File

@@ -9,6 +9,7 @@
#include "common.h"
typedef struct {
int c;
const System *system;
const Spline_data *spline_data;
} Integrand_params;
@@ -20,7 +21,7 @@ inline double xy_to_b(double x, double y) {
return b;
}
double integrand_r(double *xy, size_t dim, void *integrand_params_void) {
double integrand(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
Integrand_params *integrand_params = (Integrand_params *)integrand_params_void;
@@ -29,34 +30,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]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[0];
}
double integrand_g(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
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]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[1];
}
double integrand_b(double *xy, size_t dim, void *integrand_params_void) {
double cotpsi = 1.0/hypot(xy[0], xy[1]);
if (cotpsi >= cotpsi_max) return 0;
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]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+POW2(cotpsi));
return rgb[2];
return rgb[integrand_params->c]*fabs(gsl_spline_eval_deriv(integrand_params->spline_data->spline, cotpsi, integrand_params->spline_data->acc))*POW2(sin(theta))*POW2(1+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){
@@ -73,23 +47,14 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline
double color;
double err;
Integrand_params integrand_params = {system, &spline_data};
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++) {
switch (c) {
case 0:
F.f = integrand_r;
break;
case 1:
F.f = integrand_g;
break;
case 2:
F.f = integrand_b;
break;
}
integrand_params.c = c;
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);
@@ -132,7 +97,7 @@ int render(System *system, double *buffer, int pixel_render_max, double pixel_re
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++) {