commit 1a1eb57276d2d57455b34cde392dc6013935f95a
parent 1b508a1927856bc55a8aa51411881b58c30ee52f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 9 Oct 2015 15:15:19 +0200
Begin the estimation of the radiative properties
Diffstat:
2 files changed, 189 insertions(+), 32 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -58,12 +58,19 @@ struct ssp_rng;
struct sschiff_refractive_index { double real, imag; };
+enum sschiff_material_property {
+ SSCHIFF_MATERIAL_MEDIUM_REFRACTIVE_INDEX,
+ SSCHIFF_MATERIAL_RELATIVE_REAL_REFRACTIVE_INDEX,
+ SSCHIFF_MATERIAL_RELATIVE_IMAGINARY_REFRACTIVE_INDEX
+};
+
struct sschiff_material {
/* Retrieve the refractive complex index of the material */
- void (*get_refractive_index)
+ void (*get_property)
(void* material,
const double wavelength,
- struct sschiff_refractive_index* value);
+ const enum sschiff_material_property property,
+ double* value);
void* material; /* Opaque user defined representation of the material */
};
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -26,6 +26,8 @@
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL license and that you accept its terms. */
+#define _POSIX_C_SOURCE 200112L
+
#include "sschiff.h"
#include "sschiff_device.h"
@@ -33,12 +35,15 @@
#include <rsys/float3.h>
#include <rsys/image.h>
#include <rsys/logger.h>
+#include <rsys/math.h>
#include <rsys/mem_allocator.h>
#include <rsys/stretchy_array.h>
#include <star/s3d.h>
#include <star/ssp.h>
+#include <math.h>
+
#define NDIRS 1 /* # Sampled directions */
struct sschiff_estimator {
@@ -49,6 +54,12 @@ struct sschiff_estimator {
ref_T ref;
};
+struct accum {
+ double weight;
+ double sqr_weight;
+ size_t naccums;
+};
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -99,7 +110,7 @@ build_transform
{
ASSERT(pos && basis && transform);
- /* The linear part of the transformamatrix is the inverse of the
+ /* 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]);
@@ -111,8 +122,130 @@ build_transform
transform[11] = -f3_dot(basis+6, pos);
}
+static res_T
+compute_path_length
+ (struct s3d_scene* scn,
+ const struct s3d_hit* first_hit,
+ const float ray_org[3],
+ const float ray_dir[3],
+ double* length)
+{
+ float range[2] = { 0, FLT_MAX };
+ struct s3d_primitive prim_from;
+ struct s3d_hit hit;
+ double len = 0;
+ ASSERT(scn && first_hit && ray_org && ray_dir && !S3D_HIT_NONE(first_hit));
+
+ range[0] = first_hit->distance;
+ prim_from = first_hit->prim;
+ do {
+ do {
+ range[0] = nextafterf(range[0], range[1]);
+ S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
+ } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from)); /* Self-intersection */
+
+ if(S3D_HIT_NONE(&hit)) return RES_BAD_ARG;
+
+ len += hit.distance;
+ range[0] = hit.distance;
+ prim_from = hit.prim;
+
+ do {
+ range[0] = nextafterf(range[0], range[1]);
+ S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
+ } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from)); /* Self-intersection */
+
+ range[0] = hit.distance;
+ prim_from = hit.prim;
+
+ } while(!S3D_HIT_NONE(&hit));
+
+ *length = len * 1.e-6/* Micron to meter */;
+ return RES_OK;
+}
+
+static res_T
+compute_monte_carlo_weight
+ (struct s3d_scene* scn,
+ struct ssp_rng* rng,
+ struct sschiff_material* mtl,
+ const double wavelength,
+ 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 */
+ struct accum* accum)
+{
+ struct s3d_hit hit;
+ double sample[2];
+ double weight;
+ 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];
+ float x[3], y[3];
+ float rcp_pdf;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && mtl && basis && pos && lower && upper && accum);
+
+ f2_sub(plane_size, upper, lower);
+ rcp_pdf = plane_size[0] * plane_size[1];
+
+ /* Define the projection axis */
+ 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]));
+
+ /* Uniformly sample a position onto the projection plane and use it as the
+ * origin of the ray to trace */
+ sample[0] = ssp_rng_uniform_double(rng, -1.0, 1.0);
+ sample[1] = ssp_rng_uniform_double(rng, -1.0, 1.0);
+ f3_mulf(x, axis_x, (float)sample[0]);
+ f3_mulf(y, axis_y, (float)sample[1]);
+ f3_add(ray_org, f3_add(ray_org, x, y), org);
+
+ S3D(scene_trace_ray(scn, ray_org, axis_z, range, &hit));
+ if(S3D_HIT_NONE(&hit)) {
+ weight = 0.0;
+ } else {
+ double n_r, k_r;
+ double n_e, k_e, lambda_e;
+ double path_length;
+
+ res = compute_path_length(scn, &hit, ray_org, axis_z, &path_length);
+ if(res != RES_OK) goto error;
+
+ /* Retrieve material properties */
+ mtl->get_property(mtl->material, wavelength,
+ SSCHIFF_MATERIAL_MEDIUM_REFRACTIVE_INDEX, &n_e);
+ mtl->get_property(mtl->material, wavelength,
+ SSCHIFF_MATERIAL_RELATIVE_REAL_REFRACTIVE_INDEX, &n_r);
+ mtl->get_property(mtl->material, wavelength,
+ SSCHIFF_MATERIAL_RELATIVE_IMAGINARY_REFRACTIVE_INDEX, &k_r);
+
+ lambda_e = wavelength / n_e;
+ k_e = 2.0 * PI / lambda_e;
+
+ weight = 2.0 * rcp_pdf
+ * (1.0 - exp(-k_e*k_r*path_length) * cos(k_e*(n_r - 1.0)*path_length));
+ }
+
+ accum->weight += weight;
+ accum->sqr_weight += weight * weight;
+ accum->naccums += 1;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/* Cast rays into the RT volume */
-static void
+static INLINE void
draw_thumbnail
(struct s3d_scene* scn,
const float basis[9], /* World space basis of the RT volume */
@@ -186,18 +319,23 @@ draw_thumbnail
}
static res_T
-radiative_path_length
+radiative_properties
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_integrator_desc* desc,
- const int istep)
+ const double wavelength,
+ const int istep,
+ struct accum* accum)
{
struct s3d_shape* shape = NULL;
struct s3d_scene* scene = NULL;
struct sschiff_material material = SSCHIFF_NULL_MATERIAL;
+ float lower[3], upper[3];
+ float aabb_pt[8][3];
+ float centroid[3];
int idir;
res_T res = RES_OK;
- ASSERT(dev && rng && desc);
+ ASSERT(dev && rng && desc && accum);
res = s3d_shape_create_mesh(dev->s3d, &shape);
if(res != RES_OK) {
@@ -223,44 +361,43 @@ radiative_path_length
}
S3D(scene_begin_session(scene, S3D_TRACE));
+ S3D(scene_get_aabb(scene, 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]);
+
+ f3_mulf(centroid, f3_add(centroid, lower, upper), 0.5f);
+
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[12];
float pt[8][3];
int i;
ssp_ran_sphere_uniform(rng, dir);
- 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);
build_basis(dir, basis);
build_transform(centroid, basis, transform);
/* Transform the AABB from world to footprint space */
FOR_EACH(i, 0, 8) {
- f3_add(pt[i], f33_mulf3(pt[i], transform, pt[i]), transform + 9);
+ f3_add(pt[i], f33_mulf3(pt[i], transform, aabb_pt[i]), transform + 9);
}
/* Compute the AABB of the footprint */
@@ -270,9 +407,20 @@ radiative_path_length
f3_max(upper, upper, pt[i]);
}
+#if 0
+ (void)istep;
+ res = compute_monte_carlo_weight
+ (scene, rng, &material, wavelength, basis, centroid, lower, upper, accum);
+ if(res != RES_OK) goto error;
+#else
/* Compute an image of the shape transformed in the footprint space */
- sprintf(thumb_name, "%d_%d.ppm", istep, idir);
- draw_thumbnail(scene, basis, centroid, lower, upper, thumb_name);
+ {
+ char thumb_name[64];
+ (void)wavelength;
+ sprintf(thumb_name, "%d_%d.ppm", istep, idir);
+ draw_thumbnail(scene, basis, centroid, lower, upper, thumb_name);
+ }
+#endif
}
S3D(scene_end_session(scene));
@@ -347,6 +495,7 @@ sschiff_integrate
struct sschiff_estimator** out_estimator)
{
struct sschiff_estimator* estimator = NULL;
+ struct accum accum;
int i;
res_T res = RES_OK;
@@ -361,8 +510,9 @@ sschiff_integrate
estimator->nsteps = steps_count;
estimator->wavelength = wavelength;
+ memset(&accum, 0, sizeof(struct accum));
FOR_EACH(i, 0, (int)steps_count) {
- res = radiative_path_length(dev, rng, desc, i);
+ res = radiative_properties(dev, rng, desc, wavelength, i, &accum);
if(res != RES_OK) goto error;
}
