htrdr_slab.c (5238B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "core/htrdr.h" 25 #include "core/htrdr_log.h" 26 #include "core/htrdr_slab.h" 27 28 #include <rsys/cstr.h> 29 #include <math.h> 30 31 /******************************************************************************* 32 * Local function 33 ******************************************************************************/ 34 res_T 35 htrdr_slab_trace_ray 36 (struct htrdr* htrdr, 37 const double org[3], 38 const double dir[3], 39 const double range[2], 40 const double cell_low[3], 41 const double cell_upp[3], 42 htrdr_trace_cell_T trace_cell, 43 const size_t max_steps, 44 void* trace_cell_context) 45 { 46 double pos[2]; 47 double org_cs[3]; /* Origin of the ray transformed in local cell space */ 48 double cell_low_ws[3]; /* Cell lower bound in world space */ 49 double cell_upp_ws[3]; /* Cell upper bound in world space */ 50 double cell_sz[3]; /* Size of a cell */ 51 double t_max[3], t_delta[2], t_min_z; 52 size_t istep; 53 int64_t xy[2]; /* 2D index of the repeated cell */ 54 int incr[2]; /* Index increment */ 55 res_T res = RES_OK; 56 ASSERT(htrdr && org && dir && range && cell_low && cell_upp && trace_cell); 57 ASSERT(range[0] < range[1]); 58 59 /* Check that the ray intersects the slab */ 60 t_min_z = (cell_low[2] - org[2]) / dir[2]; 61 t_max[2] = (cell_upp[2] - org[2]) / dir[2]; 62 if(t_min_z > t_max[2]) SWAP(double, t_min_z, t_max[2]); 63 t_min_z = MMAX(t_min_z, range[0]); 64 t_max[2] = MMIN(t_max[2], range[1]); 65 if(t_min_z > t_max[2]) return RES_OK; 66 67 /* Compute the size of a cell */ 68 cell_sz[0] = cell_upp[0] - cell_low[0]; 69 cell_sz[1] = cell_upp[1] - cell_low[1]; 70 cell_sz[2] = cell_upp[2] - cell_low[2]; 71 72 /* Define the 2D index of the current cell. (0,0) is the index of the 73 * non duplicated cell */ 74 pos[0] = org[0] + t_min_z*dir[0]; 75 pos[1] = org[1] + t_min_z*dir[1]; 76 xy[0] = (int64_t)floor((pos[0] - cell_low[0]) / cell_sz[0]); 77 xy[1] = (int64_t)floor((pos[1] - cell_low[1]) / cell_sz[1]); 78 79 /* Define the 2D index increment wrt dir sign */ 80 incr[0] = dir[0] < 0 ? -1 : 1; 81 incr[1] = dir[1] < 0 ? -1 : 1; 82 83 /* Compute the world space AABB of the repeated cell currently hit */ 84 cell_low_ws[0] = cell_low[0] + (double)xy[0]*cell_sz[0]; 85 cell_low_ws[1] = cell_low[1] + (double)xy[1]*cell_sz[1]; 86 cell_low_ws[2] = cell_low[2]; 87 cell_upp_ws[0] = cell_low_ws[0] + cell_sz[0]; 88 cell_upp_ws[1] = cell_low_ws[1] + cell_sz[1]; 89 cell_upp_ws[2] = cell_upp[2]; 90 91 /* Compute the max ray intersection with the current cell */ 92 t_max[0] = ((dir[0]<0 ? cell_low_ws[0] : cell_upp_ws[0]) - org[0]) / dir[0]; 93 t_max[1] = ((dir[1]<0 ? cell_low_ws[1] : cell_upp_ws[1]) - org[1]) / dir[1]; 94 /*t_max[2] = ((dir[2]<0 ? cell_low_ws[2] : cell_upp_ws[2]) - org[2]) / dir[2];*/ 95 ASSERT(t_max[0] >= 0 && t_max[1] >= 0 && t_max[2] >= 0); 96 97 /* Compute the distance along the ray to traverse in order to move of a 98 * distance equal to the cloud size along the X and Y axis */ 99 t_delta[0] = (dir[0]<0 ? -cell_sz[0] : cell_sz[0]) / dir[0]; 100 t_delta[1] = (dir[1]<0 ? -cell_sz[1] : cell_sz[1]) / dir[1]; 101 ASSERT(t_delta[0] >= 0 && t_delta[1] >= 0); 102 103 org_cs[2] = org[2]; 104 FOR_EACH(istep, 0, max_steps) { 105 int iaxis; 106 int hit; 107 108 /* Transform the ray origin in the local cell space */ 109 org_cs[0] = org[0] - (double)xy[0]*cell_sz[0]; 110 org_cs[1] = org[1] - (double)xy[1]*cell_sz[1]; 111 112 res = trace_cell(org_cs, dir, range, trace_cell_context, &hit); 113 if(res != RES_OK) { 114 htrdr_log_err(htrdr, 115 "%s: could not trace the ray in the repeated cells -- %s.\n", 116 FUNC_NAME, res_to_cstr(res)); 117 goto error; 118 } 119 if(hit) goto exit; 120 121 /* Define the next axis to traverse */ 122 iaxis = t_max[0] < t_max[1] 123 ? (t_max[0] < t_max[2] ? 0 : 2) 124 : (t_max[1] < t_max[2] ? 1 : 2); 125 126 if(iaxis == 2) break; /* The ray traverse the slab */ 127 128 if(t_max[iaxis] >= range[1]) break; /* Out of bound */ 129 130 t_max[iaxis] += t_delta[iaxis]; 131 132 /* Define the 2D index of the next traversed cloud */ 133 xy[iaxis] += incr[iaxis]; 134 } 135 136 exit: 137 return res; 138 error: 139 goto exit; 140 } 141 142