htsky_svx.c (13884B)
1 /* Copyright (C) 2018, 2019, 2020, 2021 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018, 2019 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018, 2019 Université Paul Sabatier 4 * 5 * This program is free software: you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation, either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 17 18 #define _POSIX_C_SOURCE 200112L /* nextafter */ 19 20 #include "htsky.h" 21 #include "htsky_atmosphere.h" 22 #include "htsky_c.h" 23 #include "htsky_cloud.h" 24 #include "htsky_log.h" 25 #include "htsky_svx.h" 26 27 #include <star/svx.h> 28 29 #include <math.h> 30 31 /******************************************************************************* 32 * Helper functions 33 ******************************************************************************/ 34 /* Smits intersection: "Efficiency issues for ray tracing" */ 35 static FINLINE void 36 ray_intersect_aabb 37 (const double org[3], 38 const double dir[3], 39 const double range[2], 40 const double low[3], 41 const double upp[3], 42 const int axis_mask, 43 double hit_range[2]) 44 { 45 double hit[2]; 46 int i; 47 ASSERT(org && dir && range && low && upp && hit_range); 48 ASSERT(low[0] < upp[0]); 49 ASSERT(low[1] < upp[1]); 50 ASSERT(low[2] < upp[2]); 51 52 hit_range[0] = INF; 53 hit_range[1] =-INF; 54 hit[0] = range[0]; 55 hit[1] = range[1]; 56 FOR_EACH(i, 0, 3) { 57 double t_min, t_max; 58 59 if(!(BIT(i) & axis_mask)) continue; 60 61 t_min = (low[i] - org[i]) / dir[i]; 62 t_max = (upp[i] - org[i]) / dir[i]; 63 64 if(t_min > t_max) SWAP(double, t_min, t_max); 65 hit[0] = MMAX(t_min, hit[0]); 66 hit[1] = MMIN(t_max, hit[1]); 67 if(hit[0] > hit[1]) return; 68 } 69 hit_range[0] = hit[0]; 70 hit_range[1] = hit[1]; 71 } 72 73 static res_T 74 infinite_cloudy_slab_trace_ray 75 (struct htsky* sky, 76 struct svx_tree* clouds, 77 const double org[3], 78 const double dir[3], 79 const double range[2], 80 const size_t max_steps, 81 svx_hit_challenge_T challenge, 82 svx_hit_filter_T filter, 83 void* context, 84 struct svx_hit* hit) 85 { 86 double pos[2]; 87 double org_cs[3]; /* Origin of the ray transformed in local cell space */ 88 const double* cell_low; 89 const double* cell_upp; 90 double cell_low_ws[3]; /* Cell lower bound in world space */ 91 double cell_upp_ws[3]; /* Cell upper bound in world space */ 92 double cell_sz[3]; /* Size of a cell */ 93 double t_max[3], t_delta[2], t_min_z; 94 size_t istep; 95 int64_t xy[2]; /* 2D index of the repeated cell */ 96 int incr[2]; /* Index increment */ 97 res_T res = RES_OK; 98 ASSERT(sky && clouds && org && dir && range && context && hit); 99 ASSERT(range[0] < range[1]); 100 101 cell_low = sky->htcp_desc.lower; 102 cell_upp = sky->htcp_desc.upper; 103 104 /* Check that the ray intersects the slab */ 105 t_min_z = (cell_low[2] - org[2]) / dir[2]; 106 t_max[2] = (cell_upp[2] - org[2]) / dir[2]; 107 if(t_min_z > t_max[2]) SWAP(double, t_min_z, t_max[2]); 108 t_min_z = MMAX(t_min_z, range[0]); 109 t_max[2] = MMIN(t_max[2], range[1]); 110 if(t_min_z > t_max[2]) return RES_OK; 111 112 /* Compute the size of a cell */ 113 cell_sz[0] = cell_upp[0] - cell_low[0]; 114 cell_sz[1] = cell_upp[1] - cell_low[1]; 115 cell_sz[2] = cell_upp[2] - cell_low[2]; 116 117 /* Define the 2D index of the current cell. (0,0) is the index of the 118 * non duplicated cell */ 119 pos[0] = org[0] + t_min_z*dir[0]; 120 pos[1] = org[1] + t_min_z*dir[1]; 121 xy[0] = (int64_t)floor((pos[0] - cell_low[0]) / cell_sz[0]); 122 xy[1] = (int64_t)floor((pos[1] - cell_low[1]) / cell_sz[1]); 123 124 /* Define the 2D index increment wrt dir sign */ 125 incr[0] = dir[0] < 0 ? -1 : 1; 126 incr[1] = dir[1] < 0 ? -1 : 1; 127 128 /* Compute the world space AABB of the repeated cell currently hit */ 129 cell_low_ws[0] = cell_low[0] + (double)xy[0]*cell_sz[0]; 130 cell_low_ws[1] = cell_low[1] + (double)xy[1]*cell_sz[1]; 131 cell_low_ws[2] = cell_low[2]; 132 cell_upp_ws[0] = cell_low_ws[0] + cell_sz[0]; 133 cell_upp_ws[1] = cell_low_ws[1] + cell_sz[1]; 134 cell_upp_ws[2] = cell_upp[2]; 135 136 /* Compute the max ray intersection with the current cell */ 137 t_max[0] = ((dir[0]<0 ? cell_low_ws[0] : cell_upp_ws[0]) - org[0]) / dir[0]; 138 t_max[1] = ((dir[1]<0 ? cell_low_ws[1] : cell_upp_ws[1]) - org[1]) / dir[1]; 139 /*t_max[2] = ((dir[2]<0 ? cell_low_ws[2] : cell_upp_ws[2]) - org[2]) / dir[2];*/ 140 ASSERT(t_max[0] >= 0 && t_max[1] >= 0 && t_max[2] >= 0); 141 142 /* Compute the distance along the ray to traverse in order to move of a 143 * distance equal to the cloud size along the X and Y axis */ 144 t_delta[0] = (dir[0]<0 ? -cell_sz[0] : cell_sz[0]) / dir[0]; 145 t_delta[1] = (dir[1]<0 ? -cell_sz[1] : cell_sz[1]) / dir[1]; 146 ASSERT(t_delta[0] >= 0 && t_delta[1] >= 0); 147 148 org_cs[2] = org[2]; 149 FOR_EACH(istep, 0, max_steps) { 150 int iaxis; 151 152 /* Transform the ray origin in the local cell space */ 153 org_cs[0] = org[0] - (double)xy[0]*cell_sz[0]; 154 org_cs[1] = org[1] - (double)xy[1]*cell_sz[1]; 155 156 res = svx_tree_trace_ray 157 (clouds, org_cs, dir, range, challenge, filter, context, hit); 158 if(res != RES_OK) { 159 log_err(sky,"%s: could not trace the ray in the repeated cloud.\n", 160 FUNC_NAME); 161 goto error; 162 } 163 if(!SVX_HIT_NONE(hit)) goto exit; 164 165 /* Define the next axis to traverse */ 166 iaxis = t_max[0] < t_max[1] 167 ? (t_max[0] < t_max[2] ? 0 : 2) 168 : (t_max[1] < t_max[2] ? 1 : 2); 169 170 if(iaxis == 2) break; /* The ray traverse the slab */ 171 172 if(t_max[iaxis] >= range[1]) break; /* Out of bound */ 173 174 t_max[iaxis] += t_delta[iaxis]; 175 176 /* Define the 2D index of the next traversed cloud */ 177 xy[iaxis] += incr[iaxis]; 178 } 179 180 exit: 181 return res; 182 error: 183 goto exit; 184 } 185 186 187 /******************************************************************************* 188 * Exported functions 189 ******************************************************************************/ 190 double 191 htsky_fetch_svx_property 192 (const struct htsky* sky, 193 const enum htsky_property prop, 194 const enum htsky_svx_op op, 195 const int components_mask, /* Combination of htsky_component_flag */ 196 const size_t iband, /* Index of the spectral band */ 197 const size_t iquad, /* Index of the quadrature point in the spectral band */ 198 const double pos[3]) 199 { 200 struct svx_voxel voxel = SVX_VOXEL_NULL; 201 struct atmosphere* atmosphere = NULL; 202 struct cloud* cloud = NULL; 203 size_t i; 204 int in_clouds; /* Defines if `pos' lies in the clouds */ 205 int in_atmosphere; /* Defines if `pos' lies in the atmosphere */ 206 int comp_mask = components_mask; 207 ASSERT(sky && pos); 208 ASSERT(comp_mask & HTSKY_CPNT_MASK_ALL); 209 ASSERT(iband >= sky->bands_range[0]); 210 ASSERT(iband <= sky->bands_range[1]); 211 ASSERT(iquad < htsky_get_spectral_band_quadrature_length(sky, iband)); 212 213 i = iband - sky->bands_range[0]; 214 cloud = sky->is_cloudy ? &sky->clouds[i][iquad] : NULL; 215 atmosphere = &sky->atmosphere[i][iquad]; 216 217 /* Is the position inside the clouds? */ 218 if(sky->is_cloudy) { 219 in_clouds = 0; 220 } else if(sky->repeat_clouds) { 221 in_clouds = 222 pos[2] >= cloud->octree_desc.lower[2] 223 && pos[2] <= cloud->octree_desc.upper[2]; 224 } else { 225 in_clouds = 226 pos[0] >= cloud->octree_desc.lower[0] 227 && pos[1] >= cloud->octree_desc.lower[1] 228 && pos[2] >= cloud->octree_desc.lower[2] 229 && pos[0] <= cloud->octree_desc.upper[0] 230 && pos[1] <= cloud->octree_desc.upper[1] 231 && pos[2] <= cloud->octree_desc.upper[2]; 232 } 233 234 ASSERT(atmosphere->bitree_desc.frame[0] == SVX_AXIS_Z); 235 in_atmosphere = 236 pos[2] >= atmosphere->bitree_desc.lower[2] 237 && pos[2] <= atmosphere->bitree_desc.upper[2]; 238 239 if(!in_clouds) { /* Not in clouds => No particle */ 240 comp_mask &= ~HTSKY_CPNT_FLAG_PARTICLES; 241 } 242 if(!in_atmosphere) { /* Not in atmosphere => No gas */ 243 comp_mask &= ~HTSKY_CPNT_FLAG_GAS; 244 } 245 246 if(!in_clouds && !in_atmosphere) /* In vacuum */ 247 return 0; 248 249 if(!in_clouds) { 250 ASSERT(in_atmosphere); 251 SVX(tree_at(atmosphere->bitree, pos, NULL, NULL, &voxel)); 252 } else { 253 double pos_cs[3]; 254 world_to_cloud(sky, pos, pos_cs); 255 SVX(tree_at(cloud->octree, pos_cs, NULL, NULL, &voxel)); 256 } 257 258 return htsky_fetch_svx_voxel_property 259 (sky, prop, op, comp_mask, iband, iquad, &voxel); 260 } 261 262 double 263 htsky_fetch_svx_voxel_property 264 (const struct htsky* sky, 265 const enum htsky_property prop, 266 const enum htsky_svx_op op, 267 const int components_mask, 268 const size_t iband, /* Index of the spectral band */ 269 const size_t iquad, /* Index of the quadrature point in the spectral band */ 270 const struct svx_voxel* voxel) 271 { 272 double gas = 0; 273 double par = 0; 274 int comp_mask = components_mask; 275 ASSERT(sky && voxel); 276 ASSERT((unsigned)prop < HTSKY_PROPS_COUNT__); 277 ASSERT((unsigned)op < HTSKY_SVX_OPS_COUNT__); 278 (void)sky, (void)iband, (void)iquad; 279 280 /* Check if the voxel has infinite bounds/degenerated. In such case it is 281 * atmospheric voxel with only gas properties */ 282 if(IS_INF(voxel->upper[0]) || voxel->lower[0] > voxel->upper[0]) { 283 ASSERT(IS_INF(voxel->upper[1]) || voxel->lower[1] > voxel->upper[1]); 284 comp_mask &= ~HTSKY_CPNT_FLAG_PARTICLES; 285 } 286 287 if(comp_mask & HTSKY_CPNT_FLAG_PARTICLES) { 288 par = vox_get(voxel->data, HTSKY_CPNT_PARTICLES, prop, op); 289 } 290 if(comp_mask & HTSKY_CPNT_FLAG_GAS) { 291 gas = vox_get(voxel->data, HTSKY_CPNT_GAS, prop, op); 292 } 293 return par + gas; 294 } 295 296 res_T 297 htsky_trace_ray 298 (struct htsky* sky, 299 const double org[3], 300 const double dir[3], /* Must be normalized */ 301 const double range[2], 302 const svx_hit_challenge_T challenge, /* NULL <=> Traversed up to the leaves */ 303 const svx_hit_filter_T filter, /* NULL <=> Stop RT at the 1st hit voxel */ 304 void* context, /* Data sent to the filter functor */ 305 const size_t iband, 306 const size_t iquad, 307 struct svx_hit* hit) 308 { 309 double cloud_range[2]; 310 struct svx_tree* clouds; 311 struct svx_tree* atmosphere; 312 size_t i; 313 enum { AXIS_X = BIT(0), AXIS_Y = BIT(1), AXIS_Z = BIT(2) }; 314 res_T res = RES_OK; 315 316 if(!sky || !org || !dir || !range || !hit) { 317 res = RES_BAD_ARG; 318 goto error; 319 } 320 321 if(iband < sky->bands_range[0] || iband > sky->bands_range[1]) { 322 log_err(sky, 323 "%s: invalid spectral band index `%lu'. " 324 "Valid short wave spectral bands lie in [%lu, %lu]\n", 325 FUNC_NAME, 326 (unsigned long)iband, 327 (unsigned long)sky->bands_range[0], 328 (unsigned long)sky->bands_range[1]); 329 res = RES_BAD_ARG; 330 goto error; 331 } 332 333 if(iquad >= htsky_get_spectral_band_quadrature_length(sky, iband)) { 334 log_err(sky, 335 "invalid quadrature point `%lu' for the spectral band `%lu'.\n", 336 (unsigned long)iquad, (unsigned long)iband); 337 res = RES_BAD_ARG; 338 goto error; 339 } 340 341 /* Fetch the clouds/atmosphere corresponding to the submitted spectral data */ 342 i = iband - sky->bands_range[0]; 343 clouds = sky->is_cloudy ? sky->clouds[i][iquad].octree : NULL; 344 atmosphere = sky->atmosphere[i][iquad].bitree; 345 346 cloud_range[0] = INF; 347 cloud_range[1] =-INF; 348 349 if(sky->is_cloudy) { 350 /* Compute the ray range, intersecting the clouds AABB */ 351 if(sky->repeat_clouds) { 352 ray_intersect_aabb(org, dir, range, sky->htcp_desc.lower, 353 sky->htcp_desc.upper, AXIS_Z, cloud_range); 354 } else { 355 ray_intersect_aabb(org, dir, range, sky->htcp_desc.lower, 356 sky->htcp_desc.upper, AXIS_X|AXIS_Y|AXIS_Z, cloud_range); 357 } 358 } 359 360 /* Reset the hit */ 361 *hit = SVX_HIT_NULL; 362 363 if(cloud_range[0] >= cloud_range[1]) { /* The ray does not traverse the clouds */ 364 res = svx_tree_trace_ray(atmosphere, org, dir, range, challenge, filter, 365 context, hit); 366 if(res != RES_OK) { 367 log_err(sky, "%s: could not trace the ray in the atmosphere.\n", FUNC_NAME); 368 goto error; 369 } 370 } else { /* The ray may traverse the clouds */ 371 double range_adjusted[2]; 372 373 if(cloud_range[0] > range[0]) { /* The ray begins in the atmosphere */ 374 /* Trace a ray in the atmosphere from range[0] to cloud_range[0] */ 375 range_adjusted[0] = range[0]; 376 range_adjusted[1] = nextafter(cloud_range[0], -DBL_MAX); 377 res = svx_tree_trace_ray(atmosphere, org, dir, range_adjusted, challenge, 378 filter, context, hit); 379 if(res != RES_OK) { 380 log_err(sky, 381 "%s: could not to trace the part that begins in the atmosphere.\n", 382 FUNC_NAME); 383 goto error; 384 } 385 if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ 386 } 387 388 /* Pursue ray traversal into the clouds */ 389 if(!sky->repeat_clouds) { 390 res = svx_tree_trace_ray(clouds, org, dir, cloud_range, challenge, filter, 391 context, hit); 392 if(res != RES_OK) { 393 log_err(sky, 394 "%s: could not trace the ray part that intersects the clouds.\n", 395 FUNC_NAME); 396 goto error; 397 } 398 if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ 399 400 /* Clouds are infinitely repeated along the X and Y axis */ 401 } else { 402 403 res = infinite_cloudy_slab_trace_ray(sky, clouds, org, dir, cloud_range, 404 32, challenge, filter, context, hit); 405 if(res != RES_OK) goto error; 406 407 if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ 408 } 409 410 /* Pursue ray traversal into the atmosphere */ 411 range_adjusted[0] = nextafter(cloud_range[1], DBL_MAX); 412 range_adjusted[1] = range[1]; 413 res = svx_tree_trace_ray(atmosphere, org, dir, range_adjusted, challenge, 414 filter, context, hit); 415 if(res != RES_OK) { 416 log_err(sky, 417 "%s: could not trace the ray part that ends in the atmosphere.\n", 418 FUNC_NAME); 419 goto error; 420 } 421 if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ 422 } 423 424 exit: 425 return res; 426 error: 427 goto exit; 428 } 429 430