htrdr_ran_wlen_discrete.c (9375B)
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 #define _POSIX_C_SOURCE 200112L /* nextafter */ 25 26 #include "core/htrdr_c.h" 27 #include "core/htrdr_log.h" 28 #include "core/htrdr_ran_wlen_discrete.h" 29 30 #include <rsys/algorithm.h> 31 #include <rsys/dynamic_array_double.h> 32 #include <rsys/mem_allocator.h> 33 #include <rsys/ref_count.h> 34 35 #include <math.h> 36 37 struct htrdr_ran_wlen_discrete { 38 struct darray_double cumul; 39 struct darray_double proba; 40 struct darray_double radia; /* In W/m²/sr/m */ 41 struct darray_double wlens; /* In nm */ 42 double range[2]; /* Boundaries of the spectral integration interval in nm */ 43 size_t nbands; /* #bands */ 44 45 struct htrdr* htrdr; 46 ref_T ref; 47 }; 48 49 /******************************************************************************* 50 * Helper functions 51 ******************************************************************************/ 52 static INLINE res_T 53 check_htrdr_ran_wlen_discrete_create_args 54 (const struct htrdr_ran_wlen_discrete_create_args* args) 55 { 56 if(!args) return RES_BAD_ARG; 57 58 /* Invalid number of wavelength */ 59 if(!args->nwavelengths) 60 return RES_BAD_ARG; 61 62 /* Invalid functor */ 63 if(!args->get) 64 return RES_BAD_ARG; 65 66 return RES_OK; 67 } 68 69 static res_T 70 setup_per_wlen_radiance 71 (struct htrdr_ran_wlen_discrete* ran, 72 const struct htrdr_ran_wlen_discrete_create_args* args) 73 { 74 double* wlens = NULL; 75 double* radia = NULL; 76 size_t iwlen = 0; 77 res_T res = RES_OK; 78 ASSERT(ran && args); 79 80 res = darray_double_resize(&ran->wlens, args->nwavelengths); 81 if(res != RES_OK) { 82 htrdr_log_err(ran->htrdr, 83 "Error allocating discrete wavelength distribution wavelength list.\n"); 84 goto error; 85 } 86 res = darray_double_resize(&ran->radia, args->nwavelengths); 87 if(res != RES_OK) { 88 htrdr_log_err(ran->htrdr, 89 "Error allocating discrete wavelength distribution radiance list.\n"); 90 goto error; 91 } 92 93 wlens = darray_double_data_get(&ran->wlens); 94 radia = darray_double_data_get(&ran->radia); 95 96 /* Store the discrete values */ 97 FOR_EACH(iwlen, 0, args->nwavelengths) { 98 args->get(args->context, iwlen, wlens+iwlen, radia+iwlen); 99 100 if(iwlen > 0 && wlens[iwlen] <= wlens[iwlen-1]) { 101 htrdr_log_err(ran->htrdr, 102 "Failed to calculate discrete luminance distribution probabilities. " 103 "Wavelengths are not sorted in ascending order.\n"); 104 res = RES_BAD_ARG; 105 goto error; 106 } 107 } 108 109 /* Setup the spectral range */ 110 ran->range[0] = wlens[0]; 111 ran->range[1] = wlens[args->nwavelengths-1]; 112 113 exit: 114 return res; 115 error: 116 darray_double_clear(&ran->wlens); 117 darray_double_clear(&ran->radia); 118 goto exit; 119 } 120 121 static res_T 122 setup_distribution 123 (struct htrdr_ran_wlen_discrete* ran, 124 const struct htrdr_ran_wlen_discrete_create_args* args) 125 { 126 double* cumul = NULL; 127 double* proba = NULL; 128 double sum = 0; 129 size_t iband; 130 res_T res = RES_OK; 131 ASSERT(ran && check_htrdr_ran_wlen_discrete_create_args(args) == RES_OK); 132 ASSERT(ran->nbands >= 1); /* At least one band */ 133 134 res = darray_double_resize(&ran->cumul, ran->nbands); 135 if(res != RES_OK) { 136 htrdr_log_err(ran->htrdr, 137 "Error allocating the cumulative discrete wavelength distribution.\n"); 138 goto error; 139 } 140 res = darray_double_resize(&ran->proba, ran->nbands); 141 if(res != RES_OK) { 142 htrdr_log_err(ran->htrdr, 143 "Error allocating the discrete wavelength distribution probabilities.\n"); 144 goto error; 145 } 146 147 cumul = darray_double_data_get(&ran->cumul); 148 proba = darray_double_data_get(&ran->proba); 149 150 /* Compute the unormalized probabilities to sample a band */ 151 FOR_EACH(iband, 0, ran->nbands) { 152 const size_t iw0 = iband+0; 153 const size_t iw1 = iband+1; 154 double w0, L0; 155 double w1, L1; 156 double area; 157 158 args->get(args->context, iw0, &w0, &L0); 159 args->get(args->context, iw1, &w1, &L1); 160 ASSERT(w0 < w1); 161 162 area = (L0 + L1) * (w1-w0) * 0.5; 163 164 proba[iband] = area; 165 sum += area; 166 } 167 168 htrdr_log(ran->htrdr, "Discrete radiance integral = %g W/m²/sr\n", sum); 169 170 /* Normalize the probabilities and setup the cumulative */ 171 FOR_EACH(iband, 0, ran->nbands) { 172 proba[iband] /= sum; 173 if(iband == 0) { 174 cumul[iband] = proba[iband]; 175 } else { 176 cumul[iband] = proba[iband] + cumul[iband-1]; 177 ASSERT(cumul[iband] >= cumul[iband-1]); 178 } 179 } 180 ASSERT(eq_eps(cumul[ran->nbands-1], 1, 1e-6)); 181 cumul[ran->nbands-1] = 1.0; /* Fix numerical imprecision */ 182 183 exit: 184 return res; 185 error: 186 darray_double_clear(&ran->proba); 187 darray_double_clear(&ran->cumul); 188 goto exit; 189 } 190 191 static void 192 release_discrete(ref_T* ref) 193 { 194 struct htrdr_ran_wlen_discrete* ran = NULL; 195 struct htrdr* htrdr = NULL; 196 ASSERT(ref); 197 ran = CONTAINER_OF(ref, struct htrdr_ran_wlen_discrete, ref); 198 darray_double_release(&ran->cumul); 199 darray_double_release(&ran->proba); 200 darray_double_release(&ran->radia); 201 darray_double_release(&ran->wlens); 202 htrdr = ran->htrdr; 203 MEM_RM(htrdr_get_allocator(htrdr), ran); 204 htrdr_ref_put(htrdr); 205 } 206 207 /******************************************************************************* 208 * Exported functions 209 ******************************************************************************/ 210 res_T 211 htrdr_ran_wlen_discrete_create 212 (struct htrdr* htrdr, 213 const struct htrdr_ran_wlen_discrete_create_args* args, 214 struct htrdr_ran_wlen_discrete** out_ran) 215 { 216 struct htrdr_ran_wlen_discrete* ran = NULL; 217 res_T res = RES_OK; 218 ASSERT(htrdr && args && out_ran); 219 ASSERT(check_htrdr_ran_wlen_discrete_create_args(args) == RES_OK); 220 221 ran = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*ran)); 222 if(!ran) { 223 res = RES_MEM_ERR; 224 htrdr_log_err(htrdr, 225 "%s: error allocating discrete wavelength distribution\n", 226 FUNC_NAME); 227 goto error; 228 } 229 230 ref_init(&ran->ref); 231 darray_double_init(htrdr_get_allocator(htrdr), &ran->cumul); 232 darray_double_init(htrdr_get_allocator(htrdr), &ran->proba); 233 darray_double_init(htrdr_get_allocator(htrdr), &ran->radia); 234 darray_double_init(htrdr_get_allocator(htrdr), &ran->wlens); 235 htrdr_ref_get(htrdr); 236 ran->htrdr = htrdr; 237 238 ran->nbands = args->nwavelengths - 1; 239 240 res = setup_per_wlen_radiance(ran, args); 241 if(res != RES_OK) goto error; 242 if(ran->nbands != 0) { 243 res = setup_distribution(ran, args); 244 if(res != RES_OK) goto error; 245 } 246 247 exit: 248 *out_ran = ran; 249 return res; 250 error: 251 if(ran) { htrdr_ran_wlen_discrete_ref_put(ran); ran = NULL; } 252 goto exit; 253 } 254 255 void 256 htrdr_ran_wlen_discrete_ref_get(struct htrdr_ran_wlen_discrete* ran) 257 { 258 ASSERT(ran); 259 ref_get(&ran->ref); 260 } 261 262 void 263 htrdr_ran_wlen_discrete_ref_put(struct htrdr_ran_wlen_discrete* ran) 264 { 265 ASSERT(ran); 266 ref_put(&ran->ref, release_discrete); 267 } 268 269 double 270 htrdr_ran_wlen_discrete_sample 271 (struct htrdr_ran_wlen_discrete* ran, 272 const double r0, 273 const double r1, 274 double* pdf) /* In nm⁻¹ */ 275 { 276 double lambda = 0; 277 278 ASSERT(ran); 279 ASSERT(0 <= r0 && r0 < 1); 280 ASSERT(0 <= r1 && r1 < 1); 281 282 if(ran->nbands == 0) { 283 ASSERT(darray_double_size_get(&ran->wlens) == 1); 284 lambda = darray_double_cdata_get(&ran->wlens)[0]; 285 if(pdf) *pdf = 1; 286 287 } else { 288 double* find = NULL; 289 const double* proba = NULL; 290 const double* cumul = NULL; 291 const double* wlens = NULL; 292 double w0, w1; /* Wavelengths of the sampled band in nm */ 293 size_t iband = 0; 294 double r0_next = nextafter(r0, DBL_MAX); 295 ASSERT(0 <= r0 && r0 < 1); 296 ASSERT(0 <= r1 && r1 < 1); 297 298 cumul = darray_double_cdata_get(&ran->cumul); 299 proba = darray_double_cdata_get(&ran->proba); 300 wlens = darray_double_cdata_get(&ran->wlens); 301 302 /* Sample a band. Use r0_next rather than r0 to find the first entry that is 303 * not less than *or equal* to r0 */ 304 find = search_lower_bound 305 (&r0_next, cumul, ran->nbands, sizeof(double), cmp_dbl); 306 ASSERT(find); 307 308 iband = (size_t)(find - cumul); 309 ASSERT(iband < ran->nbands); 310 ASSERT(cumul[iband] > r0 && (!iband || cumul[iband-1] <= r0)); 311 312 /* Retrieve the boundaries of the sampled band */ 313 w0 = wlens[iband+0]; 314 w1 = wlens[iband+1]; 315 316 /* Uniformely sample the wavelength in [w0, w1[ */ 317 lambda = w0 + r1 * (w1 - w0); 318 319 if(pdf) { 320 const double pdf_wlen = 1.f / (w1-w0); 321 *pdf = proba[iband] * pdf_wlen; 322 } 323 } 324 return lambda; 325 }