test_smc_light_path.c (13258B)
1 /* Copyright (C) 2015-2018, 2021-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 "smc.h" 17 #include "test_smc_utils.h" 18 19 #include <star/s3d.h> 20 #include <star/ssp.h> 21 22 #include <rsys/float2.h> 23 #include <rsys/float3.h> 24 #include <rsys/image.h> 25 #include <rsys/math.h> 26 #include <rsys/mem_allocator.h> 27 28 #define IMG_WIDTH 640 29 #define IMG_HEIGHT 480 30 #define LIGHT_PATH_DEPTH 16 31 #define ALBEDO 0.6f 32 #define EPSILON 1.e-1f 33 #define GAMMA 2.2 34 35 struct cbox_desc{ 36 const float* vertices; 37 const unsigned* indices; 38 }; 39 40 /******************************************************************************* 41 * Box walls data 42 ******************************************************************************/ 43 static const float cbox_walls[] = { 44 552.f, 0.f, 0.f, 45 0.f, 0.f, 0.f, 46 0.f, 559.f, 0.f, 47 552.f, 559.f, 0.f, 48 552.f, 0.f, 548.f, 49 0.f, 0.f, 548.f, 50 0.f, 559.f, 548.f, 51 552.f, 559.f, 548.f 52 }; 53 54 const unsigned cbox_walls_ids[] = { 55 0, 1, 2, 2, 3, 0, /* Bottom */ 56 4, 5, 6, 6, 7, 4, /* Top */ 57 1, 2, 6, 6, 5, 1, /* Left */ 58 0, 3, 7, 7, 4, 0, /* Right */ 59 2, 3, 7, 7, 6, 2 /* Back */ 60 }; 61 62 static const unsigned cbox_walls_ntris = sizeof(cbox_walls_ids)/(3*sizeof(unsigned)); 63 static const unsigned cbox_walls_nverts = sizeof(cbox_walls)/(3*sizeof(float)); 64 static struct cbox_desc cbox_walls_desc = { cbox_walls, cbox_walls_ids }; 65 66 /******************************************************************************* 67 * Short/tall blocks data 68 ******************************************************************************/ 69 static const float cbox_short_block[] = { 70 130.f, 65.f, 0.f, 71 82.f, 225.f, 0.f, 72 240.f, 272.f, 0.f, 73 290.f, 114.f, 0.f, 74 130.f, 65.f, 165.f, 75 82.f, 225.f, 165.f, 76 240.f, 272.f, 165.f, 77 290.f, 114.f, 165.f 78 }; 79 80 static const float cbox_tall_block[] = { 81 423.0f, 247.0f, 0.f, 82 265.0f, 296.0f, 0.f, 83 314.0f, 456.0f, 0.f, 84 472.0f, 406.0f, 0.f, 85 423.0f, 247.0f, 330.f, 86 265.0f, 296.0f, 330.f, 87 314.0f, 456.0f, 330.f, 88 472.0f, 406.0f, 330.f 89 }; 90 91 static const unsigned cbox_block_ids[] = { 92 4, 5, 6, 6, 7, 4, 93 1, 2, 6, 6, 5, 1, 94 0, 3, 7, 7, 4, 0, 95 2, 3, 7, 7, 6, 2, 96 0, 1, 5, 5, 4, 0 97 }; 98 99 static const unsigned cbox_block_ntris = sizeof(cbox_block_ids)/(3*sizeof(unsigned)); 100 101 static const unsigned cbox_short_block_nverts = 102 sizeof(cbox_short_block)/(3*sizeof(float)); 103 static struct cbox_desc cbox_short_block_desc = 104 { cbox_short_block, cbox_block_ids }; 105 106 static const unsigned cbox_tall_block_nverts = 107 sizeof(cbox_tall_block)/(3*sizeof(float)); 108 static struct cbox_desc cbox_tall_block_desc = 109 { cbox_tall_block, cbox_block_ids }; 110 111 /******************************************************************************* 112 * Cornell box callbacks 113 ******************************************************************************/ 114 static INLINE void 115 cbox_get_ids(const unsigned itri, unsigned ids[3], void* data) 116 { 117 const unsigned id = itri * 3; 118 struct cbox_desc* desc = data; 119 CHK(desc != NULL); 120 ids[0] = desc->indices[id + 0]; 121 ids[1] = desc->indices[id + 1]; 122 ids[2] = desc->indices[id + 2]; 123 } 124 125 static INLINE void 126 cbox_get_position(const unsigned ivert, float position[3], void* data) 127 { 128 struct cbox_desc* desc = data; 129 CHK(desc != NULL); 130 position[0] = desc->vertices[ivert*3 + 0]; 131 position[1] = desc->vertices[ivert*3 + 1]; 132 position[2] = desc->vertices[ivert*3 + 2]; 133 } 134 135 /******************************************************************************* 136 * Camera 137 ******************************************************************************/ 138 struct camera { 139 float pos[3]; 140 float x[3], y[3], z[3]; /* Basis */ 141 }; 142 143 static void 144 camera_init(struct camera* cam) 145 { 146 const float pos[3] = { 178.f, -1000.f, 273.f }; 147 const float tgt[3] = { 178.f, 0.f, 273.f }; 148 const float up[3] = { 0.f, 0.f, 1.f }; 149 const float proj_ratio = (float)IMG_WIDTH/(float)IMG_HEIGHT; 150 const float fov_x = (float)PI * 0.25f; 151 float f = 0.f; 152 ASSERT(cam); 153 154 f3_set(cam->pos, pos); 155 f = f3_normalize(cam->z, f3_sub(cam->z, tgt, pos)); CHK(f != 0); 156 f = f3_normalize(cam->x, f3_cross(cam->x, cam->z, up)); CHK(f != 0); 157 f = f3_normalize(cam->y, f3_cross(cam->y, cam->z, cam->x)); CHK(f != 0); 158 f3_divf(cam->z, cam->z, (float)tan(fov_x*0.5f)); 159 f3_divf(cam->y, cam->y, proj_ratio); 160 } 161 162 static void 163 camera_ray 164 (const struct camera* cam, 165 const float pixel[2], 166 float org[3], 167 float dir[3]) 168 { 169 float x[3], y[3], f; 170 ASSERT(cam && pixel && org && dir); 171 172 f3_mulf(x, cam->x, pixel[0]*2.f - 1.f); 173 f3_mulf(y, cam->y, pixel[1]*2.f - 1.f); 174 f3_add(dir, f3_add(dir, x, y), cam->z); 175 f = f3_normalize(dir, dir); CHK(f != 0); 176 f3_set(org, cam->pos); 177 } 178 179 /******************************************************************************* 180 * Integrator 181 ******************************************************************************/ 182 struct integrator_context { 183 struct s3d_scene_view* view; 184 struct camera* cam; 185 float pixel_size[2]; /* Normalized pixel size */ 186 size_t ipixel[2]; /* Image space pixel coordinates */ 187 }; 188 189 static float 190 direct_lighting 191 (struct s3d_scene_view* view, const float pos[3], const float N[3]) 192 { 193 const float light_pos[3] = { 200.f, 200.f, 400.f }; 194 const float flux = 60.0; /* Radiant flux in watt */ 195 const float intensity = (float)(flux/(4.0*PI)); /*Radiant intensity in W/sr*/ 196 struct s3d_hit hit; 197 float len; 198 float wi[3]; 199 float range[2]; 200 201 CHK(view != NULL); 202 CHK(pos != NULL); 203 CHK(N != NULL); 204 CHK(f3_is_normalized(N) == 1); 205 206 f3_sub(wi, light_pos, pos); 207 len = f3_normalize(wi, wi); 208 209 /* Trace shadow ray */ 210 range[0] = EPSILON; 211 range[1] = len; 212 CHK(s3d_scene_view_trace_ray(view, pos, wi, range, NULL, &hit) == RES_OK); 213 if(!S3D_HIT_NONE(&hit)) return 0.f; /* Light is occluded */ 214 215 len *= 0.01f; /* Transform len from centimer to meter */ 216 return 217 intensity / (len*len) /* radiance in W/(sr.m^2) */ 218 * (float)(ALBEDO / PI) /* BRDF */ 219 * f3_dot(wi, N); /* cos(wi,N) */ 220 } 221 222 static float 223 skydome_lighting(const float wi[3]) 224 { 225 const float ground_irradiance = 0.01f; 226 const float sky_irradiance = 12.0f; 227 float u; 228 229 CHK(wi != NULL); 230 u = CLAMP(wi[2], 0.f, 0.05f) * 1.f; 231 return u * sky_irradiance + (1.f - u) * ground_irradiance; 232 } 233 234 static res_T 235 light_path_integrator 236 (void* value, 237 struct ssp_rng* rng, 238 const unsigned ithread, 239 const uint64_t irealisation, 240 void* data) 241 { 242 struct integrator_context* ctx = data; 243 float pix_lower[2], pix_upper[2]; /* Pixel AABB */ 244 float pix_samp[2]; /* Pixel sample */ 245 float ray_org[3], ray_dir[3], ray_range[2]; 246 float L = 0.f; 247 float throughput = 1.f; 248 int idepth; 249 250 (void)ithread, (void)irealisation; 251 252 CHK(value != NULL); 253 CHK(rng != NULL); 254 CHK(data != NULL); 255 256 /* Compute the pixel bound */ 257 pix_lower[0] = (float)ctx->ipixel[0] / (float)IMG_WIDTH; 258 pix_lower[1] = (float)ctx->ipixel[1] / (float)IMG_HEIGHT; 259 f2_add(pix_upper, pix_lower, ctx->pixel_size); 260 261 /* Randomly sample the pixel quad */ 262 pix_samp[0] = (float)ssp_rng_uniform_double(rng, pix_lower[0], pix_upper[0]); 263 pix_samp[1] = (float)ssp_rng_uniform_double(rng, pix_lower[1], pix_upper[1]); 264 265 /* Build a camera ray */ 266 camera_ray(ctx->cam, pix_samp, ray_org, ray_dir); 267 ray_range[0] = 0.f; 268 ray_range[1] = FLT_MAX; 269 270 FOR_EACH(idepth, 0, LIGHT_PATH_DEPTH) { 271 struct s3d_hit hit = S3D_HIT_NULL; 272 float cos_theta; 273 float pdf; 274 float pos[3]; 275 float N[3] = {0, 0, 0}; 276 277 CHK(s3d_scene_view_trace_ray 278 (ctx->view, ray_org, ray_dir, ray_range, NULL, &hit) == RES_OK); 279 280 if(S3D_HIT_NONE(&hit)) { /* Skydome lighting */ 281 L += throughput * skydome_lighting(ray_dir); 282 break; 283 } 284 285 f3_normalize(N, hit.normal); 286 cos_theta = f3_dot(N, ray_dir); 287 if(cos_theta >= 0.f) f3_minus(N, N); 288 f3_add(pos, f3_mulf(pos, ray_dir, hit.distance), ray_org); 289 290 /* Direct lighting */ 291 L += throughput * direct_lighting(ctx->view, pos, N); 292 293 /* New ray */ 294 ssp_ran_hemisphere_cos_float(rng, N, ray_dir, &pdf); 295 cos_theta = f3_dot(N, ray_dir); 296 f3_set(ray_org, pos); 297 ray_range[0] = EPSILON; 298 throughput *= (float)(ALBEDO / PI) / pdf * cos_theta; 299 } 300 *(float*)value = L; 301 return RES_OK; 302 } 303 304 /******************************************************************************* 305 * Application 306 ******************************************************************************/ 307 int 308 main(int argc, char** argv) 309 { 310 struct mem_allocator allocator; 311 struct image img; 312 FILE* fp = NULL; 313 struct integrator_context* contexts; 314 struct s3d_device* dev; 315 struct s3d_scene* scn; 316 struct s3d_scene_view* view; 317 struct s3d_shape* shape; 318 struct s3d_vertex_data attrib; 319 struct smc_device_create_args args = SMC_DEVICE_CREATE_ARGS_DEFAULT; 320 struct smc_device* smc; 321 struct smc_integrator integrator = SMC_INTEGRATOR_NULL; 322 struct smc_estimator** estimators; 323 struct camera cam; 324 size_t ix, iy; 325 float pos[3]; 326 (void)argc, (void)argv; 327 328 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 329 330 CHK(image_init(&allocator, &img) == RES_OK); 331 CHK(image_setup 332 (&img, IMG_WIDTH, IMG_HEIGHT, 3*IMG_WIDTH, IMAGE_RGB8, NULL) == RES_OK); 333 334 CHK(s3d_device_create(NULL, &allocator, 0, &dev) == RES_OK); 335 CHK(s3d_scene_create(dev, &scn) == RES_OK); 336 337 attrib.usage = S3D_POSITION; 338 attrib.type = S3D_FLOAT3; 339 attrib.get = cbox_get_position; 340 341 CHK(s3d_shape_create_mesh(dev, &shape) == RES_OK); 342 CHK(s3d_mesh_setup_indexed_vertices(shape, cbox_walls_ntris, cbox_get_ids, 343 cbox_walls_nverts, &attrib, 1, &cbox_walls_desc) == RES_OK); 344 CHK(s3d_scene_attach_shape(scn, shape) == RES_OK); 345 CHK(s3d_shape_ref_put(shape) == RES_OK); 346 347 CHK(s3d_shape_create_mesh(dev, &shape) == RES_OK); 348 CHK(s3d_mesh_setup_indexed_vertices(shape, cbox_block_ntris, cbox_get_ids, 349 cbox_short_block_nverts, &attrib, 1, &cbox_short_block_desc) == RES_OK); 350 CHK(s3d_scene_attach_shape(scn, shape) == RES_OK); 351 CHK(s3d_shape_ref_put(shape) == RES_OK); 352 353 CHK(s3d_shape_create_mesh(dev, &shape) == RES_OK); 354 CHK(s3d_mesh_setup_indexed_vertices(shape, cbox_block_ntris, cbox_get_ids, 355 cbox_tall_block_nverts, &attrib, 1, &cbox_tall_block_desc) == RES_OK); 356 CHK(s3d_scene_attach_shape(scn, shape) == RES_OK); 357 CHK(s3d_shape_ref_put(shape) == RES_OK); 358 359 CHK(s3d_scene_instantiate(scn, &shape) == RES_OK); 360 CHK(s3d_scene_ref_put(scn) == RES_OK); 361 CHK(s3d_instance_set_position(shape, f3(pos, -100.f, 0.f, -2.f)) == RES_OK); 362 363 CHK(s3d_scene_create(dev, &scn) == RES_OK); 364 CHK(s3d_scene_attach_shape(scn, shape) == RES_OK); 365 CHK(s3d_shape_ref_put(shape) == RES_OK); 366 367 CHK(s3d_scene_view_create(scn, S3D_TRACE, &view) == RES_OK); 368 369 args.allocator = &allocator; 370 args.verbose = 1; 371 CHK(smc_device_create(&args, &smc) == RES_OK); 372 373 integrator.integrand = light_path_integrator; 374 integrator.type = &smc_float; 375 integrator.max_realisations = 64; 376 377 contexts = MEM_CALLOC 378 (&allocator, IMG_WIDTH*IMG_HEIGHT, sizeof(struct integrator_context)); 379 CHK(contexts != NULL); 380 estimators = MEM_CALLOC 381 (&allocator, IMG_WIDTH*IMG_HEIGHT, sizeof(struct smc_estimator*)); 382 CHK(estimators != NULL); 383 384 camera_init(&cam); 385 386 FOR_EACH(iy, 0, IMG_HEIGHT) { 387 FOR_EACH(ix, 0, IMG_WIDTH) { 388 const size_t ictx = iy * IMG_WIDTH + ix; 389 contexts[ictx].view = view; 390 contexts[ictx].cam = &cam; 391 contexts[ictx].pixel_size[0] = 1.f / (float)IMG_WIDTH; 392 contexts[ictx].pixel_size[1] = 1.f / (float)IMG_HEIGHT; 393 contexts[ictx].ipixel[0] = ix; 394 contexts[ictx].ipixel[1] = iy; 395 }} 396 397 CHK(smc_solve_N(smc, &integrator, IMG_WIDTH * IMG_HEIGHT, contexts, 398 sizeof(struct integrator_context), estimators) == RES_OK); 399 400 FOR_EACH(iy, 0, IMG_HEIGHT) { 401 FOR_EACH(ix, 0, IMG_WIDTH) { 402 const size_t iestimator = (iy*IMG_WIDTH + ix); 403 uint8_t* pix = (uint8_t*)(img.pixels + iy*img.pitch + ix*3/*RGB*/); 404 struct smc_estimator_status status; 405 float col; 406 unsigned char colu; 407 408 CHK(smc_estimator_get_status(estimators[iestimator], &status) == RES_OK); 409 col = (float)pow(SMC_FLOAT(status.E), 1.0/GAMMA); /* Gamma correction */ 410 colu = (uint8_t)(CLAMP(col, 0.f, 1.f) * 255.f); /* Float to U8 */ 411 pix[0] = pix[1] = pix[2] = colu; 412 CHK(smc_estimator_ref_put(estimators[iestimator]) == RES_OK); 413 }} 414 415 CHK(fp = fopen("image.ppm", "w")); 416 image_write_ppm_stream(&img, 0, fp); 417 CHK(fclose(fp) == 0); 418 419 CHK(image_release(&img) == RES_OK); 420 MEM_RM(&allocator, contexts); 421 MEM_RM(&allocator, estimators); 422 423 CHK(s3d_scene_view_ref_put(view) == RES_OK); 424 425 CHK(s3d_device_ref_put(dev) == RES_OK); 426 CHK(s3d_scene_ref_put(scn) == RES_OK); 427 CHK(smc_device_ref_put(smc) == RES_OK); 428 429 check_memory_allocator(&allocator); 430 mem_shutdown_proxy_allocator(&allocator); 431 CHK(mem_allocated_size() == 0); 432 return 0; 433 }