ssf_phase_rdgfa.c (11743B)
1 /* Copyright (C) 2016-2018, 2021-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "ssf.h" 17 #include "ssf_phase_c.h" 18 19 #include <star/ssp.h> 20 21 #include <rsys/algorithm.h> 22 #include <rsys/dynamic_array_double.h> 23 24 #if defined(SSF_USE_SIMD_128) || defined(SSF_USE_SIMD_256) 25 /* Helper macro */ 26 #define USE_SIMD 27 #endif 28 29 #ifdef USE_SIMD 30 #include <rsimd/math.h> 31 #include <rsimd/rsimd.h> 32 33 /* Generate the struct darray_simdf dynamic array */ 34 #define DARRAY_NAME simdf 35 #define DARRAY_DATA float 36 #ifdef SSF_USE_SIMD_256 37 #define DARRAY_ALIGNMENT 64 38 #else 39 #define DARRAY_ALIGNMENT 16 40 #endif 41 #include <rsys/dynamic_array.h> 42 #endif 43 44 #define EXP1 2.7182818284590452354 45 46 struct rdgfa { 47 double wavelength; /* In nm */ 48 double fractal_dimension; /* No unit */ 49 double gyration_radius; /* In nm */ 50 51 double rcp_normalize_factor; /* Reciprocal normalization factor of the CDF */ 52 53 /* Discretized cumulative (#entries = nintervals) */ 54 struct darray_double cdf; 55 56 #ifdef USE_SIMD 57 struct darray_simdf f_list; 58 struct darray_simdf d_omega_list; 59 #endif 60 61 /* Precomputed values */ 62 double Rg2; /* gyration_radius^2 */ 63 double cst_4pi_div_lambda; /* (4*PI)/wavelength */ 64 double cst_3Df_div_2E; /* 3*Df/(2*exp(1)) */ 65 double g; /* g function */ 66 67 /* #intervals used to discretize the [0,pi[ angular domain */ 68 size_t nintervals; 69 70 /* Length of an angular interval */ 71 double dtheta; /* In rad */ 72 }; 73 74 /******************************************************************************* 75 * Helper functions 76 ******************************************************************************/ 77 static INLINE int 78 check_phase_rdgfa_setup_args(const struct ssf_phase_rdgfa_setup_args* args) 79 { 80 return args 81 && args->wavelength > 0 82 && args->fractal_dimension > 0 83 && args->gyration_radius > 0 84 && args->nintervals != 0; 85 } 86 87 static INLINE int 88 cmp_dbl(const void* a, const void* b) 89 { 90 const double key = *((double*)a); 91 const double val = *((double*)b); 92 if(key < val) { 93 return -1; 94 } else if(key > val) { 95 return 1; 96 } else { 97 return 0; 98 } 99 } 100 101 static INLINE double 102 eval_f(struct rdgfa* rdgfa, const double theta) 103 { 104 double Df, Rg2, q, q2Rg2, f; 105 ASSERT(rdgfa); 106 107 /* Input arguments */ 108 Df = rdgfa->fractal_dimension; 109 110 /* Fech precomputed constants */ 111 Rg2 = rdgfa->Rg2; 112 113 /* Precompute values */ 114 q = rdgfa->cst_4pi_div_lambda * sin(theta*0.5); 115 q2Rg2 = q*q*Rg2; 116 117 /* Evaluate f(theta) */ 118 f = q2Rg2 < 1.5*Df 119 ? exp(-1.0/3.0 * q2Rg2) 120 : pow(rdgfa->cst_3Df_div_2E * 1/q2Rg2, Df*0.5); 121 return f; 122 } 123 124 /* Not normalized */ 125 static INLINE double 126 eval2 127 (struct rdgfa* rdgfa, 128 const double theta, 129 const double cos_theta) 130 { 131 double f, cos2_theta, phase; 132 ASSERT(rdgfa && eq_eps(cos_theta, cos(theta), 1.e-6)); 133 134 /* Precompute values */ 135 cos2_theta = cos_theta * cos_theta; 136 137 /* Evaluate phase(theta) */ 138 f = eval_f(rdgfa, theta); 139 phase = 3.0/(16*PI) * f / rdgfa->g * (1 + cos2_theta); 140 return phase; 141 } 142 143 /* Not normalized */ 144 static FINLINE double 145 eval(struct rdgfa* rdgfa, const double theta) 146 { 147 return eval2(rdgfa, theta, cos(theta)); 148 } 149 150 static INLINE res_T 151 compute_cumulative(struct rdgfa* rdgfa) 152 { 153 double* cdf = NULL; 154 double f1; 155 double theta1; 156 size_t i; 157 res_T res = RES_OK; 158 ASSERT(rdgfa); 159 160 /* Allocate the cumulative array */ 161 res = darray_double_resize(&rdgfa->cdf, rdgfa->nintervals); 162 if(res != RES_OK) goto error; 163 cdf = darray_double_data_get(&rdgfa->cdf); 164 165 /* Compute the angular step for the angular domain */ 166 rdgfa->dtheta = PI / (double)rdgfa->nintervals; 167 168 theta1 = 0; 169 f1 = eval(rdgfa, 0); 170 FOR_EACH(i, 0, rdgfa->nintervals) { 171 /* Compute the upper bound of the current angular range */ 172 const double theta2 = theta1 + rdgfa->dtheta; 173 174 /* Compute the (unormalized) cumulative for current interval */ 175 const double delta_omega = 2*PI*sin((theta1+theta2)*0.5)*rdgfa->dtheta; 176 const double f2 = eval(rdgfa, theta2); 177 const double tmp = (f1 + f2) * 0.5 * delta_omega; 178 cdf[i] = (i == 0 ? tmp : tmp + cdf[i-1]); 179 180 /* Go to the next interval */ 181 f1 = f2; 182 theta1 = theta2; 183 } 184 185 /* Save the normamlization factor */ 186 rdgfa->rcp_normalize_factor = 1.0 / cdf[rdgfa->nintervals-1]; 187 188 /* Finally normalize the CDF */ 189 FOR_EACH(i, 0, rdgfa->nintervals) { 190 cdf[i] *= rdgfa->rcp_normalize_factor; 191 } 192 193 exit: 194 return res; 195 error: 196 darray_double_clear(&rdgfa->cdf); 197 goto exit; 198 } 199 200 /* Generate the simd functions if required */ 201 #ifdef SSF_USE_SIMD_128 202 #define SIMD_WIDTH__ 4 203 #include "ssf_phase_rdgfa_simdX.h" 204 #endif 205 #ifdef SSF_USE_SIMD_256 206 #define SIMD_WIDTH__ 8 207 #include "ssf_phase_rdgfa_simdX.h" 208 #endif 209 210 /******************************************************************************* 211 * Private functions 212 ******************************************************************************/ 213 static res_T 214 rdgfa_init(struct mem_allocator* allocator, void* phase) 215 { 216 struct rdgfa* rdgfa = phase; 217 ASSERT(phase); 218 memset(rdgfa, 0, sizeof(*rdgfa)); 219 darray_double_init(allocator, &rdgfa->cdf); 220 #ifdef USE_SIMD 221 darray_simdf_init(allocator, &rdgfa->f_list); 222 darray_simdf_init(allocator, &rdgfa->d_omega_list); 223 #endif /* USE_SIMD */ 224 return RES_OK; 225 } 226 227 static void 228 rdgfa_release(void* phase) 229 { 230 struct rdgfa* rdgfa = phase; 231 ASSERT(phase); 232 darray_double_release(&rdgfa->cdf); 233 #ifdef USE_SIMD 234 darray_simdf_release(&rdgfa->f_list); 235 darray_simdf_release(&rdgfa->d_omega_list); 236 #endif /* USE_SIMD */ 237 } 238 239 static double 240 rdgfa_eval(void* data, const double wo[3], const double wi[3]) 241 { 242 const struct rdgfa* rdgfa = data; 243 double cos_theta, theta; 244 double v[3]; 245 ASSERT(d3_is_normalized(wo) && d3_is_normalized(wi)); 246 247 d3_minus(v, wo); 248 cos_theta = d3_dot(v, wi); 249 theta = acos(cos_theta); 250 return eval2(data, theta, cos_theta) * rdgfa->rcp_normalize_factor; 251 } 252 253 static void 254 rdgfa_sample 255 (void* data, 256 struct ssp_rng* rng, 257 const double wo[3], 258 double wi[3], 259 double* pdf) 260 { 261 struct rdgfa* rdgfa = data; 262 const double* cdf = NULL; 263 double* find = NULL; 264 double frame[9]; 265 double w[3]; 266 double thetas[2]; 267 double theta; 268 double sin_theta; 269 double cos_theta; 270 double phi; 271 double r; 272 double u; 273 size_t n; 274 size_t i; 275 ASSERT(data && rng && wo && wi); 276 277 /* Fetch the CDF array and its number of entries */ 278 cdf = darray_double_cdata_get(&rdgfa->cdf); 279 n = darray_double_size_get(&rdgfa->cdf); 280 281 /* Sample the CDF */ 282 r = ssp_rng_canonical(rng); 283 find = search_lower_bound(&r, cdf, n, sizeof(double), cmp_dbl); 284 ASSERT(find && (*find) >= r); 285 i = (size_t)(find - cdf); 286 287 /* Compute the angular range into which the sample lies */ 288 thetas[0] = rdgfa->dtheta * (double)(i + 0); 289 thetas[1] = rdgfa->dtheta * (double)(i + 1); 290 291 /* Compute the sampled theta angle by linearly interpolate it from the 292 * boundaries of its interval */ 293 if(i == 0) { 294 u = r / cdf[0]; 295 } else { 296 u = (r - cdf[i-1]) / (cdf[i] - cdf[i-1]); 297 } 298 ASSERT(0 <= u && u < 1); 299 theta = u * (thetas[1] - thetas[0]) + thetas[0]; 300 301 /* Uniformly sample a phi angle in [0, 2PI[ */ 302 phi = ssp_rng_uniform_double(rng, 0, 2*PI); 303 sin_theta = sin(theta); 304 cos_theta = cos(theta); 305 306 /* Compute the cartesian coordinates of the sampled direction in the _local_ 307 * phase function space */ 308 wi[0] = cos(phi) * sin_theta; 309 wi[1] = sin(phi) * sin_theta; 310 wi[2] = cos_theta; 311 312 /* Compute the transformation matrix from local phase function to world 313 * space. Note that by convention, in Star-SF the directions point outward 314 * the scattering position. Revert 'wo' to match the convention used by the 315 * previous computations */ 316 d33_basis(frame, d3_minus(w, wo)); 317 d33_muld3(wi, frame, wi); 318 319 ASSERT(eq_eps(d3_dot(wi, w), cos_theta, fabs(cos_theta*1.e-6))); 320 321 if(pdf) *pdf = rdgfa_eval(rdgfa, wo, wi); 322 } 323 324 /******************************************************************************* 325 * Exported symbols 326 ******************************************************************************/ 327 const struct ssf_phase_type ssf_phase_rdgfa = { 328 rdgfa_init, 329 rdgfa_release, 330 rdgfa_sample, 331 rdgfa_eval, 332 rdgfa_eval, 333 sizeof(struct rdgfa), 334 ALIGNOF(struct rdgfa) 335 }; 336 337 res_T 338 ssf_phase_rdgfa_setup 339 (struct ssf_phase* phase, 340 const struct ssf_phase_rdgfa_setup_args* args) 341 { 342 struct rdgfa* rdgfa = NULL; 343 double lambda, Df, Rg, k, k2; 344 res_T res = RES_OK; 345 346 if(!phase 347 || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) 348 || !check_phase_rdgfa_setup_args(args)) { 349 res = RES_BAD_ARG; 350 goto error; 351 } 352 353 rdgfa = phase->data; 354 rdgfa->wavelength = args->wavelength; 355 rdgfa->fractal_dimension = args->fractal_dimension; 356 rdgfa->gyration_radius = args->gyration_radius; 357 rdgfa->nintervals = args->nintervals; 358 359 /* Fetch input data */ 360 lambda = rdgfa->wavelength; 361 Df = rdgfa->fractal_dimension; 362 Rg = rdgfa->gyration_radius; 363 364 /* Precompute constants */ 365 rdgfa->Rg2 = Rg*Rg; /* gyration_radius^2 */ 366 rdgfa->cst_4pi_div_lambda = 4*PI / lambda; 367 rdgfa->cst_3Df_div_2E = 3*Df/(2.0*EXP1); 368 369 /* Precompute the function g */ 370 k = (2.0*PI) / lambda; 371 k2 = k*k; 372 rdgfa->g = pow(1 + 4*k2*rdgfa->Rg2/(3*Df), -Df*0.5); 373 374 /* Precompute the phase function cumulative */ 375 switch(args->simd) { 376 case SSF_SIMD_NONE: 377 res = compute_cumulative(rdgfa); 378 if(res != RES_OK) goto error; 379 break; 380 case SSF_SIMD_128: 381 #ifdef SSF_USE_SIMD_128 382 res = compute_cumulative_simd4(rdgfa); 383 if(res != RES_OK) goto error; 384 #else 385 res = RES_BAD_OP; 386 goto error; 387 #endif 388 break; 389 case SSF_SIMD_256: 390 #ifdef SSF_USE_SIMD_256 391 res = compute_cumulative_simd8(rdgfa); 392 if(res != RES_OK) goto error; 393 #else 394 res = RES_BAD_OP; 395 goto error; 396 #endif 397 break; 398 default: FATAL("Unreachable code.\n"); break; 399 } 400 401 exit: 402 return res; 403 error: 404 goto exit; 405 } 406 407 res_T 408 ssf_phase_rdgfa_get_desc 409 (const struct ssf_phase* phase, 410 struct ssf_phase_rdgfa_desc* desc) 411 { 412 struct rdgfa* rdgfa = NULL; 413 res_T res = RES_OK; 414 415 if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) || !desc) { 416 res = RES_BAD_ARG; 417 goto error; 418 } 419 420 rdgfa = phase->data; 421 desc->wavelength = rdgfa->wavelength; 422 desc->gyration_radius = rdgfa->gyration_radius; 423 desc->fractal_dimension = rdgfa->fractal_dimension; 424 desc->normalization_factor = 1.0/rdgfa->rcp_normalize_factor; 425 desc->nintervals = rdgfa->nintervals; 426 427 exit: 428 return res; 429 error: 430 goto exit; 431 } 432 433 res_T 434 ssf_phase_rdgfa_get_interval 435 (const struct ssf_phase* phase, 436 const size_t interval_id, /* In [0, #intervals[ */ 437 struct ssf_phase_rdgfa_interval* interval) 438 { 439 struct rdgfa* rdgfa = NULL; 440 res_T res = RES_OK; 441 442 if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) || !interval) { 443 res = RES_BAD_ARG; 444 goto error; 445 } 446 447 rdgfa = phase->data; 448 if(interval_id >= darray_double_size_get(&rdgfa->cdf)) { 449 res = RES_BAD_ARG; 450 goto error; 451 } 452 453 interval->range[0] = rdgfa->dtheta * (double)(interval_id + 0); 454 interval->range[1] = rdgfa->dtheta * (double)(interval_id + 1); 455 interval->cumulative = darray_double_cdata_get(&rdgfa->cdf)[interval_id]; 456 457 exit: 458 return res; 459 error: 460 goto exit; 461 } 462