commit 1a9645c524411ace416ae5ea63b585cc32cb09c0
parent b3e48bafd23ff9eef9e329e4eb87a51b54358031
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 9 Mar 2017 15:54:39 +0100
Merge remote-tracking branch 'origin/develop' into feature-path_tracing
Diffstat:
6 files changed, 654 insertions(+), 45 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -168,6 +168,7 @@ if(NOT NO_TEST)
new_test(test_ssol_solver4)
new_test(test_ssol_solver5)
new_test(test_ssol_solver6)
+ new_test(test_ssol_solver7)
new_test(test_ssol_sun)
build_test(test_ssol_draw)
diff --git a/src/ssol.h b/src/ssol.h
@@ -93,6 +93,7 @@ enum ssol_parametrization_type {
enum ssol_quadric_type {
SSOL_QUADRIC_PLANE,
SSOL_QUADRIC_PARABOL,
+ SSOL_QUADRIC_HYPERBOL,
SSOL_QUADRIC_PARABOLIC_CYLINDER,
SSOL_QUADRIC_TYPE_COUNT__
};
@@ -157,6 +158,15 @@ struct ssol_quadric_parabol {
static const struct ssol_quadric_parabol SSOL_QUADRIC_PARABOL_NULL =
SSOL_QUADRIC_PARABOL_NULL__;
+struct ssol_quadric_hyperbol {
+ double img_focal, real_focal;
+ /* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0
+ * with a^2 = f - f^2; b = f -1/2; f = real_focal / (img_focal + real_focal) */
+};
+#define SSOL_QUADRIC_HYPERBOL_NULL__ { -1.0 , -1.0 }
+static const struct ssol_quadric_hyperbol SSOL_QUADRIC_HYPERBOL_NULL =
+SSOL_QUADRIC_HYPERBOL_NULL__;
+
struct ssol_quadric_parabolic_cylinder {
double focal; /* Define y^2 - 4 focal z = 0 */
};
@@ -169,6 +179,7 @@ struct ssol_quadric {
union {
struct ssol_quadric_plane plane;
struct ssol_quadric_parabol parabol;
+ struct ssol_quadric_hyperbol hyperbol;
struct ssol_quadric_parabolic_cylinder parabolic_cylinder;
} data;
diff --git a/src/ssol_shape.c b/src/ssol_shape.c
@@ -44,7 +44,7 @@ struct mesh_context {
struct quadric_mesh_context {
const double* coords;
const size_t* ids;
- double focal; /* Use by parabol and parabolic cylinder quadrics */
+ const struct priv_quadric_data* quadric;
const double* transform; /* 3x4 column major matrix */
};
@@ -64,6 +64,12 @@ check_parabol(const struct ssol_quadric_parabol* parabol)
}
static INLINE int
+check_hyperbol(const struct ssol_quadric_hyperbol* hyperbol)
+{
+ return hyperbol && hyperbol->img_focal > 0 && hyperbol->real_focal > 0;
+}
+
+static INLINE int
check_parabolic_cylinder
(const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder)
{
@@ -80,6 +86,8 @@ check_quadric(const struct ssol_quadric* quadric)
return check_plane(&quadric->data.plane);
case SSOL_QUADRIC_PARABOL:
return check_parabol(&quadric->data.parabol);
+ case SSOL_QUADRIC_HYPERBOL:
+ return check_hyperbol(&quadric->data.hyperbol);
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
return check_parabolic_cylinder(&quadric->data.parabolic_cylinder);
default: return 0;
@@ -183,6 +191,33 @@ quadric_mesh_plane_get_pos(const unsigned ivert, float pos[3], void* ctx)
f3_set_d3(pos, p);
}
+static FINLINE double
+hyperbol_z
+ (const double p[2],
+ const struct priv_hyperbol_data* hyperbol)
+{
+ const double z0 = hyperbol->g_2 + hyperbol->abs_b;
+ const double r2 = p[0] * p[0] + p[1] * p[1];
+ return hyperbol->abs_b * sqrt(1 + r2 * hyperbol->_1_a2) + hyperbol->g_2 - z0;
+}
+
+static FINLINE double
+parabol_z
+ (const double p[2],
+ const struct priv_parabol_data* parabol)
+{
+ const double r2 = p[0] * p[0] + p[1] * p[1];
+ return r2 * parabol->_1_4f;
+}
+
+static FINLINE double
+parabolic_cylinder_z
+ (const double p[2],
+ const struct priv_pcylinder_data* pcyl)
+{
+ return (p[1] * p[1]) * pcyl->_1_4f;
+}
+
static void
quadric_mesh_parabol_get_pos(const unsigned ivert, float pos[3], void* ctx)
{
@@ -192,7 +227,25 @@ quadric_mesh_parabol_get_pos(const unsigned ivert, float pos[3], void* ctx)
ASSERT(pos && ctx);
p[0] = msh->coords[i+0];
p[1] = msh->coords[i+1];
- p[2] = (p[0]*p[0] + p[1]*p[1]) / (4.0*msh->focal);
+ p[2] = parabol_z(p, &msh->quadric->data.parabol);
+
+ /* Transform the position in object space */
+ d33_muld3(p, msh->transform, p);
+ d3_add(p, p, msh->transform+9);
+
+ f3_set_d3(pos, p);
+}
+
+static void
+quadric_mesh_hyperbol_get_pos(const unsigned ivert, float pos[3], void* ctx)
+{
+ const size_t i = ivert * 2/*#coords per vertex*/;
+ const struct quadric_mesh_context* msh = ctx;
+ double p[3]; /* Temporary quadric space position */
+ ASSERT(pos && ctx);
+ p[0] = msh->coords[i+0];
+ p[1] = msh->coords[i+1];
+ p[2] = hyperbol_z(p, &msh->quadric->data.hyperbol);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -211,7 +264,7 @@ quadric_mesh_parabolic_cylinder_get_pos
ASSERT(pos && ctx);
p[0] = msh->coords[i+0];
p[1] = msh->coords[i+1];
- p[2] = ((p[1]*p[1]) / (4.0*msh->focal));
+ p[2] = parabolic_cylinder_z(p, &msh->quadric->data.pcylinder);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -457,10 +510,10 @@ mesh_compute_area
* quadric equation */
static res_T
quadric_setup_s3d_shape_rt
- (const struct ssol_quadric* quadric,
+ (const struct ssol_shape* shape,
const struct darray_double* coords,
const struct darray_size_t* ids,
- struct s3d_shape* shape,
+ struct s3d_shape* s3dshape,
double* rt_area)
{
struct quadric_mesh_context ctx;
@@ -468,7 +521,7 @@ quadric_setup_s3d_shape_rt
unsigned nverts;
unsigned ntris;
res_T res;
- ASSERT(quadric && coords && ids && shape);
+ ASSERT(shape && coords && ids && s3dshape && rt_area);
ASSERT(darray_double_size_get(coords)%2 == 0);
ASSERT(darray_size_t_size_get(ids)%3 == 0);
ASSERT(darray_double_size_get(coords)/2 <= UINT_MAX);
@@ -478,17 +531,20 @@ quadric_setup_s3d_shape_rt
ntris = (unsigned)darray_size_t_size_get(ids) / 3/*#ids per triangle*/;
ctx.coords = darray_double_cdata_get(coords);
ctx.ids = darray_size_t_cdata_get(ids);
- ctx.transform = quadric->transform;
+ ctx.transform = shape->quadric.transform;
vdata.usage = S3D_POSITION;
vdata.type = S3D_FLOAT3;
- switch(quadric->type) {
+ vdata.get = NULL;
+ ctx.quadric = &shape->priv_quadric;
+ switch (shape->quadric.type) {
case SSOL_QUADRIC_PARABOL:
- ctx.focal = quadric->data.parabol.focal;
vdata.get = quadric_mesh_parabol_get_pos;
break;
+ case SSOL_QUADRIC_HYPERBOL:
+ vdata.get = quadric_mesh_hyperbol_get_pos;
+ break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
- ctx.focal = quadric->data.parabolic_cylinder.focal;
vdata.get = quadric_mesh_parabolic_cylinder_get_pos;
break;
case SSOL_QUADRIC_PLANE:
@@ -498,9 +554,10 @@ quadric_setup_s3d_shape_rt
}
res = s3d_mesh_setup_indexed_vertices
- (shape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx);
+ (s3dshape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx);
if(res != RES_OK) return res;
+ ASSERT(vdata.get);
*rt_area = mesh_compute_area
(ntris, quadric_mesh_get_ids, nverts, vdata.get, &ctx);
return RES_OK;
@@ -536,9 +593,9 @@ quadric_setup_s3d_shape_samp
vdata.usage = S3D_POSITION;
vdata.type = S3D_FLOAT3;
vdata.get = quadric_mesh_plane_get_pos;
- res = s3d_mesh_setup_indexed_vertices
+ res = s3d_mesh_setup_indexed_vertices
(shape, ntris, quadric_mesh_get_ids, nverts, &vdata, 1, &ctx);
- if (res != RES_OK) return res;
+ if(res != RES_OK) return res;
*samp_area = mesh_compute_area
(ntris, quadric_mesh_get_ids, nverts, quadric_mesh_plane_get_pos, &ctx);
return RES_OK;
@@ -634,7 +691,7 @@ quadric_plane_gradient_local(double grad[3])
static FINLINE void
quadric_parabol_gradient_local
- (const struct ssol_quadric_parabol* quad,
+ (const struct priv_parabol_data* quad,
const double pt[3],
double grad[3])
{
@@ -645,8 +702,23 @@ quadric_parabol_gradient_local
}
static FINLINE void
+quadric_hyperbol_gradient_local
+ (const struct priv_hyperbol_data* quad,
+ const double pt[3],
+ double grad[3])
+{
+ ASSERT(quad && pt && grad);
+ {
+ const double z0 = quad->g_2 + quad->abs_b;
+ grad[0] = pt[0];
+ grad[1] = pt[1];
+ grad[2] = -(pt[2] + z0 - quad->g_2) * quad->_a2_b2;
+ }
+}
+
+static FINLINE void
quadric_parabolic_cylinder_gradient_local
- (const struct ssol_quadric_parabolic_cylinder* quad,
+ (const struct priv_pcylinder_data* quad,
const double pt[3],
double grad[3])
{
@@ -681,7 +753,7 @@ quadric_plane_intersect_local
static FINLINE int
quadric_parabol_intersect_local
- (const struct ssol_quadric_parabol* quad,
+ (const struct priv_parabol_data* quad,
const double org[3],
const double dir[3],
const double hint,
@@ -705,8 +777,37 @@ quadric_parabol_intersect_local
}
static FINLINE int
+quadric_hyperbol_intersect_local
+ (const struct priv_hyperbol_data* quad,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double hit_pt[3],
+ double grad[3],
+ double* dist)
+{
+ double dst;
+ const double b2 = quad->abs_b * quad->abs_b;
+ const double b2_a2 = b2 * quad->_1_a2;
+ const double z0 = quad->g_2 + quad->abs_b;
+ const double a =
+ b2_a2 * (dir[0] * dir[0] + dir[1] * dir[1]) - dir[2] * dir[2];
+ const double b =
+ 2 * (b2_a2 * (org[0] * dir[0] + org[1] * dir[1]) - (org[2] + z0 - quad->g_2) * dir[2]);
+ const double c = b2_a2 * (org[0] * org[0] + org[1] * org[1]) + b2
+ - (org[2] + z0 - quad->g_2) * (org[2] + z0 - quad->g_2);
+ const int sol = quadric_solve_second(a, b, c, hint, &dst);
+
+ if(!sol) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst));
+ quadric_hyperbol_gradient_local(quad, hit_pt, grad);
+ *dist = dst;
+ return 1;
+}
+
+static FINLINE int
quadric_parabolic_cylinder_intersect_local
- (const struct ssol_quadric_parabolic_cylinder* quad,
+ (const struct priv_pcylinder_data* quad,
const double org[3],
const double dir[3],
const double hint,
@@ -735,17 +836,15 @@ punched_shape_set_z_local(const struct ssol_shape* shape, double pt[3])
case SSOL_QUADRIC_PLANE:
pt[2] = 0;
break;
- case SSOL_QUADRIC_PARABOLIC_CYLINDER: {
- const struct ssol_quadric_parabolic_cylinder* quad
- = &shape->quadric.data.parabolic_cylinder;
- pt[2] = (pt[1] * pt[1]) / (4.0 * quad->focal);
+ case SSOL_QUADRIC_PARABOLIC_CYLINDER:
+ pt[2] = parabolic_cylinder_z(pt, &shape->priv_quadric.data.pcylinder);
break;
- }
- case SSOL_QUADRIC_PARABOL: {
- const struct ssol_quadric_parabol* quad = &shape->quadric.data.parabol;
- pt[2] = (pt[0] * pt[0] + pt[1] * pt[1]) / (4.0 * quad->focal);
+ case SSOL_QUADRIC_PARABOL:
+ pt[2] = parabol_z(pt, &shape->priv_quadric.data.parabol);
+ break;
+ case SSOL_QUADRIC_HYPERBOL:
+ pt[2] = hyperbol_z(pt, &shape->priv_quadric.data.hyperbol);
break;
- }
default: FATAL("Unreachable code\n"); break;
}
}
@@ -764,11 +863,15 @@ punched_shape_set_normal_local
break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
quadric_parabolic_cylinder_gradient_local
- (&shape->quadric.data.parabolic_cylinder, pt, normal);
+ (&shape->priv_quadric.data.pcylinder, pt, normal);
break;
case SSOL_QUADRIC_PARABOL: {
quadric_parabol_gradient_local
- (&shape->quadric.data.parabol, pt, normal);
+ (&shape->priv_quadric.data.parabol, pt, normal);
+ break;
+ case SSOL_QUADRIC_HYPERBOL:
+ quadric_hyperbol_gradient_local
+ (&shape->priv_quadric.data.hyperbol, pt, normal);
break;
}
default: FATAL("Unreachable code\n"); break;
@@ -797,11 +900,15 @@ punched_shape_intersect_local
break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
hit = quadric_parabolic_cylinder_intersect_local
- (&shape->quadric.data.parabolic_cylinder, org, dir, hint, pt, N, dist);
+ (&shape->priv_quadric.data.pcylinder, org, dir, hint, pt, N, dist);
break;
case SSOL_QUADRIC_PARABOL:
hit = quadric_parabol_intersect_local
- (&shape->quadric.data.parabol, org, dir, hint, pt, N, dist);
+ (&shape->priv_quadric.data.parabol, org, dir, hint, pt, N, dist);
+ break;
+ case SSOL_QUADRIC_HYPERBOL:
+ hit = quadric_hyperbol_intersect_local
+ (&shape->priv_quadric.data.hyperbol, org, dir, hint, pt, N, dist);
break;
default: FATAL("Unreachable code\n"); break;
}
@@ -906,7 +1013,7 @@ punched_shape_trace_ray
valid = punched_shape_intersect_local
(shape, org_local, dir_local, hint_dst, hit_local, N_local, &dst);
if(!valid) return INF;
-
+
/* Transform the quadric normal in world space */
d33_muld3(N_quadric, R_invtrans, N_local);
d3_normalize(N_quadric, N_quadric);
@@ -1067,21 +1174,43 @@ ssol_punched_surface_setup
nslices = 1;
break;
case SSOL_QUADRIC_PARABOL: {
- double z[2];
- z[0] = (lower[0] * lower[0] + lower[1] * lower[1])
- / (4.0 * psurf->quadric->data.parabol.focal);
- z[1] = (upper[0] * upper[0] + upper[1] * upper[1])
- / (4.0 * psurf->quadric->data.parabol.focal);
- nslices = MMIN(50, (size_t)(3 + sqrt(MMAX(z[0], z[1])) * 6));
+ const struct ssol_quadric_parabol* parabol
+ = &psurf->quadric->data.parabol;
+ struct priv_parabol_data* data = &shape->priv_quadric.data.parabol;
+ double max_z;
+ data->focal = parabol->focal;
+ data->_1_4f = 1 / (4.0 * parabol->focal);
+ max_z = MMAX(parabol_z(lower, data), parabol_z(upper, data));
+ nslices = MMIN(50, (size_t) (3 + sqrt(max_z) * 6));
+ break;
+ }
+ case SSOL_QUADRIC_HYPERBOL: {
+ const struct ssol_quadric_hyperbol* hyperbol =
+ &psurf->quadric->data.hyperbol;
+ struct priv_hyperbol_data* data = &shape->priv_quadric.data.hyperbol;
+ /* re-dimensionalize */
+ const double g = hyperbol->real_focal + hyperbol->img_focal;
+ const double f = hyperbol->real_focal / g;
+ const double a2 = g * g * (f - f * f);
+ double max_z;
+ data->g_2 = g * 0.5;
+ data->abs_b = g * fabs(f - 0.5);
+ data->_a2_b2 = a2 / (data->abs_b * data->abs_b);
+ data->_1_a2 = 1 / a2;
+ max_z = MMAX(hyperbol_z(lower, data), hyperbol_z(upper, data));
+ nslices = MMIN(50, (size_t) (3 + sqrt(max_z) * 6));
break;
}
case SSOL_QUADRIC_PARABOLIC_CYLINDER: {
- double z[2];
- z[0] = (lower[1] * lower[1]) /
- (4.0 * psurf->quadric->data.parabolic_cylinder.focal);
- z[1] = (upper[1] * upper[1]) /
- (4.0 * psurf->quadric->data.parabolic_cylinder.focal);
- nslices = MMIN(50, (size_t)(3 + sqrt(MMAX(z[0], z[1])) * 6));
+ const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder
+ = &psurf->quadric->data.parabolic_cylinder;
+ struct priv_pcylinder_data* data = &shape->priv_quadric.data.pcylinder;
+ double max_z;
+ data->focal = psurf->quadric->data.parabolic_cylinder.focal;
+ data->_1_4f = 1 / (4.0 * parabolic_cylinder->focal);
+ max_z = MMAX(parabolic_cylinder_z(lower, data),
+ parabolic_cylinder_z(upper, data));
+ nslices = MMIN(50, (size_t) (3 + sqrt(max_z) * 6));
break;
}
default: FATAL("Unreachable code\n"); break;
@@ -1096,7 +1225,7 @@ ssol_punched_surface_setup
/* Setup the Star-3D shape to ray-trace */
res = quadric_setup_s3d_shape_rt
- (psurf->quadric, &coords, &ids, shape->shape_rt, &shape->shape_rt_area);
+ (shape, &coords, &ids, shape->shape_rt, &shape->shape_rt_area);
if(res != RES_OK) goto error;
/* Setup the Star-3D shape to sample */
@@ -1168,7 +1297,8 @@ ssol_mesh_setup
res = s3d_mesh_setup_indexed_vertices
(shape->shape_rt, ntris, get_indices, nverts, attrs, nattribs, data);
if(res != RES_OK) goto error;
- shape->shape_rt_area = mesh_compute_area(ntris, get_indices, nverts, get_position, data);
+ shape->shape_rt_area =
+ mesh_compute_area(ntris, get_indices, nverts, get_position, data);
/* The Star-3D shape to sample is the same of the one to ray-traced */
res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp);
diff --git a/src/ssol_shape_c.h b/src/ssol_shape_c.h
@@ -26,11 +26,38 @@ enum shape_type {
SHAPE_TYPES_COUNT__
};
+struct priv_parabol_data {
+ double focal;
+ double _1_4f;
+};
+
+struct priv_hyperbol_data {
+ double g_2;
+ double _a2_b2;
+ double _1_a2;
+ double abs_b;
+};
+
+struct priv_pcylinder_data {
+ double focal;
+ double _1_4f;
+};
+
+struct priv_quadric_data {
+ enum ssol_quadric_type type;
+ union {
+ struct priv_hyperbol_data hyperbol;
+ struct priv_parabol_data parabol;
+ struct priv_pcylinder_data pcylinder;
+ } data;
+};
+
struct ssol_shape {
enum shape_type type;
struct s3d_shape* shape_rt; /* Star-3D shape to ray-trace */
struct s3d_shape* shape_samp; /* Star-3D shape to sample */
+ struct priv_quadric_data priv_quadric;
struct ssol_quadric quadric;
double shape_rt_area, shape_samp_area;
diff --git a/src/test_ssol_shape.c b/src/test_ssol_shape.c
@@ -200,6 +200,19 @@ main(int argc, char** argv)
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
quadric.data.parabol.focal = 1;
+ quadric.type = SSOL_QUADRIC_HYPERBOL;
+ quadric.data.hyperbol.real_focal = 1;
+ quadric.data.hyperbol.img_focal = 1;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+
+ quadric.data.hyperbol.real_focal = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ quadric.data.hyperbol.real_focal = 1;
+
+ quadric.data.hyperbol.img_focal = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ quadric.data.hyperbol.img_focal = 1;
+
quadric.type = SSOL_QUADRIC_PARABOLIC_CYLINDER;
quadric.data.parabolic_cylinder.focal = 1;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
diff --git a/src/test_ssol_solver7.c b/src/test_ssol_solver7.c
@@ -0,0 +1,427 @@
+/* Copyright (C) CNRS 2016-2017
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "ssol.h"
+#include "test_ssol_utils.h"
+
+#define REFLECTIVITY 0
+#include "test_ssol_materials.h"
+
+#define TARGET_SZ 2
+#define PLANE_NAME TARGET
+#define HALF_X (TARGET_SZ/2)
+#define HALF_Y (TARGET_SZ/2)
+STATIC_ASSERT((HALF_X * 2 == TARGET_SZ), ONLY_ENVEN_VALUES_FOR_SZ);
+#include "test_ssol_rect_geometry.h"
+
+#define HELIOSTAT_SZ 4
+#define POLYGON_NAME HELIOSTAT
+#define HALF_X (HELIOSTAT_SZ/2)
+#define HALF_Y (HELIOSTAT_SZ/2)
+STATIC_ASSERT(HALF_X * 2 == HELIOSTAT_SZ, ONLY_ENVEN_VALUES_FOR_SZ);
+#include "test_ssol_rect2D_geometry.h"
+
+#define HYPERBOL_SZ 24
+#define POLYGON_NAME HYPERBOL
+#define HALF_X (HYPERBOL_SZ/2)
+#define HALF_Y (HYPERBOL_SZ/2)
+STATIC_ASSERT(HALF_X * 2 == HYPERBOL_SZ, ONLY_ENVEN_VALUES_FOR_SZ);
+#include "test_ssol_rect2D_geometry.h"
+
+#include <rsys/double33.h>
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+static void
+get_wlen(const size_t i, double* wlen, double* data, void* ctx)
+{
+ double wavelengths[3] = { 1, 2, 3 };
+ double intensities[3] = { 1, 0.8, 1 };
+ CHECK(i < 3, 1);
+ (void)ctx;
+ *wlen = wavelengths[i];
+ *data = intensities[i];
+}
+
+#include <rsys/mem_allocator.h>
+#include <rsys/image.h>
+#include "ssol_scene_c.h"
+#include "ssol_device_c.h"
+
+#define SCREEN_GAMMA 2.2
+
+#define HEIGHT 800
+#define WIDTH 1000
+
+/* Assume that the pixel format of the src is DOUBLE3 in gray scale while the
+* pixel format of dst is UBYTE */
+static void
+tone_map(const double* src, unsigned char* dst, const size_t count)
+{
+ size_t i;
+ ASSERT(src && dst && count);
+ FOR_EACH(i, 0, count) {
+ double val;
+ val = pow(src[i * 3/*#channels*/], 1 / SCREEN_GAMMA);/* Gamma correction */
+ val = CLAMP(val, 0, 1);
+ dst[i] = (unsigned char) ((val * 255) + 0.5/*round*/);
+ }
+}
+
+static void
+tone_map_image(const struct ssol_image* img, unsigned char* dst)
+{
+ struct ssol_image_layout layout;
+ size_t irow = 0;
+ void* mem;
+ ASSERT(img && dst);
+
+ SSOL(image_get_layout(img, &layout));
+ ASSERT(layout.pixel_format == SSOL_PIXEL_DOUBLE3);
+
+ SSOL(image_map(img, &mem));
+ FOR_EACH(irow, 0, layout.height) {
+ const void* src_row = ((char*) mem) + layout.offset + irow * layout.row_pitch;
+ unsigned char* dst_row = dst + irow * layout.width;
+ tone_map(src_row, dst_row, layout.width);
+ }
+}
+
+static res_T
+setup_framebuffer(struct ssol_device* dev, struct ssol_image** framebuffer)
+{
+ struct ssol_image* fbuf = NULL;
+ res_T res = RES_OK;
+ ASSERT(dev && framebuffer);
+
+ res = ssol_image_create(dev, &fbuf);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not create the rendering framebuffer.\n");
+ goto error;
+ }
+
+ res = ssol_image_setup(fbuf, WIDTH, HEIGHT, SSOL_PIXEL_DOUBLE3);
+ if (res != RES_OK) {
+ fprintf(stderr,
+ "Could not set the framebuffer definition to %dx%d.\n",
+ WIDTH, HEIGHT);
+ goto error;
+ }
+
+exit:
+ *framebuffer = fbuf;
+ return res;
+
+error:
+ if (fbuf) {
+ SSOL(image_ref_put(fbuf));
+ fbuf = NULL;
+ }
+ goto exit;
+}
+
+static res_T
+setup_camera(struct ssol_device* dev, struct ssol_camera** camera)
+{
+ struct ssol_camera* cam = NULL;
+ double proj_ratio = 0;
+ res_T res = RES_OK;
+ double pos[3] = { 50, -120, 50 }, tgt[3] = { 50, 0, 50 }, up[3] = { 0, 0, 1 };
+ ASSERT(dev && camera);
+
+ res = ssol_camera_create(dev, &cam);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not create the rendering camera.\n");
+ goto error;
+ }
+
+#define FOV 60.0
+
+ proj_ratio = WIDTH / HEIGHT;
+ res = ssol_camera_set_proj_ratio(cam, proj_ratio);
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid image ratio '%g'.\n", proj_ratio);
+ goto error;
+ }
+
+ res = ssol_camera_set_fov(cam, MDEG2RAD(FOV));
+ if (res != RES_OK) {
+ fprintf(stderr, "Invalid horizontal field of view '%g' degrees.\n", FOV);
+ goto error;
+ }
+
+ res = ssol_camera_look_at(cam, pos, tgt, up);
+ if (res != RES_OK) {
+ fprintf(stderr,
+ "Invalid camera point of view:\n"
+ " position = %g %g %g\n"
+ " target = %g %g %g\n"
+ " up = %g %g %g\n",
+ SPLIT3(pos),
+ SPLIT3(tgt),
+ SPLIT3(up));
+ goto error;
+ }
+
+exit:
+ *camera = cam;
+ return res;
+error:
+ if (cam) {
+ SSOL(camera_ref_put(cam));
+ cam = NULL;
+ }
+ goto exit;
+}
+
+/* TODO Remove this dead code or move and refactor it in the test utilities */
+static INLINE res_T
+solstice_draw(struct ssol_scene* scene, const char* name)
+{
+ struct ssol_image_layout layout;
+ unsigned char* ubytes = NULL;
+ struct ssol_image* framebuffer;
+ struct ssol_camera* camera;
+ FILE* output; /* Output stream */
+ res_T res = RES_OK;
+ ASSERT(scene && name);
+
+ output = fopen(name, "w");
+ if (!output) return RES_IO_ERR;
+
+ res = setup_framebuffer(scene->dev, &framebuffer);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not setup the framebuffer.\n");
+ goto error;
+ }
+
+ SSOL(image_get_layout(framebuffer, &layout));
+ ubytes = MEM_ALLOC(scene->dev->allocator, layout.width*layout.height);
+ if (!ubytes) {
+ fprintf(stderr, "Could not allocate the 8-bits image buffer.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ res = setup_camera(scene->dev, &camera);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not setup the camera.\n");
+ goto error;
+ }
+
+ res = ssol_draw(scene, camera, layout.width,
+ layout.height, ssol_image_write, framebuffer);
+ if (res != RES_OK) {
+ fprintf(stderr, "Rendering error\n");
+ goto error;
+ }
+
+ tone_map_image(framebuffer, ubytes);
+ res = image_ppm_write_stream(output, (int) layout.width, (int) layout.height, 1, ubytes);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not write the rendered image to the output stream.\n");
+ goto error;
+ }
+
+exit:
+ if (output) fclose(output);
+ if (ubytes) MEM_RM(scene->dev->allocator, ubytes);
+ return res;
+error:
+ goto exit;
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct ssol_device* dev;
+ struct ssp_rng* rng;
+ struct ssol_scene* scene;
+ struct ssol_shape* square;
+ struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ };
+ struct ssol_carving carving1 = SSOL_CARVING_NULL;
+ struct ssol_carving carving2 = SSOL_CARVING_NULL;
+ struct ssol_material *m_mtl, *bck_mtl, *v_mtl;
+ struct ssol_mirror_shader m_shader = SSOL_MIRROR_SHADER_NULL;
+ struct ssol_matte_shader bck_shader = SSOL_MATTE_SHADER_NULL;
+ struct ssol_object* t_object;
+ struct ssol_instance* target;
+ struct ssol_sun* sun;
+ struct ssol_spectrum* spectrum;
+ struct ssol_estimator* estimator;
+ struct ssol_mc_global mc_global;
+ struct ssol_mc_receiver mc_rcv;
+ double dir[3];
+ double transform[12]; /* 3x4 column major matrix */
+ FILE* tmp;
+ /* primary is a parabol */
+ struct ssol_quadric quadric1 = SSOL_QUADRIC_DEFAULT;
+ struct ssol_punched_surface punched1 = SSOL_PUNCHED_SURFACE_NULL;
+ struct ssol_object* m_object1;
+ struct ssol_shape* quad_square1;
+ struct ssol_instance* heliostat;
+ /* secondary is an hyperbol */
+ struct ssol_quadric quadric2 = SSOL_QUADRIC_DEFAULT;
+ struct ssol_punched_surface punched2 = SSOL_PUNCHED_SURFACE_NULL;
+ struct ssol_object* m_object2;
+ struct ssol_shape* quad_square2;
+ struct ssol_instance* secondary;
+ (void) argc, (void) argv;
+
+#define H 10
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssol_device_create
+ (NULL, &allocator, 1, 0, &dev), RES_OK);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+ CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK);
+ CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK);
+ CHECK(ssol_sun_create_directional(dev, &sun), RES_OK);
+ CHECK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)), RES_OK);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
+ CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+ CHECK(ssol_scene_create(dev, &scene), RES_OK);
+ CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
+
+ /* create scene content */
+
+ CHECK(ssol_shape_create_mesh(dev, &square), RES_OK);
+ attribs[0].usage = SSOL_POSITION;
+ attribs[0].get = get_position;
+ CHECK(ssol_mesh_setup(square, TARGET_NTRIS__, get_ids,
+ TARGET_NVERTS__, attribs, 1, (void*) &TARGET_DESC__), RES_OK);
+
+ CHECK(ssol_material_create_mirror(dev, &m_mtl), RES_OK);
+ m_shader.normal = get_shader_normal;
+ m_shader.reflectivity = get_shader_reflectivity;
+ m_shader.roughness = get_shader_roughness;
+ CHECK(ssol_mirror_set_shader(m_mtl, &m_shader), RES_OK);
+ CHECK(ssol_material_create_matte(dev, &bck_mtl), RES_OK);
+ bck_shader.normal = get_shader_normal;
+ bck_shader.reflectivity = get_shader_reflectivity_2;
+ CHECK(ssol_matte_set_shader(bck_mtl, &bck_shader), RES_OK);
+ CHECK(ssol_material_create_virtual(dev, &v_mtl), RES_OK);
+
+ carving1.get = get_polygon_vertices;
+ carving1.operation = SSOL_AND;
+ carving1.nb_vertices = HELIOSTAT_NVERTS__;
+ carving1.context = &HELIOSTAT_EDGES__;
+
+ CHECK(ssol_shape_create_punched_surface(dev, &quad_square1), RES_OK);
+ quadric1.type = SSOL_QUADRIC_PARABOL;
+ quadric1.data.parabol.focal = 10 * sqrt(2) * H;
+ punched1.nb_carvings = 1;
+ punched1.quadric = &quadric1;
+ punched1.carvings = &carving1;
+ CHECK(ssol_punched_surface_setup(quad_square1, &punched1), RES_OK);
+ CHECK(ssol_object_create(dev, &m_object1), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(m_object1, quad_square1, m_mtl, m_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(m_object1, &heliostat), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, heliostat), RES_OK);
+ d33_rotation_yaw(transform, -0.25 * PI);
+ d3_splat(transform + 9, 0);
+ transform[9] = 10 * H; /* target the img focal point of the hyperbol */
+ CHECK(ssol_instance_set_transform(heliostat, transform), RES_OK);
+
+ carving2.get = get_polygon_vertices;
+ carving2.operation = SSOL_AND;
+ carving2.nb_vertices = HYPERBOL_NVERTS__;
+ carving2.context = &HYPERBOL_EDGES__;
+
+ CHECK(ssol_shape_create_punched_surface(dev, &quad_square2), RES_OK);
+ quadric2.type = SSOL_QUADRIC_HYPERBOL;
+ quadric2.data.hyperbol.real_focal = 9 * H;
+ quadric2.data.hyperbol.img_focal = 1 * H;
+ punched2.nb_carvings = 1;
+ punched2.quadric = &quadric2;
+ punched2.carvings = &carving2;
+ CHECK(ssol_punched_surface_setup(quad_square2, &punched2), RES_OK);
+ CHECK(ssol_object_create(dev, &m_object2), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(m_object2, quad_square2, m_mtl, v_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(m_object2, &secondary), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, secondary), RES_OK);
+ d33_set_identity(transform);
+ d3_splat(transform + 9, 0);
+ transform[11] = 9 * H;
+ CHECK(ssol_instance_set_transform(secondary, transform), RES_OK);
+ CHECK(ssol_instance_sample(secondary, 0), RES_OK);
+
+ CHECK(ssol_object_create(dev, &t_object), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(t_object, square, bck_mtl, v_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(t_object, &target), RES_OK);
+ CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 0), RES_OK);
+ CHECK(ssol_instance_sample(target, 0), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
+ d33_set_identity(transform);
+ d3_splat(transform + 9, 0);
+ CHECK(ssol_instance_set_transform(target, transform), RES_OK);
+
+#define N__ 10000
+#define DNI_cos (1000 * cos(0))
+#define TOTAL (HELIOSTAT_SZ * HELIOSTAT_SZ * DNI_cos)
+#define GET_MC_RCV ssol_estimator_get_mc_receiver
+
+ NCHECK(tmp = tmpfile(), 0);
+ CHECK(ssol_solve(scene, rng, N__, tmp, &estimator), RES_OK);
+ CHECK(fclose(tmp), 0);
+
+ /*CHECK(solstice_draw(scene, "full_scene.ppm"), RES_OK);*/
+
+ printf("Total = %g\n", TOTAL);
+ CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
+ printf("Shadows = %g +/- %g\n",
+ mc_global.shadowed.E, mc_global.shadowed.SE);
+ CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
+ printf("Missing = %g +/- %g\n",
+ mc_global.missing.E, mc_global.missing.SE);
+ CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1);
+
+ CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
+ printf("Ir(target1) = %g +/- %g\n",
+ mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
+ CHECK(eq_eps(mc_rcv.integrated_irradiance.E, TOTAL, TOTAL * 1e-4), 1);
+ CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, 0, 1e-4), 1);
+
+ /* Free data */
+ CHECK(ssol_instance_ref_put(heliostat), RES_OK);
+ CHECK(ssol_instance_ref_put(secondary), RES_OK);
+ CHECK(ssol_instance_ref_put(target), RES_OK);
+ CHECK(ssol_object_ref_put(m_object1), RES_OK);
+ CHECK(ssol_object_ref_put(m_object2), RES_OK);
+ CHECK(ssol_object_ref_put(t_object), RES_OK);
+ CHECK(ssol_shape_ref_put(square), RES_OK);
+ CHECK(ssol_shape_ref_put(quad_square1), RES_OK);
+ CHECK(ssol_shape_ref_put(quad_square2), RES_OK);
+ CHECK(ssol_material_ref_put(m_mtl), RES_OK);
+ CHECK(ssol_material_ref_put(bck_mtl), RES_OK);
+ CHECK(ssol_material_ref_put(v_mtl), RES_OK);
+ CHECK(ssol_estimator_ref_put(estimator), RES_OK);
+ CHECK(ssol_device_ref_put(dev), RES_OK);
+ CHECK(ssol_scene_ref_put(scene), RES_OK);
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+ CHECK(ssol_spectrum_ref_put(spectrum), RES_OK);
+ CHECK(ssol_sun_ref_put(sun), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+
+ return 0;
+}
+
