star-schiff

Library for estimating radiative properties
git clone git://git.meso-star.com/star-schiff.git
Log | Files | Refs | README | LICENSE

commit ccdf7c495e787245e14b31873c024555741e7c0b
parent cb59b68ec96f35bae58e60379c51767150b3e65c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  7 Oct 2015 09:39:22 +0200

Draft of the projective footprint sampling.

Build the transformation matrix from world space to footprint space with
respect to a sampled direction. Transform the Scene AABB in footprint
space to define the ray sampling domain. Finally ray-cast rays in
footprint space to generate an image of the geometry footprint.

Diffstat:
Msrc/salgae_estimator.c | 153++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 152 insertions(+), 1 deletion(-)

diff --git a/src/salgae_estimator.c b/src/salgae_estimator.c @@ -29,8 +29,12 @@ #include "salgae.h" #include "salgae_device.h" +#include <rsys/float2.h> +#include <rsys/float3.h> +#include <rsys/image.h> #include <rsys/logger.h> #include <rsys/mem_allocator.h> +#include <rsys/stretchy_array.h> #include <star/s3d.h> #include <star/ssp.h> @@ -48,6 +52,108 @@ struct salgae_estimator { /******************************************************************************* * Helper functions ******************************************************************************/ +static FINLINE uint16_t +morton2D_decode(const uint32_t u32) +{ + uint32_t x = u32 & 0x55555555; + x = (x | (x >> 1)) & 0x33333333; + x = (x | (x >> 2)) & 0x0F0F0F0F; + x = (x | (x >> 4)) & 0x00FF00FF; + x = (x | (x >> 8)) & 0x0000FFFF; + return (uint16_t)x; +} + +static FINLINE void +build_transform(const float pos[3], const float basis[9], float transform[12]) +{ + ASSERT(pos && basis && transform); + + /* `transform' is the matrix from world to footprint space while basis is a + * 3x3 matrix from local to world space. The linear part of the transform + * matrix is thus the inverse of the `basis' matrix. Since basis encodes an + * orthonormal transformation, a simple transposition is sufficient to + * inverse it */ + f3(transform + 0, basis[0], basis[3], basis[6]); + f3(transform + 3, basis[1], basis[4], basis[8]); + f3(transform + 6, basis[2], basis[5], basis[9]); + /* Affine part of the transform matrix */ + transform[9] = -f3_dot(basis+0, pos); + transform[10] = -f3_dot(basis+3, pos); + transform[11] = -f3_dot(basis+6, pos); +} + +static void +trace_rays + (struct s3d_scene* scn, + const float basis[9], + const float pos[3], + const float lower[3], /* Transform space lower boundary of the RT volume */ + const float upper[3], /* Transform space upper boundary of the RT volume */ + const char* name) +{ + +#define DEFINITION 128 + STATIC_ASSERT(IS_POW2(DEFINITION), Invalid_thumbnail_definition); + + float axis_x[3], axis_y[3], axis_z[3]; + float plane_size[2]; + float pixel_size[2]; + float ray_org[3]; + float org[3]; + float depth; + float range[2] = {0, FLT_MAX}; + uint32_t npixels = (uint32_t)(DEFINITION * DEFINITION); + uint32_t ipixel; + uint16_t ipixel_x, ipixel_y; + unsigned char* img = NULL; + ASSERT(basis && pos && lower && upper); + + img = sa_add(img, DEFINITION*DEFINITION); + ASSERT(img != NULL); + + f2_sub(plane_size, upper, lower); + f3_mulf(axis_x, basis + 0, plane_size[0] * 0.5f); + f3_mulf(axis_y, basis + 3, plane_size[1] * 0.5f); + f3_set(axis_z, basis + 6); + + depth = lower[2] + 1.e-6f; /* Ensure to be outside the RT volume */ + f3_add(org, pos, f3_mulf(org, basis + 9, depth)); + + pixel_size[0] = 1.f / (float)DEFINITION; + pixel_size[1] = 1.f / (float)DEFINITION; + + FOR_EACH(ipixel, 0, npixels) { + struct s3d_hit hit; + float sample[2], x[3], y[3]; + size_t idst; + + ipixel_x = morton2D_decode(ipixel); + ipixel_y = morton2D_decode(ipixel>>1); + sample[0] = ((float)ipixel_x + 0.5f) * pixel_size[0]; + sample[1] = ((float)ipixel_y + 0.5f) * pixel_size[1]; + + f3_mulf(x, axis_x, sample[0]*2.f - 1.f); + f3_mulf(y, axis_y, sample[1]*2.f - 1.f); + f3_add(ray_org, f3_add(ray_org, x, y), org); + + S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit)); + + idst = (size_t)(ipixel_y * DEFINITION + ipixel_x); + if(S3D_HIT_NONE(&hit)) { + img[idst] = 0; + } else { + float N[3], cosine; + f3_normalize(N, hit.normal); + if(0 > (cosine = f3_dot(N, axis_z))) { + cosine = f3_dot(f3_minus(N, N), axis_z); + } + img[idst] = (unsigned char)(cosine*255.f); + } + } + image_ppm_write(name, DEFINITION, DEFINITION, 1, img); + sa_release(img); +} + static res_T radiative_path_length (struct salgae_device* dev, @@ -86,9 +192,54 @@ radiative_path_length S3D(scene_begin_session(scene, S3D_TRACE)); FOR_EACH(i, 0, NDIRS) { + char thumb_name[64]; float dir[4]; + float centroid[3]; + float basis[9]; + float lower[3], upper[3]; + float transform[16]; + float pt[6][3]; + int i; + ssp_ran_sphere_uniform(rng, dir); - /* TODO */ + S3D(scene_get_aabb(scene, lower, upper)); + + /* AABB vertex layout + * 6-------7 + * /' /| + * 2-------3 | + * | 4.....|.5 + * |, |/ + * 0-------1 */ + f3(pt[0], lower[0], lower[1], lower[2]); + f3(pt[1], upper[0], lower[1], lower[2]); + f3(pt[2], lower[0], upper[1], lower[2]); + f3(pt[3], upper[0], upper[1], lower[2]); + f3(pt[4], lower[0], lower[1], upper[2]); + f3(pt[5], upper[0], lower[1], upper[2]); + f3(pt[6], lower[0], upper[1], upper[2]); + f3(pt[7], upper[0], upper[1], upper[2]); + + /* Build the transformation matrix from world to footprint space. Use the + * AABB centroid as the origin of the footprint space. */ + f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f); + f33_basis(basis, dir); + build_transform(centroid, basis, transform); + + /* Transform the AABB from world to footprint space */ + FOR_EACH(i, 0, 6) { + f3_add(pt[i], f33_mulf3(pt[i], transform, pt[i]), transform + 9); + } + + /* Compute the AABB of the footprint */ + f3_splat(lower, FLT_MAX); + f3_splat(upper,-FLT_MAX); + FOR_EACH(i, 0, 6) { + f3_min(lower, lower, pt[i]); + f3_max(upper, upper, pt[i]); + } + sprintf(thumb_name, "thumb_%d.ppm", i); + trace_rays(scene, basis, centroid, lower, upper, thumb_name); } S3D(scene_end_session(scene));