summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/sphere.c85
1 files changed, 27 insertions, 58 deletions
diff --git a/libdimension/sphere.c b/libdimension/sphere.c
index 8f1e1a7..468014e 100644
--- a/libdimension/sphere.c
+++ b/libdimension/sphere.c
@@ -59,74 +59,43 @@ dmnsn_sphere_inside_fn(const dmnsn_object *sphere, dmnsn_vector point)
return point.x*point.x + point.y*point.y + point.z*point.z < 1.0;
}
-/// Sphere bounding callback.
-static dmnsn_bounding_box
-dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+/// Helper for sphere bounding box calculation.
+static inline double
+dmnsn_implicit_dot(double row[4])
{
- // Get a tight bound using the conic representation of a sphere:
- //
- // S = [ 1 0 0 0 ]
- // [ 0 1 0 0 ]
- // [ 0 0 1 0 ]
- // [ 0 0 0 -1 ].
- //
- // The surface is defined by
- // p^T * S * p = 0,
- // and the tangent planes are defined by
- // q * S^-1 * q^T = 0.
- // Note that S = S^-1.
- //
- // The symmetric matrix R, defined by
- // R = M * S^-1 * M^T,
- // characterizes the tangent planes. Specifically,
- // min.x = (R[0,3] - sqrt(R[0,3]^2 - R[0,0]*R[3,3]))/R[3,3]
- // max.x = (R[0,3] + sqrt(R[0,3]^2 - R[0,0]*R[3,3]))/R[3,3]
- // min.y = (R[1,3] - sqrt(R[1,3]^2 - R[1,1]*R[3,3]))/R[3,3]
- // max.y = (R[1,3] + sqrt(R[1,3]^2 - R[1,1]*R[3,3]))/R[3,3]
- // min.z = (R[2,3] - sqrt(R[2,3]^2 - R[2,2]*R[3,3]))/R[3,3]
- // max.z = (R[2,3] + sqrt(R[2,3]^2 - R[2,2]*R[3,3]))/R[3,3]
- //
- // Unfortunately, we can't use dmnsn_matrix because the matrices are not
- // affine
-
- // MS = M * S^-1 = M * S
- // Last row is [ 0 0 0 -1 ] implicitly
- double MS[3][4];
+ double ret = 0.0;
for (int i = 0; i < 3; ++i) {
- for (int j = 0; j < 3; ++j) {
- MS[i][j] = trans.n[i][j];
- }
- MS[i][3] = -trans.n[i][3];
+ ret += row[i]*row[i];
}
+ return ret;
+}
- // R = MS * M^T
- // We only compute the upper triangular portion
- // R[3][3] is implicitly -1
- double R[4][4];
- for (int i = 0; i < 3; ++i) {
- for (int j = i; j < 3; ++j) {
- R[i][j] = 0.0;
- for (int k = 0; k < 4; ++k) {
- R[i][j] += MS[i][k]*trans.n[j][k];
- }
- }
- R[i][3] = MS[i][3];
- }
+/// Sphere bounding callback.
+static dmnsn_bounding_box
+dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+ // Get a tight bound using the quadric representation of a sphere. For
+ // details, see
+ // http://tavianator.com/2014/06/exact-bounding-boxes-for-spheres-ellipsoids
dmnsn_bounding_box box;
- double dx = sqrt(R[0][3]*R[0][3] + R[0][0]);
- box.min.x = -R[0][3] - dx;
- box.max.x = -R[0][3] + dx;
+ double cx = trans.n[0][3];
+ double dx = sqrt(dmnsn_implicit_dot(trans.n[0]));
+ box.min.x = cx - dx;
+ box.max.x = cx + dx;
- double dy = sqrt(R[1][3]*R[1][3] + R[1][1]);
- box.min.y = -R[1][3] - dy;
- box.max.y = -R[1][3] + dy;
+ double cy = trans.n[1][3];
+ double dy = sqrt(dmnsn_implicit_dot(trans.n[1]));
+ box.min.y = cy - dy;
+ box.max.y = cy + dy;
- double dz = sqrt(R[2][3]*R[2][3] + R[2][2]);
- box.min.z = -R[2][3] - dz;
- box.max.z = -R[2][3] + dz;
+ double cz = trans.n[2][3];
+ double dz = sqrt(dmnsn_implicit_dot(trans.n[2]));
+ box.min.z = cz - dz;
+ box.max.z = cz + dz;
+ printf(DMNSN_BOUNDING_BOX_FORMAT "\n", DMNSN_BOUNDING_BOX_PRINTF(box));
return box;
}