star-schiff

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

commit 613498cb1e8a4edfbcaa9e03d8152f0aab99f34b
parent ccdf7c495e787245e14b31873c024555741e7c0b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  7 Oct 2015 14:51:24 +0200

Begin the tests on the Star-Algae estimator

Diffstat:
Mcmake/CMakeLists.txt | 7++++---
Msrc/salgae.h | 6+++++-
Msrc/salgae_estimator.c | 66++++++++++++++++++++++++++++++++++++++++++++++++------------------
Asrc/test_salgae_estimator.c | 176+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 233 insertions(+), 22 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -38,9 +38,9 @@ option(NO_TEST "Disable the test" OFF) ################################################################################ find_package(RCMake 0.2.1 REQUIRED) find_package(RSys 0.3 REQUIRED) -find_package(StarSP 0.2 REQUIRED) -find_package(StarMC 0.2 REQUIRED) -find_package(Star3D 0.2 REQUIRED) +find_package(StarSP 0.3 REQUIRED) +find_package(StarMC 0.3 REQUIRED) +find_package(Star3D 0.3 REQUIRED) include_directories( ${RSys_INCLUDE_DIR} @@ -93,6 +93,7 @@ if(NOT NO_TEST) add_test(${_name} ${_name}) endfunction(new_test) new_test(test_salgae_device) + new_test(test_salgae_estimator) endif() ################################################################################ diff --git a/src/salgae.h b/src/salgae.h @@ -77,7 +77,11 @@ struct salgae_integrator_desc { struct s3d_shape* shape, /* Sampled shape */ void* sampler_context); void* sampler_context; - const double scattering_angles_count; /* # of scattering angles to handle */ + size_t scattering_angles_count; /* # of scattering angles to handle */ +}; + +static const struct salgae_integrator_desc SALGAE_NULL_INTEGRATOR_DESC = { + NULL, NULL, 0 }; /* Result of the Monte Carlo estimation */ diff --git a/src/salgae_estimator.c b/src/salgae_estimator.c @@ -63,6 +63,32 @@ morton2D_decode(const uint32_t u32) return (uint16_t)x; } +static void +build_basis(const float dir[3], float basis[9]) +{ + float dir_abs[3]; + float axis_min[3]; + int i; + ASSERT(dir && basis); + + /* 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 */ + 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); + else f3(axis_min, 0.f, 0.f, 1.f); + } else { + if(dir_abs[1] < dir_abs[2]) f3(axis_min, 0.f, 1.f, 0.f); + else f3(axis_min, 0.f, 0.f, 1.f); + } + + f3_normalize(basis + 0, f3_cross(basis + 0, dir, axis_min)); + f3_normalize(basis + 3, f3_cross(basis + 3, dir, basis + 0)); + f3_set(basis + 6, dir); +} + static FINLINE void build_transform(const float pos[3], const float basis[9], float transform[12]) { @@ -74,10 +100,10 @@ build_transform(const float pos[3], const float basis[9], float transform[12]) * 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]); + 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[9 ] = -f3_dot(basis+0, pos); transform[10] = -f3_dot(basis+3, pos); transform[11] = -f3_dot(basis+6, pos); } @@ -91,7 +117,6 @@ trace_rays 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); @@ -116,11 +141,11 @@ trace_rays 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)); + depth = lower[2] + 1.e-3f; /* Ensure to be outside the RT volume */ + f3_add(org, pos, f3_mulf(org, axis_z, depth)); - pixel_size[0] = 1.f / (float)DEFINITION; - pixel_size[1] = 1.f / (float)DEFINITION; + pixel_size[0] = 1.f/(float)DEFINITION; + pixel_size[1] = 1.f/(float)DEFINITION; FOR_EACH(ipixel, 0, npixels) { struct s3d_hit hit; @@ -142,7 +167,8 @@ trace_rays if(S3D_HIT_NONE(&hit)) { img[idst] = 0; } else { - float N[3], cosine; + float N[3] = {0.f,0.f,0.f}; + float cosine; f3_normalize(N, hit.normal); if(0 > (cosine = f3_dot(N, axis_z))) { cosine = f3_dot(f3_minus(N, N), axis_z); @@ -152,6 +178,7 @@ trace_rays } image_ppm_write(name, DEFINITION, DEFINITION, 1, img); sa_release(img); +#undef DEFINITION } static res_T @@ -163,7 +190,7 @@ radiative_path_length struct s3d_shape* shape = NULL; struct s3d_scene* scene = NULL; struct salgae_material material = SALGAE_NULL_MATERIAL; - size_t i; + int idir; res_T res = RES_OK; ASSERT(dev && rng && desc); @@ -191,14 +218,14 @@ radiative_path_length } S3D(scene_begin_session(scene, S3D_TRACE)); - FOR_EACH(i, 0, NDIRS) { + FOR_EACH(idir, 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]; + float transform[12]; + float pt[8][3]; int i; ssp_ran_sphere_uniform(rng, dir); @@ -223,22 +250,22 @@ radiative_path_length /* 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_basis(dir, basis); build_transform(centroid, basis, transform); /* Transform the AABB from world to footprint space */ - FOR_EACH(i, 0, 6) { + FOR_EACH(i, 0, 8) { 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) { + FOR_EACH(i, 0, 8) { f3_min(lower, lower, pt[i]); f3_max(upper, upper, pt[i]); } - sprintf(thumb_name, "thumb_%d.ppm", i); + sprintf(thumb_name, "thumb_%d.ppm", idir); trace_rays(scene, basis, centroid, lower, upper, thumb_name); } S3D(scene_end_session(scene)); @@ -263,9 +290,12 @@ static void estimator_release(ref_T* ref) { struct salgae_estimator* estimator; + struct salgae_device* dev; ASSERT(ref); estimator = CONTAINER_OF(ref, struct salgae_estimator, ref); - MEM_RM(estimator->dev->allocator, estimator); + dev = estimator->dev; + MEM_RM(dev->allocator, estimator); + SALGAE(device_ref_put(dev)); } static res_T diff --git a/src/test_salgae_estimator.c b/src/test_salgae_estimator.c @@ -0,0 +1,176 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is governed by the CeCILL license under French law and + * abiding by the rules of distribution of free software. You can use, + * modify and/or redistribute the software under the terms of the CeCILL + * license as circulated by CEA, CNRS and INRIA at the following URL + * "http://www.cecill.info". + * + * As a counterpart to the access to the source code and rights to copy, + * modify and redistribute granted by the license, users are provided only + * with a limited warranty and the software's author, the holder of the + * economic rights, and the successive licensors have only limited + * liability. + * + * In this respect, the user's attention is drawn to the risks associated + * with loading, using, modifying and/or developing or reproducing the + * software by the user in light of its specific status of free software, + * that may mean that it is complicated to manipulate, and that also + * therefore means that it is reserved for developers and experienced + * professionals having in-depth computer knowledge. Users are therefore + * encouraged to load and test the software's suitability as regards their + * requirements in conditions enabling the security of their systems and/or + * data to be ensured and, more generally, to use and operate it in the + * same conditions as regards security. + * + * The fact that you are presently reading this means that you have had + * knowledge of the CeCILL license and that you accept its terms. */ + +#include "salgae.h" +#include "test_salgae_utils.h" + +#include <rsys/float3.h> + +#include <star/s3d.h> +#include <star/ssp.h> + +static const float verts[] = { + /* Box */ + 552.f, 0.f, 0.f, + 0.f, 0.f, 0.f, + 0.f, 559.f, 0.f, + 552.f, 559.f, 0.f, + 552.f, 0.f, 548.f, + 0.f, 0.f, 548.f, + 0.f, 559.f, 548.f, + 552.f, 559.f, 548.f, + /* Short block */ + 130.f, 65.f, 0.f, + 82.f, 225.f, 0.f, + 240.f, 272.f, 0.f, + 290.f, 114.f, 0.f, + 130.f, 65.f, 165.f, + 82.f, 225.f, 165.f, + 240.f, 272.f, 165.f, + 290.f, 114.f, 165.f, + /* Tall block */ + 423.0f, 247.0f, 0.f, + 265.0f, 296.0f, 0.f, + 314.0f, 456.0f, 0.f, + 472.0f, 406.0f, 0.f, + 423.0f, 247.0f, 330.f, + 265.0f, 296.0f, 330.f, + 314.0f, 456.0f, 330.f, + 472.0f, 406.0f, 330.f +}; +static const unsigned nverts = (unsigned)(sizeof(verts) / (sizeof(float)*3)); +static const unsigned ids[] = { + /* Box */ + 0, 1, 2, 2, 3, 0, + 4, 5, 6, 6, 7, 4, + 1, 2, 6, 6, 5, 1, + 0, 3, 7, 7, 4, 0, + 2, 3, 7, 7, 6, 2, + /* Short block */ + 12, 13, 14, 14, 15, 12, + 9, 10, 14, 14, 13, 9, + 8, 11, 15, 15, 12, 8, + 10, 11, 15, 15, 14, 10, + 8, 9, 13, 13, 12, 8, + /* Tall block */ + 20, 21, 22, 22, 23, 20, + 17, 18, 22, 22, 21, 17, + 16, 19, 23, 23, 20, 16, + 18, 19, 23, 23, 22, 18, + 16, 17, 21, 21, 20, 16 +}; +static const unsigned nids = (unsigned)(sizeof(ids)/sizeof(uint32_t)); + +static void +get_ids(const unsigned itri, unsigned indices[3], void* data) +{ + const unsigned id = itri * 3; + (void)data; + indices[0] = ids[id + 0]; + indices[1] = ids[id + 1]; + indices[2] = ids[id + 2]; +} + +static void +get_pos(const unsigned ivert, float vertex[3], void* data) +{ + (void)data; + f3_set(vertex, verts + (ivert*3)); +} + +static res_T +sample_geometry + (struct ssp_rng* rng, + struct salgae_material* mtl, + struct s3d_shape* shape, + void* sampler_context) +{ + struct s3d_vertex_data attrib; + (void)sampler_context; + NCHECK(rng, NULL); + NCHECK(mtl, NULL); + NCHECK(shape, NULL); + + attrib.usage = S3D_POSITION; + attrib.type = S3D_FLOAT3; + attrib.get = get_pos; + + *mtl = SALGAE_NULL_MATERIAL; + return s3d_mesh_setup_indexed_vertices + (shape, nids/3, get_ids, nverts, &attrib, 1, NULL); +} + +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct salgae_device* dev; + struct salgae_estimator* estimator; + struct salgae_integrator_desc desc = SALGAE_NULL_INTEGRATOR_DESC; + struct ssp_rng* rng; + (void)argc, (void)argv; + + mem_init_proxy_allocator(&allocator, &mem_default_allocator); + + CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK); + CHECK(salgae_device_create(NULL, &allocator, 1, NULL, &dev), RES_OK); + + desc.sample_geometry = sample_geometry; + desc.scattering_angles_count = 1; + CHECK(salgae_integrate(NULL, NULL, NULL, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(dev, NULL, NULL, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(NULL, rng, NULL, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(dev, rng, NULL, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(NULL, NULL, &desc, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(dev, NULL, &desc, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(NULL, rng, &desc, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(dev, rng, &desc, 0, 0, NULL), RES_BAD_ARG); + CHECK(salgae_integrate(NULL, NULL, NULL, 0, 1, &estimator), RES_BAD_ARG); + CHECK(salgae_integrate(dev, NULL, NULL, 0, 1, &estimator), RES_BAD_ARG); + CHECK(salgae_integrate(NULL, rng, NULL, 0, 1, &estimator), RES_BAD_ARG); + CHECK(salgae_integrate(dev, rng, NULL, 0, 1, &estimator), RES_BAD_ARG); + 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, 1, &estimator), RES_OK); + + CHECK(salgae_estimator_ref_get(NULL), RES_BAD_ARG); + CHECK(salgae_estimator_ref_get(estimator), RES_OK); + CHECK(salgae_estimator_ref_put(NULL), RES_BAD_ARG); + CHECK(salgae_estimator_ref_put(estimator), RES_OK); + CHECK(salgae_estimator_ref_put(estimator), RES_OK); + + CHECK(salgae_device_ref_put(dev), RES_OK); + CHECK(ssp_rng_ref_put(rng), RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHECK(mem_allocated_size(), 0); + return 0; +} +