htsky

Load and structure a vertically stratified atmosphere
git clone git://git.meso-star.fr/htsky.git
Log | Files | Refs | README | LICENSE

htsky_svx.h (5000B)


      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 #ifndef HTSKY_SVX_H
     19 #define HTSKY_SVX_H
     20 
     21 #include "htsky.h"
     22 
     23 #include <rsys/math.h>
     24 
     25 /*
     26  * SVX Memory layout
     27  * -----------------
     28  *
     29  * For each SVX voxel, the data of the optical property are stored
     30  * linearly as N single precision floating point data, with N computed as
     31  * bellow:
     32  *
     33  *  N = HTSKY_PROPS_COUNT__   #optical properties per voxel
     34  *    * HTSKY_SVX_OPS_COUNT__ #supported operations on each properties
     35  *    * HTSKY_CPNTS_COUNT__;  #components on which properties are defined
     36  *
     37  * In a given voxel, the index `id' in [0, N-1] corresponding to the optical
     38  * property `enum htrdr_sky_property P' of the component `enum
     39  * htrdr_sky_component C' according to the operation `enum htsky_svx_op O' is
     40  * then computed as bellow:
     41  *
     42  *  id = C * NFLOATS_PER_CPNT + P * HTSKY_SVX_OPS_COUNT__ + O;
     43  *  NFLOATS_PER_CPNT = HTSKY_SVX_OPS_COUNT__ * HTSKY_PROPS_COUNT__;
     44  */
     45 
     46 /* Constant defining the number of floating point data per component */
     47 #define NFLOATS_PER_CPNT (HTSKY_SVX_OPS_COUNT__ * HTSKY_PROPS_COUNT__)
     48 
     49 /* Constant defining the overall number of floating point data of a voxel */
     50 #define NFLOATS_PER_VOXEL (NFLOATS_PER_CPNT * HTSKY_CPNTS_COUNT__)
     51 
     52 /* Context used to build the SVX hierarchical data structures */
     53 struct build_tree_context {
     54   const struct htsky* sky;
     55   double vxsz[3];
     56   double tau_threshold; /* Threshold criteria for the merge process */
     57   size_t iband; /* Index of the band that overlaps the CIE XYZ color space */
     58   size_t quadrature_range[2]; /* Range of quadrature point indices to handle */
     59 };
     60 static const struct build_tree_context BUILD_TREE_CONTEXT_NULL;
     61 
     62 static FINLINE float
     63 vox_get
     64   (const float* data,
     65    const enum htsky_component cpnt,
     66    const enum htsky_property prop,
     67    const enum htsky_svx_op op)
     68 {
     69   ASSERT(data);
     70   return data[cpnt*NFLOATS_PER_CPNT+ prop*HTSKY_SVX_OPS_COUNT__ + op];
     71 }
     72 
     73 static FINLINE void
     74 vox_set
     75   (float* data,
     76    const enum htsky_component cpnt,
     77    const enum htsky_property prop,
     78    const enum htsky_svx_op op,
     79    const float val)
     80 {
     81   ASSERT(data);
     82   data[cpnt*NFLOATS_PER_CPNT+ prop*HTSKY_SVX_OPS_COUNT__ + op] = val;
     83 }
     84 
     85 static INLINE void
     86 vox_merge_component
     87   (float* vox_out,
     88    const enum htsky_component cpnt,
     89    const float* voxs[],
     90    const size_t nvoxs)
     91 {
     92   float ka_min = FLT_MAX;
     93   float ka_max =-FLT_MAX;
     94   float ks_min = FLT_MAX;
     95   float ks_max =-FLT_MAX;
     96   float kext_min = FLT_MAX;
     97   float kext_max =-FLT_MAX;
     98   size_t ivox;
     99   ASSERT(vox_out && voxs && nvoxs);
    100 
    101   FOR_EACH(ivox, 0, nvoxs) {
    102     const float* vox = voxs[ivox];
    103     ka_min = MMIN(ka_min, vox_get(vox, cpnt, HTSKY_Ka, HTSKY_SVX_MIN));
    104     ka_max = MMAX(ka_max, vox_get(vox, cpnt, HTSKY_Ka, HTSKY_SVX_MAX));
    105     ks_min = MMIN(ks_min, vox_get(vox, cpnt, HTSKY_Ks, HTSKY_SVX_MIN));
    106     ks_max = MMAX(ks_max, vox_get(vox, cpnt, HTSKY_Ks, HTSKY_SVX_MAX));
    107     kext_min = MMIN(kext_min, vox_get(vox, cpnt, HTSKY_Kext, HTSKY_SVX_MIN));
    108     kext_max = MMAX(kext_max, vox_get(vox, cpnt, HTSKY_Kext, HTSKY_SVX_MAX));
    109   }
    110 
    111   vox_set(vox_out, cpnt, HTSKY_Ka, HTSKY_SVX_MIN, ka_min);
    112   vox_set(vox_out, cpnt, HTSKY_Ka, HTSKY_SVX_MAX, ka_max);
    113   vox_set(vox_out, cpnt, HTSKY_Ks, HTSKY_SVX_MIN, ks_min);
    114   vox_set(vox_out, cpnt, HTSKY_Ks, HTSKY_SVX_MAX, ks_max);
    115   vox_set(vox_out, cpnt, HTSKY_Kext, HTSKY_SVX_MIN, kext_min);
    116   vox_set(vox_out, cpnt, HTSKY_Kext, HTSKY_SVX_MAX, kext_max);
    117 }
    118 
    119 static INLINE int
    120 vox_challenge_merge_component
    121   (const enum htsky_component comp,
    122    const struct svx_voxel voxels[],
    123    const size_t nvoxs,
    124    struct build_tree_context* ctx)
    125 {
    126   double lower_z = DBL_MAX;
    127   double upper_z =-DBL_MAX;
    128   double dst;
    129   float kext_min = FLT_MAX;
    130   float kext_max =-FLT_MAX;
    131   size_t ivox;
    132   ASSERT(voxels && nvoxs && ctx);
    133 
    134   FOR_EACH(ivox, 0, nvoxs) {
    135     const float* vox = voxels[ivox].data;
    136     kext_min = MMIN(kext_min, vox_get(vox, comp, HTSKY_Kext, HTSKY_SVX_MIN));
    137     kext_max = MMAX(kext_max, vox_get(vox, comp, HTSKY_Kext, HTSKY_SVX_MAX));
    138     lower_z = MMIN(voxels[ivox].lower[2], lower_z);
    139     upper_z = MMAX(voxels[ivox].upper[2], upper_z);
    140   }
    141   dst = upper_z - lower_z;
    142   return (kext_max - kext_min)*dst <= ctx->tau_threshold;
    143 }
    144 
    145 #endif /* HTSKY_SVX_H */
    146