test_suvm_primitive_intersection.c (9593B)
1 /* Copyright (C) 2020-2023 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "suvm.h" 17 #include "test_suvm_ball.h" 18 #include "test_suvm_box.h" 19 #include "test_suvm_utils.h" 20 21 #include <rsys/float3.h> 22 #include <string.h> 23 24 struct mesh { 25 const double* vertices; /* List of double[3] */ 26 size_t nvertices; 27 const size_t* tetrahedra; /* List of size_t[4] */ 28 size_t ntetrahedra; 29 }; 30 static const struct mesh MESH_NULL; 31 32 /******************************************************************************* 33 * Helper functions 34 *************v*****************************************************************/ 35 static void 36 get_indices(const size_t itetra, size_t ids[4], void* ctx) 37 { 38 struct mesh* msh = ctx; 39 CHK(itetra < msh->ntetrahedra); 40 CHK(ids != NULL); 41 ids[0] = msh->tetrahedra[itetra*4+0]; 42 ids[1] = msh->tetrahedra[itetra*4+1]; 43 ids[2] = msh->tetrahedra[itetra*4+2]; 44 ids[3] = msh->tetrahedra[itetra*4+3]; 45 } 46 47 static void 48 get_position(const size_t ivert, double pos[3], void* ctx) 49 { 50 struct mesh* msh = ctx; 51 CHK(ctx != NULL); 52 CHK(pos != NULL); 53 CHK(ivert < msh->nvertices); 54 pos[0] = msh->vertices[ivert*3+0]; 55 pos[1] = msh->vertices[ivert*3+1]; 56 pos[2] = msh->vertices[ivert*3+2]; 57 } 58 59 static void 60 write_voxels 61 (FILE* stream, 62 const int* vxls, 63 const size_t def[3], 64 const double vxl_sz[3]) 65 { 66 size_t ix, iy, iz; 67 size_t ivxl; 68 size_t i; 69 CHK(stream && vxls && def && def[0] && def[1] && def[2]); 70 71 fprintf(stream, "# vtk DataFile Version 2.0\n"); 72 fprintf(stream, "nothing\n"); 73 fprintf(stream, "ASCII\n"); 74 75 fprintf(stream, "DATASET RECTILINEAR_GRID\n"); 76 fprintf(stream, "DIMENSIONS %lu %lu %lu\n", def[0]+1, def[1]+1, def[2]+1); 77 fprintf(stream, "X_COORDINATES %lu float\n", def[0]+1); 78 FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0]); 79 fprintf(stream, "\n"); 80 fprintf(stream, "Y_COORDINATES %lu float\n", def[1]+1); 81 FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1]); 82 fprintf(stream, "\n"); 83 fprintf(stream, "Z_COORDINATES %lu float\n", def[2]+1); 84 FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2]); 85 fprintf(stream, "\n"); 86 87 fprintf(stream, "CELL_DATA %lu\n", def[0]*def[1]*def[2]); 88 fprintf(stream, "SCALARS intersect_type int 1\n"); 89 fprintf(stream, "LOOKUP_TABLE default\n"); 90 91 ivxl = 0; 92 FOR_EACH(iz, 0, def[2]) { 93 FOR_EACH(iy, 0, def[1]) { 94 FOR_EACH(ix, 0, def[0]) { 95 fprintf(stream, "%i\n", vxls[ivxl]); 96 ++ivxl; 97 } 98 } 99 } 100 } 101 102 static void 103 voxelise_volume 104 (const struct suvm_volume* vol, 105 int* vxls, 106 const size_t def[3]) 107 { 108 double low[3]; 109 double upp[3]; 110 double vxl_sz[3]; 111 size_t iprim; 112 size_t nprims; 113 114 CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK); 115 CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK); 116 117 memset(vxls, 0, sizeof(*vxls)*def[0]*def[1]*def[2]); 118 119 vxl_sz[0] = (upp[0] - low[0]) / (double)def[0]; 120 vxl_sz[1] = (upp[1] - low[1]) / (double)def[1]; 121 vxl_sz[2] = (upp[2] - low[2]) / (double)def[2]; 122 123 FOR_EACH(iprim, 0, nprims) { 124 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 125 struct suvm_polyhedron poly; 126 size_t ivxl_low[3]; 127 size_t ivxl_upp[3]; 128 size_t ivxl[3]; 129 size_t i = 0; 130 131 CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK); 132 CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK); 133 134 /* Transform the polyhedron AABB in voxel space */ 135 ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]); 136 ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]); 137 ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]); 138 ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]); 139 ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]); 140 ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]); 141 CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]); 142 CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]); 143 CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]); 144 145 FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) { 146 float vxl_low[3]; 147 float vxl_upp[3]; 148 vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]); 149 vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2]; 150 FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) { 151 vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]); 152 vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1]; 153 FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) { 154 vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]); 155 vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0]; 156 157 i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1]; 158 vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp); 159 } 160 } 161 } 162 } 163 } 164 165 static void 166 test_mesh_voxelization 167 (struct suvm_device* dev, 168 struct mem_allocator* allocator, 169 struct mesh* msh, 170 const size_t def[3], 171 const char* filename) /* NULL <=> do not write the result onto disk */ 172 { 173 struct suvm_volume* vol = NULL; 174 struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 175 double low[3], upp[3]; 176 double vxl_sz[3]; 177 int* vxls = NULL; 178 179 args.ntetrahedra = msh->ntetrahedra; 180 args.nvertices = msh->nvertices; 181 args.get_indices = get_indices; 182 args.get_position = get_position; 183 args.context = msh; 184 185 CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK); 186 CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK); 187 vxl_sz[0] = (upp[0] - low[0]) / (double)def[0]; 188 vxl_sz[1] = (upp[1] - low[1]) / (double)def[1]; 189 vxl_sz[2] = (upp[2] - low[2]) / (double)def[2]; 190 191 CHK(allocator != NULL); 192 vxls = MEM_CALLOC(allocator, def[0]*def[1]*def[2], sizeof(*vxls)); 193 CHK(vxls != NULL); 194 195 voxelise_volume(vol, vxls, def); 196 197 if(filename) { 198 FILE* fp = fopen(filename, "w"); 199 CHK(fp != NULL); 200 write_voxels(fp, vxls, def, vxl_sz); 201 CHK(fclose(fp) == 0); 202 } 203 204 MEM_RM(allocator, vxls); 205 CHK(suvm_volume_ref_put(vol) == RES_OK); 206 } 207 208 209 static void 210 test_tetra_aabb_intersection 211 (struct suvm_device* dev, 212 struct mem_allocator* allocator) 213 { 214 const double vertices[] = { 215 0.0, 0.0, 0.0, 216 0.0, 0.0, 1.0, 217 1.0, 0.0, 1.0, 218 0.0, 1.0, 1.0 219 }; 220 const size_t tetra[] = { 0, 1, 2, 3 }; 221 struct mesh msh = MESH_NULL; 222 struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 223 struct suvm_polyhedron poly; 224 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 225 struct suvm_volume* vol= NULL; 226 float low[3]; 227 float upp[3]; 228 double vxl_sz[3]; 229 const size_t def[3] = {16, 16, 16}; 230 size_t nprims; 231 int* vxls = NULL; 232 (void)vxl_sz; 233 234 msh.vertices = vertices; 235 msh.nvertices = 4; 236 msh.tetrahedra = tetra; 237 msh.ntetrahedra = 1; 238 239 args.ntetrahedra = 1; 240 args.nvertices = 4; 241 args.get_indices = get_indices; 242 args.get_position = get_position; 243 args.context = &msh; 244 245 CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK); 246 247 CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK); 248 CHK(nprims == 1); 249 250 CHK(suvm_volume_get_primitive(vol, 0, &prim) == RES_OK); 251 CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK); 252 253 f3(low, 0.f, 0.f, 0.f); 254 f3(upp, 1.f, 1.f, 1.f); 255 CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_IS_INCLUDED); 256 257 f3(low, 0.f, 0.f, 0.5f); 258 f3(upp, 0.5f, 0.5f, 1.f); 259 CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_PARTIAL); 260 261 f3(low, 0.f, 0.f, 0.66667f /* ~1-1/3 */); 262 f3(upp, 0.33332f/*~1/3*/, 0.33332f/*~1-1/3*/, 1.f); 263 CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_INCLUDE); 264 265 f3(low, 0.33334f/*~1.3*/, 0.33334f/* ~1/3 */, 0); 266 f3(upp, 1.f, 1.f, 0.66665f /*~ 1-1/3 */); 267 CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_NONE); 268 269 CHK(allocator != NULL); 270 vxls = MEM_CALLOC(allocator, def[0]*def[1]*def[2], sizeof(*vxls)); 271 CHK(vxls != NULL); 272 273 voxelise_volume(vol, vxls, def); 274 275 MEM_RM(allocator, vxls); 276 CHK(suvm_volume_ref_put(vol) == RES_OK); 277 } 278 279 /******************************************************************************* 280 * Main function 281 ******************************************************************************/ 282 int 283 main(int argc, char** argv) 284 { 285 struct mesh msh = MESH_NULL; 286 struct suvm_device* dev = NULL; 287 size_t def[3] = {64, 64, 64}; 288 (void)argc, (void)argv; 289 290 CHK(suvm_device_create(NULL, &mem_default_allocator, 1, &dev) == RES_OK); 291 292 test_tetra_aabb_intersection(dev, &mem_default_allocator); 293 294 msh.vertices = box_vertices; 295 msh.nvertices = box_nverts; 296 msh.tetrahedra = box_indices; 297 msh.ntetrahedra = box_ntetras; 298 test_mesh_voxelization(dev, &mem_default_allocator, &msh, def, "box.vtk"); 299 300 msh.vertices = ball_vertices; 301 msh.nvertices = ball_nvertices; 302 msh.tetrahedra = ball_tetrahedra; 303 msh.ntetrahedra = ball_ntetrahedra; 304 test_mesh_voxelization(dev, &mem_default_allocator, &msh, def, "ball.vtk"); 305 306 CHK(suvm_device_ref_put(dev) == RES_OK); 307 check_memory_allocator(&mem_default_allocator); 308 CHK(mem_allocated_size() == 0); 309 return 0; 310 }