sphere: Use tightest possible bounding boxes.

3 files changed, 136 insertions, 6 deletions

diff --git a/dimension/tests/Makefile.am b/dimension/tests/Makefile.am
index 684b2aa..b5b6882 100644
--- a/dimension/tests/Makefile.am
+++ b/dimension/tests/Makefile.am
@@ -1,5 +1,5 @@
 ###########################################################################
-## Copyright (C) 2009-2011 Tavian Barnes <tavianator@tavianator.com>     ##
+## Copyright (C) 2009-2014 Tavian Barnes <tavianator@tavianator.com>     ##
 ##                                                                       ##
 ## This file is part of The Dimension Build Suite.                       ##
 ##                                                                       ##
@@ -17,8 +17,9 @@
 ## along with this program.  If not, see <http://www.gnu.org/licenses/>. ##
 ###########################################################################
 
-TESTS = demo.dmnsn                                                             \
-        cube.dmnsn
+TESTS = cube.dmnsn                                                             \
+        demo.dmnsn                                                             \
+        ellipsoid.dmnsn
 TEST_EXTENSIONS = .dmnsn
 DMNSN_LOG_COMPILER = $(top_srcdir)/dimension/dimension
 AM_DMNSN_LOG_FLAGS = --strict -v
diff --git a/dimension/tests/ellipsoid.dmnsn b/dimension/tests/ellipsoid.dmnsn
new file mode 100644
index 0000000..9b5e08e
--- /dev/null
+++ b/dimension/tests/ellipsoid.dmnsn
@@ -0,0 +1,67 @@
+#########################################################################
+# Copyright (C) 2010-2014 Tavian Barnes <tavianator@tavianator.com>     #
+#                                                                       #
+# This file is part of The Dimension Test Suite.                        #
+#                                                                       #
+# The Dimension Test Suite 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.                       #
+#                                                                       #
+# The Dimension Test Suite 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/>. #
+#########################################################################
+
+camera = PerspectiveCamera(location = -7*Z, look_at = 0)
+
+background = 0.5*sRGB(0.73, 0.90, 0.97)
+
+lights.append(PointLight(location = (0, 7, -7), color = White))
+
+objects.append(
+  Plane(
+    normal = Y, distance = -4,
+
+    texture = Texture(
+      pigment = sRGB(0.73, 0.90, 0.97),
+      finish = Ambient(sRGB(0.5)),
+    )
+  )
+)
+
+objects.append(
+  Sphere(
+    center = 0,
+    radius = 1,
+
+    pigment = White
+  )
+  .translate(-3*X)
+)
+
+objects.append(
+  Sphere(
+    center = 0,
+    radius = 1,
+
+    pigment = White
+  )
+  .scale(1.25, 0.75, 0.75)
+)
+
+objects.append(
+  Sphere(
+    center = 0,
+    radius = 1,
+
+    pigment = White
+  )
+  .scale(1.25, 0.75, 0.75)
+  .rotate(45*Z)
+  .translate(3*X)
+)
diff --git a/libdimension/sphere.c b/libdimension/sphere.c
index 5555499..8f1e1a7 100644
--- a/libdimension/sphere.c
+++ b/libdimension/sphere.c
@@ -63,9 +63,71 @@ dmnsn_sphere_inside_fn(const dmnsn_object *sphere, dmnsn_vector point)
 static dmnsn_bounding_box
 dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
 {
-  // TODO: tighter bound
-  dmnsn_bounding_box box = dmnsn_symmetric_bounding_box(dmnsn_new_vector(1.0, 1.0, 1.0));
-  return dmnsn_transform_bounding_box(trans, box);
+  // 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];
+  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];
+  }
+
+  // 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];
+  }
+
+  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 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 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;
+
+  return box;
 }
 
 /// Sphere vtable.