From: Amos Robinson Date: Mon, 3 Dec 2012 07:45:51 +0000 (+1100) Subject: dph-examples: pluecker. add matrix/vector mult for rotation X-Git-Tag: ghc-7.8.1-release~40 X-Git-Url: http://git.haskell.org/packages/dph.git/commitdiff_plain/1de8c9da821df59bec37145b10d735467ba16f4d dph-examples: pluecker. add matrix/vector mult for rotation --- diff --git a/dph-examples/examples/spectral/Pluecker/Common.hs b/dph-examples/examples/spectral/Pluecker/Common.hs index 599b564..e5cdac6 100644 --- a/dph-examples/examples/spectral/Pluecker/Common.hs +++ b/dph-examples/examples/spectral/Pluecker/Common.hs @@ -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) diff --git a/dph-examples/examples/spectral/Pluecker/CommonVectorised.hs b/dph-examples/examples/spectral/Pluecker/CommonVectorised.hs index 3cabe2a..d0a6465 100644 --- a/dph-examples/examples/spectral/Pluecker/CommonVectorised.hs +++ b/dph-examples/examples/spectral/Pluecker/CommonVectorised.hs @@ -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) diff --git a/dph-examples/examples/spectral/Pluecker/MainGloss.hs b/dph-examples/examples/spectral/Pluecker/MainGloss.hs index 9845406..f962d34 100644 --- a/dph-examples/examples/spectral/Pluecker/MainGloss.hs +++ b/dph-examples/examples/spectral/Pluecker/MainGloss.hs @@ -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) diff --git a/dph-examples/examples/spectral/Pluecker/Solver.hs b/dph-examples/examples/spectral/Pluecker/Solver.hs index 0a17a8a..df88e25 100644 --- a/dph-examples/examples/spectral/Pluecker/Solver.hs +++ b/dph-examples/examples/spectral/Pluecker/Solver.hs @@ -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)) diff --git a/dph-examples/examples/spectral/Pluecker/Vector.hs b/dph-examples/examples/spectral/Pluecker/Vector.hs index d600ff2..b196723 100644 --- a/dph-examples/examples/spectral/Pluecker/Vector.hs +++ b/dph-examples/examples/spectral/Pluecker/Vector.hs @@ -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) diff --git a/dph-examples/examples/spectral/Pluecker/Vectorised.hs b/dph-examples/examples/spectral/Pluecker/Vectorised.hs index 29d4d77..3033bc9 100644 --- a/dph-examples/examples/spectral/Pluecker/Vectorised.hs +++ b/dph-examples/examples/spectral/Pluecker/Vectorised.hs @@ -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)