test_svx_octree_trace_ray.c (14217B)
1 /* Copyright (C) 2018, 2020-2025 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018 Université Paul Sabatier 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* nextafter function */ 18 19 #include "svx.h" 20 #include "test_svx_utils.h" 21 22 #include <math.h> 23 24 #include <rsys/clock_time.h> 25 #include <rsys/double2.h> 26 #include <rsys/double3.h> 27 #include <rsys/image.h> 28 #include <rsys/math.h> 29 #include <rsys/morton.h> 30 31 struct scene { 32 double origin[3]; 33 double vxsz[3]; 34 double sphere_pos[3]; 35 double sphere_radius; 36 }; 37 38 static char 39 sphere_intersect_aabb 40 (const double pos[3], /* Sphere position */ 41 const double radius, /* Sphere radius */ 42 const double low[3], /* AABB lower bound */ 43 const double upp[3]) /* AABB upper bound */ 44 { 45 double v[3]; 46 double sqr_dst; 47 48 CHK(pos != NULL); 49 CHK(low != NULL); 50 CHK(upp != NULL); 51 CHK(radius > 0); 52 53 v[0] = pos[0]<low[0] ? low[0]-pos[0] : (pos[0]>upp[0] ? pos[0]-upp[0] : 0); 54 v[1] = pos[1]<low[1] ? low[1]-pos[1] : (pos[1]>upp[1] ? pos[1]-upp[1] : 0); 55 v[2] = pos[2]<low[2] ? low[2]-pos[2] : (pos[2]>upp[2] ? pos[2]-upp[2] : 0); 56 57 sqr_dst = v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; 58 59 return sqr_dst <= radius*radius; 60 } 61 62 static void 63 voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* ctx) 64 { 65 const struct scene* scn = ctx; 66 char* val = dst; 67 double low[3]; 68 double upp[3]; 69 uint32_t ui3[3]; 70 CHK(xyz != NULL); 71 CHK(dst != NULL); 72 CHK(ctx != NULL); 73 74 /* Compute the AABB of the voxel */ 75 low[0] = (double)xyz[0] * scn->vxsz[0] + scn->origin[0]; 76 low[1] = (double)xyz[1] * scn->vxsz[1] + scn->origin[1]; 77 low[2] = (double)xyz[2] * scn->vxsz[2] + scn->origin[2]; 78 upp[0] = low[0] + scn->vxsz[0]; 79 upp[1] = low[1] + scn->vxsz[1]; 80 upp[2] = low[2] + scn->vxsz[2]; 81 82 ui3[0] = (uint32_t)xyz[0]; 83 ui3[1] = (uint32_t)xyz[1]; 84 ui3[2] = (uint32_t)xyz[2]; 85 CHK(mcode == morton_xyz_encode_u21(ui3)); 86 87 /* Binary octree, i.e. it stores if the voxel intersect the sphere or not */ 88 *val = sphere_intersect_aabb(scn->sphere_pos, scn->sphere_radius, low, upp); 89 } 90 91 static void 92 voxels_merge(void* dst, const void* voxels[], const size_t nvoxels, void* ctx) 93 { 94 size_t ivoxel; 95 char tmp = 0; 96 char* val = dst; 97 98 CHK(dst && voxels && nvoxels && ctx); 99 100 for(ivoxel=0; !tmp && ivoxel<nvoxels; ++ivoxel) { 101 const char* voxel_data = voxels[ivoxel]; 102 tmp = *voxel_data; 103 } 104 *val = tmp; 105 } 106 107 static int 108 voxels_challenge_merge 109 (const struct svx_voxel voxels[], const size_t nvoxels, void* ctx) 110 { 111 size_t ivoxel; 112 int merge = 1; 113 char ref; 114 115 CHK(voxels && nvoxels && ctx); 116 117 ref = *(char*)(voxels[0].data); 118 119 for(ivoxel=1; merge && ivoxel<nvoxels; ++ivoxel) { 120 const char* voxel_data = voxels[ivoxel].data; 121 merge = (ref == *voxel_data); 122 } 123 return merge; 124 } 125 126 static INLINE void 127 write_scalars 128 (const struct svx_voxel* leaf, 129 const size_t ileaf, 130 void* context) 131 { 132 FILE* stream = context; 133 (void)ileaf; 134 CHK(stream != NULL); 135 CHK(leaf != NULL); 136 fprintf(stream, "%d\n", *(char*)leaf->data); 137 } 138 139 static int 140 hit_filter 141 (const struct svx_hit* hit, 142 const double ray_org[3], 143 const double ray_dir[3], 144 const double ray_range[3], 145 void* context) 146 { 147 const struct ray* ray = context; 148 149 CHK(hit && ray_org && ray_dir && ray_range && context); 150 CHK(d3_eq(ray->org, ray_org)); 151 CHK(d3_eq(ray->dir, ray_dir)); 152 CHK(d2_eq(ray->range, ray_range)); 153 CHK(!SVX_HIT_NONE(hit)); 154 155 return *((char*)hit->voxel.data) == 0; 156 } 157 158 static int 159 hit_filter2 160 (const struct svx_hit* hit, 161 const double ray_org[3], 162 const double ray_dir[3], 163 const double ray_range[3], 164 void* context) 165 { 166 const struct svx_voxel* voxel = context; 167 168 CHK(hit && ray_org && ray_dir && ray_range && context); 169 CHK(!SVX_HIT_NONE(hit)); 170 return SVX_VOXEL_EQ(&hit->voxel, voxel) 171 || *((char*)hit->voxel.data) == 0; 172 } 173 174 static int 175 hit_filter3 176 (const struct svx_hit* hit, 177 const double ray_org[3], 178 const double ray_dir[3], 179 const double ray_range[3], 180 void* context) 181 { 182 int* accum = context; 183 CHK(hit && ray_org && ray_dir && ray_range && context); 184 CHK(!SVX_HIT_NONE(hit)); 185 *accum += *(char*)hit->voxel.data; 186 return 1; 187 } 188 189 static int 190 hit_challenge 191 (const struct svx_hit* hit, 192 const double ray_org[3], 193 const double ray_dir[3], 194 const double ray_range[2], 195 void* context) 196 { 197 (void)context; 198 CHK(hit && ray_org && ray_dir && ray_range); 199 CHK(!SVX_HIT_NONE(hit)); 200 return 1; 201 } 202 203 static int 204 hit_challenge_pass_through 205 (const struct svx_hit* hit, 206 const double ray_org[3], 207 const double ray_dir[3], 208 const double ray_range[2], 209 void* context) 210 { 211 const struct ray* ray = context; 212 CHK(hit && ray_org && ray_dir && ray_range); 213 CHK(!SVX_HIT_NONE(hit)); 214 CHK(d3_eq(ray->org, ray_org)); 215 CHK(d3_eq(ray->dir, ray_dir)); 216 CHK(d2_eq(ray->range, ray_range)); 217 return 0; 218 } 219 220 static int 221 hit_challenge_root 222 (const struct svx_hit* hit, 223 const double ray_org[3], 224 const double ray_dir[3], 225 const double ray_range[2], 226 void* context) 227 { 228 (void)hit, (void)ray_org, (void)ray_dir, (void)ray_range, (void)context; 229 return hit->voxel.depth == 0; 230 } 231 232 static void 233 draw_image(struct image* img, struct svx_tree* oct, const struct scene* scn) 234 { 235 char buf[32]; 236 struct time t0, t1; 237 struct camera cam; 238 unsigned char* pixels = NULL; 239 const size_t width = 512; 240 const size_t height = 512; 241 const double pos[3] = { 0,-1.5, 0}; 242 const double tgt[3] = { 0, 0, 0}; 243 const double up[3] = { 0, 0, 1}; 244 double pix[2]; 245 size_t ix, iy; 246 247 CHK(oct && img); 248 camera_init(&cam, pos, tgt, up, (double)width/(double)height); 249 250 CHK(image_setup(img, width, height, sizeof_image_format(IMAGE_RGB8)*width, 251 IMAGE_RGB8, NULL) == RES_OK); 252 pixels = (unsigned char*)img->pixels; 253 254 time_current(&t0); 255 FOR_EACH(iy, 0, height) { 256 pix[1] = (double)iy / (double)height; 257 FOR_EACH(ix, 0, width) { 258 struct svx_hit hit; 259 struct ray r; 260 size_t ipix = (iy*width + ix)*3/*#channels per pixel*/; 261 pix[0] = (double)ix / (double)width; 262 camera_ray(&cam, pix, &r); 263 264 CHK(svx_tree_trace_ray(oct, r.org, r.dir, r.range, 265 hit_challenge_pass_through, hit_filter, &r, &hit) == RES_OK); 266 if(SVX_HIT_NONE(&hit)) { 267 pixels[ipix+0] = 0; 268 pixels[ipix+1] = 0; 269 pixels[ipix+2] = 0; 270 } else { 271 double N[3]; 272 N[0] = (r.org[0] + hit.distance[0] * r.dir[0]) - scn->sphere_pos[0]; 273 N[1] = (r.org[1] + hit.distance[0] * r.dir[1]) - scn->sphere_pos[1]; 274 N[2] = (r.org[2] + hit.distance[0] * r.dir[2]) - scn->sphere_pos[2]; 275 CHK(d3_normalize(N, N) != 0); 276 N[0] = fabs(N[0]); 277 N[1] = fabs(N[1]); 278 N[2] = fabs(N[2]); 279 pixels[ipix+0] = (unsigned char)(N[0] * 255.0); 280 pixels[ipix+1] = (unsigned char)(N[1] * 255.0); 281 pixels[ipix+2] = (unsigned char)(N[2] * 255.0); 282 } 283 } 284 } 285 time_sub(&t0, time_current(&t1), &t0); 286 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 287 fprintf(stderr, "Render time: %s\n", buf); 288 } 289 290 int 291 main(int argc, char** argv) 292 { 293 struct scene scn; 294 struct ray r; 295 struct image img, img2; 296 FILE* stream = NULL; 297 struct svx_device* dev = NULL; 298 struct svx_tree* oct = NULL; 299 struct svx_tree* oct2 = NULL; 300 struct svx_tree_desc tree_desc = SVX_TREE_DESC_NULL; 301 struct svx_voxel_desc voxel_desc = SVX_VOXEL_DESC_NULL; 302 struct svx_hit hit = SVX_HIT_NULL; 303 struct svx_hit hit2 = SVX_HIT_NULL; 304 const double lower[3] = {-1,-1,-1}; 305 const double upper[3] = { 1, 1, 1}; 306 const size_t def[3] = {32,32,32}; 307 double scnsz[3]; 308 double vxsz; 309 double dst; 310 int accum; 311 (void)argc, (void)argv; 312 313 CHK(svx_device_create(NULL, NULL, 1, &dev) == RES_OK); 314 315 scnsz[0] = upper[0] - lower[0]; 316 scnsz[1] = upper[1] - lower[1]; 317 scnsz[2] = upper[2] - lower[2]; 318 319 scn.origin[0] = lower[0]; 320 scn.origin[1] = lower[1]; 321 scn.origin[2] = lower[2]; 322 scn.vxsz[0] = scnsz[0] / (double)def[0]; 323 scn.vxsz[1] = scnsz[1] / (double)def[1]; 324 scn.vxsz[2] = scnsz[2] / (double)def[2]; 325 scn.sphere_pos[0] = lower[0] + scnsz[0] * 0.5; 326 scn.sphere_pos[1] = lower[1] + scnsz[1] * 0.5; 327 scn.sphere_pos[2] = lower[2] + scnsz[2] * 0.5; 328 scn.sphere_radius = MMIN(MMIN(scnsz[0], scnsz[1]), scnsz[2]) * 0.25; 329 330 voxel_desc.get = voxel_get; 331 voxel_desc.merge = voxels_merge; 332 voxel_desc.challenge_merge = voxels_challenge_merge; 333 voxel_desc.context = &scn; 334 voxel_desc.size = sizeof(char); 335 336 CHK(svx_octree_create(dev, lower, upper, def, &voxel_desc, &oct) == RES_OK); 337 338 /* Duplicate the octree through serialization */ 339 CHK(stream = tmpfile()); 340 CHK(svx_tree_write(oct, stream) == RES_OK); 341 CHK(svx_tree_create_from_stream(dev, stream, &oct2) == RES_BAD_ARG); 342 rewind(stream); 343 CHK(svx_tree_create_from_stream(dev, stream, &oct2) == RES_OK); 344 fclose(stream); 345 346 /*dump_data(stdout, oct, CHAR__, 1, write_scalars);*/ 347 348 #define RT svx_tree_trace_ray 349 d3(r.org, -5,-5, 0); 350 d3(r.dir, 0, 1, 0); 351 d2(r.range, 0, INF); 352 353 CHK(RT(NULL, r.org, r.dir, r.range, NULL, NULL, NULL, &hit) == RES_BAD_ARG); 354 CHK(RT(oct, NULL, r.dir, r.range, NULL, NULL, NULL, &hit) == RES_BAD_ARG); 355 CHK(RT(oct, r.org, NULL, r.range, NULL, NULL, NULL, &hit) == RES_BAD_ARG); 356 CHK(RT(oct, r.org, r.dir, NULL, NULL, NULL, NULL, &hit) == RES_BAD_ARG); 357 CHK(RT(oct, r.org, r.dir, r.range, NULL, NULL, NULL, NULL) == RES_BAD_ARG); 358 359 CHK(RT(oct, r.org, r.dir, r.range, NULL, NULL, NULL, &hit) == RES_OK); 360 CHK(SVX_HIT_NONE(&hit)); 361 362 r.org[0] = 0; 363 CHK(RT(oct, r.org, r.dir, r.range, NULL, NULL, NULL, &hit) == RES_OK); 364 CHK(!SVX_HIT_NONE(&hit)); 365 CHK(eq_eps(hit.distance[0], hit.voxel.lower[1] - r.org[1], 1.e-6)); 366 CHK(eq_eps(hit.distance[1], hit.voxel.upper[1] - r.org[1], 1.e-6)); 367 CHK(hit.voxel.is_leaf); 368 CHK(*(char*)hit.voxel.data == 0); 369 370 /* Use challenge functor to intersect voxels that are not leaves */ 371 CHK(svx_tree_get_desc(oct, &tree_desc) == RES_OK); 372 CHK(RT(oct, r.org, r.dir, r.range, hit_challenge, NULL, NULL, &hit2) == RES_OK); 373 CHK(!SVX_HIT_NONE(&hit2)); 374 CHK(!SVX_VOXEL_EQ(&hit.voxel, &hit2.voxel)); 375 CHK(hit2.voxel.is_leaf == 0); 376 CHK(*(char*)hit2.voxel.data == 1); 377 CHK(eq_eps(hit.distance[0], hit2.distance[0], 1.e-6)); 378 379 /* Use filter function to discard leaves with a value == 0 */ 380 CHK(RT(oct, r.org, r.dir, r.range, hit_challenge_pass_through, hit_filter, 381 &r, &hit) == RES_OK); 382 CHK(!SVX_HIT_NONE(&hit)); 383 CHK(eq_eps(hit.distance[0], hit.voxel.lower[1] - r.org[1], 1.e-6)); 384 CHK(eq_eps(hit.distance[1], hit.voxel.upper[1] - r.org[1], 1.e-6)); 385 CHK(hit.voxel.is_leaf); 386 CHK(*(char*)hit.voxel.data == 1); 387 388 /* Use the ray range to avoid the intersection with the closest voxel */ 389 r.range[1] = nextafter(hit.distance[0], -1); 390 CHK(RT(oct, r.org, r.dir, r.range, hit_challenge_pass_through, hit_filter, 391 &r, &hit) == RES_OK); 392 CHK(SVX_HIT_NONE(&hit)); 393 394 r.range[1] = INF; 395 CHK(RT(oct, r.org, r.dir, r.range, NULL, hit_filter, &r, &hit) == RES_OK); 396 CHK(!SVX_HIT_NONE(&hit)); 397 CHK(*(char*)hit.voxel.data == 1); 398 CHK(hit.voxel.is_leaf); 399 400 /* Use the ray range to discard the closest voxel */ 401 dst = hit.voxel.upper[1] - r.org[1]; 402 r.range[0] = nextafter(dst, DBL_MAX); 403 CHK(RT(oct, r.org, r.dir, r.range, NULL, hit_filter, &r, &hit2) == RES_OK); 404 CHK(!SVX_HIT_NONE(&hit2)); 405 CHK(!SVX_VOXEL_EQ(&hit.voxel, &hit2.voxel)); 406 vxsz = hit2.voxel.upper[1] - hit2.voxel.lower[1]; 407 CHK(eq_eps(hit2.distance[0], r.range[0], 1.e-6)); 408 CHK(eq_eps(hit2.distance[1], dst + vxsz, 1.e-6)); 409 CHK(*(char*)hit.voxel.data == 1); 410 CHK(hit.voxel.is_leaf); 411 412 /* Adjust the ray range to hit the interior of the closest voxel */ 413 r.range[0] = nextafter(hit.distance[0], DBL_MAX); 414 CHK(RT(oct, r.org, r.dir, r.range, NULL, hit_filter, &r, &hit2) == RES_OK); 415 CHK(!SVX_HIT_NONE(&hit2)); 416 CHK(SVX_VOXEL_EQ(&hit.voxel, &hit2.voxel)); 417 CHK(eq_eps(hit2.distance[0], r.range[0], 1.e-6)); 418 CHK(eq_eps(hit2.distance[1], dst, 1.e-6)); 419 CHK(*(char*)hit.voxel.data == 1); 420 CHK(hit.voxel.is_leaf); 421 422 /* Discard the closest voxel with the filter function */ 423 r.range[0] = 0; 424 CHK(RT(oct, r.org, r.dir, r.range, NULL, hit_filter2, &hit.voxel, &hit2) == RES_OK); 425 CHK(eq_eps(hit2.distance[0], hit2.voxel.lower[1] - r.org[1], 1.e-6)); 426 CHK(eq_eps(hit2.distance[1], hit2.voxel.upper[1] - r.org[1], 1.e-6)); 427 CHK(!SVX_HIT_NONE(&hit2)); 428 CHK(!SVX_VOXEL_EQ(&hit.voxel, &hit2.voxel)); 429 CHK(eq_eps(hit2.distance[0], dst, 1.e-6)); 430 CHK(*(char*)hit.voxel.data == 1); 431 CHK(hit.voxel.is_leaf); 432 433 /* Use the filter functor to accumulate the leaves */ 434 accum = 0; 435 CHK(RT(oct, r.org, r.dir, r.range, NULL, hit_filter3, &accum, &hit) == RES_OK); 436 CHK(SVX_HIT_NONE(&hit)); 437 CHK(accum != 0); 438 439 /* Check the root node challenge */ 440 CHK(RT(oct, r.org, r.dir, r.range, hit_challenge_root, NULL, NULL, &hit) 441 == RES_OK); 442 CHK(!SVX_HIT_NONE(&hit)); 443 CHK(d3_eq_eps(hit.voxel.lower, lower, 1.e-6)); 444 CHK(d3_eq_eps(hit.voxel.upper, upper, 1.e-6)); 445 CHK(hit.voxel.depth == 0); 446 CHK(hit.voxel.is_leaf == 0); 447 CHK(*(char*)hit.voxel.data == 1); 448 CHK(eq_eps(hit.distance[0], hit.voxel.lower[1] - r.org[1], 1.e-6)); 449 CHK(eq_eps(hit.distance[1], hit.voxel.upper[1] - r.org[1], 1.e-6)); 450 451 image_init(NULL, &img); 452 image_init(NULL, &img2); 453 draw_image(&img, oct, &scn); 454 draw_image(&img2, oct2, &scn); 455 456 /* Check that using oct or oct2 produces effectively the same image */ 457 check_img_eq(&img, &img2); 458 459 /* Write the drawn image on stdout */ 460 CHK(image_write_ppm_stream(&img, 0/*binary*/, stdout) == RES_OK); 461 image_release(&img); 462 image_release(&img2); 463 464 CHK(svx_tree_ref_put(oct) == RES_OK); 465 CHK(svx_tree_ref_put(oct2) == RES_OK); 466 CHK(svx_device_ref_put(dev) == RES_OK); 467 CHK(mem_allocated_size() == 0); 468 return 0; 469 } 470