commit e64e4ec43361c69ecf845d60d384edf1061b969c
parent 28309351246550297945b6966634b23bb85f22e6
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 11 Apr 2017 10:57:55 +0200
Add analytic shapes to the solver, starting with cylinders.
Diffstat:
9 files changed, 535 insertions(+), 111 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -26,6 +26,7 @@ option(NO_TEST "Do not build tests" OFF)
find_package(RCMake 0.2.3 REQUIRED)
find_package(RSys 0.4 REQUIRED)
find_package(Star3D 0.4 REQUIRED)
+find_package(Star3DUT REQUIRED)
find_package(StarCPR REQUIRED)
find_package(StarSF 0.1 REQUIRED)
find_package(StarSP 0.4 REQUIRED)
@@ -38,11 +39,12 @@ include(rcmake_runtime)
include_directories(
${RSys_INCLUDE_DIR}
${Star3D_INCLUDE_DIR}
+ ${Star3DUT_INCLUDE_DIR}
${StarCPR_INCLUDE_DIR}
${StarSF_INCLUDE_DIR}
${StarSP_INCLUDE_DIR})
-rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarCPR StarSF StarSP)
+rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D Star3DUT StarCPR StarSF StarSP)
################################################################################
# Configure and define targets
@@ -104,7 +106,7 @@ rcmake_prepend_path(SSOL_FILES_INC_API ${SSOL_SOURCE_DIR})
rcmake_prepend_path(SSOL_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(solstice-solver SHARED ${SSOL_FILES_SRC} ${SSOL_FILES_INC} ${SSOL_FILES_INC_API})
-target_link_libraries(solstice-solver RSys Star3D StarCPR StarSF StarSP)
+target_link_libraries(solstice-solver RSys Star3D Star3DUT StarCPR StarSF StarSP)
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(solstice-solver m)
diff --git a/src/ssol.h b/src/ssol.h
@@ -203,7 +203,7 @@ struct ssol_quadric {
/* 3x4 column major transformation of the quadric in object space */
double transform[12];
- /* Hint on the how to discretised */
+ /* Hint on how to discretise */
size_t slices_count_hint;
};
@@ -237,6 +237,41 @@ struct ssol_punched_surface {
static const struct ssol_punched_surface SSOL_PUNCHED_SURFACE_NULL =
SSOL_PUNCHED_SURFACE_NULL__;
+enum ssol_analytic_type {
+ SSOL_ANALYTIC_CYLINDER,
+ SSOL_ANALYTIC_SURFACE_TYPES_COUNT__
+};
+
+struct ssol_analytic_cylinder {
+ double radius;
+ double height;
+ unsigned nslices;
+ unsigned nstacks;
+};
+
+#define SSOL_ANALYTIC_CYLINDER_NULL__ { -1, -1, 0, 0 }
+static const struct ssol_analytic_cylinder SSOL_ANALYTIC_CYLINDER_NULL =
+SSOL_ANALYTIC_CYLINDER_NULL__;
+
+struct ssol_analytic_surface {
+ enum ssol_analytic_type type;
+ union {
+ struct ssol_analytic_cylinder cylinder;
+ } data;
+
+ /* 3x4 column major transformation of the quadric in object space */
+ double transform[12];
+};
+
+#define SSOL_ANALYTIC_SURFACE_NULL__ { \
+ SSOL_ANALYTIC_SURFACE_TYPES_COUNT__, \
+ SSOL_ANALYTIC_CYLINDER_NULL__, \
+ {1,0,0, 0,1,0, 0,0,1, 0,0,0}, \
+}
+
+static const struct ssol_analytic_surface SSOL_ANALYTIC_SURFACE_NULL =
+SSOL_ANALYTIC_SURFACE_NULL__;
+
struct ssol_medium {
double absorptivity;
double refractive_index;
@@ -395,6 +430,14 @@ struct ssol_mc_receiver {
static const struct ssol_mc_receiver SSOL_MC_RECEIVER_NULL =
SSOL_MC_RECEIVER_NULL__;
+#define MC_RCV_NONE__ { \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ 0, NULL, NULL \
+}
+
struct ssol_mc_shape {
/* Internal data */
size_t N__;
@@ -635,6 +678,11 @@ ssol_shape_create_punched_surface
struct ssol_shape** shape);
SSOL_API res_T
+ssol_shape_create_analytic_surface
+(struct ssol_device* dev,
+ struct ssol_shape** shape);
+
+SSOL_API res_T
ssol_shape_ref_get
(struct ssol_shape* shape);
@@ -684,6 +732,11 @@ ssol_mesh_setup
const unsigned nattribs,
void* data);
+SSOL_API res_T
+ssol_analytic_surface_setup
+ (struct ssol_shape* shape,
+ const struct ssol_analytic_surface* analytic_surface);
+
/*******************************************************************************
* Material API - Define the surfacic (e.g.: BRDF) as well as the volumic
* (e.g.: refractive index) properties of a geometry.
diff --git a/src/ssol_draw_draft.c b/src/ssol_draw_draft.c
@@ -65,8 +65,13 @@ Li
/* Retrieve and normalized the hit normal */
switch(sshape->shape->type) {
- case SHAPE_MESH: f3_normalize(N, hit.normal); break;
- case SHAPE_PUNCHED: f3_normalize(N, f3_set_d3(N, ray_data.N)); break;
+ case SHAPE_MESH:
+ f3_normalize(N, hit.normal);
+ break;
+ case SHAPE_PUNCHED:
+ case SHAPE_ANALYTIC:
+ f3_normalize(N, f3_set_d3(N, ray_data.N));
+ break;
default: FATAL("Unreachable code"); break;
}
ASSERT(f3_is_normalized(N));
diff --git a/src/ssol_draw_pt.c b/src/ssol_draw_pt.c
@@ -188,8 +188,13 @@ Li(struct ssol_scene* scn,
/* Retrieve and normalized the hit normal */
switch(sshape->shape->type) {
- case SHAPE_MESH: d3_normalize(N, d3_set_f3(N, hit.normal)); break;
- case SHAPE_PUNCHED: d3_normalize(N, ray_data.N); break;
+ case SHAPE_MESH:
+ d3_normalize(N, d3_set_f3(N, hit.normal));
+ break;
+ case SHAPE_PUNCHED:
+ case SHAPE_ANALYTIC:
+ d3_normalize(N, ray_data.N);
+ break;
default: FATAL("Unreachable code"); break;
}
diff --git a/src/ssol_scene.c b/src/ssol_scene.c
@@ -478,10 +478,10 @@ hit_filter_function
/* Project the hit position into the punched shape */
d3_set_f3(dir, dirf);
d3_set_f3(org, orgf);
- dst = punched_shape_trace_ray(sshape->shape, inst->transform, org, dir,
- hit->distance, N);
+ dst = shape_trace_ray(sshape->shape, inst->transform, org, dir,
+ hit->distance, N, punched_shape_intersect_local);
if(dst >= FLT_MAX) {
- /* No projection is found => the ray does not intersect the quadric */
+ /* No projection found => the ray does not intersect the quadric */
return 1;
}
if((float)dst <= rdata->range_min) {
@@ -499,6 +499,17 @@ hit_filter_function
return 1;
}
break;
+ case SHAPE_ANALYTIC:
+ d3_set_f3(dir, dirf);
+ d3_set_f3(org, orgf);
+ dst = shape_trace_ray(sshape->shape, inst->transform, org, dir,
+ hit->distance, N, analytic_intersect_local);
+ if(dst >= FLT_MAX) {
+ /* No projection found => the ray does not intersect the analytic surface */
+ return 1;
+ }
+ hit_side = d3_dot(dir, N) < 0 ? SSOL_FRONT : SSOL_BACK;
+ break;
default: FATAL("Unreachable code.\n"); break;
}
ASSERT(hit_side == SSOL_BACK || hit_side == SSOL_FRONT);
@@ -513,8 +524,8 @@ hit_filter_function
if((inst->receiver_mask & (int)hit_side) == 0) return 1;
}
- /* Save the nearest intersected quadric point */
- if(sshape->shape->type == SHAPE_PUNCHED && rdata->dst >= dst) {
+ /* Save the nearest intersected quadric/analytic point */
+ if(sshape->shape->type != SHAPE_MESH && rdata->dst >= dst) {
d3_set(rdata->N, N);
rdata->dst = dst;
}
diff --git a/src/ssol_shape.c b/src/ssol_shape.c
@@ -32,10 +32,16 @@
#include <rsys/math.h>
#include <star/scpr.h>
+#include <star/s3dut.h>
#include <limits.h> /* UINT_MAX constant */
#include <math.h> /* copysign function */
+struct mesh_ctx_s3dut {
+ struct s3dut_mesh_data data;
+ double transform[12];
+};
+
struct mesh_context {
const double* coords;
const size_t* ids;
@@ -44,7 +50,8 @@ struct mesh_context {
struct quadric_mesh_context {
const double* coords;
const size_t* ids;
- const union priv_quadric_data* quadric;
+ const union private_data* data;
+ union private_type private_type;
const double* transform; /* 3x4 column major matrix */
};
@@ -57,32 +64,32 @@ struct get_ctx {
/*******************************************************************************
* Helper functions
******************************************************************************/
-static INLINE int
+static FINLINE int
check_plane(const struct ssol_quadric_plane* plane)
{
return plane != NULL;
}
-static INLINE int
+static FINLINE int
check_parabol(const struct ssol_quadric_parabol* parabol)
{
return parabol && parabol->focal > 0;
}
-static INLINE int
+static FINLINE int
check_hyperbol(const struct ssol_quadric_hyperbol* hyperbol)
{
return hyperbol && hyperbol->img_focal > 0 && hyperbol->real_focal > 0;
}
-static INLINE int
+static FINLINE int
check_parabolic_cylinder
(const struct ssol_quadric_parabolic_cylinder* parabolic_cylinder)
{
return parabolic_cylinder && parabolic_cylinder->focal > 0;
}
-static INLINE int
+static FINLINE int
check_hemisphere(const struct ssol_quadric_hemisphere* hemisphere)
{
return hemisphere && hemisphere->radius > 0;
@@ -138,6 +145,28 @@ check_punched_surface(const struct ssol_punched_surface* punched_surface)
return 1;
}
+static FINLINE int
+check_cylinder(const struct ssol_analytic_cylinder* cylinder)
+{
+ return cylinder != NULL
+ && cylinder->radius > 0 && cylinder->height > 0
+ && cylinder->nslices > 2 && cylinder->nstacks > 0;
+}
+
+static INLINE int
+check_analytic_surface(const struct ssol_analytic_surface* analytic_surface)
+{
+ if(!analytic_surface
+ || analytic_surface->type >= SSOL_ANALYTIC_SURFACE_TYPES_COUNT__)
+ return 0;
+
+ switch(analytic_surface->type) {
+ case SSOL_ANALYTIC_CYLINDER:
+ return check_cylinder(&analytic_surface->data.cylinder);
+ default: return 0;
+ }
+}
+
static INLINE int
check_shape(const struct ssol_shape* shape)
{
@@ -255,7 +284,7 @@ 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] = parabol_z(p, &msh->quadric->parabol);
+ p[2] = parabol_z(p, &msh->data->parabol);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -273,7 +302,7 @@ quadric_mesh_hyperbol_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] = hyperbol_z(p, &msh->quadric->hyperbol);
+ p[2] = hyperbol_z(p, &msh->data->hyperbol);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -292,7 +321,7 @@ quadric_mesh_parabolic_cylinder_get_pos
ASSERT(pos && ctx);
p[0] = msh->coords[i+0];
p[1] = msh->coords[i+1];
- p[2] = parabolic_cylinder_z(p, &msh->quadric->pcylinder);
+ p[2] = parabolic_cylinder_z(p, &msh->data->pcylinder);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -311,7 +340,7 @@ quadric_mesh_hemisphere_get_pos
ASSERT(pos && ctx);
p[0] = msh->coords[i + 0];
p[1] = msh->coords[i + 1];
- p[2] = hemisphere_z(p, &msh->quadric->hemisphere);
+ p[2] = hemisphere_z(p, &msh->data->hemisphere);
/* Transform the position in object space */
d33_muld3(p, msh->transform, p);
@@ -578,13 +607,14 @@ 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 = shape->quadric.transform;
+ ctx.transform = shape->transform;
vdata.usage = S3D_POSITION;
vdata.type = S3D_FLOAT3;
vdata.get = NULL;
- ctx.quadric = &shape->priv_quadric;
- switch (shape->quadric.type) {
+ ctx.data = &shape->private_data;
+ ctx.private_type = shape->private_type;
+ switch (shape->private_type.quadric) {
case SSOL_QUADRIC_PARABOL:
vdata.get = quadric_mesh_parabol_get_pos;
break;
@@ -698,36 +728,42 @@ error:
}
/* Solve a 2nd degree equation. "hint" is used to select among the 2 solutions
- * (if applies) the selected solution is then the closest to hint */
+ * (if applies) the selected solution is then the closest to hint ans is
+ * returned in dist[0].
+ * If there is a second solution, it is returned in dist[1].
+ * Returns the number of roots. */
static int
-quadric_solve_second
+solve_second
(const double a,
const double b,
const double c,
const double hint,
- double* dist)
+ double dist[2])
{
- double t = -1;
ASSERT(dist && hint >= 0);
-
if(a == 0) {
- if(b != 0) t = -c / b; /* Degenerated case: 1st degree only */
+ if(b != 0) {
+ dist[0] = -c / b; /* Degenerated case: 1st degree only */
+ return 1;
+ }
+ return 0; /* 0 distance determined */
} else { /* Standard case: 2nd degree */
const double delta = b*b - 4*a*c;
if(delta == 0) {
- t = -b / (2*a);
+ dist[0] = -b / (2*a);
+ return 1;
} else {
const double sqrt_delta = sqrt(delta);
/* Precise formula */
const double t1 = (-b - copysign(sqrt_delta, b)) / (2*a);
const double t2 = c / (a*t1);
- /* Choose the closest value to hint */
- t = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2;
+ /* dist[0] is the closest value to hint */
+ dist[0] = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2;
+ dist[1] = fabs(t1 - hint) < fabs(t2 - hint) ? t2 : t1;
+ return 2;
}
}
- *dist = t;
- return t >= 0;
}
static FINLINE void
@@ -790,6 +826,17 @@ quadric_hemisphere_gradient_local
grad[2] = quad->radius - pt[2];
}
+static FINLINE void
+analytic_cylinder_gradient_local
+ (const double pt[3],
+ double grad[3])
+{
+ ASSERT(pt && grad);
+ grad[0] = pt[0];
+ grad[1] = pt[1];
+ grad[2] = 0;
+}
+
static FINLINE int
quadric_plane_intersect_local
(const double org[3],
@@ -803,13 +850,13 @@ quadric_plane_intersect_local
const double a = 0;
const double b = dir[2];
const double c = org[2];
- double dst;
- int sol = quadric_solve_second(a, b, c, hint, &dst);
+ double dst[2];
+ const int n = solve_second(a, b, c, hint, dst);
- if(!sol) return 0;
- d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst));
+ if(!n) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst));
quadric_plane_gradient_local(grad);
- *dist = dst;
+ *dist = *dst;
return 1;
}
@@ -824,17 +871,17 @@ quadric_parabol_intersect_local
double* dist) /* in/out: */
{
/* Define x^2 + y^2 - 4*focal*z = 0 */
- double dst;
+ double dst[2];
const double a = dir[0] * dir[0] + dir[1] * dir[1];
const double b =
2 * org[0] * dir[0] + 2 * org[1] * dir[1] - 4 * quad->focal * dir[2];
const double c = org[0] * org[0] + org[1] * org[1] - 4 * quad->focal * org[2];
- const int sol = quadric_solve_second(a, b, c, hint, &dst);
+ const int n = solve_second(a, b, c, hint, dst);
- if(!sol) return 0;
- d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst));
+ if(!n) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst));
quadric_parabol_gradient_local(quad, hit_pt, grad);
- *dist = dst;
+ *dist = *dst;
return 1;
}
@@ -848,7 +895,7 @@ quadric_hyperbol_intersect_local
double grad[3],
double* dist)
{
- double dst;
+ double dst[2];
const double b2 = quad->abs_b * quad->abs_b;
const double b2_a2 = b2 * quad->one_over_a_square;
const double z0 = quad->g_square + quad->abs_b;
@@ -859,12 +906,12 @@ quadric_hyperbol_intersect_local
- (org[2] + z0 - quad->g_square) * dir[2]);
const double c = b2_a2 * (org[0] * org[0] + org[1] * org[1]) + b2
- (org[2] + z0 - quad->g_square) * (org[2] + z0 - quad->g_square);
- const int sol = quadric_solve_second(a, b, c, hint, &dst);
+ const int n = solve_second(a, b, c, hint, dst);
- if(!sol) return 0;
- d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst));
+ if(!n) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst));
quadric_hyperbol_gradient_local(quad, hit_pt, grad);
- *dist = dst;
+ *dist = *dst;
return 1;
}
@@ -878,19 +925,19 @@ quadric_hemisphere_intersect_local
double grad[3],
double* dist)
{
- double dst;
+ double dst[2];
double z0 = -quad->radius;
const double a = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2];
const double b = 2 * (org[0] * dir[0] + org[1] * dir[1] + org[2] * dir[2] + z0 * dir[2]);
const double c =
org[0] * org[0] + org[1] * org[1] + org[2] * org[2] - quad->sqr_radius
+ 2 * z0 * org[2] + z0 * z0;
- const int sol = quadric_solve_second(a, b, c, hint, &dst);
+ const int n = solve_second(a, b, c, hint, dst);
- if(!sol) return 0;
- d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst));
+ if(!n) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst));
quadric_hemisphere_gradient_local(quad, hit_pt, grad);
- *dist = dst;
+ *dist = *dst;
return 1;
}
@@ -904,15 +951,49 @@ quadric_parabolic_cylinder_intersect_local
double grad[3],
double* dist)
{
- /* Define y^2 - 4 focal z = 0 */
+ double dst[2];
const double a = dir[1] * dir[1];
const double b = 2 * org[1] * dir[1] - 4 * quad->focal * dir[2];
const double c = org[1] * org[1] - 4 * quad->focal * org[2];
- const int sol = quadric_solve_second(a, b, c, hint, dist);
+ const int n = solve_second(a, b, c, hint, dst);
- if(!sol) return 0;
- d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dist));
+ if(!n) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, *dst));
quadric_parabolic_cylinder_gradient_local(quad, hit_pt, grad);
+ *dist = *dst;
+ return 1;
+}
+
+static FINLINE int
+analytic_cylinder_intersect_local
+ (const struct priv_analytic_cylinder* analytic,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double hit_pt[3],
+ double grad[3],
+ double* dist) /* in/out: */
+{
+ double dst[2];
+ const double a = dir[0] * dir[0] + dir[1] * dir[1];
+ const double b = 2 * (org[0] * dir[0] + org[1] * dir[1]);
+ const double c = org[0] * org[0] + org[1] * org[1]
+ - analytic->radius * analytic->radius;
+ const int n = solve_second(a, b, c, hint, dst);
+
+ if(!n) return 0;
+ /* We just solved intersection with an infinite cylinder.
+ * Now we filter according to cylinder's height. */
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst[0]));
+ if(fabs(hit_pt[2]) > 0.5 * analytic->height) {
+ if(n == 1) return 0;
+ d3_add(hit_pt, org, d3_muld(hit_pt, dir, dst[1]));
+ if(fabs(hit_pt[2]) > 0.5 * analytic->height) return 0;
+ *dist = dst[1];
+ } else {
+ *dist = dst[0];
+ }
+ analytic_cylinder_gradient_local(hit_pt, grad);
return 1;
}
@@ -921,21 +1002,21 @@ punched_shape_set_z_local(const struct ssol_shape* shape, double pt[3])
{
ASSERT(shape && pt);
ASSERT(shape->type == SHAPE_PUNCHED);
- switch (shape->quadric.type) {
+ switch (shape->private_type.quadric) {
case SSOL_QUADRIC_PLANE:
pt[2] = 0;
break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
- pt[2] = parabolic_cylinder_z(pt, &shape->priv_quadric.pcylinder);
+ pt[2] = parabolic_cylinder_z(pt, &shape->private_data.pcylinder);
break;
case SSOL_QUADRIC_PARABOL:
- pt[2] = parabol_z(pt, &shape->priv_quadric.parabol);
+ pt[2] = parabol_z(pt, &shape->private_data.parabol);
break;
case SSOL_QUADRIC_HYPERBOL:
- pt[2] = hyperbol_z(pt, &shape->priv_quadric.hyperbol);
+ pt[2] = hyperbol_z(pt, &shape->private_data.hyperbol);
break;
case SSOL_QUADRIC_HEMISPHERE:
- pt[2] = hemisphere_z(pt, &shape->priv_quadric.hemisphere);
+ pt[2] = hemisphere_z(pt, &shape->private_data.hemisphere);
break;
default: FATAL("Unreachable code\n"); break;
}
@@ -949,31 +1030,31 @@ punched_shape_set_normal_local
{
ASSERT(shape && pt);
ASSERT(shape->type == SHAPE_PUNCHED);
- switch (shape->quadric.type) {
+ switch (shape->private_type.quadric) {
case SSOL_QUADRIC_PLANE:
quadric_plane_gradient_local(normal);
break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
quadric_parabolic_cylinder_gradient_local
- (&shape->priv_quadric.pcylinder, pt, normal);
+ (&shape->private_data.pcylinder, pt, normal);
break;
case SSOL_QUADRIC_PARABOL:
quadric_parabol_gradient_local
- (&shape->priv_quadric.parabol, pt, normal);
+ (&shape->private_data.parabol, pt, normal);
break;
case SSOL_QUADRIC_HYPERBOL:
quadric_hyperbol_gradient_local
- (&shape->priv_quadric.hyperbol, pt, normal);
+ (&shape->private_data.hyperbol, pt, normal);
break;
case SSOL_QUADRIC_HEMISPHERE:
quadric_hemisphere_gradient_local
- (&shape->priv_quadric.hemisphere, pt, normal);
+ (&shape->private_data.hemisphere, pt, normal);
break;
default: FATAL("Unreachable code\n"); break;
}
}
-static FINLINE int
+int
punched_shape_intersect_local
(const struct ssol_shape* shape,
const double org[3],
@@ -989,31 +1070,56 @@ punched_shape_intersect_local
ASSERT(dir[0] || dir[1] || dir[2]);
/* Hits on quadrics must be recomputed more accurately */
- switch (shape->quadric.type) {
+ switch (shape->private_type.quadric) {
case SSOL_QUADRIC_PLANE:
hit = quadric_plane_intersect_local(org, dir, hint, pt, N, dist);
break;
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
hit = quadric_parabolic_cylinder_intersect_local
- (&shape->priv_quadric.pcylinder, org, dir, hint, pt, N, dist);
+ (&shape->private_data.pcylinder, org, dir, hint, pt, N, dist);
break;
case SSOL_QUADRIC_PARABOL:
hit = quadric_parabol_intersect_local
- (&shape->priv_quadric.parabol, org, dir, hint, pt, N, dist);
+ (&shape->private_data.parabol, org, dir, hint, pt, N, dist);
break;
case SSOL_QUADRIC_HYPERBOL:
hit = quadric_hyperbol_intersect_local
- (&shape->priv_quadric.hyperbol, org, dir, hint, pt, N, dist);
+ (&shape->private_data.hyperbol, org, dir, hint, pt, N, dist);
break;
case SSOL_QUADRIC_HEMISPHERE:
hit = quadric_hemisphere_intersect_local
- (&shape->priv_quadric.hemisphere, org, dir, hint, pt, N, dist);
+ (&shape->private_data.hemisphere, org, dir, hint, pt, N, dist);
break;
default: FATAL("Unreachable code\n"); break;
}
return hit;
}
+int
+analytic_intersect_local
+ (const struct ssol_shape* shape,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double pt[3],
+ double N[3],
+ double* dist)
+{
+ int hit;
+ ASSERT(shape && org && dir && hint >= 0 && pt && N && dist);
+ ASSERT(shape->type == SHAPE_ANALYTIC);
+ ASSERT(dir[0] || dir[1] || dir[2]);
+
+ switch (shape->private_type.analytic) {
+ case SSOL_ANALYTIC_CYLINDER:
+ hit = analytic_cylinder_intersect_local
+ (&shape->private_data.cylinder, org, dir, hint, pt, N, dist);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return hit;
+}
+
static void
shape_release(ref_T* ref)
{
@@ -1080,7 +1186,7 @@ priv_hemisphere_data_setup
static INLINE void
priv_quadric_data_setup
- (union priv_quadric_data* priv_data,
+ (union private_data* priv_data,
const struct ssol_quadric* quadric)
{
ASSERT(priv_data && quadric);
@@ -1109,7 +1215,7 @@ priv_quadric_data_setup
static INLINE size_t
priv_quadric_data_compute_slices_count
(const enum ssol_quadric_type type,
- const union priv_quadric_data* priv_data,
+ const union private_data* priv_data,
const double lower[2],
const double upper[2])
{
@@ -1143,11 +1249,11 @@ priv_quadric_data_compute_slices_count
double r;
d2_set(u, upper);
r = d2_dot(u, u);
- if (r > priv_data->hemisphere.sqr_radius)
+ if(r > priv_data->hemisphere.sqr_radius)
d2_muld(u, u, sqrt(priv_data->hemisphere.sqr_radius / r));
d2_set(l, lower);
r = d2_dot(l, l);
- if (r > priv_data->hemisphere.sqr_radius)
+ if(r > priv_data->hemisphere.sqr_radius)
d2_muld(l, l, sqrt(priv_data->hemisphere.sqr_radius / r));
max_z = MMAX
(hemisphere_z(l, &priv_data->hemisphere),
@@ -1170,6 +1276,120 @@ get_circular(const size_t ivert, double position[2], void* ctx)
position[1] = data->radius * sin(a);
}
+static void
+mesh_ctx_s3dut_get_ids(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ const struct mesh_ctx_s3dut* mesh = ctx;
+ ASSERT(ids && ctx && itri < mesh->data.nprimitives);
+ ASSERT(mesh->data.indices[itri * 3 + 0] <= UINT_MAX);
+ ASSERT(mesh->data.indices[itri * 3 + 1] <= UINT_MAX);
+ ASSERT(mesh->data.indices[itri * 3 + 2] <= UINT_MAX);
+ ids[0] = (unsigned)mesh->data.indices[itri * 3 + 0];
+ ids[1] = (unsigned)mesh->data.indices[itri * 3 + 1];
+ ids[2] = (unsigned)mesh->data.indices[itri * 3 + 2];
+}
+
+static void
+mesh_ctx_s3dut_get_pos(const unsigned ivert, float pos[3], void* ctx)
+{
+ const struct mesh_ctx_s3dut* mesh = ctx;
+ double tmp[3];
+ ASSERT(pos && ctx && ivert < mesh->data.nvertices);
+ d3_set(tmp, mesh->data.positions + ivert * 3);
+ d33_muld3(tmp, mesh->transform, tmp);
+ d3_add(tmp, mesh->transform + 9, tmp);
+ f3_set_d3(pos, tmp);
+}
+
+/* Setup the Star-3D shape of the analytic surface to ray-trace */
+static res_T
+analytic_setup_s3d_shape_rt
+ (struct ssol_shape* shape,
+ double* rt_area)
+{
+ struct s3dut_mesh* mesh = NULL;
+ struct mesh_ctx_s3dut mesh_ctx;
+ struct s3d_vertex_data attrs;
+ res_T res;
+
+ switch(shape->private_type.analytic) {
+ case SSOL_ANALYTIC_CYLINDER:
+ ASSERT(shape->private_data.cylinder.nslices < UINT_MAX);
+ res = s3dut_create_cylinder(
+ shape->dev->allocator,
+ shape->private_data.cylinder.radius,
+ shape->private_data.cylinder.height,
+ shape->private_data.cylinder.nslices,
+ shape->private_data.cylinder.nstacks,
+ &mesh);
+ if(res != RES_OK) {
+ fprintf(stderr, "Could not create the cylinder 3D data.\n");
+ goto error;
+ }
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ S3DUT(mesh_get_data(mesh, &mesh_ctx.data));
+ ASSERT(mesh_ctx.data.nprimitives <= UINT_MAX);
+ ASSERT(mesh_ctx.data.nvertices <= UINT_MAX);
+ d33_set(mesh_ctx.transform, shape->transform);
+ d3_set(mesh_ctx.transform + 9, shape->transform + 9);
+
+ attrs.usage = SSOL_TO_S3D_POSITION;
+ attrs.type = S3D_FLOAT3;
+ attrs.get = mesh_ctx_s3dut_get_pos;
+
+ res = s3d_mesh_setup_indexed_vertices
+ (shape->shape_rt, (unsigned)mesh_ctx.data.nprimitives,
+ mesh_ctx_s3dut_get_ids, (unsigned)mesh_ctx.data.nvertices,
+ &attrs, 1, &mesh_ctx);
+ if(res != RES_OK) {
+ fprintf(stderr, "Could not setup a Solstice Solver mesh shape.\n");
+ goto error;
+ }
+
+ /* FIXME: use mesh area or analytic area??? */
+ *rt_area = mesh_compute_area
+ ((unsigned)mesh_ctx.data.nprimitives, mesh_ctx_s3dut_get_ids,
+ (unsigned)mesh_ctx.data.nvertices, attrs.get, &mesh_ctx);
+ *rt_area = 2 * PI *
+ shape->private_data.cylinder.radius *
+ shape->private_data.cylinder.height;
+end:
+ if(mesh) s3dut_mesh_ref_put(mesh);
+ return RES_OK;
+error:
+ goto end;
+}
+
+static INLINE void
+priv_analytic_cylinder_data_setup
+ (struct priv_analytic_cylinder* priv_data,
+ const struct ssol_analytic_cylinder* analytic)
+{
+ ASSERT(priv_data && analytic);
+ priv_data->radius = analytic->radius;
+ priv_data->height = analytic->height;
+ priv_data->nslices = analytic->nslices;
+ priv_data->nstacks = analytic->nstacks;
+}
+
+static INLINE void
+priv_analytic_data_setup
+ (union private_data* priv_data,
+ const struct ssol_analytic_surface* analytic)
+{
+ ASSERT(priv_data && analytic);
+ switch (analytic->type) {
+ case SSOL_ANALYTIC_CYLINDER:
+ priv_analytic_cylinder_data_setup
+ (&priv_data->cylinder, &analytic->data.cylinder);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -1191,8 +1411,8 @@ punched_shape_project_point
ASSERT(shape->type == SHAPE_PUNCHED);
/* Compute world<->quadric space transformations */
- d33_muld33(R, transform, shape->quadric.transform);
- d33_muld3(T, transform, shape->quadric.transform+9);
+ d33_muld33(R, transform, shape->transform);
+ d33_muld3(T, transform, shape->transform+9);
d3_add(T, T, transform + 9);
d33_invtrans(R_invtrans, R);
d3_minus(T_inv, T);
@@ -1215,13 +1435,14 @@ punched_shape_project_point
}
double
-punched_shape_trace_ray
+shape_trace_ray
(struct ssol_shape* shape,
const double transform[12], /* Shape to world space transformation */
const double org[3], /* World space position near of the ray origin */
const double dir[3], /* World space ray direction */
const double hint_dst, /* Hint on the hit distance */
- double N_quadric[3]) /* World space normal onto the quadric */
+ double N_shape[3], /* World space normal onto the shape */
+ intersect_local_fn local) /* the intersection function for this shape */
{
double R[9]; /* Quadric to world rotation matrix */
double R_invtrans[9]; /* Inverse transpose of R */
@@ -1233,12 +1454,11 @@ punched_shape_trace_ray
double N_local[3];
double dst; /* Hit distance */
int valid;
- ASSERT(shape && transform && org && N_quadric);
- ASSERT(shape->type == SHAPE_PUNCHED);
+ ASSERT(shape && transform && org && shape);
/* Compute world<->quadric space transformations */
- d33_muld33(R, transform, shape->quadric.transform);
- d33_muld3(T, transform, shape->quadric.transform+9);
+ d33_muld33(R, transform, shape->transform);
+ d33_muld3(T, transform, shape->transform+9);
d3_add(T, T, transform + 9);
d33_invtrans(R_invtrans, R);
d3_minus(T_inv, T);
@@ -1250,14 +1470,14 @@ punched_shape_trace_ray
/* Transform dir in quadric space */
d3_muld33(dir_local, dir, R_invtrans);
- /* Project pos_local onto the quadric and compute its associated normal */
- valid = punched_shape_intersect_local
+ /* Project pos_local onto the shape and compute its associated normal */
+ valid = 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);
+ /* Transform the shape normal in world space */
+ d33_muld3(N_shape, R_invtrans, N_local);
+ d3_normalize(N_shape, N_shape);
return dst;
}
@@ -1310,6 +1530,14 @@ ssol_shape_create_punched_surface
}
res_T
+ssol_shape_create_analytic_surface
+ (struct ssol_device* dev,
+ struct ssol_shape** out_shape)
+{
+ return shape_create(dev, out_shape, SHAPE_ANALYTIC);
+}
+
+res_T
ssol_shape_ref_get(struct ssol_shape* shape)
{
if(!shape) return RES_BAD_ARG;
@@ -1350,11 +1578,11 @@ ssol_shape_get_vertex_attrib
/* Transform the fetch attrib */
if(shape->type == SHAPE_PUNCHED) {
if(usage == SSOL_POSITION) {
- d33_muld3(value, shape->quadric.transform, value);
- d3_add(value, shape->quadric.transform + 9, value);
+ d33_muld3(value, shape->transform, value);
+ d3_add(value, shape->transform + 9, value);
} else if(usage == SSOL_NORMAL) {
double R_invtrans[9];
- d33_invtrans(R_invtrans, shape->quadric.transform);
+ d33_invtrans(R_invtrans, shape->transform);
d33_muld3(value, R_invtrans, value);
}
}
@@ -1399,7 +1627,9 @@ ssol_punched_surface_setup
}
/* Save quadric for further object instancing */
- shape->quadric = *psurf->quadric;
+ d33_set(shape->transform, psurf->quadric->transform);
+ d3_set(shape->transform+9, psurf->quadric->transform+9);
+ shape->private_type.quadric = psurf->quadric->type;
if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) {
size_t c;
@@ -1420,7 +1650,7 @@ ssol_punched_surface_setup
tmp_carvings[psurf->nb_carvings].nb_vertices = 1024;
ctx.nbvert = tmp_carvings[psurf->nb_carvings].nb_vertices;
ctx.two_pi_over_nbvert = 2 * PI / (double)ctx.nbvert;
- ctx.radius = shape->quadric.data.hemisphere.radius;
+ ctx.radius = psurf->quadric->data.hemisphere.radius;
tmp_carvings_count = psurf->nb_carvings + 1;
}
else {
@@ -1438,14 +1668,14 @@ ssol_punched_surface_setup
}
/* Setup internal data */
- priv_quadric_data_setup(&shape->priv_quadric, psurf->quadric);
+ priv_quadric_data_setup(&shape->private_data, psurf->quadric);
/* Define the #slices of the discretized quadric */
if(psurf->quadric->slices_count_hint != SIZE_MAX) {
nslices = psurf->quadric->slices_count_hint;
} else {
nslices = priv_quadric_data_compute_slices_count
- (shape->quadric.type, &shape->priv_quadric, lower, upper);
+ (shape->private_type.quadric, &shape->private_data, lower, upper);
}
res = build_triangulated_plane(&coords, &ids, lower, upper, nslices);
@@ -1535,7 +1765,7 @@ ssol_mesh_setup
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 */
+ /* The Star-3D shape to sample is the same that the one to ray-trace */
res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp);
if(res != RES_OK) goto error;
shape->shape_samp_area = shape->shape_rt_area;
@@ -1546,3 +1776,38 @@ error:
goto exit;
}
+res_T
+ssol_analytic_surface_setup
+ (struct ssol_shape* shape,
+ const struct ssol_analytic_surface* asurf)
+{
+ res_T res = RES_OK;
+
+ if(!check_shape(shape)
+ || !check_analytic_surface(asurf)
+ || shape->type != SHAPE_ANALYTIC) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Save analytic surface for further object instancing */
+ d33_set(shape->transform, asurf->transform);
+ d3_set(shape->transform + 9, asurf->transform + 9);
+ shape->private_type.analytic = asurf->type;
+
+ /* Setup internal data */
+ priv_analytic_data_setup(&shape->private_data, asurf);
+
+ /* Setup the Star-3D shape to ray-trace */
+ res = analytic_setup_s3d_shape_rt(shape, &shape->shape_rt_area);
+ if(res != RES_OK) goto error;
+
+ /* The Star-3D shape to sample is the same that the one to ray-trace */
+ res = s3d_mesh_copy(shape->shape_rt, shape->shape_samp);
+ if(res != RES_OK) goto error;
+ shape->shape_samp_area = shape->shape_rt_area;
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/ssol_shape_c.h b/src/ssol_shape_c.h
@@ -23,6 +23,7 @@
enum shape_type {
SHAPE_MESH,
SHAPE_PUNCHED,
+ SHAPE_ANALYTIC,
SHAPE_TYPES_COUNT__
};
@@ -48,11 +49,24 @@ struct priv_hemisphere_data {
double sqr_radius;
};
-union priv_quadric_data {
+struct priv_analytic_cylinder {
+ double radius;
+ double height;
+ unsigned nslices;
+ unsigned nstacks;
+};
+
+union private_data {
struct priv_hyperbol_data hyperbol;
struct priv_parabol_data parabol;
struct priv_pcylinder_data pcylinder;
struct priv_hemisphere_data hemisphere;
+ struct priv_analytic_cylinder cylinder;
+};
+
+union private_type {
+ enum ssol_quadric_type quadric;
+ enum ssol_analytic_type analytic;
};
struct ssol_shape {
@@ -60,14 +74,24 @@ struct ssol_shape {
struct s3d_shape* shape_rt; /* Star-3D shape to ray-trace */
struct s3d_shape* shape_samp; /* Star-3D shape to sample */
- union priv_quadric_data priv_quadric;
- struct ssol_quadric quadric;
+ union private_data private_data;
+ union private_type private_type;
+ double transform[12];
double shape_rt_area, shape_samp_area;
struct ssol_device* dev;
ref_T ref;
};
+typedef int(*intersect_local_fn)
+ (const struct ssol_shape* shape,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double pt[3],
+ double N[3],
+ double* dist);
+
/* Project pos onto the punched surface and retrieve its associated normal */
extern LOCAL_SYM void
punched_shape_project_point
@@ -97,5 +121,36 @@ shape_fetched_raw_vertex_attrib
const enum ssol_attrib_usage usage,
double value[]);
+/* Compute ray/punched shape intersection */
+extern LOCAL_SYM int punched_shape_intersect_local
+ (const struct ssol_shape* shape,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double pt[3],
+ double N[3],
+ double* dist);
+
+/* Compute ray/analytic shape intersection */
+extern LOCAL_SYM int analytic_intersect_local
+ (const struct ssol_shape* shape,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double pt[3],
+ double N[3],
+ double* dist);
+
+/* Compute ray/shape intersection */
+extern LOCAL_SYM double
+shape_trace_ray
+ (struct ssol_shape* shape,
+ const double transform[12], /* Shape to world space transformation */
+ const double org[3], /* World space position near of the ray origin */
+ const double dir[3], /* World space ray direction */
+ const double hint_dst, /* Hint on the hit distance */
+ double N_quadric[3], /* World space normal onto the quadric */
+ intersect_local_fn local);
+
#endif /* SSOL_SHAPE_C_H */
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -327,6 +327,7 @@ point_update_from_hit
d3_set_f3(pt->pos, tmpf);
break;
case SHAPE_PUNCHED:
+ case SHAPE_ANALYTIC:
d3_normalize(pt->N, rdata->N);
d3_muld(tmp, pt->dir, rdata->dst);
d3_add(pt->pos, pt->pos, tmp);
diff --git a/src/test_ssol_shape.c b/src/test_ssol_shape.c
@@ -36,6 +36,7 @@ main(int argc, char** argv)
-1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 0.f, -2.f
};
const size_t npolygon_verts = sizeof(polygon)/sizeof(double[2]);
+ struct ssol_analytic_surface analytic = SSOL_ANALYTIC_SURFACE_NULL__;
double val[3];
unsigned ids[3];
unsigned i, n;
@@ -151,13 +152,7 @@ main(int argc, char** argv)
CHECK(ssol_shape_create_punched_surface(dev, NULL), RES_BAD_ARG);
CHECK(ssol_shape_create_punched_surface(NULL, &shape), RES_BAD_ARG);
CHECK(ssol_shape_create_punched_surface(dev, &shape), RES_OK);
-
- CHECK(ssol_shape_ref_get(NULL), RES_BAD_ARG);
- CHECK(ssol_shape_ref_get(shape), RES_OK);
-
- CHECK(ssol_shape_ref_put(NULL), RES_BAD_ARG);
- CHECK(ssol_shape_ref_put(shape), RES_OK);
-
+
carving.get = get_polygon_vertices;
carving.operation = SSOL_AND;
carving.nb_vertices = npolygon_verts;
@@ -230,6 +225,38 @@ main(int argc, char** argv)
quadric.data.hemisphere.radius = 10;
CHECK(ssol_shape_ref_put(shape), RES_OK);
+
+ CHECK(ssol_shape_create_analytic_surface(NULL, &shape), RES_BAD_ARG);
+ CHECK(ssol_shape_create_analytic_surface(dev, NULL), RES_BAD_ARG);
+ CHECK(ssol_shape_create_analytic_surface(dev, &shape), RES_OK);
+
+ analytic.type = SSOL_ANALYTIC_CYLINDER;
+ analytic.data.cylinder.height = 10;
+ analytic.data.cylinder.radius = 10;
+ analytic.data.cylinder.nslices = 10;
+ analytic.data.cylinder.nstacks = 10;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_OK);
+
+ analytic.data.cylinder.height = 0;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_BAD_ARG);
+ analytic.data.cylinder.height = 10;
+
+ analytic.data.cylinder.radius = 0;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_BAD_ARG);
+ analytic.data.cylinder.radius = 10;
+
+ analytic.data.cylinder.nslices = 0;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_BAD_ARG);
+ analytic.data.cylinder.nslices = 10;
+
+ analytic.data.cylinder.nstacks = 0;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_BAD_ARG);
+ analytic.data.cylinder.nstacks = 10;
+
+ analytic.type = (enum ssol_analytic_type)999;
+ CHECK(ssol_analytic_surface_setup(shape, &analytic), RES_BAD_ARG);
+
+ CHECK(ssol_shape_ref_put(shape), RES_OK);
CHECK(ssol_device_ref_put(dev), RES_OK);
check_memory_allocator(&allocator);
