schiff

Estimate the radiative properties of soft particless
git clone git://git.meso-star.com/schiff.git
Log | Files | Refs | README | LICENSE

commit cf5b83c4776a3160e169fe288ebb213e8f9a80e4
parent 665a94251e0a4723d7a7cb47d6e984687f12ced0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon,  2 Nov 2015 15:29:53 +0100

Begin the implementation of the thumbnail dump

Diffstat:
Mcmake/CMakeLists.txt | 6++++--
Asrc/schiff_thumbnail.c | 236+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/schiff_thumbnail.h | 42++++++++++++++++++++++++++++++++++++++++++
3 files changed, 282 insertions(+), 2 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -68,12 +68,14 @@ set(SCHIFF_FILES_SRC schiff_args.c schiff_geometry.c schiff_mesh.c - schiff_optical_properties.c) + schiff_optical_properties.c + schiff_thumbnail.c) set(SCHIFF_FILES_INC schiff_args.h schiff_geometry.h schiff_mesh.h - schiff_optical_properties.h) + schiff_optical_properties.h + schiff_thumbnail.h) set(SCHIFF_FILES_DOC COPYING.fr COPYING.en README.md) # Prepend each file in the `SCHIFF_FILES_<SRC|INC>' list by the absolute diff --git a/src/schiff_thumbnail.c b/src/schiff_thumbnail.c @@ -0,0 +1,236 @@ +/* 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 "schiff_thumbnail.h" + +#include <rsys/float2.h> +#include <rsys/float3.h> +#include <rsys/image.h> +#include <rsys/stretchy_array.h> + +#include <star/s3d.h> +#include <star/ssp.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static 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; +} + +/* 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 && f3_is_normalized(dir)); + + /* Define which axis direction is minimal and use its unit vector to compute + * 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); + 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_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], /* Position of the origin */ + const float basis[9], /* Column major */ + float transform[12]) /* Column major */ +{ + ASSERT(pos && basis && transform); + + /* The linear part of the transform matrix 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); +} + +static void +draw_thumbnail + (struct s3d_scene* scn, + const size_t definition, /* Thumbnail definition */ + 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) +{ + const float range[2] = { 0, FLT_MAX }; + float axis_x[3], axis_y[3], axis_z[3]; + float plane_size[2]; + float ray_org[3]; + float org[3]; + const size_t definition_adjusted = round_up_pow2(definition); + const float pixel_size = 1.f/(float)definition; + uint32_t npixels = (uint32_t)(definition_adjusted * definition_adjusted); + uint32_t mcode; + unsigned char* img = NULL; + ASSERT(definition && basis && pos && lower && upper && name); + + /* 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); + + /* Compute the origin of the projection plane */ + f3_add(org, pos, f3_mulf(org, axis_z, lower[2])); + + FOR_EACH(mcode, 0, npixels) { + struct s3d_hit hit; + size_t ipixel; + float sample[2], x[3], y[3]; + uint16_t ipixel_x, ipixel_y; + + /* Compute the sample position in [0, 1)^2 onto the projection plane */ + ipixel_x = morton2D_decode(mcode); + if(ipixel_x > definition) continue; /* Discard out of bound pixel */ + ipixel_y = morton2D_decode(mcode>>1); + if(ipixel_y > definition) continue; /* Discard out of bound pixel */ + 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)); + + /* Simple shading */ + ipixel = (size_t)(ipixel_y * definition + ipixel_x); + if(S3D_HIT_NONE(&hit)) { + img[ipixel] = 0; + } else { + 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); + } + img[ipixel] = (unsigned char)(cosine*255.f + 0.5f/*Round*/); + } + } + image_ppm_write(name, (int)definition, (int)definition, 1, img); + sa_release(img); +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +void +schiff_thumbnail + (struct s3d_scene* scn, + const size_t definition) +{ + struct ssp_rng* rng = NULL; + float lower[3], upper[3]; + float aabb_pt[8][3]; + float basis[9]; + float transform[12]; + float centroid[3]; + float dir[4]; + int i; + ASSERT(scn); + + SSP(rng_create(NULL, &ssp_rng_threefry, &rng)); + S3D(scene_get_aabb(scn, lower, upper)); + + /* AABB vertex layout + * 6-------7 + * /' /| + * 2-------3 | + * | 4.....|.5 + * |, |/ + * 0-------1 */ + f3(aabb_pt[0], lower[0], lower[1], lower[2]); + f3(aabb_pt[1], upper[0], lower[1], lower[2]); + f3(aabb_pt[2], lower[0], upper[1], lower[2]); + f3(aabb_pt[3], upper[0], upper[1], lower[2]); + f3(aabb_pt[4], lower[0], lower[1], upper[2]); + f3(aabb_pt[5], upper[0], lower[1], upper[2]); + f3(aabb_pt[6], lower[0], upper[1], upper[2]); + f3(aabb_pt[7], upper[0], upper[1], upper[2]); + + ssp_ran_sphere_uniform(rng, dir); + + /* Build the transformation matrix from world to camera space. Use the AABB + * centroid as the origin of the camera space. */ + build_basis(dir, basis); + build_transform(centroid, basis, transform); + + /* Transform the AABB from world to camera space */ + f3_splat(lower, FLT_MAX); f3_splat(upper,-FLT_MAX); + FOR_EACH(i, 0, 8) { + f3_add(aabb_pt[i], f33_mulf3(aabb_pt[i], transform, aabb_pt[i]), transform+9); + f3_min(lower, lower, aabb_pt[i]); + f3_max(upper, upper, aabb_pt[i]); + } + + /* Ensure an aspect ratio of 1 */ + if(lower[0] > lower[1]) lower[1] = lower[0]; + else lower[0] = lower[1]; + if(upper[0] > upper[1]) upper[1] = upper[0]; + else upper[0] = upper[1]; + + draw_thumbnail(scn, definition, basis, centroid, lower, upper, "pouet.ppm"); + SSP(rng_ref_put(rng)); +} + diff --git a/src/schiff_thumbnail.h b/src/schiff_thumbnail.h @@ -0,0 +1,42 @@ +/* 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. */ + +#ifndef SCHIFF_THUMBNAIL_H +#define SCHIFF_THUMBNAIL_H + +#include <rsys/rsys.h> + +struct s3d_scene; + +extern LOCAL_SYM void +schiff_thumbnail + (struct s3d_scene* scn, + const size_t definition); + +#endif /* SCHIFF_THUMBNAIL_H */ +