commit 1732a933072aea3d2557ee88bc38f172f36b7e90
parent 4fbad2c17c7074c1ba25b5c0f41a3b65c6d9abd3
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 21 Apr 2017 15:23:49 +0200
Merge remote-tracking branch 'origin/feature_new_shapes' into develop
Diffstat:
8 files changed, 485 insertions(+), 118 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 0.2 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
@@ -107,6 +107,7 @@ enum ssol_quadric_type {
SSOL_QUADRIC_PARABOL,
SSOL_QUADRIC_HYPERBOL,
SSOL_QUADRIC_PARABOLIC_CYLINDER,
+ SSOL_QUADRIC_HEMISPHERE,
SSOL_QUADRIC_TYPE_COUNT__
};
@@ -171,7 +172,7 @@ static const struct ssol_quadric_parabol SSOL_QUADRIC_PARABOL_NULL =
SSOL_QUADRIC_PARABOL_NULL__;
struct ssol_quadric_hyperbol {
- /* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0
+ /* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0; z > 0
* with a^2 = f - f^2; b = f -1/2; f = real_focal / (img_focal + real_focal) */
double img_focal, real_focal;
};
@@ -186,6 +187,14 @@ struct ssol_quadric_parabolic_cylinder {
static const struct ssol_quadric_parabolic_cylinder
SSOL_QUADRIC_PARABOLIC_CYLINDER_NULL = SSOL_QUADRIC_PARABOLIC_CYLINDER_NULL__;
+struct ssol_quadric_hemisphere {
+ /* Define x^2 + y^2 + (z-radius)^2 - radius^2 = 0 with z <= r */
+ double radius;
+};
+#define SSOL_QUADRIC_HEMISPHERE_NULL__ { -1.0 }
+static const struct ssol_quadric_hemisphere SSOL_QUADRIC_HEMISPHERE_NULL =
+SSOL_QUADRIC_HEMISPHERE_NULL__;
+
struct ssol_quadric {
enum ssol_quadric_type type;
union {
@@ -193,12 +202,13 @@ struct ssol_quadric {
struct ssol_quadric_parabol parabol;
struct ssol_quadric_hyperbol hyperbol;
struct ssol_quadric_parabolic_cylinder parabolic_cylinder;
+ struct ssol_quadric_hemisphere hemisphere;
} data;
/* 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;
};
@@ -401,6 +411,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__;
diff --git a/src/ssol_draw_draft.c b/src/ssol_draw_draft.c
@@ -72,6 +72,7 @@ Li
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;
+ break;
default: FATAL("Unreachable code"); break;
}
diff --git a/src/ssol_draw_pt.c b/src/ssol_draw_pt.c
@@ -188,8 +188,12 @@ 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:
+ 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) {
@@ -514,7 +514,7 @@ hit_filter_function
}
/* Save the nearest intersected quadric point */
- if(sshape->shape->type == SHAPE_PUNCHED && rdata->dst >= dst) {
+ 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,40 +50,53 @@ 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 */
double lower[2];
double upper[2];
};
+struct get_ctx {
+ size_t nbvert;
+ double two_pi_over_nbvert;
+ double radius;
+};
+
/*******************************************************************************
* 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 FINLINE int
+check_hemisphere(const struct ssol_quadric_hemisphere* hemisphere)
+{
+ return hemisphere && hemisphere->radius > 0;
+}
+
static INLINE int
check_quadric(const struct ssol_quadric* quadric)
{
@@ -92,6 +111,8 @@ check_quadric(const struct ssol_quadric* quadric)
return check_hyperbol(&quadric->data.hyperbol);
case SSOL_QUADRIC_PARABOLIC_CYLINDER:
return check_parabolic_cylinder(&quadric->data.parabolic_cylinder);
+ case SSOL_QUADRIC_HEMISPHERE:
+ return check_hemisphere(&quadric->data.hemisphere);
default: return 0;
}
}
@@ -110,8 +131,9 @@ check_punched_surface(const struct ssol_punched_surface* punched_surface)
size_t i;
if(!punched_surface
- || punched_surface->nb_carvings == 0
- || !punched_surface->carvings
+ || (punched_surface->nb_carvings == 0
+ && punched_surface->quadric->type != SSOL_QUADRIC_HEMISPHERE)
+ || (punched_surface->nb_carvings && !punched_surface->carvings)
|| !punched_surface->quadric
|| !check_quadric(punched_surface->quadric))
return 0;
@@ -235,6 +257,18 @@ parabolic_cylinder_z
return (p[1] * p[1]) * pcyl->one_over_4focal;
}
+static FINLINE double
+hemisphere_z
+ (const double p[2],
+ const struct priv_hemisphere_data* hemisphere)
+{
+ const double r2 = p[0] * p[0] + p[1] * p[1];
+ const double z2 = hemisphere->sqr_radius - r2;
+ /* manage numerical unaccuracy */
+ ASSERT(z2 >= -hemisphere->sqr_radius * FLT_EPSILON);
+ return (z2 > 0) ? -sqrt(z2) + hemisphere->radius : hemisphere->radius;
+}
+
static void
quadric_mesh_parabol_get_pos(const unsigned ivert, float pos[3], void* ctx)
{
@@ -244,7 +278,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);
@@ -262,7 +296,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);
@@ -281,7 +315,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);
@@ -290,6 +324,25 @@ quadric_mesh_parabolic_cylinder_get_pos
f3_set_d3(pos, p);
}
+static void
+quadric_mesh_hemisphere_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] = hemisphere_z(p, &msh->data->hemisphere);
+
+ /* 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 FINLINE int
aabb_is_degenerated(const double lower[2], const double upper[2])
{
@@ -324,6 +377,92 @@ carvings_compute_aabb
}
}
+static void
+carvings_compute_radius
+ (const struct ssol_carving* carvings,
+ const size_t ncarvings,
+ double* radius)
+{
+ size_t icarving;
+ double r2 = - DBL_MAX;
+ ASSERT(carvings && radius);
+
+ if(!ncarvings) {
+ *radius = DBL_MAX;
+ return;
+ }
+
+ FOR_EACH(icarving, 0, ncarvings) {
+ size_t ivert;
+ FOR_EACH(ivert, 0, carvings[icarving].nb_vertices) {
+ double pos[2];
+ /* Discard the polygons to subtract */
+ if (carvings[icarving].operation == SSOL_SUB) continue;
+
+ carvings[icarving].get(ivert, pos, carvings[icarving].context);
+ r2 = MMAX(r2, d2_dot(pos, pos));
+ }
+ }
+
+ *radius = sqrt(r2);
+}
+
+static res_T
+build_triangulated_disk
+ (struct darray_double* coords,
+ struct darray_size_t* ids,
+ const double radius,
+ const size_t nsteps)
+{
+ struct s3dut_mesh_data data;
+ struct s3dut_mesh* mesh = NULL;
+ double *c_ptr = NULL;
+ size_t* i_ptr = NULL;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(coords && ids && nsteps && radius > 0);
+ ASSERT(nsteps < UINT_MAX);
+
+ s3dut_create_hemisphere
+ (coords->allocator, radius, (unsigned)nsteps, (unsigned)nsteps, &mesh);
+ if (res != RES_OK) {
+ fprintf(stderr, "Could not create the hemisphere 3D data.\n");
+ goto error;
+ }
+
+ S3DUT(mesh_get_data(mesh, &data));
+ if (!data.nprimitives || !data.nvertices) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ darray_double_clear(coords);
+ darray_size_t_clear(ids);
+
+ /* Reserve the memory space for the plane vertices */
+ res = darray_double_resize(coords, data.nvertices * 2/*#coords per vertex*/);
+ if (res != RES_OK) goto error;
+
+ /* Reserve the memory space for the plane indices */
+ res = darray_size_t_resize(ids, data.nprimitives * 3/*#ids per triangle*/);
+ if (res != RES_OK) goto error;
+
+ c_ptr = darray_double_data_get(coords);
+ FOR_EACH(i, 0, data.nvertices) {
+ d2_set(c_ptr + i * 2, data.positions + i * 3); /* don't get z */
+ }
+ i_ptr = darray_size_t_data_get(ids);
+ FOR_EACH(i, 0, data.nprimitives * 3) i_ptr[i] = data.indices[i];
+
+exit:
+ if(mesh) S3DUT(mesh_ref_put(mesh));
+ return res;
+error:
+ darray_double_clear(coords);
+ darray_size_t_clear(ids);
+ goto exit;
+}
+
static res_T
build_triangulated_plane
(struct darray_double* coords,
@@ -339,7 +478,7 @@ build_triangulated_plane
double size_min;
double delta;
res_T res = RES_OK;
- ASSERT(coords && lower && upper && nsteps);
+ ASSERT(coords && ids && lower && upper && nsteps);
ASSERT(!aabb_is_degenerated(lower, upper));
darray_double_clear(coords);
@@ -411,7 +550,7 @@ error:
}
static res_T
-clip_triangulated_plane
+clip_triangulated_sheet
(struct darray_double* coords,
struct darray_size_t* ids,
struct scpr_mesh* mesh,
@@ -551,7 +690,7 @@ 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;
d2_set(ctx.lower, lower);
d2_set(ctx.upper, upper);
@@ -563,8 +702,9 @@ quadric_setup_s3d_shape_rt
vdata[1].type = S3D_FLOAT2;
vdata[1].get = quadric_mesh_get_uv;
- 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[0].get = quadric_mesh_parabol_get_pos;
break;
@@ -577,6 +717,9 @@ quadric_setup_s3d_shape_rt
case SSOL_QUADRIC_PLANE:
vdata[0].get = quadric_mesh_plane_get_pos;
break;
+ case SSOL_QUADRIC_HEMISPHERE:
+ vdata[0].get = quadric_mesh_hemisphere_get_pos;
+ break;
default: FATAL("Unreachable code.\n"); break;
}
@@ -684,36 +827,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
@@ -764,6 +913,18 @@ quadric_parabolic_cylinder_gradient_local
grad[2] = 2 * quad->focal;
}
+static FINLINE void
+quadric_hemisphere_gradient_local
+ (const struct priv_hemisphere_data* quad,
+ const double pt[3],
+ double grad[3])
+{
+ ASSERT(pt && grad);
+ grad[0] = -pt[0];
+ grad[1] = -pt[1];
+ grad[2] = quad->radius - pt[2];
+}
+
static FINLINE int
quadric_plane_intersect_local
(const double org[3],
@@ -777,13 +938,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;
}
@@ -798,17 +959,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;
}
@@ -822,7 +983,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;
@@ -833,12 +994,38 @@ 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;
+}
+
+static FINLINE int
+quadric_hemisphere_intersect_local
+ (const struct priv_hemisphere_data* quad,
+ const double org[3],
+ const double dir[3],
+ const double hint,
+ double hit_pt[3],
+ double grad[3],
+ double* dist)
+{
+ 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 n = solve_second(a, b, c, hint, 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;
return 1;
}
@@ -852,15 +1039,16 @@ 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;
}
@@ -869,18 +1057,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->private_data.hemisphere);
break;
default: FATAL("Unreachable code\n"); break;
}
@@ -894,28 +1085,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: {
+ 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->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],
@@ -931,21 +1125,25 @@ 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->private_data.hemisphere, org, dir, hint, pt, N, dist);
break;
default: FATAL("Unreachable code\n"); break;
}
@@ -1006,9 +1204,19 @@ priv_parabolic_cylinder_data_setup
data->one_over_4focal = 1 / (4.0 * parabolic_cylinder->focal);
}
+static FINLINE void
+priv_hemisphere_data_setup
+ (struct priv_hemisphere_data* data,
+ const struct ssol_quadric_hemisphere* hemisphere)
+{
+ ASSERT(data && hemisphere);
+ data->radius = hemisphere->radius;
+ data->sqr_radius = hemisphere->radius * hemisphere->radius;
+}
+
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);
@@ -1026,14 +1234,18 @@ priv_quadric_data_setup
priv_parabolic_cylinder_data_setup
(&priv_data->pcylinder, &quadric->data.parabolic_cylinder);
break;
+ case SSOL_QUADRIC_HEMISPHERE:
+ priv_hemisphere_data_setup
+ (&priv_data->hemisphere, &quadric->data.hemisphere);
+ break;
default: FATAL("Unreachable code\n"); break;
}
}
static INLINE size_t
-priv_quadric_data_compute_slices_count
+priv_quadric_data_compute_slices_count_aabb
(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])
{
@@ -1066,6 +1278,56 @@ priv_quadric_data_compute_slices_count
return nslices;
}
+static INLINE size_t
+priv_quadric_data_compute_slices_count_radius
+ (const enum ssol_quadric_type type,
+ const union private_data* priv_data,
+ const double radius)
+{
+ size_t nslices;
+ double pt[2];
+ double max_z = -DBL_MAX;
+
+ ASSERT(priv_data && radius > 0);
+ switch (type) {
+ case SSOL_QUADRIC_HEMISPHERE:
+ d2(pt, 0, radius);
+ max_z = MMAX(max_z, hemisphere_z(pt, &priv_data->hemisphere));
+
+ nslices = MMIN(50, (size_t)(3 + sqrt(max_z) * 6));
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ return nslices;
+}
+
+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);
+}
+
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -1087,8 +1349,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);
@@ -1111,13 +1373,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 */
@@ -1129,12 +1392,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 && N_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);
@@ -1146,14 +1408,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;
}
@@ -1246,11 +1508,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);
}
}
@@ -1278,6 +1540,7 @@ ssol_punched_surface_setup
const struct ssol_punched_surface* psurf)
{
double lower[2], upper[2]; /* Carvings AABB */
+ double radius;
struct darray_double coords;
struct darray_size_t ids;
size_t nslices;
@@ -1294,34 +1557,55 @@ ssol_punched_surface_setup
}
/* Save quadric for further object instancing */
- shape->quadric = *psurf->quadric;
-
- carvings_compute_aabb(psurf->carvings, psurf->nb_carvings, lower, upper);
- if(aabb_is_degenerated(lower, upper)) {
- log_error(shape->dev,
- "%s: infinite or null punched surface.\n",
- FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
+ 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) {
+ carvings_compute_radius(psurf->carvings, psurf->nb_carvings, &radius);
+ radius = MMIN(radius, psurf->quadric->data.hemisphere.radius);
+ lower[0] = lower[1] = -radius;
+ upper[0] = upper[1] = +radius;
+ } else {
+ carvings_compute_aabb(psurf->carvings, psurf->nb_carvings, lower, upper);
+ if(aabb_is_degenerated(lower, upper)) {
+ log_error(shape->dev,
+ "%s: infinite or null punched surface.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
}
/* 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);
+ if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) {
+ nslices = priv_quadric_data_compute_slices_count_radius
+ (shape->private_type.quadric, &shape->private_data, radius);
+ }
+ else {
+ nslices = priv_quadric_data_compute_slices_count_aabb
+ (shape->private_type.quadric, &shape->private_data, lower, upper);
+ }
}
- res = build_triangulated_plane(&coords, &ids, lower, upper, nslices);
- if(res != RES_OK) goto error;
-
- res = clip_triangulated_plane
- (&coords, &ids, shape->dev->scpr_mesh, psurf->carvings, psurf->nb_carvings);
+ if(psurf->quadric->type == SSOL_QUADRIC_HEMISPHERE) {
+ res = build_triangulated_disk(&coords, &ids, radius, nslices);
+ } else {
+ res = build_triangulated_plane(&coords, &ids, lower, upper, nslices);
+ }
if(res != RES_OK) goto error;
+
+ if(psurf->nb_carvings) {
+ res = clip_triangulated_sheet
+ (&coords, &ids, shape->dev->scpr_mesh, psurf->carvings, psurf->nb_carvings);
+ if(res != RES_OK) goto error;
+ }
/* Setup the Star-3D shape to ray-trace */
res = quadric_setup_s3d_shape_rt(shape, &coords, &ids, lower, upper,
@@ -1400,7 +1684,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;
@@ -1410,4 +1694,3 @@ exit:
error:
goto exit;
}
-
diff --git a/src/ssol_shape_c.h b/src/ssol_shape_c.h
@@ -43,10 +43,20 @@ struct priv_pcylinder_data {
double one_over_4focal;
};
-union priv_quadric_data {
+struct priv_hemisphere_data {
+ double radius;
+ double sqr_radius;
+};
+
+union private_data {
struct priv_hyperbol_data hyperbol;
struct priv_parabol_data parabol;
struct priv_pcylinder_data pcylinder;
+ struct priv_hemisphere_data hemisphere;
+};
+
+union private_type {
+ enum ssol_quadric_type quadric;
};
struct ssol_shape {
@@ -54,14 +64,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
@@ -91,5 +111,26 @@ 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/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/test_ssol_shape.c b/src/test_ssol_shape.c
@@ -151,13 +151,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;
@@ -196,6 +190,10 @@ main(int argc, char** argv)
quadric.data.parabol.focal = 1;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+ punched_surface.nb_carvings = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ punched_surface.nb_carvings = 1;
+
quadric.data.parabol.focal = 0;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
quadric.data.parabol.focal = 1;
@@ -205,6 +203,10 @@ main(int argc, char** argv)
quadric.data.hyperbol.img_focal = 1;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+ punched_surface.nb_carvings = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ punched_surface.nb_carvings = 1;
+
quadric.data.hyperbol.real_focal = 0;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
quadric.data.hyperbol.real_focal = 1;
@@ -217,10 +219,26 @@ main(int argc, char** argv)
quadric.data.parabolic_cylinder.focal = 1;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+ punched_surface.nb_carvings = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ punched_surface.nb_carvings = 1;
+
quadric.data.parabolic_cylinder.focal = 0;
CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
quadric.data.parabolic_cylinder.focal = 1;
+ quadric.type = SSOL_QUADRIC_HEMISPHERE;
+ quadric.data.hemisphere.radius = 10;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+
+ punched_surface.nb_carvings = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_OK);
+ punched_surface.nb_carvings = 1;
+
+ quadric.data.hemisphere.radius = 0;
+ CHECK(ssol_punched_surface_setup(shape, &punched_surface), RES_BAD_ARG);
+ quadric.data.hemisphere.radius = 10;
+
CHECK(ssol_shape_ref_put(shape), RES_OK);
CHECK(ssol_device_ref_put(dev), RES_OK);
