dph-examples: pluecker. add matrix/vector mult for rotation
authorAmos Robinson <amos.robinson@gmail.com>
Mon, 3 Dec 2012 07:45:51 +0000 (18:45 +1100)
committerAmos Robinson <amos.robinson@gmail.com>
Mon, 3 Dec 2012 07:45:51 +0000 (18:45 +1100)
dph-examples/examples/spectral/Pluecker/Common.hs
dph-examples/examples/spectral/Pluecker/CommonVectorised.hs
dph-examples/examples/spectral/Pluecker/MainGloss.hs
dph-examples/examples/spectral/Pluecker/Solver.hs
dph-examples/examples/spectral/Pluecker/Vector.hs
dph-examples/examples/spectral/Pluecker/Vectorised.hs

index 599b564..e5cdac6 100644 (file)
@@ -18,7 +18,6 @@ type Triangle = (Vec3, Vec3, Vec3)
 -- Normal and distance of normal from (0,0,0)
 type Plane = (Vec3,Double)
 
-
 {-# INLINE dot #-}
 dot :: Vec3 -> Vec3 -> Double
 dot (u,v,w) (x,y,z) = u*x + v*y + w*z
@@ -96,3 +95,23 @@ lineOnPlane (p,q) (n,d)
 lineOnTriangle :: Line -> Triangle -> Double
 lineOnTriangle p t
  = lineOnPlane p (planeOfTriangle t)
+
+-- matrix mult
+{-# INLINE mvecmul #-}
+mvecmul :: (Vec3,Vec3,Vec3) -> Vec3 -> Vec3
+mvecmul
+   ((xx, xy, xz)
+   ,(yx, yy, yz)
+   ,(zx, zy, zz))
+    (x,   y,  z)
+ = ( x * xx + y * yx + z * zx
+   , x * xy + y * yy + z * zz
+   , x * xz + y * yz + z * zz)
+
+{-# INLINE rotate #-}
+rotate :: Vec3 -> Double -> Vec3
+rotate v d
+ =(mvecmul ((   cos d, 0, sin d)
+           ,(       0, 1,     0)
+           ,(-(sin d), 0, cos d))
+           (v `vsub` (0,0,5))) `vadd` (0,0,15)
index 3cabe2a..d0a6465 100644 (file)
@@ -107,3 +107,22 @@ lineOnTriangle :: Line -> Triangle -> Double
 lineOnTriangle p t
  = lineOnPlane p (planeOfTriangle t)
 
+-- matrix mult
+{-# INLINE mvecmul #-}
+mvecmul :: (Vec3,Vec3,Vec3) -> Vec3 -> Vec3
+mvecmul
+   ((xx, xy, xz)
+   ,(yx, yy, yz)
+   ,(zx, zy, zz))
+    (x,   y,  z)
+ = ( x * xx + y * yx + z * zx
+   , x * xy + y * yy + z * zz
+   , x * xz + y * yz + z * zz)
+
+{-# INLINE rotate #-}
+rotate :: Vec3 -> Double -> Vec3
+rotate v d
+ =(mvecmul ((   cos d, 0, sin d)
+           ,(       0, 1,     0)
+           ,(negate (sin d), 0, cos d))
+           (v `vsub` (0,0,5))) `vadd` (0,0,15)
index 9845406..f962d34 100644 (file)
@@ -21,8 +21,8 @@ model numverts numtris
    return (vs,ts)
  where
    mkvec _ = do
-     x <- randomRIO (-30,30)
-     y <- randomRIO (-30,30)
+     x <- randomRIO (-10,10)
+     y <- randomRIO (-10,10)
      z <- randomRIO (0,10)
      return (x,y,z)
    mktri _ = do
@@ -67,13 +67,13 @@ mainGloss
         
 mainGloss v t numr solver windowSize
  = let  mkray _
-         = do   x <- randomRIO (-10,10)
-                y <- randomRIO (-10,10)
+         = do   x <- randomRIO (-3,3)
+                y <- randomRIO (-3,3)
                 return (x,y,1)
             
-        draw _
+        draw time
          = do   rays <- VU.generateM numr mkray
-                let pts = solver v t rays
+                let pts = solver v t rays (realToFrac time)
                 return $ Pictures $ map drawPoint $ VU.toList pts
 
    in   animateIO
@@ -85,11 +85,15 @@ mainGloss v t numr solver windowSize
 
 
 
-pointSize = 20
+pointSize = 15
 
 drawPoint :: (Vec3, Double) -> Picture
 drawPoint ((x,y,z), d)
-       = Translate (realToFrac x * 50) (realToFrac y * 50) 
-    $ Color (makeColor d' d' d' 1)
-       $ ThickCircle (pointSize / 2) pointSize
- where d' = 1 - (realToFrac $ d / 10)
+       = Translate (realToFrac x * 200) (realToFrac y * 200) 
+    $ Color (makeColor d' d' d' 0.5)
+    $ Polygon [ (-pointSize, -pointSize)
+              , ( pointSize, -pointSize)
+              , ( pointSize,  pointSize)
+              , (-pointSize,  pointSize) ]
+       -- $ ThickCircle (pointSize / 2) pointSize
+ where d' = 1 - (realToFrac $ d / 20)
index 0a17a8a..df88e25 100644 (file)
@@ -11,13 +11,13 @@ import qualified Data.Array.Parallel        as P
 import qualified Data.Array.Parallel.PArray    as P
 
 
-type Solver = VU.Vector Vec3 -> VU.Vector (Int,Int,Int) -> VU.Vector Vec3 -> VU.Vector (Vec3, Double)
+type Solver = VU.Vector Vec3 -> VU.Vector (Int,Int,Int) -> VU.Vector Vec3 -> Double -> VU.Vector (Vec3, Double)
 
 solvers :: [(String,Solver)]
 solvers =
  [("vectorised", solverPA)
  ,("vector",     SV.solveV)]
 
-solverPA verts tris rays
+solverPA verts tris rays time
  = let fv a = P.fromVector (VU.convert a)
-   in  VU.convert (P.toVector (SPA.solvePA (fv verts) (fv tris) (fv rays)))
+   in  VU.convert (P.toVector (SPA.solvePA (fv verts) (fv tris) (fv rays) time))
index d600ff2..b196723 100644 (file)
@@ -12,29 +12,35 @@ solveV
     :: Vector Vec3             -- ^ vertices of the surface
     -> Vector (Int,Int,Int)    -- ^ triangles, each 3 vertex indices
     -> Vector Vec3             -- ^ rays to cast
+    -> Double                  -- ^ time
     -> Vector (Vec3,Double)    -- ^ rays and their distance
-solveV vertices triangles rays
+solveV vertices triangles rays time
  = map cast' rays
  where
-  cast' = cast vertices triangles
+  cast' = cast vertices triangles time
 
 
 cast 
     :: Vector Vec3          -- ^ vertices of the surface
     -> Vector (Int,Int,Int) -- ^ triangles, each 3 vertex indices
+    -> Double               -- ^ time
     -> Vec3                 -- ^ ray
     -> (Vec3,Double)
-cast vertices triangles ray
+cast vertices triangles time ray
  = let r' = ((0,0,0), ray)
        pl = plueckerOfLine r'
        mi = minimum
-          $ map (\t -> check r' pl $ tri vertices t) triangles
+          $ filter (>0)
+          $ map (\t -> check r' pl $ tri vertices t time) triangles
    in  (ray, mi)
 
 check r pl t
   | inside pl t = lineOnTriangle r t
   | otherwise = 1e100
 
-tri :: Vector Vec3 -> (Int,Int,Int) -> Triangle
-tri v (a,b,c) = (v ! a, v ! b, v ! c)
+tri :: Vector Vec3 -> (Int,Int,Int) -> Double -> Triangle
+tri v (a,b,c) time = (get a, get b, get c)
+ where
+  {-# INLINE get #-}
+  get i = rotate (v `unsafeIndex` i) (time / 4)
 
index 29d4d77..3033bc9 100644 (file)
@@ -19,36 +19,42 @@ solvePA
     :: PArray Vec3           -- ^ vertices of the surface
     -> PArray (Int,Int,Int)  -- ^ triangles, each 3 vertex indices
     -> PArray Vec3           -- ^ rays to cast
+    -> Double                -- ^ time, used for rotation
     -> PArray (Vec3,Double)  -- ^ rays and their distance
-solvePA vertices triangles rays
- = toPArrayP (solveV (fromPArrayP vertices) (fromPArrayP triangles) (fromPArrayP rays))
+solvePA vertices triangles rays time
+ = toPArrayP (solveV (fromPArrayP vertices) (fromPArrayP triangles) (fromPArrayP rays) time)
 
 solveV 
     :: [:Vec3:]             -- ^ vertices of the surface
     -> [:(Int,Int,Int):]    -- ^ triangles, each 3 vertex indices
     -> [:Vec3:]             -- ^ rays to cast
+    -> Double               -- ^ time
     -> [:(Vec3,Double):]    -- ^ rays and their distance
-solveV vertices triangles rays
+solveV vertices triangles rays time
  = mapP cast' rays
  where
-  cast' = cast vertices triangles
+  cast' = cast vertices triangles time
 
 
 cast 
     :: [:Vec3:]          -- ^ vertices of the surface
     -> [:(Int,Int,Int):] -- ^ triangles, each 3 vertex indices
+    -> Double            -- ^ time
     -> Vec3              -- ^ ray
     -> (Vec3,Double)
-cast vertices triangles ray
+cast vertices triangles time ray
  = let r' = ((0,0,0), ray)
        pl = plueckerOfLine r'
        mi = minimumP
-            (mapP (\t -> check r' pl (tri vertices t)) triangles)
+            (mapP (\t -> check r' pl (tri vertices t time)) triangles)
    in  (ray, mi)
 
 check r pl t
   | inside pl t = lineOnTriangle r t
   | otherwise = 1e100
 
-tri :: [:Vec3:] -> (Int,Int,Int) -> Triangle
-tri v (a,b,c) = (v !: a, v !: b, v !: c)
+tri :: [:Vec3:] -> (Int,Int,Int) -> Double -> Triangle
+tri v (a,b,c) time = (get a, get b, get c)
+ where
+  {-# INLINE get #-}
+  get i = rotate (v !: i) (time / 4)