htsky_dump_cloud_vtk.c (9095B)
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 #include "htsky.h" 19 #include "htsky_c.h" 20 #include "htsky_log.h" 21 #include "htsky_cloud.h" 22 23 #include <high_tune/htgop.h> 24 25 #include <rsys/dynamic_array_double.h> 26 #include <rsys/dynamic_array_size_t.h> 27 #include <rsys/hash_table.h> 28 29 #include <star/svx.h> 30 31 struct vertex { 32 double x; 33 double y; 34 double z; 35 }; 36 37 static char 38 vertex_eq(const struct vertex* v0, const struct vertex* v1) 39 { 40 return eq_eps(v0->x, v1->x, 1.e-6) 41 && eq_eps(v0->y, v1->y, 1.e-6) 42 && eq_eps(v0->z, v1->z, 1.e-6); 43 } 44 45 /* Generate the htable_vertex data structure */ 46 #define HTABLE_NAME vertex 47 #define HTABLE_KEY struct vertex 48 #define HTABLE_DATA size_t 49 #define HTABLE_KEY_FUNCTOR_EQ vertex_eq 50 #include <rsys/hash_table.h> 51 52 /* Temporary data structure used to dump the octree data into a VTK file */ 53 struct octree_data { 54 struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */ 55 struct darray_double vertices; /* Array of unique vertices */ 56 struct darray_double data; /* List of registered leaf data */ 57 struct darray_size_t cells; /* Ids of the cell vertices */ 58 size_t iband; /* Index of the band that overlaps the CIE XYZ color space */ 59 size_t iquad; /* Index of the quadrature point into the band */ 60 const struct htsky* sky; 61 }; 62 63 /******************************************************************************* 64 * Helper functions 65 ******************************************************************************/ 66 static INLINE void 67 octree_data_init 68 (const struct htsky* sky, 69 const size_t iband, 70 const size_t iquad, 71 struct octree_data* data) 72 { 73 ASSERT(data); 74 ASSERT(iband >= sky->bands_range[0]); 75 ASSERT(iband <= sky->bands_range[1]); 76 (void)iquad; 77 htable_vertex_init(sky->allocator, &data->vertex2id); 78 darray_double_init(sky->allocator, &data->vertices); 79 darray_double_init(sky->allocator, &data->data); 80 darray_size_t_init(sky->allocator, &data->cells); 81 data->sky = sky; 82 data->iband = iband; 83 data->iquad = iquad; 84 } 85 86 static INLINE void 87 octree_data_release(struct octree_data* data) 88 { 89 ASSERT(data); 90 htable_vertex_release(&data->vertex2id); 91 darray_double_release(&data->vertices); 92 darray_double_release(&data->data); 93 darray_size_t_release(&data->cells); 94 } 95 96 static INLINE void 97 register_leaf 98 (const struct svx_voxel* leaf, 99 const size_t ileaf, 100 void* context) 101 { 102 struct octree_data* ctx = context; 103 struct vertex v[8]; 104 double kext_min; 105 double kext_max; 106 int i; 107 ASSERT(leaf && ctx); 108 (void)ileaf; 109 110 /* Compute the leaf vertices */ 111 v[0].x = leaf->lower[0]; v[0].y = leaf->lower[1]; v[0].z = leaf->lower[2]; 112 v[1].x = leaf->upper[0]; v[1].y = leaf->lower[1]; v[1].z = leaf->lower[2]; 113 v[2].x = leaf->lower[0]; v[2].y = leaf->upper[1]; v[2].z = leaf->lower[2]; 114 v[3].x = leaf->upper[0]; v[3].y = leaf->upper[1]; v[3].z = leaf->lower[2]; 115 v[4].x = leaf->lower[0]; v[4].y = leaf->lower[1]; v[4].z = leaf->upper[2]; 116 v[5].x = leaf->upper[0]; v[5].y = leaf->lower[1]; v[5].z = leaf->upper[2]; 117 v[6].x = leaf->lower[0]; v[6].y = leaf->upper[1]; v[6].z = leaf->upper[2]; 118 v[7].x = leaf->upper[0]; v[7].y = leaf->upper[1]; v[7].z = leaf->upper[2]; 119 120 FOR_EACH(i, 0, 8) { 121 size_t *pid = htable_vertex_find(&ctx->vertex2id, v+i); 122 size_t id; 123 if(pid) { 124 id = *pid; 125 } else { /* Register the leaf vertex */ 126 id = darray_double_size_get(&ctx->vertices)/3; 127 CHK(RES_OK == htable_vertex_set(&ctx->vertex2id, v+i, &id)); 128 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].x)); 129 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].y)); 130 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].z)); 131 } 132 /* Add the vertex id to the leaf cell */ 133 CHK(RES_OK == darray_size_t_push_back(&ctx->cells, &id)); 134 } 135 136 /* Register the leaf data */ 137 kext_max = htsky_fetch_svx_voxel_property(ctx->sky, HTSKY_Kext, 138 HTSKY_SVX_MAX, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, leaf); 139 kext_min = htsky_fetch_svx_voxel_property(ctx->sky, HTSKY_Kext, 140 HTSKY_SVX_MIN, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, leaf); 141 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min)); 142 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max)); 143 } 144 145 /******************************************************************************* 146 * Exported functions 147 ******************************************************************************/ 148 res_T 149 htsky_dump_cloud_vtk 150 (const struct htsky* sky, 151 const size_t iband, /* Index of the spectral band */ 152 const size_t iquad, /* Index of the quadrature point */ 153 FILE* stream) 154 { 155 const struct cloud* cloud; 156 struct htgop_spectral_interval specint; 157 struct octree_data data; 158 const double* leaf_data; 159 size_t nvertices; 160 size_t ncells; 161 size_t i; 162 res_T res = RES_OK; 163 164 if(!sky || !stream) { 165 res = RES_BAD_ARG; 166 goto error; 167 } 168 169 if(iband < sky->bands_range[0] || iband > sky->bands_range[1]) { 170 log_err(sky, 171 "%s: invalid spectral band index `%lu'. " 172 "Valid short wave spectral bands lie in [%lu, %lu]\n", 173 FUNC_NAME, 174 (unsigned long)iband, 175 (unsigned long)sky->bands_range[0], 176 (unsigned long)sky->bands_range[1]); 177 res = RES_BAD_ARG; 178 goto error; 179 } 180 181 if(iquad >= htsky_get_spectral_band_quadrature_length(sky, iband)) { 182 log_err(sky, 183 "invalid quadrature point `%lu' for the spectral band `%lu'.\n", 184 (unsigned long)iquad, (unsigned long)iband); 185 res = RES_BAD_ARG; 186 goto error; 187 } 188 189 if(!sky->is_cloudy) { 190 log_warn(sky, "%s: the sky has no cloud.\n", FUNC_NAME); 191 return RES_OK; 192 } 193 194 i = iband - sky->bands_range[0]; 195 196 octree_data_init(sky, iband, iquad, &data); 197 cloud = &sky->clouds[i][iquad]; 198 199 ASSERT(cloud->octree_desc.type == SVX_OCTREE); 200 201 /* Register leaf data */ 202 SVX(tree_for_each_leaf(cloud->octree, register_leaf, &data)); 203 nvertices = darray_double_size_get(&data.vertices) / 3/*#coords per vertex*/; 204 ncells = darray_size_t_size_get(&data.cells)/8/*#ids per cell*/; 205 ASSERT(ncells == cloud->octree_desc.nleaves); 206 207 /* Fetch the spectral interval descriptor */ 208 switch(sky->spectral_type) { 209 case HTSKY_SPECTRAL_LW: 210 HTGOP(get_lw_spectral_interval(sky->htgop, iband, &specint)); 211 break; 212 case HTSKY_SPECTRAL_SW: 213 HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint)); 214 break; 215 default: FATAL("Unreachable code.\n"); break; 216 } 217 218 /* Write headers */ 219 fprintf(stream, "# vtk DataFile Version 2.0\n"); 220 fprintf(stream, "Clouds optical properties in [%g, %g] nanometers\n", 221 wavenumber_to_wavelength(specint.wave_numbers[1]), 222 wavenumber_to_wavelength(specint.wave_numbers[0])); 223 fprintf(stream, "ASCII\n"); 224 fprintf(stream, "DATASET UNSTRUCTURED_GRID\n"); 225 226 /* Write vertex coordinates */ 227 fprintf(stream, "POINTS %lu float\n", (unsigned long)nvertices); 228 FOR_EACH(i, 0, nvertices) { 229 fprintf(stream, "%g %g %g\n", 230 SPLIT3(darray_double_cdata_get(&data.vertices) + i*3)); 231 } 232 233 /* Write the cells */ 234 fprintf(stream, "CELLS %lu %lu\n", 235 (unsigned long)ncells, 236 (unsigned long)(ncells*(8/*#verts per cell*/ + 1/*1st field of a cell*/))); 237 FOR_EACH(i, 0, ncells) { 238 fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n", 239 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+0], 240 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+1], 241 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+2], 242 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+3], 243 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+4], 244 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+5], 245 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+6], 246 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+7]); 247 } 248 249 /* Write the cell type */ 250 fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)ncells); 251 FOR_EACH(i, 0, ncells) fprintf(stream, "11\n"); 252 253 /* Write the cell data */ 254 leaf_data = darray_double_cdata_get(&data.data); 255 fprintf(stream, "CELL_DATA %lu\n", (unsigned long)ncells); 256 fprintf(stream, "SCALARS Kext double 2\n"); 257 fprintf(stream, "LOOKUP_TABLE default\n"); 258 FOR_EACH(i, 0, ncells) { 259 fprintf(stream, "%g %g\n", leaf_data[i*2+0], leaf_data[i*2+1]); 260 } 261 octree_data_release(&data); 262 263 exit: 264 return res; 265 error: 266 goto exit; 267 } 268