generalize Snell's law for ND universes; fixes #16

scenes/4d_fresnel.json

@@ -0,0 +1,86 @@
+{
+  "Universe4": {
+    "camera": {
+      "FreeCamera4": []
+    },
+    "entities": [
+      {
+        "Entity4Impl::new": [
+          {
+            "Sphere4::new": [
+              {
+                "Point4::new": [
+                  10,
+                  0,
+                  0,
+                  0
+                ]
+              },
+              3
+            ]
+          },
+          {
+            "Vacuum4::new": []
+          },
+          {
+            "ComposableSurface4": {
+              "reflection_ratio": {
+                "reflection_ratio_fresnel_4": [
+                  1.458,
+                  1
+                ]
+              },
+              "reflection_direction": {
+                "reflection_direction_specular_4": []
+              },
+              "threshold_direction": {
+                "threshold_direction_snell_4": [
+                    1.458
+                ]
+              },
+              "surface_color": {
+                "surface_color_uniform_4": [
+                  {
+                    "Rgba::new": [
+                      0,
+                      0,
+                      0,
+                      0
+                    ]
+                  }
+                ]
+              }
+            }
+          }
+        ]
+      },
+      {
+        "Void4::new_with_vacuum": []
+      }
+    ],
+    "background": {
+      "MappedTextureImpl4::new": [
+        {
+          "uv_derank_4": [
+            {
+              "uv_sphere_3": [
+                {
+                  "Point3::new": [
+                    0,
+                    0,
+                    0
+                  ]
+                }
+              ]
+            }
+          ]
+        },
+        {
+          "texture_image": [
+            "./resources/pixelcg_uv.jpg"
+          ]
+        }
+      ]
+    }
+  }
+}

src/scene.rs

@@ -1242,6 +1242,14 @@ impl Parser {
                 }
             }
 
+            add_deserializer! {
+                "threshold_direction_snell_4";
+                [refractive_index: F]
+                -> Box<ThresholdDirectionProvider<F, Point4<F>, Vector4<F>>> {
+                    threshold_direction_snell(refractive_index)
+                }
+            }
+
             add_deserializer! {
                 "threshold_direction_identity_3";
                 -> Box<ThresholdDirectionProvider<F, Point3<F>, Vector3<F>>> {

src/universe/d3/entity/surface.rs

@@ -1,10 +1,6 @@
 use num::traits::NumCast;
 use rand::StdRng;
 use rand::Rand;
-use num::One;
-use na;
-use na::Rotate;
-use na::UnitQuaternion;
 use na::Cast;
 use na::Point2;
 use na::Point3;
@@ -16,9 +12,7 @@ use palette;
 use palette::Hsv;
 use palette::RgbHue;
 use util::CustomFloat;
-use util::AngleBetween;
 use universe::entity::surface::Surface;
-use universe::entity::surface::ThresholdDirectionProvider;
 use universe::entity::surface::UVFn;
 use universe::entity::surface::SurfaceColorProvider;
 use universe::entity::shape::TracingContext;
@@ -43,28 +37,6 @@ pub fn surface_color_perlin_hue<F: CustomFloat>
     })
 }
 
-// Possibly generalize for n-dimensional spaces
-// http://math.stackexchange.com/questions/1402362/rotation-in-4d
-pub fn threshold_direction_snell<F: CustomFloat>
-    (refractive_index: F)
-    -> Box<ThresholdDirectionProvider<F, Point3<F>, Vector3<F>>>
-{
-    Box::new(move |context: &TracingContext<F, Point3<F>, Vector3<F>>| {
-        let normal = -context.intersection_normal_closer;
-        let axis = na::cross(&context.intersection.direction, &normal);
-        let from_theta = context.intersection.direction.angle_between(&normal);
-        let refractive_index_modifier = if context.exiting {
-            refractive_index
-        } else {
-            <F as One>::one() / refractive_index
-        };
-        let to_theta = (refractive_index_modifier * from_theta.sin()).asin();
-        let quaternion = UnitQuaternion::new(axis * (from_theta - to_theta));
-
-        quaternion.rotate(&context.intersection.direction)
-    })
-}
-
 pub fn surface_color_perlin_hue_seed<F: CustomFloat>
     (seed: u32,
      size: F,

src/universe/entity/surface.rs

@@ -17,7 +17,6 @@ use util::CustomPoint;
 use util::CustomVector;
 use num::Zero;
 use num::One;
-use na;
 use na::Cast;
 use na::Point2;
 use palette::Rgb;
@@ -222,7 +221,7 @@ pub fn reflection_ratio_fresnel<F: CustomFloat, P: CustomPoint<F, V>, V: CustomV
         };
         let to_theta = ((from_index / to_index) * from_theta.sin()).asin();
 
-        let ratio = if to_theta.is_nan() {
+        if to_theta.is_nan() {
             <F as One>::one()
         } else {
             // s-polarized light
@@ -237,9 +236,7 @@ pub fn reflection_ratio_fresnel<F: CustomFloat, P: CustomPoint<F, V>, V: CustomV
 
             // to get the reflectance of unpolarised light, we take the average
             (reflectance_s + reflectance_p) / (<F as One>::one() + <F as One>::one())
-        };
-
-        ratio
+        }
     })
 }
 
@@ -265,6 +262,28 @@ pub fn threshold_direction_identity<F: CustomFloat, P: CustomPoint<F, V>, V: Cus
     })
 }
 
+pub fn threshold_direction_snell<F: CustomFloat, P: CustomPoint<F, V>, V: CustomVector<F, P>>
+    (refractive_index: F)
+    -> Box<ThresholdDirectionProvider<F, P, V>>
+{
+    Box::new(move |context: &TracingContext<F, P, V>| {
+        let normal = -context.intersection_normal_closer;
+        let from_theta = context.intersection.direction.angle_between(&normal);
+        let refractive_index_modifier = if context.exiting {
+            refractive_index
+        } else {
+            <F as One>::one() / refractive_index
+        };
+        let to_theta = (refractive_index_modifier * from_theta.sin()).asin();
+        let angle_delta = to_theta - from_theta;
+        let mut data = [context.intersection.direction];
+
+        normal.general_rotation(&context.intersection.direction, angle_delta, &mut data);
+
+        data[0]
+    })
+}
+
 pub type BlendFunction<F> = Fn(Rgba<F>, Rgba<F>) -> Rgba<F>;
 pub type PaletteBlendFunction<C: Blend<Color=C> + ComponentWise> =
     Fn(PreAlpha<C, C::Scalar>, PreAlpha<C, C::Scalar>) -> PreAlpha<C, C::Scalar>;

src/util.rs

@@ -39,6 +39,9 @@ use num::traits::ToPrimitive;
 use palette;
 use image::Rgba;
 use na;
+use na::Eye;
+use na::Column;
+use na::Transpose;
 use na::BaseFloat;
 use na::Cast;
 use na::ApproxEq;
@@ -61,14 +64,19 @@ use na::Mean;
 use na::Translate;
 use na::Point2;
 use na::Vector2;
+use na::Matrix2;
 use na::Point3;
 use na::Vector3;
+use na::Matrix3;
 use na::Point4;
 use na::Vector4;
+use na::Matrix4;
 use na::Point5;
 use na::Vector5;
+use na::Matrix5;
 use na::Point6;
 use na::Vector6;
+use na::Matrix6;
 use core::iter::FromIterator;
 use core::marker::Reflect;
 use core::ops::DerefMut;
@@ -532,6 +540,7 @@ pub trait CustomPoint<F: CustomFloat, V: CustomVector<F, Self>>:
 
 pub trait CustomVector<F: CustomFloat, P: CustomPoint<F, Self>>:
 // Rotate<O> +
+    GeneralRotation<F> +
     AngleBetween<F> +
     PartialOrder +
     Div<F, Output=Self> +
@@ -618,8 +627,15 @@ pub trait Derank {
     fn derank(&self) -> Self::Type;
 }
 
+/// Use the Gram-Schmidt process to orthonormalize this vector in relation to the other
+/// and return the rotation matrix.
+// TODO: Not very ergonomic having to pass in a slice.
+pub trait GeneralRotation<F: CustomFloat>: Sized {
+    fn general_rotation(&self, other: &Self, angle: F, vectors_to_rotate: &mut [Self]);
+}
+
 macro_rules! dimension {
-    ($point:ident, $vector:ident) => {
+    ($point:ident, $vector:ident, $matrix:ident, $dimension:expr) => {
         impl<F: CustomFloat> CustomPoint<F, $vector<F>> for $point<F> {}
         impl<F: CustomFloat> CustomVector<F, $point<F>> for $vector<F> {}
 
@@ -638,14 +654,51 @@ macro_rules! dimension {
                 *self.as_mut() = *coords.as_ref();
             }
         }
+
+        impl<F: CustomFloat> GeneralRotation<F> for $vector<F> {
+            fn general_rotation(&self, other: &Self, angle: F, vectors_to_rotate: &mut [$vector<F>]) {
+                let mut original = $matrix::new_identity($dimension);
+
+                // Set the rotation plane
+                original.set_column(0, *self);
+                original.set_column(1, *other);
+
+                // Orthonormalize the rotation plane vectors using the Gram-Schmidt process
+                // in order to get the translation matrix
+                let mut result = original;
+
+                {
+                    result.set_column(0, original.column(0));
+
+                    for i in 1 .. $dimension {
+                        for j in 0 .. i {
+                            let original_column = original.column(i);
+                            original.set_column(i, original_column - result.column(j) * result.column(j).dot(&original_column));
+                        }
+
+                        result.set_column(i, original.column(i).normalize());
+                    }
+                }
+
+                // Create the rotation matrix that rotates translated vectors
+                let mut rotation_matrix = $matrix::new_identity($dimension);
+                rotation_matrix.m11 = angle.cos(); rotation_matrix.m12 = -angle.sin();
+                rotation_matrix.m21 = angle.sin(); rotation_matrix.m22 = angle.cos();
+                let result = result * (rotation_matrix * result.transpose());
+
+                for vector in vectors_to_rotate {
+                    *vector = result * *vector;
+                }
+            }
+        }
     }
 }
 
-dimension!(Point2, Vector2);
-dimension!(Point3, Vector3);
-dimension!(Point4, Vector4);
-dimension!(Point5, Vector5);
-dimension!(Point6, Vector6);
+dimension!(Point2, Vector2, Matrix2, 2);
+dimension!(Point3, Vector3, Matrix3, 3);
+dimension!(Point4, Vector4, Matrix4, 4);
+dimension!(Point5, Vector5, Matrix5, 5);
+dimension!(Point6, Vector6, Matrix6, 6);
 
 impl<F: CustomFloat> Derank for Point4<F> {
     type Type = Point3<F>;