rnatm_octrees_storage.c (19392B)
1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace 3 * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris 4 * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2022, 2023, 2025 Observatoire de Paris 6 * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne 7 * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin 8 * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier 9 * 10 * This program is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU General Public License as published by 12 * the Free Software Foundation, either version 3 of the License, or 13 * (at your option) any later version. 14 * 15 * This program is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU General Public License for more details. 19 * 20 * You should have received a copy of the GNU General Public License 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 22 23 #define _POSIX_C_SOURCE 200112L /* nextafter */ 24 25 #include "rnatm_c.h" 26 #include "rnatm_log.h" 27 #include "rnatm_octrees_storage.h" 28 29 #include <star/sars.h> 30 #include <star/sck.h> 31 #include <star/suvm.h> 32 #include <star/svx.h> 33 34 #include <rsys/cstr.h> 35 36 #include <math.h> /* nextafter */ 37 38 /******************************************************************************* 39 * Helper functions 40 ******************************************************************************/ 41 static int 42 cmp_hash256(const void* a, const void* b) 43 { 44 const char* hash0 = *(const hash256_T*)a; 45 const char* hash1 = *(const hash256_T*)b; 46 size_t i; 47 FOR_EACH(i, 0, sizeof(hash256_T)) { 48 if(hash0[i] < hash1[i]) { 49 return -1; 50 } else if(hash0[i] > hash1[i]) { 51 return +1; 52 } 53 } 54 return 0; 55 } 56 57 static INLINE res_T 58 compute_gas_mesh_signature(const struct gas* gas, hash256_T hash) 59 { 60 struct suvm_mesh_desc mesh; 61 res_T res = RES_OK; 62 ASSERT(gas && hash); 63 64 res = suvm_volume_get_mesh_desc(gas->volume, &mesh); 65 if(res != RES_OK) goto error; 66 res = suvm_mesh_desc_compute_hash(&mesh, hash); 67 if(res != RES_OK) goto error; 68 69 exit: 70 return res; 71 error: 72 goto exit; 73 } 74 75 static INLINE res_T 76 compute_aerosol_mesh_signature(const struct aerosol* aerosol, hash256_T hash) 77 { 78 struct suvm_mesh_desc mesh; 79 res_T res = RES_OK; 80 ASSERT(aerosol &&hash); 81 82 res = suvm_volume_get_mesh_desc(aerosol->volume, &mesh); 83 if(res != RES_OK) goto error; 84 res = suvm_mesh_desc_compute_hash(&mesh, hash); 85 if(res != RES_OK) goto error; 86 87 exit: 88 return res; 89 error: 90 goto exit; 91 } 92 93 static res_T 94 compute_mesh_signature(const struct rnatm* atm, hash256_T mesh_signature) 95 { 96 hash256_T gas_signature; 97 hash256_T* aerosol_signatures = NULL; 98 size_t naerosols; 99 size_t i; 100 res_T res = RES_OK; 101 ASSERT(atm && mesh_signature); 102 103 /* Compute the mesh signature of the gas */ 104 res = compute_gas_mesh_signature(&atm->gas, gas_signature); 105 if(res != RES_OK) goto error; 106 107 naerosols = darray_aerosol_size_get(&atm->aerosols); 108 109 if(!naerosols) { 110 memcpy(mesh_signature, gas_signature, sizeof(hash256_T)); 111 } else { 112 struct sha256_ctx ctx; 113 114 aerosol_signatures = MEM_CALLOC(atm->allocator, naerosols, sizeof(hash256_T)); 115 if(!aerosol_signatures) { res = RES_MEM_ERR; goto error; } 116 117 FOR_EACH(i, 0, naerosols) { 118 const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+i; 119 120 res = compute_aerosol_mesh_signature(aerosol, aerosol_signatures[i]); 121 if(res != RES_OK) goto error; 122 } 123 124 /* Pay attention to the order of how aerosol hashes are registered. 125 * Therefore, the same set of aerosols will have the same signature (i.e. the 126 * same hash), regardless of the order in which the aerosols were listed in 127 * the input arguments */ 128 qsort(aerosol_signatures, naerosols, sizeof(hash256_T), cmp_hash256); 129 130 /* Build the signature for the dataset */ 131 sha256_ctx_init(&ctx); 132 sha256_ctx_update(&ctx, gas_signature, sizeof(hash256_T)); 133 FOR_EACH(i, 0, naerosols) { 134 sha256_ctx_update(&ctx, aerosol_signatures[i], sizeof(hash256_T)); 135 } 136 sha256_ctx_finalize(&ctx, mesh_signature); 137 } 138 139 exit: 140 if(aerosol_signatures) MEM_RM(atm->allocator, aerosol_signatures); 141 return res; 142 error: 143 log_err(atm, "mesh signature calculation error -- %s\n", res_to_cstr(res)); 144 goto exit; 145 } 146 147 static res_T 148 compute_aerosol_sars_signature 149 (const struct aerosol* aerosol, 150 const double range[2], /* In nm. Limits are inclusive */ 151 hash256_T hash) 152 { 153 size_t ibands[2]; 154 res_T res = RES_OK; 155 ASSERT(aerosol && range && hash && range[0] < range[1]); 156 157 res = sars_find_bands(aerosol->sars, range, ibands); 158 if(res != RES_OK) goto error; 159 160 /* No band are covered by the spectral range */ 161 if(ibands[0] > ibands[1]) { 162 memset(hash, 0, sizeof(hash256_T)); 163 164 /* Hash the data of the covered aerosol bands */ 165 } else { 166 struct sha256_ctx ctx; 167 size_t iband; 168 169 sha256_ctx_init(&ctx); 170 FOR_EACH(iband, ibands[0], ibands[1]+1) { 171 res = sars_band_compute_hash(aerosol->sars, iband, hash); 172 if(res != RES_OK) goto error; 173 sha256_ctx_update(&ctx, hash, sizeof(hash256_T)); 174 } 175 sha256_ctx_finalize(&ctx, hash); 176 } 177 178 exit: 179 return res; 180 error: 181 goto exit; 182 } 183 184 static res_T 185 compute_band_signature 186 (const struct rnatm* atm, 187 const size_t iband, 188 hash256_T band_signature) 189 { 190 hash256_T gas_signature; 191 hash256_T* aerosol_signatures = NULL; 192 size_t naerosols = 0; 193 size_t i; 194 res_T res = RES_OK; 195 ASSERT(atm && band_signature); 196 197 /* Compute the signature for the gas */ 198 res = sck_band_compute_hash(atm->gas.ck, iband, gas_signature); 199 if(res != RES_OK) goto error; 200 201 naerosols = darray_aerosol_size_get(&atm->aerosols); 202 203 if(!naerosols) { 204 memcpy(band_signature, gas_signature, sizeof(hash256_T)); 205 } else { 206 struct sha256_ctx ctx; 207 struct sck_band band = SCK_BAND_NULL; 208 double range[2]; 209 210 aerosol_signatures = MEM_CALLOC(atm->allocator, naerosols, sizeof(hash256_T)); 211 if(!aerosol_signatures) { res = RES_MEM_ERR; goto error; } 212 213 /* Retrieve the spectral range of the band */ 214 SCK(get_band(atm->gas.ck, iband, &band)); 215 range[0] = band.lower; 216 range[1] = nextafter(band.upper, DBL_MAX); /* Make limit inclusive */ 217 218 FOR_EACH(i, 0, naerosols) { 219 const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+i; 220 221 res = compute_aerosol_sars_signature(aerosol, range, aerosol_signatures[i]); 222 if(res != RES_OK) goto error; 223 } 224 225 /* Pay attention to the order of how aerosol hashes are registered. 226 * Therefore, the same set of aerosols will have the same signature (i.e. the 227 * same hash), regardless of the order in which the aerosols were listed in 228 * the input arguments */ 229 qsort(aerosol_signatures, naerosols, sizeof(hash256_T), cmp_hash256); 230 231 /* Build the signature for the dataset */ 232 sha256_ctx_init(&ctx); 233 sha256_ctx_update(&ctx, gas_signature, sizeof(hash256_T)); 234 FOR_EACH(i, 0, naerosols) { 235 sha256_ctx_update(&ctx, aerosol_signatures[i], sizeof(hash256_T)); 236 } 237 sha256_ctx_finalize(&ctx, band_signature); 238 } 239 240 exit: 241 if(aerosol_signatures) MEM_RM(atm->allocator, aerosol_signatures); 242 return res; 243 error: 244 log_err(atm, "signature calculation error for the band %lu -- %s\n", 245 (unsigned long)iband, res_to_cstr(res)); 246 goto exit; 247 } 248 249 /******************************************************************************* 250 * Local functions 251 ******************************************************************************/ 252 res_T 253 octrees_storage_desc_init 254 (const struct rnatm* atm, 255 struct octrees_storage_desc* desc) 256 { 257 size_t iband; 258 size_t nbands; 259 res_T res = RES_OK; 260 ASSERT(atm && desc); 261 262 *desc = OCTREES_STORAGE_DESC_NULL; 263 264 desc->version = OCTREES_STORAGE_VERSION; 265 desc->ibands[0] = atm->ibands[0]; 266 desc->ibands[1] = atm->ibands[1]; 267 desc->optical_thickness = atm->optical_thickness; 268 desc->grid_definition[0] = atm->grid_definition[0]; 269 desc->grid_definition[1] = atm->grid_definition[1]; 270 desc->grid_definition[2] = atm->grid_definition[2]; 271 desc->allocator = atm->allocator; 272 273 log_info(atm, "sign octree data\n"); 274 275 /* Calculate the signature of volumetric meshes */ 276 res = compute_mesh_signature(atm, desc->mesh_signature); 277 if(res != RES_OK) goto error; 278 279 /* Allocate the list of band signatures */ 280 nbands = atm->ibands[1] - atm->ibands[0] + 1; 281 desc->band_signatures = MEM_CALLOC(desc->allocator, nbands, sizeof(hash256_T)); 282 if(!desc->band_signatures) { 283 log_err(atm, "error allocating signatures per band\n"); 284 res = RES_MEM_ERR; 285 goto error; 286 } 287 288 /* Calculate the signature of spectral bands */ 289 FOR_EACH(iband, atm->ibands[0], atm->ibands[1]+1) { 290 const size_t i = iband - atm->ibands[0]; 291 res = compute_band_signature(atm, iband, desc->band_signatures[i]); 292 if(res != RES_OK) goto error; 293 } 294 295 exit: 296 return res; 297 error: 298 goto exit; 299 } 300 301 res_T 302 octrees_storage_desc_init_from_stream 303 (const struct rnatm* atm, 304 struct octrees_storage_desc* desc, 305 FILE* stream) 306 { 307 size_t iband; 308 size_t nbands; 309 res_T res = RES_OK; 310 ASSERT(atm && desc && stream); 311 312 *desc = OCTREES_STORAGE_DESC_NULL; 313 314 desc->allocator = atm->allocator; 315 316 #define READ(Var, N) { \ 317 if(fread((Var), sizeof(*(Var)), (N), stream) != (N)) { \ 318 if(feof(stream)) { \ 319 res = RES_BAD_ARG; \ 320 } else if(ferror(stream)) { \ 321 res = RES_IO_ERR; \ 322 } else { \ 323 res = RES_UNKNOWN_ERR; \ 324 } \ 325 goto error; \ 326 } \ 327 } (void)0 328 READ(&desc->version, 1); 329 if(desc->version != OCTREES_STORAGE_VERSION) { /* Basic versioning */ 330 res = RES_BAD_ARG; 331 goto error; 332 } 333 READ(desc->ibands, 2); 334 READ(&desc->optical_thickness, 1); 335 READ(desc->grid_definition, 3); 336 READ(desc->mesh_signature, sizeof(hash256_T)); 337 338 /* Allocate the list of band signatures */ 339 nbands = desc->ibands[1] - desc->ibands[0] + 1; 340 desc->band_signatures = MEM_CALLOC(desc->allocator, nbands, sizeof(hash256_T)); 341 if(!desc->band_signatures) { 342 log_err(atm, "error allocating signatures per band\n"); 343 res = RES_MEM_ERR; 344 goto error; 345 } 346 347 FOR_EACH(iband, 0, nbands) { 348 READ(desc->band_signatures[iband], sizeof(hash256_T)); 349 } 350 #undef READ 351 352 exit: 353 return res; 354 error: 355 log_err(atm, "unable to read tree storage descriptor -- %s\n", 356 res_to_cstr(res)); 357 goto exit; 358 } 359 360 void 361 octrees_storage_desc_release(struct octrees_storage_desc* desc) 362 { 363 ASSERT(desc); 364 if(desc->band_signatures) MEM_RM(desc->allocator, desc->band_signatures); 365 } 366 367 res_T 368 write_octrees_storage_desc 369 (const struct rnatm* atm, 370 const struct octrees_storage_desc* desc, 371 FILE* stream) 372 { 373 size_t iband; 374 size_t nbands; 375 res_T res = RES_OK; 376 ASSERT(atm && desc && stream); 377 378 nbands = desc->ibands[1] - desc->ibands[0] + 1; 379 380 #define WRITE(Var, N) { \ 381 if(fwrite((Var), sizeof(*(Var)), (N), stream) != (N)) { \ 382 res = RES_IO_ERR; \ 383 goto error; \ 384 } \ 385 } (void)0 386 WRITE(&desc->version, 1); 387 WRITE(desc->ibands, 2); 388 WRITE(&desc->optical_thickness, 1); 389 WRITE(desc->grid_definition, 3); 390 WRITE(desc->mesh_signature, sizeof(hash256_T)); 391 FOR_EACH(iband, 0, nbands) { 392 WRITE(desc->band_signatures[iband], sizeof(hash256_T)); 393 } 394 #undef WRITE 395 396 exit: 397 return res; 398 error: 399 log_err(atm, "unable to write tree storage descriptor\n"); 400 goto exit; 401 } 402 403 res_T 404 check_octrees_storage_compatibility 405 (const struct octrees_storage_desc* reference, 406 const struct octrees_storage_desc* challenge) 407 { 408 size_t iband_challenge; 409 ASSERT(reference && challenge); 410 411 /* Check version equality */ 412 if(reference->version != challenge->version) 413 return RES_BAD_ARG; 414 415 /* Verify that the bands in the repository to be challenged are included in 416 * the bands in the reference repository */ 417 if(reference->ibands[0] > challenge->ibands[0] 418 || reference->ibands[1] < challenge->ibands[1]) 419 return RES_BAD_ARG; 420 421 /* Check the equality of optical thicknesses */ 422 if(reference->optical_thickness != challenge->optical_thickness) 423 return RES_BAD_ARG; 424 425 /* Check the equality of grid definitions */ 426 if(reference->grid_definition[0] != challenge->grid_definition[0] 427 || reference->grid_definition[1] != challenge->grid_definition[1] 428 || reference->grid_definition[2] != challenge->grid_definition[2]) 429 return RES_BAD_ARG; 430 431 /* Check the equality of meshes */ 432 if(!hash256_eq(reference->mesh_signature, challenge->mesh_signature)) 433 return RES_BAD_ARG; 434 435 /* Check the equality of the per band data */ 436 FOR_EACH(iband_challenge, challenge->ibands[0], challenge->ibands[1]+1) { 437 const size_t ichallenge = iband_challenge - challenge->ibands[0]; 438 const size_t ireference = iband_challenge - reference->ibands[0]; 439 hash256_T* band_challenge = challenge->band_signatures+ichallenge; 440 hash256_T* band_reference = reference->band_signatures+ireference; 441 if(!hash256_eq(*band_challenge, *band_reference)) 442 return RES_BAD_ARG; 443 } 444 445 return RES_OK; 446 } 447 448 res_T 449 reserve_octrees_storage_toc(const struct rnatm* atm, FILE* stream) 450 { 451 size_t noctrees = 0; 452 size_t sizeof_toc = 0; 453 int err; 454 res_T res = RES_OK; 455 ASSERT(atm && stream); 456 457 noctrees = darray_accel_struct_size_get(&atm->accel_structs); 458 459 /* Compute the size of TOC */ 460 sizeof_toc = 461 sizeof(size_t)/*#octrees*/ 462 + noctrees * (sizeof(size_t)/*iband*/+sizeof(size_t)/*iquad*/+sizeof(fpos_t)); 463 ASSERT(sizeof_toc < LONG_MAX); 464 465 /* Reserve the space for the table of contents */ 466 err = fseek(stream, (long)sizeof_toc, SEEK_CUR); 467 if(err != 0) { res = RES_IO_ERR; goto error; } 468 469 exit: 470 return res; 471 error: 472 log_err(atm, "error reserving octree storage toc -- %s\n", 473 strerror(errno)); 474 goto exit; 475 } 476 477 res_T 478 write_octrees_storage_toc(const struct rnatm* atm, FILE* stream) 479 { 480 size_t noctrees = 0; 481 size_t ioctree = 0; 482 res_T res = RES_OK; 483 ASSERT(atm && stream); 484 485 noctrees = darray_accel_struct_size_get(&atm->accel_structs); 486 487 /* Write the number of voxelized octrees and their corresponding offset */ 488 #define WRITE(Var) { \ 489 if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \ 490 res = RES_IO_ERR; \ 491 goto error; \ 492 } \ 493 } (void)0 494 WRITE(&noctrees); 495 FOR_EACH(ioctree, 0, noctrees) { 496 const struct accel_struct* accel_struct = 497 darray_accel_struct_cdata_get(&atm->accel_structs) + ioctree; 498 WRITE(&accel_struct->iband); 499 WRITE(&accel_struct->iquad_pt); 500 WRITE(&accel_struct->fpos); 501 } 502 #undef WRITE 503 504 exit: 505 return res; 506 error: 507 log_err(atm, "error writing octree storage table of content -- %s\n", 508 res_to_cstr(res)); 509 goto exit; 510 } 511 512 res_T 513 read_octrees_storage_toc(struct rnatm* atm, FILE* stream) 514 { 515 struct accel_struct* accel_structs = NULL; 516 size_t ioctree = 0; 517 size_t noctrees = 0; 518 size_t ioctree_stream = 0; 519 size_t noctrees_stream = 0; 520 size_t iband = 0; 521 size_t iquad = 0; 522 res_T res = RES_OK; 523 ASSERT(atm && stream); 524 525 noctrees = darray_accel_struct_size_get(&atm->accel_structs); 526 527 #define READ(Var) { \ 528 if(fread((Var), sizeof(*(Var)), (1), stream) != (1)) { \ 529 if(feof(stream)) { \ 530 res = RES_BAD_ARG; \ 531 } else if(ferror(stream)) { \ 532 res = RES_IO_ERR; \ 533 } else { \ 534 res = RES_UNKNOWN_ERR; \ 535 } \ 536 goto error; \ 537 } \ 538 } (void)0 539 540 accel_structs = darray_accel_struct_data_get(&atm->accel_structs); 541 noctrees = darray_accel_struct_size_get(&atm->accel_structs); 542 543 /* Read the number of octrees stored in the stream */ 544 READ(&noctrees_stream); 545 if(noctrees_stream < noctrees) { 546 res = RES_BAD_ARG; 547 goto error; 548 } 549 550 /* Look for the first octree to load for the current atmosphere */ 551 FOR_EACH(ioctree_stream, 0, noctrees_stream) { 552 READ(&iband); 553 READ(&iquad); 554 READ(&accel_structs[0].fpos); 555 556 if(accel_structs[0].iband == iband 557 && accel_structs[0].iquad_pt == iquad) 558 break; 559 } 560 561 /* Cannot find the first octree */ 562 if(ioctree_stream >= noctrees_stream) { 563 res = RES_BAD_ARG; 564 goto error; 565 } 566 567 /* Load the remaining offsets */ 568 FOR_EACH(ioctree, 1, noctrees) { 569 570 /* No sufficient octrees in the stream */ 571 if(++ioctree_stream >= noctrees_stream) { 572 res = RES_BAD_ARG; 573 goto error; 574 } 575 576 READ(&iband); 577 READ(&iquad); 578 READ(&accel_structs[ioctree].fpos); 579 580 /* Unexpected octrees */ 581 if(accel_structs[ioctree].iband != iband 582 || accel_structs[ioctree].iquad_pt != iquad) { 583 res = RES_BAD_ARG; 584 goto error; 585 } 586 } 587 #undef READ 588 589 exit: 590 return res; 591 error: 592 log_err(atm, "error reading octree storage table of content -- %s\n", 593 res_to_cstr(res)); 594 goto exit; 595 } 596 597 extern LOCAL_SYM res_T 598 store_octrees(struct rnatm* atm, const size_t ioctrees[2], FILE* stream) 599 { 600 size_t ioctree; 601 int err; 602 res_T res = RES_OK; 603 ASSERT(atm && ioctrees && stream && ioctrees[0] <= ioctrees[1]); 604 605 FOR_EACH(ioctree, ioctrees[0], ioctrees[1]+1) { 606 struct accel_struct* accel_struct = NULL; 607 accel_struct = darray_accel_struct_data_get(&atm->accel_structs) + ioctree; 608 609 /* Save the current file offset */ 610 err = fgetpos(stream, &accel_struct->fpos); 611 if(err != 0) { 612 log_err(atm, "error retrieving octree storage position -- %s\n", 613 strerror(errno)); 614 res = RES_IO_ERR; 615 goto error; 616 } 617 618 /* Serialize the octree */ 619 res = svx_tree_write(accel_struct->octree, stream); 620 if(res != RES_OK) { 621 log_err(atm, "error writing octree %lu -- %s\n", 622 (unsigned long)ioctree, res_to_cstr(res)); 623 goto error; 624 } 625 } 626 627 exit: 628 return res; 629 error: 630 goto exit; 631 }