From 4be6add81591966cc72c414fc423515b540d7900 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 12 Dec 2025 13:47:23 -0500 Subject: [PATCH] update: multiple improvement --- src/common.h | 6 ++++-- src/init.c | 17 ++++++++--------- src/render.c | 30 +++++++++++++++++++++--------- src/write_png.c | 14 +------------- 4 files changed, 34 insertions(+), 33 deletions(-) diff --git a/src/common.h b/src/common.h index f6ed152..cf9217f 100644 --- a/src/common.h +++ b/src/common.h @@ -5,20 +5,22 @@ #include #include -#define POW2(x) (x*x) +#define POW2(x) ((x)*(x)) #define FLAT 0 #define PI 3.1415926535897932384626433832795028841971693993751058 #define Rs 1 #define M (FLAT ? 0 : 0.5*Rs) #define R0 (15*Rs) -#define f(r) ((1.0-(2*M)/((double)r))) +#define f(r) ((1.0-(2*M)/((double)(r)))) #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; diff --git a/src/init.c b/src/init.c index c705104..92408b3 100644 --- a/src/init.c +++ b/src/init.c @@ -63,8 +63,8 @@ static inline double find_root_5ord(double u1, double u2, double v1, double v2, pow(u2,2))*v1*v2); } -// \chi(cotpsi) is the totol deflection angle with parameter cotpsi=-dr/(r √f dphi) -double chi(double cotpsi){ +// \chi(tanpsi) is the totol deflection angle with parameter tanpsi=-(r √f dphi)/dr +double chi(double tanpsi){ gsl_odeiv2_system sys = {func, jac, 2, NULL}; gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-3, 1e-12, 1e-10); @@ -72,7 +72,7 @@ double chi(double cotpsi){ double y[2] = { u0, //sqrt(1.0/(b*b) - u0*u0 + 2.0*M*u0*u0*u0) - u0*sqrt(f(R0))*cotpsi + u0*sqrt(f(R0))/tanpsi }; double phi = 0; double phi_max = 8*PI; @@ -142,14 +142,13 @@ int init(double rmax, double rela_err_limit, int *size, double **x, double **y){ sample.size=1; sample.x = malloc(sizeof(double) * 1); sample.y = malloc(sizeof(double) * 1); - double cotpsi_min = 1/rmax; - sample.x[0] = cotpsi_min; - sample.y[0] = chi(cotpsi_min); - double chi_cotpsimax = chi(cotpsi_max); + sample.x[0] = tanpsi_min; + sample.y[0] = chi(tanpsi_min); + double chi_tanpsimax = chi(rmax); - refine_interval(cotpsi_min, cotpsi_max, sample.y[0], chi_cotpsimax, rela_err_limit, &sample, 1); + refine_interval(tanpsi_min, rmax, sample.y[0], chi_tanpsimax, rela_err_limit, &sample, 1); - sample_push(&sample, cotpsi_max, chi_cotpsimax); + sample_push(&sample, rmax, chi_tanpsimax); //sort diff --git a/src/render.c b/src/render.c index fdddb7d..20f3646 100644 --- a/src/render.c +++ b/src/render.c @@ -22,15 +22,17 @@ inline double xy_to_b(double x, double y) { } 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; + 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, cotpsi, integrand_params->spline_data->acc); + 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); - 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)); + //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){ @@ -39,10 +41,13 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline 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 - // χ'(ψ) = χ'(cot(ψ)) * cot'(ψ) = - χ'(cot(ψ))/\sin^2(ψ) + // χ'(ψ) = χ'(tan(ψ)) * tan'(ψ) = - χ'(tan(ψ))/cos^2(ψ) = - χ'(tan(ψ))(1+tan^2(ψ)) gsl_monte_miser_state *miser_state = gsl_monte_miser_alloc(2); double color; @@ -55,17 +60,23 @@ int MC_pixel_render(const System *system, int i, int j, const Spline_data spline //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 > pixel_render_err) { - int ncalls = 100*(err/pixel_render_err)*(err/pixel_render_err); - if (ncalls > pixel_render_max) ncalls = pixel_render_max; + 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(")\n"); + //printf("finished for (i,j)=(%d, %d)\n", i, j); return 0; } @@ -108,5 +119,6 @@ int render(System *system, double *buffer, int pixel_render_max, double pixel_re gsl_interp_accel_free(acc); gsl_rng_free(r); } + printf("rendering finished.\n"); return 0; } diff --git a/src/write_png.c b/src/write_png.c index e4fe247..aaa0e21 100644 --- a/src/write_png.c +++ b/src/write_png.c @@ -14,23 +14,11 @@ static inline double linear_to_srgb(double x){ } 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[color_index(i, j, c)]) max = buffer[color_index(i, j, c)]; - // } - // } - //} - 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++) { - buffer[color_index(i, j, c)] /= max; - buffer[color_index(i, j, c)] = linear_to_srgb(buffer[color_index(i, j, c)]); + buffer[color_index(i, j, c)] = linear_to_srgb(SCALE*buffer[color_index(i, j, c)]); } } }