commit 8387a8776375f7e6ac91d2672bffc2544307ca95
parent a17d23562d7c3d7af70c9fdb63e4cb5b1a4b7371
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 7 Oct 2016 14:43:07 +0200
Add matrix transformation to the quadric data structure
This matrix transform the quadric in object space; i.e. in object space,
the quadric is not necessarily aligned with the axis.
Diffstat:
10 files changed, 125 insertions(+), 75 deletions(-)
diff --git a/src/ssol.h b/src/ssol.h
@@ -146,10 +146,16 @@ struct ssol_quadric {
struct ssol_quadric_parabol parabol;
struct ssol_quadric_parabolic_cylinder parabolic_cylinder;
} data;
+
+ /* 3x4 column major transformation of the quadric in object space */
+ double transform[12];
};
-#define SSOL_QUADRIC_DEFAULT__ \
- {SSOL_QUADRIC_PLANE, {SSOL_QUADRIC_PLANE_DEFAULT__}}
+#define SSOL_QUADRIC_DEFAULT__ { \
+ SSOL_QUADRIC_PLANE, \
+ {SSOL_QUADRIC_PLANE_DEFAULT__}, \
+ {1,0,0, 0,1,0, 0,0,1, 0,0,0} \
+}
static const struct ssol_quadric SSOL_QUADRIC_DEFAULT = SSOL_QUADRIC_DEFAULT__;
/* Define the contour of a 2D polygon as well as the clipping operation to
diff --git a/src/ssol_scene.c b/src/ssol_scene.c
@@ -119,7 +119,7 @@ ssol_scene_attach_instance
res_T res;
if(!scene || !instance) return RES_BAD_ARG;
-
+
/* Attach the instantiated s3d shape to ray-trace to the RT scene */
res = s3d_scene_attach_shape(scene->scn_rt, instance->shape_rt);
if(res != RES_OK) return res;
@@ -358,41 +358,51 @@ hit_filter_function
shape = inst->object->shape;
seg->on_punched = (shape->type == SHAPE_PUNCHED);
switch (shape->type) {
- case SHAPE_PUNCHED: {
- /* hits on quadrics must be recomputed more accurately */
- double org_local[3], hit_pos_local[3], dir_local[3];
- const double* transform = inst->transform;
- double tr[9];
- int valid;
- d33_inverse(tr, transform);
-
- /* get org in local coordinate */
- d3_set(org_local, seg->org);
- d3_sub(org_local, org_local, transform + 9);
- d33_muld3(org_local, tr, org_local);
-
- /* get dir in local */
- d33_muld3(dir_local, tr, seg->dir);
- /* recompute hit */
- valid = punched_shape_intersect_local(shape, org_local, dir_local,
- hit->distance, hit_pos_local, seg->hit_normal, &seg->hit_distance);
- if (!valid) return 1;
- /* transform point to world */
- d33_muld3(seg->hit_pos, transform, hit_pos_local);
- d3_add(seg->hit_pos, transform + 9, seg->hit_pos);
- /* transform normal to world */
- d33_invtrans(tr, transform);
- d33_muld3(seg->hit_normal, tr, seg->hit_normal);
- break;
- }
- case SHAPE_MESH: {
- d3_set_f3(seg->hit_normal, hit->normal);
- /* use raytraced distance to fill hit_pos */
- d3_add(seg->hit_pos, seg->org, d3_muld(seg->hit_pos, seg->dir, hit->distance));
- seg->hit_distance = hit->distance;
- break;
- }
- default: FATAL("Unreachable code.\n"); break;
+ case SHAPE_MESH: {
+ d3_set_f3(seg->hit_normal, hit->normal);
+ /* use raytraced distance to fill hit_pos */
+ d3_add(seg->hit_pos, seg->org, d3_muld(seg->hit_pos, seg->dir, hit->distance));
+ seg->hit_distance = hit->distance;
+ break;
+ }
+ case SHAPE_PUNCHED: {
+ /* hits on quadrics must be recomputed more accurately */
+ double org_local[3], hit_pos_local[3], dir_local[3];
+ double R[9]; /* Rotation matrix */
+ double T[3]; /* Translation vector */
+ double R_invtrans[9]; /* Inverse transpose rotation matrix */
+ double T_inv[3]; /* Inverse of the translation vector */
+ int valid;
+
+ if(d33_is_identity(shape->quadric.transform)) {
+ d33_set(R, inst->transform);
+ d3_set (T, inst->transform+9);
+ } else {
+ d33_muld33(R, shape->quadric.transform, inst->transform);
+ d33_muld3 (T, shape->quadric.transform, inst->transform+9);
+ }
+ d33_invtrans(R_invtrans, R);
+ d3_minus(T_inv, T);
+
+ /* get org in local coordinate */
+ d3_set(org_local, seg->org);
+ d3_add(org_local, org_local, T_inv);
+ d3_muld33(org_local, org_local, R_invtrans);
+
+ /* get dir in local */
+ d3_muld33(dir_local, seg->dir, R_invtrans);
+ /* recompute hit */
+ valid = punched_shape_intersect_local(shape, org_local, dir_local,
+ hit->distance, hit_pos_local, seg->hit_normal, &seg->hit_distance);
+ if (!valid) return 1;
+ /* transform point to world */
+ d33_muld3(seg->hit_pos, R, hit_pos_local);
+ d3_add(seg->hit_pos, seg->hit_pos, T);
+ /* transform normal to world */
+ d33_muld3(seg->hit_normal, R_invtrans, seg->hit_normal);
+ break;
+ }
+ default: FATAL("Unreachable code.\n"); break;
}
d3_normalize(seg->hit_normal, seg->hit_normal);
diff --git a/src/ssol_shape.c b/src/ssol_shape.c
@@ -18,10 +18,12 @@
#include "ssol_device_c.h"
#include "ssol_shape_c.h"
-#include <rsys/double3.h>
#include <rsys/double2.h>
+#include <rsys/double3.h>
+#include <rsys/double33.h>
#include <rsys/dynamic_array_double.h>
#include <rsys/dynamic_array_size_t.h>
+#include <rsys/float3.h>
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
#include <rsys/rsys.h>
@@ -40,6 +42,7 @@ struct quadric_mesh_context {
const double* coords;
const size_t* ids;
double focal; /* Use by parabol and parabolic cylinder quadrics */
+ const double* transform; /* 3x4 column major matrix */
};
/*******************************************************************************
@@ -171,10 +174,17 @@ quadric_mesh_plane_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);
- pos[0] = (float)msh->coords[i+0];
- pos[1] = (float)msh->coords[i+1];
- pos[2] = 0.f;
+ p[0] = (float)msh->coords[i+0];
+ p[1] = (float)msh->coords[i+1];
+ p[2] = 0.f;
+
+ /* 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
@@ -182,13 +192,17 @@ quadric_mesh_parabol_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 x, y;
+ double p[3]; /* Temporary quadric space position */
ASSERT(pos && ctx);
- x = msh->coords[i+0];
- y = msh->coords[i+1];
- pos[0] = (float)x;
- pos[1] = (float)y;
- pos[2] = (float)((x*x + y*y) / (4.0*msh->focal));
+ 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);
+
+ /* 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
@@ -197,13 +211,17 @@ quadric_mesh_parabolic_cylinder_get_pos
{
const size_t i = ivert*2/*#coords per vertex*/;
const struct quadric_mesh_context* msh = ctx;
- double x, y;
+ double p[3]; /* Temporary quadric space position */
ASSERT(pos && ctx);
- x = msh->coords[i+0];
- y = msh->coords[i+1];
- pos[0] = (float)x;
- pos[1] = (float)y;
- pos[2] = (float)((y*y) / (4.0*msh->focal));
+ p[0] = msh->coords[i+0];
+ p[1] = msh->coords[i+1];
+ p[2] = ((p[1]*p[1]) / (4.0*msh->focal));
+
+ /* 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
@@ -425,6 +443,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 = quadric->transform;
vdata.usage = S3D_POSITION;
vdata.type = S3D_FLOAT3;
@@ -538,9 +557,9 @@ inject_same_sign(const double x, const double src)
return ucast.d;
}
-/* 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 positive value */
+/* 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 positive value */
static int
quadric_solve_second
(const double a,
@@ -551,11 +570,11 @@ quadric_solve_second
{
ASSERT(dist);
if (a != 0) {
- /* standard case: 2nd degree */
+ /* Standard case: 2nd degree */
const double delta = b * b - 4 * a * c;
if (delta > 0) {
const double sqrt_delta = sqrt(delta);
- /* precise formula */
+ /* Precise formula */
const double t1 = (-b - inject_same_sign(sqrt_delta, b)) / (2 * b);
const double t2 = c / (a * t1);
if (t1 < 0 && t2 < 0) return 0; /* no positive solution */
@@ -567,7 +586,7 @@ quadric_solve_second
*dist = t1; /* t1 is the only positive solution */
return 1;
}
- /* both t1 and t2 are positive: choose the closest value to hint */
+ /* Both t1 and t2 are positive: choose the closest value to hint */
*dist = fabs(t1 - hint) < fabs(t2 - hint) ? t1 : t2;
return 1;
} else if (delta == 0) {
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -502,21 +502,36 @@ sample_starting_point(struct realisation* rs)
break;
}
case SHAPE_PUNCHED: {
- const double* transform = start->instance->transform;
- double tr[9], pos_local[3];
+ struct ssol_instance* inst = start->instance;
+ double pos_local[3];
+ double R[9]; /* Rotation matrix */
+ double T[3]; /* Translation vector */
+ double R_invtrans[9]; /* Inverse transpose rotation matrix */
+ double T_inv[3]; /* Inverse of the translation vector */
+
+ if(d33_is_identity(shape->quadric.transform)) {
+ d33_set(R, inst->transform);
+ d3_set (T, inst->transform+9);
+ } else {
+ d33_muld33(R, shape->quadric.transform, inst->transform);
+ d33_muld3 (T, shape->quadric.transform, inst->transform+9);
+ }
+ d33_invtrans(R_invtrans, R);
+ d3_minus(T_inv, T);
+
/* project the sampled point on the quadric */
- d33_inverse(tr, transform);
- d3_sub(pos_local, start->pos, transform + 9);
- d33_muld3(pos_local, tr, pos_local);
+ d3_set(pos_local, start->pos);
+ d3_add(pos_local, pos_local, T_inv);
+ d3_muld33(pos_local, pos_local, R_invtrans);
+
punched_shape_set_z_local(shape, pos_local);
/* transform point to world */
- d33_muld3(start->pos, transform, pos_local);
- d3_add(start->pos, transform + 9, start->pos);
+ d33_muld3(start->pos, R, pos_local);
+ d3_add(start->pos, start->pos, T);
/* compute exact normal on the instance */
punched_shape_set_normal_local(shape, pos_local, start->rt_normal);
/* transform normal to world */
- d33_invtrans(tr, transform);
- d33_muld3(start->rt_normal, tr, start->rt_normal);
+ d33_muld3(start->rt_normal, R_invtrans, start->rt_normal);
break;
}
default: FATAL("Unreachable code.\n"); break;
diff --git a/src/test_ssol_solver2.c b/src/test_ssol_solver2.c
@@ -56,7 +56,7 @@ main(int argc, char** argv)
struct ssol_vertex_data attribs[1];
struct ssol_shape* quad_square;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
diff --git a/src/test_ssol_solver2b.c b/src/test_ssol_solver2b.c
@@ -56,7 +56,7 @@ main(int argc, char** argv)
struct ssol_shape* quad_square;
struct ssol_shape* quad_rect;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
diff --git a/src/test_ssol_solver3.c b/src/test_ssol_solver3.c
@@ -50,7 +50,7 @@ main(int argc, char** argv)
struct ssol_vertex_data attribs[1];
struct ssol_shape* quad_square;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
diff --git a/src/test_ssol_solver3N.c b/src/test_ssol_solver3N.c
@@ -94,7 +94,7 @@ main(int argc, char** argv)
struct ssol_vertex_data attribs[1];
struct ssol_shape* quad_square;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
diff --git a/src/test_ssol_solver4.c b/src/test_ssol_solver4.c
@@ -50,7 +50,7 @@ main(int argc, char** argv)
struct ssol_vertex_data attribs[1];
struct ssol_shape* quad_square;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
diff --git a/src/test_ssol_solver5.c b/src/test_ssol_solver5.c
@@ -50,7 +50,7 @@ main(int argc, char** argv)
struct ssol_vertex_data attribs[1];
struct ssol_shape* quad_square;
struct ssol_carving carving;
- struct ssol_quadric quadric;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched;
struct ssol_material* m_mtl;
struct ssol_material* v_mtl;
