star-schiff

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

commit d52298d1d16f4214dbaafcddae6a9ab0389630d0
parent 7f8e5e1d44a19afb0a0e0d00b3f5932fa5062071
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  7 Oct 2015 18:10:20 +0200

Add comments to the estimator code

Diffstat:
Msrc/salgae_estimator.c | 87++++++++++++++++++++++++++++++++++++++++++-------------------------------------
Msrc/test_salgae_estimator.c | 2+-
2 files changed, 47 insertions(+), 42 deletions(-)

diff --git a/src/salgae_estimator.c b/src/salgae_estimator.c @@ -39,7 +39,7 @@ #include <star/s3d.h> #include <star/ssp.h> -#define NDIRS 5 +#define NDIRS 1 /* # Sampled directions */ struct salgae_estimator { size_t nsteps; @@ -63,18 +63,18 @@ morton2D_decode(const uint32_t u32) return (uint16_t)x; } +/* Compute an orthonormal basis where `dir' is the Z axis */ static void build_basis(const float dir[3], float basis[9]) { float dir_abs[3]; float axis_min[3]; int i; - ASSERT(dir && basis); + ASSERT(dir && basis && f3_is_normalized(dir)); /* Define which axis direction is minimal and use its unit vector to compute - * a vector ortogonal to `dir' named `axis0'. This ensures that the unit - * vector and `dir' are not colinear and thus that their cross product is not - * a zero vector */ + * a vector ortogonal to `dir'. This ensures that the unit vector and `dir' + * are not colinear and thus that their cross product is not a zero vector */ FOR_EACH(i, 0, 3) dir_abs[i] = absf(dir[i]); if(dir_abs[0] < dir_abs[1]) { if(dir_abs[0] < dir_abs[2]) f3(axis_min, 1.f, 0.f, 0.f); @@ -85,87 +85,91 @@ build_basis(const float dir[3], float basis[9]) } f3_normalize(basis + 0, f3_cross(basis + 0, dir, axis_min)); - f3_normalize(basis + 3, f3_cross(basis + 3, dir, basis + 0)); + f3_cross(basis + 3, dir, basis + 0); f3_set(basis + 6, dir); } +/* Compute a 3x4 affine transformation from a orthonormal basis and the + * position of the origin */ static FINLINE void -build_transform(const float pos[3], const float basis[9], float transform[12]) +build_transform + (const float pos[3], /* Position of the origin */ + const float basis[9], /* Column major */ + float transform[12]) /* Column major */ { 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 */ + /* The linear part of the transformamatrix is the inverse of the + * orthonormal basis (i.e. its transpose) */ f3(transform + 0, basis[0], basis[3], basis[6]); f3(transform + 3, basis[1], basis[4], basis[7]); f3(transform + 6, basis[2], basis[5], basis[8]); + /* 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); } +/* Cast rays into the RT volume */ static void -trace_rays +draw_thumbnail (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) + const float basis[9], /* World space basis of the RT volume */ + const float pos[3], /* World space centroid of the RT volume */ + const float lower[3], /* Lower boundary of the RT volume */ + const float upper[3], /* Upper boundary of the RT volume */ + const char* name) /* Filename of the thumbnail */ { #define DEFINITION 128 STATIC_ASSERT(IS_POW2(DEFINITION), Invalid_thumbnail_definition); - + const float range[2] = { 0, FLT_MAX }; + const float pixel_size = 1.f/(float)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; + uint32_t mcode; unsigned char* img = NULL; ASSERT(basis && pos && lower && upper); + /* Allocate the thumbnail buffer */ img = sa_add(img, DEFINITION*DEFINITION); ASSERT(img != NULL); + /* Define the projection axis */ 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-3f; /* Ensure to be outside the RT volume */ - f3_add(org, pos, f3_mulf(org, axis_z, depth)); + /* Compute the origin of the projection plane */ + f3_add(org, pos, f3_mulf(org, axis_z, lower[2])); - pixel_size[0] = 1.f/(float)DEFINITION; - pixel_size[1] = 1.f/(float)DEFINITION; - - FOR_EACH(ipixel, 0, npixels) { + FOR_EACH(mcode, 0, npixels) { struct s3d_hit hit; + size_t ipixel; float sample[2], x[3], y[3]; - size_t idst; + uint16_t ipixel_x, ipixel_y; - 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]; + /* Compute the sample position in [0, 1)^2 onto the projection plane */ + ipixel_x = morton2D_decode(mcode); + ipixel_y = morton2D_decode(mcode>>1); + sample[0] = ((float)ipixel_x + 0.5f) * pixel_size; + sample[1] = ((float)ipixel_y + 0.5f) * pixel_size; + /* Compute the ray origin */ 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); + /* Simple shading */ + ipixel = (size_t)(ipixel_y * DEFINITION + ipixel_x); if(S3D_HIT_NONE(&hit)) { - img[idst] = 0; + img[ipixel] = 0; } else { float N[3] = {0.f,0.f,0.f}; float cosine; @@ -173,7 +177,7 @@ trace_rays if(0 > (cosine = f3_dot(N, axis_z))) { cosine = f3_dot(f3_minus(N, N), axis_z); } - img[idst] = (unsigned char)(cosine*255.f); + img[ipixel] = (unsigned char)(cosine*255.f + 0.5f/*Round*/); } } image_ppm_write(name, DEFINITION, DEFINITION, 1, img); @@ -260,14 +264,15 @@ radiative_path_length } /* Compute the AABB of the footprint */ - f3_splat(lower, FLT_MAX); - f3_splat(upper,-FLT_MAX); + f3_splat(lower, FLT_MAX); f3_splat(upper,-FLT_MAX); FOR_EACH(i, 0, 8) { f3_min(lower, lower, pt[i]); f3_max(upper, upper, pt[i]); } + + /* Compute an image of the shape transformed in the footprint space */ sprintf(thumb_name, "%d_%d.ppm", istep, idir); - trace_rays(scene, basis, centroid, lower, upper, thumb_name); + draw_thumbnail(scene, basis, centroid, lower, upper, thumb_name); } S3D(scene_end_session(scene)); diff --git a/src/test_salgae_estimator.c b/src/test_salgae_estimator.c @@ -184,7 +184,7 @@ main(int argc, char** argv) CHECK(salgae_integrate(NULL, NULL, &desc, 0, 1, &estimator), RES_BAD_ARG); CHECK(salgae_integrate(dev, NULL, &desc, 0, 1, &estimator), RES_BAD_ARG); CHECK(salgae_integrate(NULL, rng, &desc, 0, 1, &estimator), RES_BAD_ARG); - CHECK(salgae_integrate(dev, rng, &desc, 0, 10, &estimator), RES_OK); + CHECK(salgae_integrate(dev, rng, &desc, 0, 100, &estimator), RES_OK); CHECK(salgae_estimator_ref_get(NULL), RES_BAD_ARG); CHECK(salgae_estimator_ref_get(estimator), RES_OK);