Transform normals as pseudovectors, not vectors.

Also clarify the vector transformation API.  Instead of
dmnsn_transform_vector(), we have:
  - dmnsn_transform_point()
  - dmnsn_transform_direction()
  - dmnsn_transform_normal()

8 files changed, 71 insertions, 47 deletions

diff --git a/libdimension-python/dimension.pxd b/libdimension-python/dimension.pxd
index 42fda7e..1164be5 100644
--- a/libdimension-python/dimension.pxd
+++ b/libdimension-python/dimension.pxd
@@ -130,7 +130,9 @@ cdef extern from "../libdimension/dimension.h":
   dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix m)
 
   dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs)
-  dmnsn_vector dmnsn_transform_vector(dmnsn_matrix lhs, dmnsn_vector rhs)
+  dmnsn_vector dmnsn_transform_point(dmnsn_matrix lhs, dmnsn_vector rhs)
+  dmnsn_vector dmnsn_transform_direction(dmnsn_matrix lhs, dmnsn_vector rhs)
+  dmnsn_vector dmnsn_transform_normal(dmnsn_matrix lhs, dmnsn_vector rhs)
 
   dmnsn_matrix dmnsn_scale_matrix(dmnsn_vector s)
   dmnsn_matrix dmnsn_translation_matrix(dmnsn_vector d)
diff --git a/libdimension-python/dimension.pyx b/libdimension-python/dimension.pyx
index a1c8710..48e8de1 100644
--- a/libdimension-python/dimension.pyx
+++ b/libdimension-python/dimension.pyx
@@ -263,7 +263,7 @@ cdef class Matrix:
     if isinstance(rhs, Matrix):
       return _Matrix(dmnsn_matrix_mul(lhs._m, (<Matrix>rhs)._m))
     else:
-      return _Vector(dmnsn_transform_vector(lhs._m, (<Vector>rhs)._v))
+      return _Vector(dmnsn_transform_point(lhs._m, (<Vector>rhs)._v))
 
   def __richcmp__(Matrix lhs not None, Matrix rhs not None, int op):
     cdef double sum = 0.0
diff --git a/libdimension/bench/geometry.c b/libdimension/bench/geometry.c
index a952631..4b8c12a 100644
--- a/libdimension/bench/geometry.c
+++ b/libdimension/bench/geometry.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2009-2010 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2009-2011 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Benchmark Suite.                   *
  *                                                                       *
@@ -149,11 +149,23 @@ main(void)
   });
   printf("dmnsn_matrix_mul(): %ld\n", sandglass.grains);
 
-  /* dmnsn_transform_vector() */
+  /* dmnsn_transform_point() */
   sandglass_bench_fine(&sandglass, {
-    vector = dmnsn_transform_vector(matrix, vector);
+    vector = dmnsn_transform_point(matrix, vector);
   });
-  printf("dmnsn_transform_vector(): %ld\n", sandglass.grains);
+  printf("dmnsn_transform_point(): %ld\n", sandglass.grains);
+
+  /* dmnsn_transform_direction() */
+  sandglass_bench_fine(&sandglass, {
+    vector = dmnsn_transform_direction(matrix, vector);
+  });
+  printf("dmnsn_transform_direction(): %ld\n", sandglass.grains);
+
+  /* dmnsn_transform_normal() */
+  sandglass_bench_fine(&sandglass, {
+    vector = dmnsn_transform_normal(matrix, vector);
+  });
+  printf("dmnsn_transform_normal(): %ld\n", sandglass.grains);
 
   /* dmnsn_transform_line() */
   sandglass_bench_fine(&sandglass, {
diff --git a/libdimension/dimension/geometry.h b/libdimension/dimension/geometry.h
index 45aa6c6..d2ad362 100644
--- a/libdimension/dimension/geometry.h
+++ b/libdimension/dimension/geometry.h
@@ -375,9 +375,9 @@ dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix A);
 /** Multiply two matricies. */
 dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs);
 
-/** Transform a vector by a matrix. */
+/** Transform a point by a matrix. */
 DMNSN_INLINE dmnsn_vector
-dmnsn_transform_vector(dmnsn_matrix T, dmnsn_vector v)
+dmnsn_transform_point(dmnsn_matrix T, dmnsn_vector v)
 {
   /* 9 multiplications, 9 additions */
   dmnsn_vector r;
@@ -387,6 +387,36 @@ dmnsn_transform_vector(dmnsn_matrix T, dmnsn_vector v)
   return r;
 }
 
+/** Transform a direction by a matrix. */
+DMNSN_INLINE dmnsn_vector
+dmnsn_transform_direction(dmnsn_matrix T, dmnsn_vector v)
+{
+  /* 9 multiplications, 9 additions */
+  dmnsn_vector r;
+  r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z;
+  r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z;
+  r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z;
+  return r;
+}
+
+/**
+ * Transform a pseudovector by a matrix.
+ * @param[in] Tinv  The inverse of the transformation matrix.
+ * @param[in] v     The pseudovector to transform
+ * @return The transformed pseudovector.
+ */
+DMNSN_INLINE dmnsn_vector
+dmnsn_transform_normal(dmnsn_matrix Tinv, dmnsn_vector v)
+{
+  /* Multiply by the transpose of the inverse
+     (9 multiplications, 6 additions) */
+  dmnsn_vector r;
+  r.x = Tinv.n[0][0]*v.x + Tinv.n[1][0]*v.y + Tinv.n[2][0]*v.z;
+  r.y = Tinv.n[0][1]*v.x + Tinv.n[1][1]*v.y + Tinv.n[2][1]*v.z;
+  r.z = Tinv.n[0][2]*v.x + Tinv.n[1][2]*v.y + Tinv.n[2][2]*v.z;
+  return r;
+}
+
 /** Transform a bounding box by a matrix. */
 dmnsn_bounding_box dmnsn_transform_bounding_box(dmnsn_matrix T,
                                                 dmnsn_bounding_box box);
@@ -401,11 +431,8 @@ dmnsn_transform_line(dmnsn_matrix T, dmnsn_line l)
 {
   /* 18 multiplications, 24 additions */
   dmnsn_line ret;
-  ret.x0 = dmnsn_transform_vector(T, l.x0);
-  ret.n  = dmnsn_vector_sub(
-    dmnsn_transform_vector(T, dmnsn_vector_add(l.x0, l.n)),
-    ret.x0
-  );
+  ret.x0 = dmnsn_transform_point(T, l.x0);
+  ret.n  = dmnsn_transform_direction(T, l.n);
   return ret;
 }
 
diff --git a/libdimension/dimension/object.h b/libdimension/dimension/object.h
index 8254782..8b497ce 100644
--- a/libdimension/dimension/object.h
+++ b/libdimension/dimension/object.h
@@ -112,24 +112,6 @@ void dmnsn_delete_object(dmnsn_object *object);
 void dmnsn_initialize_object(dmnsn_object *object);
 
 /**
- * Transform a surface normal vector.
- * @param[in] trans   The transformation matrix.
- * @param[in] normal  The normal vector to transform
- * @return The transformed normal vector.
- */
-DMNSN_INLINE dmnsn_vector
-dmnsn_transform_normal(dmnsn_matrix trans, dmnsn_vector normal)
-{
-  return dmnsn_vector_normalized(
-    dmnsn_vector_sub(
-      dmnsn_transform_vector(trans, normal),
-      /* Optimized form of dmnsn_transform_vector(trans, dmnsn_zero) */
-      dmnsn_new_vector(trans.n[0][3], trans.n[1][3], trans.n[2][3])
-    )
-  );
-}
-
-/**
  * Appropriately transform a ray, then test for an intersection.
  * @param[in]  object        The object to test.
  * @param[in]  line          The ray to test.
@@ -144,9 +126,10 @@ dmnsn_object_intersection(const dmnsn_object *object, dmnsn_line line,
   intersection->object = NULL;
   if (object->intersection_fn(object, line_trans, intersection)) {
     /* Get us back into world coordinates */
-    intersection->ray    = line;
-    intersection->normal = dmnsn_transform_normal(object->trans,
-                                                  intersection->normal);
+    intersection->ray = line;
+    intersection->normal = dmnsn_vector_normalized(
+      dmnsn_transform_normal(object->trans_inv, intersection->normal)
+    );
     if (!intersection->object) {
       intersection->object = object;
     }
@@ -166,6 +149,6 @@ dmnsn_object_intersection(const dmnsn_object *object, dmnsn_line line,
 DMNSN_INLINE bool
 dmnsn_object_inside(const dmnsn_object *object, dmnsn_vector point)
 {
-  point = dmnsn_transform_vector(object->trans_inv, point);
+  point = dmnsn_transform_point(object->trans_inv, point);
   return object->inside_fn(object, point);
 }
diff --git a/libdimension/geometry.c b/libdimension/geometry.c
index 83be0f9..c6a4da2 100644
--- a/libdimension/geometry.c
+++ b/libdimension/geometry.c
@@ -113,8 +113,8 @@ dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to,
 {
   double theta1 = dmnsn_axis_angle(from, to, axis1);
   dmnsn_matrix align1 = dmnsn_rotation_matrix(dmnsn_vector_mul(theta1, axis1));
-  from  = dmnsn_transform_vector(align1, from);
-  axis2 = dmnsn_transform_vector(align1, axis2);
+  from  = dmnsn_transform_direction(align1, from);
+  axis2 = dmnsn_transform_direction(align1, axis2);
 
   double theta2 = dmnsn_axis_angle(from, to, axis2);
   dmnsn_matrix align2 = dmnsn_rotation_matrix(dmnsn_vector_mul(theta2, axis2));
@@ -381,41 +381,41 @@ dmnsn_transform_bounding_box(dmnsn_matrix trans, dmnsn_bounding_box box)
 
   dmnsn_vector corner;
   dmnsn_bounding_box ret;
-  ret.min = dmnsn_transform_vector(trans, box.min);
+  ret.min = dmnsn_transform_point(trans, box.min);
   ret.max = ret.min;
 
   corner  = dmnsn_new_vector(box.min.x, box.min.y, box.max.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.min.x, box.max.y, box.min.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.min.x, box.max.y, box.max.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.max.x, box.min.y, box.min.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.max.x, box.min.y, box.max.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.max.x, box.max.y, box.min.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
   corner  = dmnsn_new_vector(box.max.x, box.max.y, box.max.z);
-  corner  = dmnsn_transform_vector(trans, corner);
+  corner  = dmnsn_transform_point(trans, corner);
   ret.min = dmnsn_vector_min(ret.min, corner);
   ret.max = dmnsn_vector_max(ret.max, corner);
 
diff --git a/libdimension/pigment.c b/libdimension/pigment.c
index 8c1cc9a..3650ba3 100644
--- a/libdimension/pigment.c
+++ b/libdimension/pigment.c
@@ -71,7 +71,7 @@ dmnsn_color
 dmnsn_evaluate_pigment(const dmnsn_pigment *pigment, dmnsn_vector v)
 {
   if (pigment->pigment_fn) {
-    dmnsn_vector v_trans = dmnsn_transform_vector(pigment->trans_inv, v);
+    dmnsn_vector v_trans = dmnsn_transform_point(pigment->trans_inv, v);
     return pigment->pigment_fn(pigment, v_trans);
   } else {
     return pigment->quick_color;
diff --git a/libdimension/raytrace.c b/libdimension/raytrace.c
index fdfa3cb..4ee451a 100644
--- a/libdimension/raytrace.c
+++ b/libdimension/raytrace.c
@@ -139,7 +139,7 @@ dmnsn_initialize_raytrace_state(dmnsn_raytrace_state *state,
   state->interior     = intersection->object->interior;
 
   state->r = dmnsn_line_point(intersection->ray, intersection->t);
-  state->pigment_r = dmnsn_transform_vector(
+  state->pigment_r = dmnsn_transform_point(
     intersection->object->pigment_trans,
     state->r
   );