Begin porting things to vector-space

This follows the approach taken in Brent Yorgey's diagrams package,
using Conal Elliott's vector-space classes. This allows us to get rid of
the rather unidiomatic and arguably overly-polymorphic
Graphics.Implicit.SaneOperators. Since this has been removed, we can now
begin pulling out the many type annotations strewn about the codebase.

The recently introduced Criterion benchmark confirms that this has no
measurable effect on performance,

                        Before              After
Object 1/STL write      235.7 ± 9.3 ms      236.0 ± 12.8 ms
Object 2/STL write      1.54  ± 0.99 ms     1.44  ± 0.28 us

Note that I make no attempt to enforce any separate between AffineSpaces
and VectorSpaces here. This is work for a later patch.

Additionally, a small performance gain may be had by strictifying and
unpacking the fields of R3 and R2. This, however, would require that a
`data` type be defined and all packing/unpacking of vectors would need
to happen with the appropriate data constructor. Given this is even more
invasive a change than the current patch, this is left as future work.

Conflicts:
	Graphics/Implicit/Export/RayTrace.hs
	Graphics/Implicit/ObjectUtil/GetImplicit3.hs

Graphics/Implicit/Definitions.hs

@@ -1,3 +1,5 @@
+{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-}
+
 -- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
 -- Released under the GNU GPL, see LICENSE
 
@@ -7,6 +9,8 @@ module Graphics.Implicit.Definitions where
 -- we want global IO refs.
 import Data.IORef (IORef, newIORef, readIORef)
 import System.IO.Unsafe (unsafePerformIO)
+import Data.VectorSpace       
+import Control.Applicative       
 
 -- Let's make things a bit nicer. 
 -- Following math notation ℝ, ℝ², ℝ³...
@@ -16,6 +20,21 @@ type ℝ3 = (ℝ,ℝ,ℝ)
 
 type ℕ = Int
 
+-- TODO: Find a better place for this
+(⋅) :: InnerSpace a => a -> a -> Scalar a
+(⋅) = (<.>)
+
+-- TODO: Find a better way to do this?
+class ComponentWiseMultable a where
+    (⋯*) :: a -> a -> a
+    (⋯/) :: a -> a -> a
+instance ComponentWiseMultable ℝ2 where
+    (x,y) ⋯* (x',y') = (x*x', y*y')
+    (x,y) ⋯/ (x',y') = (x/x', y/y')
+instance ComponentWiseMultable ℝ3 where
+    (x,y,z) ⋯* (x',y',z') = (x*x', y*y', z*z')
+    (x,y,z) ⋯/ (x',y',z') = (x/x', y/y', z/z')
+
 -- nxn matrices
 -- eg. M2 ℝ = M₂(ℝ)
 type M2 a = ((a,a),(a,a))

Graphics/Implicit/Export/RayTrace.hs

@@ -3,14 +3,14 @@
 
 module Graphics.Implicit.Export.RayTrace where
 
-import Prelude hiding ((+),(-),(*),(/))
-import qualified Prelude as P
 import Graphics.Implicit.ObjectUtil
 import Graphics.Implicit.Definitions
-import Graphics.Implicit.SaneOperators
 import Graphics.Implicit.Export.Definitions
 import Codec.Picture
 import Control.Monad
+import Data.VectorSpace       
+import Data.AffineSpace       
+import Data.Cross
 
 import Debug.Trace
 
@@ -30,16 +30,12 @@ dynamicImage = ImageRGBA8
 
 -- Math
 
-d a b = norm (b-a)
+d a b = magnitude (b-a)
 
-instance Multiplicative ℝ Color Color where
-	s * (PixelRGBA8 a b c d) = color (s `mult` a) (s `mult` b) (s `mult` c) (d)
-		where 
-			bound = max 0 . min 254
-			mult a b = fromIntegral . round . bound $ a * (fromIntegral b :: ℝ)
-
-instance Multiplicative Color ℝ Color where
-	a * b = b * a
+s `colorMult` (PixelRGBA8 a b c d) = color (s `mult` a) (s `mult` b) (s `mult` c) d
+    where 
+        bound = max 0 . min 254
+        mult a b = fromIntegral . round . bound $ a * fromIntegral b
 
 average :: [Color] -> Color
 average l = 
@@ -57,14 +53,14 @@ average l =
 cameraRay :: Camera -> ℝ2 -> Ray
 cameraRay (Camera p vx vy f) (x,y) =
 	let
-		v  = vx ⨯ vy
-		p' = p + f*v + x*vx + y*vy
-		n  = normalized (p' - p)
+		v  = vx `cross3` vy
+		p' = p ^+^ f*^v ^+^ x*^vx ^+^ y*^vy
+		n  = normalized (p' ^-^ p)
 	in
 		Ray p' n
 
 rayFromTo :: ℝ3 -> ℝ3 -> Ray
-rayFromTo p1 p2 = Ray p1 (normalized $ p2 - p1)
+rayFromTo p1 p2 = Ray p1 (normalized $ p2 ^-^ p1)
 
 rayBounds :: Ray -> (ℝ3, ℝ3) -> ℝ2
 rayBounds ray box =
@@ -90,11 +86,11 @@ intersection r@(Ray p v) ((a, aval),b) res obj =
 			else if aval/(2::ℝ) > res then res/(2 :: ℝ)
 			else                           res/(10 :: ℝ)
 		a'  = a + step
-		a'val = obj (p + a'*v)
+		a'val = obj (p ^+^ a'*^v)
 	in if a'val < 0
 	then 
-		let a'' = refine (a,a') (\s -> obj (p+s*v))
-		in Just (p + a''*v)
+		let a'' = refine (a,a') (\s -> obj (p ^+^ s*^v))
+		in Just (p ^+^ a''*^v)
 	else if a' < b
 	then intersection r ((a',a'val), b) res obj
 	else Nothing
@@ -130,8 +126,8 @@ traceRay :: Ray -> ℝ -> (ℝ3, ℝ3) -> Scene -> Color
 traceRay ray@(Ray cameraP cameraV) step box (Scene obj objColor lights defaultColor) =
 	let
 		(a,b) = rayBounds ray box
-	in case intersection ray ((a, obj (cameraP + a*cameraV)), b) step obj of
-		Just p  -> objColor * (sum $ [0.2] ++ do
+	in case intersection ray ((a, obj (cameraP ^+^ a*^cameraV)), b) step obj of
+		Just p  -> flip colorMult objColor $ (sum $ [0.2] ++ do
 			Light lightPos lightIntensity <- lights
 			let
 				ray'@(Ray _ v) = rayFromTo p lightPos
@@ -141,20 +137,20 @@ traceRay ray@(Ray cameraP cameraV) step box (Scene obj objColor lights defaultCo
 				pval = obj p
 				step = 0.1 :: ℝ
 				dirDeriv :: ℝ3 -> ℝ
-				dirDeriv v = (obj (p + step*v) - pval)/step
+				dirDeriv v = (obj (p ^+^ step*^v) ^-^ pval)/step
 				deriv = (dirDeriv (1,0,0), dirDeriv (0,1,0), dirDeriv (0,0,1))
 				normal = normalized $ deriv
 				unitV = normalized $ v'
-				proj a b = (a⋅b)*b
+				proj a b = (a⋅b)*^b
 				dist  = d p lightPos
-				illumination = (max 0 (normal ⋅ unitV)) * lightIntensity * ((25 :: ℝ)/dist)
+				illumination = (max 0 (normal ⋅ unitV)) * lightIntensity * (25 /dist)
 				rV = 
 					let
 						normalComponent = proj v' normal
 						parComponent    = v' - normalComponent
 					in
 						normalComponent - parComponent	
-			return $ illumination*((3::ℝ) + (0.3::ℝ)*(abs $ rV ⋅ cameraV)^2)
+			return $ illumination*(3  + 0.3*(abs $ rV ⋅ cameraV)^2)
 			)
 		Nothing   -> defaultColor
 
@@ -182,16 +178,16 @@ instance DiscreteAproxable SymbolicObj3 DynamicImage where
 				in 
 					average $ [
 						traceRay 
-							(cameraRay camera ((a,b) + (( 0.25::ℝ)/w, ( 0.25::ℝ)/h)))
+							(cameraRay camera ((a,b) ^+^ ( 0.25/w, 0.25/h)))
 							2 box scene,
 						traceRay 
-							(cameraRay camera ((a,b) + ((-0.25::ℝ)/w, ( 0.25::ℝ)/h)))
+							(cameraRay camera ((a,b) ^+^ (-0.25/w, 0.25/h)))
 							0.5 box scene,
 						traceRay 
-							(cameraRay camera ((a,b) + (( 0.25::ℝ)/w,-( 0.25::ℝ)/h)))
+							(cameraRay camera ((a,b) ^+^ (0.25/w, -0.25/h)))
 							0.5 box scene,
 						traceRay 
-							(cameraRay camera ((a,b) + ((-0.25::ℝ)/w,-( 0.25::ℝ)/h)))
+							(cameraRay camera ((a,b) ^+^ (-0.25/w,-0.25/h)))
 							0.5 box scene
 						]
 
@@ -207,10 +203,9 @@ instance DiscreteAproxable SymbolicObj2 DynamicImage where
 			pixelRenderer :: Int -> Int -> Color
 			pixelRenderer a b = color
 				where
-					xy :: ℝ -> ℝ -> ℝ2
-					xy a b = (x1,y2) - (dxy-dx, dy-dxy)/(2::ℝ) + dxy*(a/w, -b/h)
+					xy a b = ((x1,y2) .-^ (dxy-dx, dy-dxy)^/2) .+^ dxy*^(a/w, -b/h)
 					s = 0.25 :: ℝ
-					(a', b') = (realToFrac a, realToFrac b) :: (ℝ,ℝ)
+					(a', b') = (realToFrac a, realToFrac b)
 					color = average [objColor $ xy a' b', objColor $ xy a' b',
 						objColor $ xy (a'+s) (b'+s),
 						objColor $ xy (a'-s) (b'-s),

Graphics/Implicit/Export/Render/HandleSquares.hs

@@ -5,9 +5,9 @@ module Graphics.Implicit.Export.Render.HandleSquares (mergedSquareTris) where
 
 import Graphics.Implicit.Definitions
 import Graphics.Implicit.Export.Render.Definitions
-import qualified Graphics.Implicit.SaneOperators as S
 import GHC.Exts (groupWith)
 import Data.List (sortBy)
+import Data.VectorSpace       
 
 -- We want small meshes. Essential to this, is getting rid of triangles.
 -- We secifically mark quads in tesselation (refer to Graphics.Implicit.
@@ -119,13 +119,13 @@ joinYaligned [] = []
 -- Reconstruct a triangle
 squareToTri (Sq (b1,b2,b3) z (x1,x2) (y1,y2)) =
 	let
-		zV = b3 S.* z
-		(x1V, x2V) = (x1 S.* b1, x2 S.* b1)
-		(y1V, y2V) = (y1 S.* b2, y2 S.* b2)
-		a = zV S.+ x1V S.+ y1V
-		b = zV S.+ x2V S.+ y1V
-		c = zV S.+ x1V S.+ y2V
-		d = zV S.+ x2V S.+ y2V
+		zV = b3 ^* z
+		(x1V, x2V) = (x1 *^ b1, x2 *^ b1)
+		(y1V, y2V) = (y1 *^ b2, y2 *^ b2)
+		a = zV ^+^ x1V ^+^ y1V
+		b = zV ^+^ x2V ^+^ y1V
+		c = zV ^+^ x1V ^+^ y2V
+		d = zV ^+^ x2V ^+^ y2V
 	in
 		[(a,b,c),(c,b,d)]
 

Graphics/Implicit/Export/Render/RefineSegs.hs

@@ -3,9 +3,8 @@
 
 module Graphics.Implicit.Export.Render.RefineSegs where
 
+import Data.VectorSpace
 import Graphics.Implicit.Definitions
-import qualified Graphics.Implicit.SaneOperators as S
-import Graphics.Implicit.SaneOperators ((⋅), (⨯), norm, normalized)
 
 -- The purpose of refine is to add detail to a polyline aproximating
 -- the boundary of an implicit function and to remove redundant points.
@@ -29,16 +28,16 @@ detail' _ _ a = a
 detail :: Int -> ℝ -> (ℝ2 -> ℝ) -> [ℝ2] -> [ℝ2]
 detail n res obj [p1@(x1,y1), p2@(x2,y2)] | n < 2 =
 	let
-		mid@(midX, midY) = (p1 S.+ p2) S./ (2 :: ℝ)
+		mid@(midX, midY) = (p1 ^+^ p2) ^/ 2
 		midval = obj mid 
 	in if abs midval < res / 40
 	then [(x1,y1), (x2,y2)]
 	else let
-		normal = (\(a,b) -> (b, -a)) $ normalized (p2 S.- p1) 
-		derivN = -(obj (mid S.- (normal S.* (midval/2))) - midval) S.* (2/midval)
+		normal = (\(a,b) -> (b, -a)) $ normalized (p2 ^-^ p1) 
+		derivN = -(obj (mid ^-^ (normal ^* (midval/2))) - midval) * (2/midval)
 	in if abs derivN > 0.5 && abs derivN < 2
 	then let
-		mid' = mid S.- (normal S.* (midval / derivN))
+		mid' = mid ^-^ (normal ^* (midval / derivN))
 	in detail (n+1) res obj [(x1,y1), mid'] 
 	   ++ tail (detail (n+1) res obj [mid', (x2,y2)] )
 	else let
@@ -64,16 +63,16 @@ simplify res = {-simplify3 . simplify2 res . -} simplify1
 
 simplify1 :: [ℝ2] -> [ℝ2]
 simplify1 (a:b:c:xs) =
-	if abs ( ((b S.- a) ⋅ (c S.- a)) - norm (b S.- a) * norm (c S.- a) ) < 0.0001
+	if abs ( ((b ^-^ a) ⋅ (c ^-^ a)) - magnitude (b ^-^ a) * magnitude (c ^-^ a) ) < 0.0001
 	then simplify1 (a:c:xs)
 	else a : simplify1 (b:c:xs)
 simplify1 a = a
 
 {-
 simplify2 :: ℝ -> [ℝ2] -> [ℝ2]
 simplify2 res [a,b,c,d] = 
-	if norm (b S.- c) < res/10
-	then [a, ((b S.+ c) S./ (2::ℝ)), d]
+	if norm (b - c) < res/10
+	then [a, ((b + c) / (2::ℝ)), d]
 	else [a,b,c,d]
 simplify2 _ a = a
 

Graphics/Implicit/Export/Render/TesselateLoops.hs

@@ -5,9 +5,8 @@ module Graphics.Implicit.Export.Render.TesselateLoops (tesselateLoop) where
 
 import Graphics.Implicit.Definitions
 import Graphics.Implicit.Export.Render.Definitions
-import qualified Graphics.Implicit.SaneOperators as S
-import Graphics.Implicit.SaneOperators ((⋅),norm,(⨯),normalized)
-import Debug.Trace
+import Data.VectorSpace
+import Data.Cross       
 
 tesselateLoop :: ℝ -> Obj3 -> [[ℝ3]] -> [TriSquare]
 
@@ -40,11 +39,11 @@ tesselateLoop res obj [as@(_:_:_:_),[_,_], bs@(_:_:_:_), [_,_] ] | length as ==
    #__#
 -}
 
-tesselateLoop res obj [[a,_],[b,_],[c,_],[d,_]] | (a S.+ c) == (b S.+ d) =
+tesselateLoop res obj [[a,_],[b,_],[c,_],[d,_]] | (a + c) == (b + d) =
 	let
-		b1 = normalized $ a S.- b
-		b2 = normalized $ c S.- b
-		b3 = b1 ⨯ b2
+		b1 = normalized $ a ^-^ b
+		b2 = normalized $ c ^-^ b
+		b3 = b1 `cross3` b2
 	in [Sq (b1,b2,b3) (a ⋅ b3) (a ⋅ b1, c ⋅ b1) (a ⋅ b2, c ⋅ b2) ]
 
 {-
@@ -53,7 +52,7 @@ tesselateLoop res obj [[a,_],[b,_],[c,_],[d,_]] | (a S.+ c) == (b S.+ d) =
    #__#      #/_#
 -}
 
-tesselateLoop res obj [[a,_],[b,_],[c,_],[d,_]] | obj ((a S.+ c) S./ (2 :: ℝ)) < res/30 =
+tesselateLoop res obj [[a,_],[b,_],[c,_],[d,_]] | obj ((a ^+^ c) / 2) < res/30 =
 	return $ Tris $ [(a,b,c),(a,c,d)]
 
 -- Fallback case: make fans
@@ -66,14 +65,14 @@ tesselateLoop res obj pathSides = return $ Tris $
 	then early_tris
 	else let
 		len = fromIntegral $ length path :: ℝ
-		mid@(midx,midy,midz) = (foldl1 (S.+) path) S./ len
+		mid@(midx,midy,midz) = (foldl1 (^+^) path) ^/ len
 		midval = obj mid
-		preNormal = foldl1 (S.+) $
-			[ a ⨯ b | (a,b) <- zip path (tail path ++ [head path]) ]
-		preNormalNorm = norm preNormal
-		normal = preNormal S./ preNormalNorm
-		deriv = (obj (mid S.+ (normal S.* (res/100)) ) - midval)/res*100
-		mid' = mid S.- normal S.* (midval/deriv)
+		preNormal = foldl1 (^+^) $
+			[ a `cross3` b | (a,b) <- zip path (tail path ++ [head path]) ]
+		preNormalNorm = magnitude preNormal
+		normal = preNormal ^/ preNormalNorm
+		deriv = (obj (mid ^+^ (normal ^* (res/100)) ) ^-^ midval)/res*100
+		mid' = mid ^-^ normal ^* (midval/deriv)
 	in if abs midval > res/50 && preNormalNorm > 0.5 && abs deriv > 0.5 
 		      && abs (deriv*midval) < 1.1*res && 5*abs (obj mid') < abs midval
 		then early_tris ++ [(a,b,mid') | (a,b) <- zip path (tail path ++ [head path]) ]
@@ -83,14 +82,14 @@ tesselateLoop res obj pathSides = return $ Tris $
 shrinkLoop :: Int -> [ℝ3] -> ℝ -> Obj3 -> ([Triangle], [ℝ3])
 
 shrinkLoop _ path@[a,b,c] res obj =
-	if   abs (obj ((a S.+ b S.+ c) S./ (3::ℝ) )) < res/50
+	if   abs (obj ((a + b + c) / 3 )) < res/50
 	then 
 		( [(a,b,c)], [])
 	else 
 		([], path)
 
 shrinkLoop n path@(a:b:c:xs) res obj | n < length path =
-	if abs (obj ((a S.+ c) S./ (2::ℝ) )) < res/50
+	if abs (obj ((a + c) / 2 )) < res/50
 	then 
 		let (tris,remainder) = shrinkLoop 0 (a:c:xs) res obj
 		in ((a,b,c):tris, remainder)

Graphics/Implicit/Export/Symbolic/Rebound2.hs

@@ -1,12 +1,12 @@
 module Graphics.Implicit.Export.Symbolic.Rebound2 (rebound2) where
 
+import Data.VectorSpace
 import Graphics.Implicit.Definitions
-import qualified Graphics.Implicit.SaneOperators as S
 
 rebound2 :: BoxedObj2 -> BoxedObj2
 rebound2 (obj, (a,b)) = 
 	let
 		d :: ℝ2
-		d = (b S.- a) S./ (10.0 :: ℝ)
+		d = (b ^-^ a) ^/ 10
 	in 
-		(obj, ((a S.- d), (b S.+ d)))
+		(obj, ((a ^-^ d), (b ^+^ d)))

Graphics/Implicit/Export/Symbolic/Rebound3.hs

@@ -1,13 +1,13 @@
 module Graphics.Implicit.Export.Symbolic.Rebound3 (rebound3) where
 
 import Graphics.Implicit.Definitions
-import qualified Graphics.Implicit.SaneOperators as S
+import Data.VectorSpace
 
 rebound3 :: BoxedObj3 -> BoxedObj3
 rebound3 (obj, (a,b)) = 
 	let
 		d :: ℝ3
-		d = (b S.- a) S./ (10.0 :: ℝ)
+		d = (b ^-^ a) ^/ 10
 	in 
-		(obj, ((a S.- d), (b S.+ d)))
+		(obj, ((a ^-^ d), (b ^+^ d)))