packages feed

vect 0.4.0 → 0.4.5

raw patch · 23 files changed

+1916/−518 lines, 23 files

Files

Data/Vect/Double/Base.hs view
@@ -1,11 +1,34 @@ {-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} -module Data.Vect.Flt.Base where+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-} +module Data.Vect.Flt.Base+  ( AbelianGroup(..) , vecSum+  , MultSemiGroup(..) , Ring , semigroupProduct+  , LeftModule(..) , RightModule(..)+  , Vector(..) , DotProd(..) , CrossProd(..)+  , normalize , distance , angle , angle'+  , UnitVector(..)+  , Pointwise(..)+  , Extend(..) , HasCoordinates(..) , Dimension(..)+  , Matrix(..) , Tensor(..) , Diagonal (..) , Determinant(..)+  , Orthogonal(..) , Projective(..) , MatrixNorms(..)+  , Vec2(..) , Vec3(..) , Vec4(..)+  , Mat2(..) , Mat3(..) , Mat4(..)+  , Ortho2 , Ortho3 , Ortho4+  , Normal2 , Normal3 , Normal4+  , Proj3 , Proj4+  , mkVec2 , mkVec3 , mkVec4+  , project , project' , projectUnsafe , flipNormal+  , householder, householderOrtho+  )+  where+ import Control.Monad import System.Random   import Foreign +-------------------------------------------------------------------------------- -- class declarations  class AbelianGroup g where@@ -20,15 +43,17 @@ vecSum :: AbelianGroup g => [g] -> g vecSum l = foldl (&+) zero l  -class (AbelianGroup r) => -      Ring r where+class MultSemiGroup r where   (.*.) :: r -> r -> r   one   :: r +class (AbelianGroup r, MultSemiGroup r) => Ring r + infixl 7 .*.  -ringProduct :: Ring r => [r] -> r-ringProduct l = foldl (.*.) one l+-- was: ringProduct :: Ring r => [r] -> r+semigroupProduct :: MultSemiGroup r => [r] -> r +semigroupProduct l = foldl (.*.) one l  class LeftModule r m where   lmul :: r -> m -> m@@ -75,7 +100,7 @@   lensqr = normsqr   dotprod :: v -> v -> Flt   normsqr v = (v &. v)  -  norm = sqrt.lensqr+  norm = sqrt . lensqr   dotprod = (&.)  infix 7 &.@@ -112,26 +137,28 @@   fromNormalRadius :: Flt -> u -> v   fromNormalRadius t n = t *& fromNormal n  --- | projects the first vector onto the direction of the second (unit) vector+-- | Projects the first vector onto the direction of the second (unit) vector project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v project' what dir = projectUnsafe what (fromNormal dir) --- | direction (second argument) is assumed to be a /unit/ vector!+-- | Direction (second argument) is assumed to be a /unit/ vector! projectUnsafe :: (Vector v, DotProd v) => v -> v -> v projectUnsafe what dir = what &- dir &* (what &. dir)  project :: (Vector v, DotProd v) => v -> v -> v project what dir = what &- dir &* ((what &. dir) / (dir &. dir)) --- | since unit vectors are not a group, we need a separate function.+-- | Since unit vectors are not a group, we need a separate function. flipNormal :: UnitVector v n => n -> n  flipNormal = toNormalUnsafe . neg . fromNormal  +-- | Cross product class CrossProd v where   crossprod :: v -> v -> v   (&^)      :: v -> v -> v   (&^) = crossprod-  + +-- | Pointwise multiplication  class Pointwise v where   pointwise :: v -> v -> v   (&!)      :: v -> v -> v@@ -166,13 +193,50 @@ "inverse is an involution"    forall m. inverse (inverse m) = m   #-}   +class Matrix m => Orthogonal m o | m->o, o->m where  +  fromOrtho     :: o -> m +  toOrthoUnsafe :: m -> o+  +class (AbelianGroup m, Matrix m) => MatrixNorms m where+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices+  operatorNorm   :: m -> Flt       -- ^ (euclidean) operator norm (not implemented yet)+  matrixDistance m n = frobeniusNorm (n &- m)+  operatorNorm = error "operatorNorm: not implemented yet"+   -- | Outer product (could be unified with Diagonal?) class Tensor t v | t->v where   outer :: v -> v -> t      class Determinant m where   det :: m -> Flt    -    ++class Dimension a where+  dim :: a -> Int+     +-- | Householder matrix, see <http://en.wikipedia.org/wiki/Householder_transformation>.  +-- In plain words, it is the reflection to the hyperplane orthogonal to the input vector.+householder :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m+householder u = idmtx &- (2 *& outer v v) +  where v = fromNormal u++householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o+householderOrtho = toOrthoUnsafe . householder++-- | \"Projective\" matrices have the following form: the top left corner+-- is an any matrix, the bottom right corner is 1, and the top-right+-- column is zero. These describe the affine orthogonal transformation of+-- the space one dimension less.+class (Vector v, Orthogonal n o, Diagonal v n) => Projective v n o m p +    | m->p, p->m, p->o, o->p, p->n, n->p, p->v, v->p, n->o, n->v, v->n where+  fromProjective     :: p -> m+  toProjectiveUnsafe :: m -> p+  orthogonal         :: o -> p+  linear             :: n -> p+  translation        :: v -> p+  scaling            :: v -> p++-------------------------------------------------------------------------------- -- Vec / Mat datatypes   data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt @@ -182,16 +246,16 @@ data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt    deriving (Read,Show) --- | these are /row/ vectors +-- | The components are /row/ vectors  data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show) data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show) data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)  -- | The assumption when dealing with these is always that they are of unit length. -- Also, interpolation works differently.-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) -newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) -newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension)   mkVec2 :: (Flt,Flt) -> Vec2 mkVec3 :: (Flt,Flt,Flt) -> Vec3@@ -201,6 +265,19 @@ mkVec3 (x,y,z)   = Vec3 x y z mkVec4 (x,y,z,w) = Vec4 x y z w +-- | Orthogonal matrices.+--+-- Note: the "Random" instances generates orthogonal matrices with determinant 1+-- (that is, orientation-preserving orthogonal transformations)!+newtype Ortho2 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)++-- | Projective matrices, encoding affine transformations in dimension one less.+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)++-------------------------------------------------------------------------------- -- Unit vectors    instance UnitVector Vec2 Normal2 where@@ -218,47 +295,163 @@   fromNormal (Normal4 v) = v    toNormalUnsafe = Normal4 -rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)-rndUnit g = +_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+_rndUnit g =    if d > 0.01     then ( v &* (1.0/d) , h )-    else rndUnit h+    else _rndUnit h   where     (v,h) = random g     d = norm v      instance Random Normal2 where-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal2 v, h)     randomR _ = random  instance Random Normal3 where-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal3 v, h)     randomR _ = random  instance Random Normal4 where-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal4 v, h)     randomR _ = random -{--instance Storable Normal2 where-  alignment _ = alignment (undefined::Vec2)-  sizeOf    _ = sizeOf    (undefined::Vec2)-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)-  poke p (Normal2 v) = poke (castPtr p) v  - -instance Storable Normal3 where-  alignment _ = alignment (undefined::Vec3)-  sizeOf    _ = sizeOf    (undefined::Vec3)-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)-  poke p (Normal3 v) = poke (castPtr p) v  +instance CrossProd Normal3 where+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w) -instance Storable Normal4 where-  alignment _ = alignment (undefined::Vec4)-  sizeOf    _ = sizeOf    (undefined::Vec4)-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)-  poke p (Normal4 v) = poke (castPtr p) v  --}+--------------------------------------------------------------------------------+-- Orthogonal matrices +instance Orthogonal Mat2 Ortho2 where+  fromOrtho (Ortho2 o) = o+  toOrthoUnsafe = Ortho2++instance Orthogonal Mat3 Ortho3 where+  fromOrtho (Ortho3 o) = o+  toOrthoUnsafe = Ortho3 ++instance Orthogonal Mat4 Ortho4 where+  fromOrtho (Ortho4 o) = o+  toOrthoUnsafe = Ortho4++------++instance Matrix Ortho2 where+  transpose (Ortho2 o) = Ortho2 (transpose o)+  idmtx = Ortho2 idmtx+  inverse = transpose++instance Matrix Ortho3 where+  transpose (Ortho3 o) = Ortho3 (transpose o)+  idmtx = Ortho3 idmtx+  inverse = transpose++instance Matrix Ortho4 where+  transpose (Ortho4 o) = Ortho4 (transpose o)+  idmtx = Ortho4 idmtx+  inverse = transpose++------++instance Random Ortho2 where+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)+  randomR _ = random++instance Random Ortho3 where+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)+  randomR _ = random++instance Random Ortho4 where+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)+  randomR _ = random++------++-- determinant will be -1+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)+_rndOrtho2 g = (h2, g1) where+  h2 = householder u2 :: Mat2 +  (u2,g1) = random g   ++-- generates a uniformly random orthogonal 3x3 matrix +-- /with determinant +1/, with respect to the Haar measure of SO3.+--+-- see Theorem 4 in:+-- Francesco Mezzadri: How to Generate Random Matrices from the Classical Compact Groups +-- Notices of the AMS, May 2007 issue+-- <http://www.ams.org/notices/200705/fea-mezzadri-web.ps>+_rndOrtho3 :: RandomGen g => g -> (Mat3, g) +_rndOrtho3 g = ( (h3 .*. m3), g2) where+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 +  h3 = householder u3 :: Mat3+  (u3,g1) = random g+  (o2,g2) = _rndOrtho2 g1++-- determinant will be -1+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) +_rndOrtho4 g = ( (h4 .*. m4), g2) where+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 +  h4 = householder u4 :: Mat4+  (u4,g1) = random g+  (o3,g2) = _rndOrtho3 g1++------++_flip1stRow2 :: Mat2 -> Mat2+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b++_flip1stRow3 :: Mat3 -> Mat3+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c++_flip1stRow4 :: Mat4 -> Mat4+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d++--------------------------------------------------------------------------------+-- projective matrices+  +instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where+  fromProjective (Proj3 m) = m+  toProjectiveUnsafe = Proj3+  orthogonal = Proj3 . extendWith 1 . fromOrtho+  linear     = Proj3 . extendWith 1+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)+  scaling     v = Proj3 $ diag (extendWith 1 v)+  +instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where+  fromProjective (Proj4 m) = m+  toProjectiveUnsafe = Proj4+  orthogonal = Proj4 . extendWith 1 . fromOrtho +  linear     = Proj4 . extendWith 1+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)+  scaling     v = Proj4 $ diag (extendWith 1 v)++instance Matrix Proj3 where+  idmtx = Proj3 idmtx+  transpose (Proj3 m) = Proj3 (transpose m)+  inverse = _invertProj3++instance Matrix Proj4 where+  idmtx = Proj4 idmtx+  transpose (Proj4 m) = Proj4 (transpose m)+  inverse = _invertProj4++_invertProj3 :: Proj3 -> Proj3+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = +  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec2) .* invm2 +    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)++-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.+-- We assume that the bottom-right corner is 1.+_invertProj4 :: Proj4 -> Proj4+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = +  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec3) .* invm3 +    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)++-------------------------------------------------------------------------------- -- Vec2 instances  instance HasCoordinates Vec2 Flt where@@ -314,7 +507,10 @@         k = sizeOf (undefined::Flt)     poke        p   x     pokeByteOff p k y-               ++instance Dimension Vec2 where dim _ = 2++--------------------------------------------------------------------------------                     -- Mat2 instances  instance HasCoordinates Mat2 Vec2 where@@ -336,19 +532,21 @@   (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)   (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)   neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  -  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)-+  zero = Mat2 zero zero  +   instance Vector Mat2 where   scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s   mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f -instance Ring Mat2 where+instance MultSemiGroup Mat2 where   (.*.) (Mat2 r1 r2) n =      let (Mat2 c1 c2) = transpose n     in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))             (Vec2 (r2 &. c1) (r2 &. c2))   one = idmtx  +instance Ring Mat2+ instance LeftModule Mat2 Vec2 where   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v)    @@ -362,11 +560,6 @@   outer (Vec2 a b) (Vec2 x y) = Mat2     (Vec2 (a*x) (a*y))     (Vec2 (b*x) (b*y))-{--  outer v w = -    let full = Mat2 (Vec2 1 1) (Vec2 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat2 where   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c @@ -393,6 +586,24 @@     poke        p   r1     pokeByteOff p k r2 +instance Random Mat2 where+  random = randomR (Mat2 v1 v1 , Mat2 v2 v2) where +    v1 = Vec2 (-1) (-1) +    v2 = Vec2   1    1+  randomR (Mat2 a b, Mat2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Mat2 x y, gen2)+          +instance Dimension Mat2 where dim _ = 2+     +instance MatrixNorms Mat2 where +  frobeniusNorm (Mat2 r1 r2) =  +    sqrt $+      normsqr r1 + +      normsqr r2+     +--------------------------------------------------------------------------------      -- Vec3 instances  instance HasCoordinates Vec3 Flt where@@ -454,7 +665,10 @@     poke        p       x     pokeByteOff p (k  ) y     pokeByteOff p (k+k) z-   ++instance Dimension Vec3 where dim _ = 3++--------------------------------------------------------------------------------    -- Mat3 instances  instance HasCoordinates Mat3 Vec3 where@@ -495,13 +709,13 @@   (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)   (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)   neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) -  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)+  zero = Mat3 zero zero zero   instance Vector Mat3 where   scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s   mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f -instance Ring Mat3 where+instance MultSemiGroup Mat3 where   (.*.) (Mat3 r1 r2 r3) n =      let (Mat3 c1 c2 c3) = transpose n     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))@@ -509,6 +723,8 @@             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))   one = idmtx  +instance Ring Mat3+ instance LeftModule Mat3 Vec3 where   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)   @@ -523,11 +739,6 @@     (Vec3 (a*x) (a*y) (a*z))     (Vec3 (b*x) (b*y) (b*z))     (Vec3 (c*x) (c*y) (c*z))-{--  outer v w = -    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat3 where   det (Mat3 r1 r2 r3) = det (r1,r2,r3)@@ -556,6 +767,26 @@     pokeByteOff p (k  ) r2     pokeByteOff p (k+k) r3 +instance Random Mat3 where+  random = randomR (Mat3 v1 v1 v1 , Mat3 v2 v2 v2) where+    v1 = Vec3 (-1) (-1) (-1)+    v2 = Vec3   1    1    1+  randomR (Mat3 a b c, Mat3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Mat3 x y z, gen3)+   +instance Dimension Mat3 where dim _ = 3+  +instance MatrixNorms Mat3 where +  frobeniusNorm (Mat3 r1 r2 r3)  = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 +    +-------------------------------------------------------------------------------- -- Vec4 instances  instance HasCoordinates Vec4 Flt where@@ -615,6 +846,9 @@     pokeByteOff p (k+k) z     pokeByteOff p (3*k) w +instance Dimension Vec4 where dim _ = 4++-------------------------------------------------------------------------------- -- Mat4 instances  instance HasCoordinates Mat4 Vec4 where@@ -642,7 +876,7 @@   scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s   mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f -instance Ring Mat4 where+instance MultSemiGroup Mat4 where   (.*.) (Mat4 r1 r2 r3 r4) n =      let (Mat4 c1 c2 c3 c4) = transpose n     in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))@@ -651,6 +885,8 @@             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))   one = idmtx  +instance Ring Mat4+ instance LeftModule Mat4 Vec4 where   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)   @@ -666,14 +902,10 @@     (Vec4 (b*x) (b*y) (b*z) (b*w))     (Vec4 (c*x) (c*y) (c*z) (c*w))     (Vec4 (d*x) (d*y) (d*z) (d*w))-{--  outer v w = -    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)-    in  (diag v) .*. full .*. (diag w)--} ---instance Determinant Mat4 where---  det (Mat4 r1 r2 r3 r4)+instance Determinant Mat4 where+  det = error "det/Mat4: not implemented yet" +  -- det (Mat4 r1 r2 r3 r4) =   {- instance Show Mat4 where@@ -701,6 +933,28 @@     pokeByteOff p (k+k) r3     pokeByteOff p (3*k) r4 +instance Random Mat4 where+  random = randomR (Mat4 v1 v1 v1 v1, Mat4 v2 v2 v2 v2) where+    v1 = Vec4 (-1) (-1) (-1) (-1)+    v2 = Vec4   1    1    1    1+  randomR (Mat4 a b c d, Mat4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Mat4 x y z w, gen4)+    +instance Dimension Mat4 where dim _ = 4+   +instance MatrixNorms Mat4 where +  frobeniusNorm (Mat4 r1 r2 r3 r4) = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 + +      normsqr r4  +    +-------------------------------------------------------------------------------- -- Extend instances  instance Extend Vec2 Vec3 where@@ -719,17 +973,20 @@   trim (Vec4 x y z _)       = Vec3 x y z  instance Extend Mat2 Mat3 where-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith w (Mat2 p q) = Mat3 (extendZero p) (extendZero q) (Vec3 0 0 w)   trim (Mat3 p q _) = Mat2 (trim p) (trim q)  instance Extend Mat2 Mat4 where-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith w (Mat2 p q) = Mat4 (extendZero p) (extendZero q) (Vec4 0 0 w 0) (Vec4 0 0 0 w)   trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)  instance Extend Mat3 Mat4 where-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith w (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) (Vec4 0 0 0 w)   trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+  +--------------------------------------------------------------------------------+
Data/Vect/Double/GramSchmidt.hs view
@@ -10,7 +10,7 @@  import Data.Vect.Flt.Base --------------------------------------------------------+--------------------------------------------------------------------------------  liftPair :: (a -> b) -> (a,a) -> (b,b) liftPair f (x,y) = (f x, f y)@@ -21,7 +21,7 @@ liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b) liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w) --------------------------------------------------------+--------------------------------------------------------------------------------      -- | produces orthogonal\/orthonormal vectors from a set of vectors     class GramSchmidt a where@@ -33,7 +33,7 @@ "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a    #-} --------------------------------------------------------+--------------------------------------------------------------------------------  instance GramSchmidt (Vec2,Vec2) where   gramSchmidt = gramSchmidtPair
Data/Vect/Double/OpenGL.hs view
@@ -3,21 +3,92 @@ -- TODO: the pointer versions of these functions should be really implemented  -- via the pointer versions of the original opengl functions... --- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.   module Data.Vect.Flt.OpenGL where  import Control.Monad import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Projective import qualified Graphics.Rendering.OpenGL as GL  import Foreign -import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)+import Graphics.Rendering.OpenGL hiding +  ( Normal3 , rotate , translate , scale+  , matrix , currentMatrix , withMatrix , multMatrix +  ) --------------------------------------------------------+-------------------------------------------------------------------------------- -{-# SPECIALISE radianToDegrees :: Float -> Float #-}+-- | There should be a big warning here about the different conventions, +-- hidden transpositions, and all the confusion this will inevitably cause...+--+-- As it stands, +--+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3+-- +-- has the same result as+--+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)+--+-- because at the interface of OpenGL and this library there is a transposition+-- to compensate for the different conventions. (This transposition is implicit+-- in the code, because the way the matrices are stored in the memory is also+-- different: OpenGL stores them column-major, and we store them row-major).++class ToOpenGLMatrix m where+  makeGLMatrix :: m -> IO (GLmatrix Flt)++class FromOpenGLMatrix m where+  peekGLMatrix :: GLmatrix Flt -> IO m+  +setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x+ +getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix++matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)++currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m+currentMatrix = matrix Nothing++multMatrix :: ToOpenGLMatrix m => m -> IO ()+multMatrix m = makeGLMatrix m >>= GL.multMatrix++instance ToOpenGLMatrix Mat4 where+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) + +instance FromOpenGLMatrix Mat4 where+  -- huh? GL.withMatrix is strange+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)+  +instance ToOpenGLMatrix Mat3 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)+ +instance ToOpenGLMatrix Mat2 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)++instance ToOpenGLMatrix Ortho4 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)++instance ToOpenGLMatrix Ortho3 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)++instance ToOpenGLMatrix Ortho2 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)++instance ToOpenGLMatrix Proj4 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)++instance ToOpenGLMatrix Proj3 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)+  +--------------------------------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float  -> Float  #-} {-# SPECIALISE radianToDegrees :: Double -> Double #-} radianToDegrees :: RealFrac a => a -> a radianToDegrees x = x * 57.295779513082322@@ -28,20 +99,61 @@ degreesToRadian x = x * 1.7453292519943295e-2  -- | The angle is in radians. (WARNING: OpenGL uses degrees!)-rotate :: Flt -> Vec3 -> IO ()-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)+glRotate :: Flt -> Vec3 -> IO ()+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z) -translate :: Vec3 -> IO ()-translate (Vec3 x y z) = GL.translate (Vector3 x y z)+glTranslate :: Vec3 -> IO ()+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z) -scale3 :: Vec3 -> IO ()-scale3 (Vec3 x y z) = GL.scale x y z+glScale3 :: Vec3 -> IO ()+glScale3 (Vec3 x y z) = GL.scale x y z -scale :: Flt -> IO ()-scale x = GL.scale x x x+glScale :: Flt -> IO ()+glScale x = GL.scale x x x --------------------------------------------------------+--------------------------------------------------------------------------------+ +-- | \"Orthogonal projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1)+orthoMatrix +  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +orthoMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2/(r-l)) 0 0 0)+  (Vec4 0 (2/(t-b)) 0 0)+  (Vec4 0 0 (-2/(f-n)) 0)+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)+  +-- | The same as "orthoMatrix", but with a different parametrization.+orthoMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f) +-- | \"Perspective projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1).+frustumMatrix+  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +frustumMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2*n/(r-l)) 0 0 0)+  (Vec4 0 (2*n/(t-b)) 0 0)+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))+  (Vec4 0 0 (-2*f*n*(f-n)) 0)+  +-- | The same as "frustumMatrix", but with a different parametrization.+frustumMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)++-------------------------------------------------------------------------------- -- Vertex instances  instance GL.Vertex Vec2 where@@ -56,17 +168,20 @@   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)   vertexv p = peek p >>= vertex    ---------------------------------------------------------+-------------------------------------------------------------------------------- -- the Normal instance -- note that there is no Normal2\/Normal4 in the OpenGL binding  instance GL.Normal Normal3 where-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normal u = GL.normal (GL.Normal3 x y z) +    where Vec3 x y z = fromNormal u    normalv p = peek p >>= normal  --------------------------------------------------------+instance GL.Normal Vec3 where+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) +  normalv p = peek p >>= normal  +-------------------------------------------------------------------------------- -- Color instances    instance GL.Color Vec3 where@@ -88,8 +203,7 @@   secondaryColorv p = peek p >>= secondaryColor -} ---------------------------------------------------------+-------------------------------------------------------------------------------- -- TexCoord instances  instance GL.TexCoord Vec2 where@@ -110,8 +224,7 @@   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)   multiTexCoordv unit p = peek p >>= multiTexCoord unit ---------------------------------------------------------    +-------------------------------------------------------------------------------- -- Vertex Attributes (experimental)  class VertexAttrib' a where@@ -130,16 +243,18 @@   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w   instance VertexAttrib' Normal2 where-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y+  vertexAttrib loc u = GL.vertexAttrib2 loc x y+    where Vec2 x y = fromNormal u   instance VertexAttrib' Normal3 where-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z+    where Vec3 x y z = fromNormal u   instance VertexAttrib' Normal4 where-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w----------------------------------------------------------+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w+    where Vec4 x y z w = fromNormal u +   +-------------------------------------------------------------------------------- -- Uniform (again, experimental)  -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; @@ -181,3 +296,5 @@   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)      #endif++
Data/Vect/Double/Util/Dim2.hs view
@@ -4,12 +4,13 @@  import Data.Vect.Flt.Base --- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@. structVec2 :: [Flt] -> [Vec2] structVec2 [] = [] structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls)  structVec2 _ = error "structVec2" +-- | The opposite of "structVec2". destructVec2 :: [Vec2] -> [Flt] destructVec2 [] = [] destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  @@ -49,10 +50,13 @@ angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt angle2' = angle2 . fromNormal --- |Rotation matrix by a given angle (in radians), counterclockwise.+-- | Rotation matrix by a given angle (in radians), counterclockwise. rotMatrix2 :: Flt -> Mat2 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a +rotMatrixOrtho2 :: Flt -> Ortho2+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2+ rotate2 :: Flt -> Vec2 -> Vec2 rotate2 a v = v .* (rotMatrix2 a)  @@ -63,3 +67,5 @@ -- |Rotates clockwise by 90 degrees. rotateCW :: Vec2 -> Vec2 rotateCW (Vec2 x y) = Vec2 y (-x)++
Data/Vect/Double/Util/Dim3.hs view
@@ -4,18 +4,26 @@  import Data.Vect.Flt.Base +--------------------------------------------------------------------------------++-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@. structVec3 :: [Flt] -> [Vec3] structVec3 [] = [] structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls)  structVec3 _ = error "structVec3" +-- | The opposite of "structVec3". destructVec3 :: [Vec3] -> [Flt] destructVec3 [] = [] destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)   +--------------------------------------------------------------------------------+ det3 :: Vec3 -> Vec3 -> Vec3 -> Flt det3 u v w = det (u,v,w) +--------------------------------------------------------------------------------+ translate3X :: Flt -> Vec3 -> Vec3 translate3Y :: Flt -> Vec3 -> Vec3 translate3Z :: Flt -> Vec3 -> Vec3@@ -43,26 +51,34 @@ rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a -rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)-         -> Normal3   -- ^ axis (should be a /unit/ vector!) -         -> Vec3      -- ^ vector-         -> Vec3      -- ^ result+--------------------------------------------------------------------------------++rotate3' {- ' CPP is sensitive to primes -} +  :: Flt       -- ^ angle (in radians)+  -> Normal3   -- ^ axis (should be a /unit/ vector!) +  -> Vec3      -- ^ vector+  -> Vec3      -- ^ result rotate3' angle axis v = v .* (rotMatrix3' axis angle) -rotate3 :: Flt    -- ^ angle (in radians)-        -> Vec3   -- ^ axis (arbitrary nonzero vector)-        -> Vec3   -- ^ vector-        -> Vec3   -- ^ result+rotate3 +  :: Flt    -- ^ angle (in radians)+  -> Vec3   -- ^ axis (arbitrary nonzero vector)+  -> Vec3   -- ^ vector+  -> Vec3   -- ^ result rotate3 angle axis v = v .* (rotMatrix3 axis angle)       --- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3 :: Vec3 -> Flt -> Mat3 rotMatrix3 v a = rotMatrix3' (mkNormal v) a --- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3+rotMatrixOrtho3 v a = toOrthoUnsafe $ rotMatrix3 v a++-- | Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3-rotMatrix3' (Normal3 v) a = -  let c = cos a+rotMatrix3' u a = +  let v = fromNormal u+      c = cos a       s = sin a       m1 = scalarMul (1-c) (outer v v)       x = _1 v@@ -73,4 +89,65 @@                 (Vec3 ( s*y) (-s*x)   c   )   in (m1 &+ m2) +rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3+rotMatrixOrtho3' u a = toOrthoUnsafe $ rotMatrix3' u a +--------------------------------------------------------------------------------++-- | Reflects a vector to an axis: that is, the result of @reflect n v@ is+-- 2\<n,v\>n - v+reflect :: Normal3 -> Vec3 -> Vec3+reflect u v = (s *& n) &- v where +  n = fromNormal u+  s = 2 * (n &. v)++reflect' :: Normal3 -> Normal3 -> Normal3+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)+  +refract :: Flt -> Normal3 -> Vec3 -> Vec3+refract eta u v = s *& fromNormal w where+  s = norm v +  w = refract' eta u (toNormalUnsafe $ v &* (1.0/s))+  +-- | Refraction.+-- First parameter (@eta@) is the relative refraction index +--+-- >        refl_inside+-- > eta = --------------+-- >        refl_outside+--+-- where \"inside\" is the direction of the second argument +-- (to vector normal to plane which models the boundary +-- between the two materials). That is, total internal reflection+-- can occur when @eta>1@.+--+-- The convention is that the origin is the point of intersection+-- of the ray and the surface, and all the vectors \"point away\"+-- from here (unlike, say, GLSL's @refract@, where the incident+-- vector \"points towards\" the material)+refract' {- ' CPP is sensitive to primes -} +  :: Flt -> Normal3 -> Normal3 -> Normal3+refract' eta u i = +  if k<0+    then reflect' u i +    else toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++-- | When total internal reflection would occur, we return "Nothing".+refractOnly' {- ' CPP is sensitive to primes -}  +  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3+refractOnly' eta u i = +  if k<0+    then Nothing +    else Just $ toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++--------------------------------------------------------------------------------
Data/Vect/Double/Util/Dim4.hs view
@@ -40,7 +40,7 @@ vec4Z = Vec4 0 0 1 0 vec4W = Vec4 0 0 0 1 ----------------------------------------------------------------------------+--------------------------------------------------------------------------------  -- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code) -- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.@@ -76,8 +76,9 @@ -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors. -- Intended for multiplication on the /right/! rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +rotMatrix4' angle (u1,u2) = m1 &+ (s *& m2) &+ m3    where+    v = fromNormal u1 ; w = fromNormal u2     c = cos angle ; s = sin angle     m1 = scalarMul (1-c) ( outer v v  &+  outer w w )     m2 = biVector4AsTensor v w
Data/Vect/Double/Util/Projective.hs view
@@ -1,6 +1,7 @@ {-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} --- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- | Classic 4x4 projective matrices, encoding the affine transformations of R^3.+-- Our convention is that they are intended for multiplication on -- the /right/, that is, they are of the form -- -- >     _____@@ -13,6 +14,8 @@ -- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices  -- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory, -- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.+--+-- Warning: The naming conventions will probably change in the future.  module Data.Vect.Flt.Util.Projective where @@ -21,68 +24,73 @@  import qualified Data.Vect.Flt.Util.Dim4 as Dim4 -class ExtendProjective v e | v->e where-  extendProj     :: v -> e-  extendProjWith :: Flt -> v -> e-  extendProj = extendProjWith 1-  -instance ExtendProjective Vec2 Vec4 where-  extendProj       (Vec2 x y) = Vec4 x y 0 1-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w-  -instance ExtendProjective Vec3 Vec4 where-  extendProj       (Vec3 x y z) = Vec4 x y z 1-  extendProjWith w (Vec3 x y z) = Vec4 x y z w--instance ExtendProjective Vec4 Vec4 where-  extendProj = id-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w+-------------------------------------------------------------------------------- -instance ExtendProjective Mat2 Mat4 where-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle -instance ExtendProjective Mat3 Mat4 where-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)+rotMatrixProj4 :: Flt -> Vec3 -> Proj4+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle -rotMatrixProj :: Flt -> Normal3 -> Mat4-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle+-- | synonym for "rotateAfterProj4"+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateProj4 = rotateAfterProj4 -rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis)  -translMatrixProj :: Vec3 -> Mat4-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m  --- | we assume that the bottom-right corner is 1.-translWithProj :: Vec3 -> Mat4 -> Mat4-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)+--------------- -scaleMatrixProj :: Vec3 -> Mat4-scaleMatrixProj v = diag $ extendProj v+--scalingUniformProj3 :: Flt -> Proj3+--scalingUniformProj3 x = scaling (Vec2 x x) -scaleMatrixUniformProj :: Flt -> Mat4-scaleMatrixUniformProj s = diag (Vec4 s s s 1)+scalingUniformProj4 :: Flt -> Proj4+scalingUniformProj4 x = scaling (Vec3 x x x) -class ProjectiveAction v where-  actProj :: v -> Mat4 -> v- -instance ProjectiveAction Vec3 where-  actProj v m = trim $ (extendProj v) .* m +-- | Equivalent to @\m -> scaling v .*. m@.+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4+scaleBeforeProj4 (Vec3 u v w) p4 =  +  toProjectiveUnsafe $ +    Mat4 (u*&a) (v*&b) (w*&c) t+  where+    Mat4 a b c t = fromProjective p4 -instance ProjectiveAction Vec4 where-  actProj v m = v .* m +-- | Equivalent to @\m -> m .*. scaling v@.+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4+scaleAfterProj4 v p4 =+  toProjectiveUnsafe $ +    Mat4 (a&!w) (b&!w) (c&!w) (t&!w)+  where+    w = extendWith 1 v+    Mat4 a b c t = fromProjective p4+    +--------------- --- | When acting on unit vectors, we ignore the translation part.-instance ProjectiveAction Normal3 where-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))+-- | Synonym for "translateAfter4"+translate4 :: Vec3 -> Proj4 -> Proj4+translate4 = translateAfter4 --- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,--- and the bottom-right corner is 1.-invertProj :: Mat4 -> Mat4-invertProj mat@(Mat4 u v w t) = -  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+-- | Equivalent to @\m -> m .*. translation v@.+translateAfter4 :: Vec3 -> Proj4 -> Proj4+translateAfter4 v p4 = +  toProjectiveUnsafe $+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)   where-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)-    +    Mat4 r1 r2 r3 r4 = fromProjective p4 ++-- | Equivalent to @\m -> translation v .*. m@.+translateBefore4 :: Vec3 -> Proj4 -> Proj4+translateBefore4 v p4 = +  toProjectiveUnsafe $ +    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) +  where +   u = v .* (trim mat :: Mat3) +   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4+   +---------------+
Data/Vect/Float/Base.hs view
@@ -1,11 +1,34 @@ {-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} -module Data.Vect.Flt.Base where+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-} +module Data.Vect.Flt.Base+  ( AbelianGroup(..) , vecSum+  , MultSemiGroup(..) , Ring , semigroupProduct+  , LeftModule(..) , RightModule(..)+  , Vector(..) , DotProd(..) , CrossProd(..)+  , normalize , distance , angle , angle'+  , UnitVector(..)+  , Pointwise(..)+  , Extend(..) , HasCoordinates(..) , Dimension(..)+  , Matrix(..) , Tensor(..) , Diagonal (..) , Determinant(..)+  , Orthogonal(..) , Projective(..) , MatrixNorms(..)+  , Vec2(..) , Vec3(..) , Vec4(..)+  , Mat2(..) , Mat3(..) , Mat4(..)+  , Ortho2 , Ortho3 , Ortho4+  , Normal2 , Normal3 , Normal4+  , Proj3 , Proj4+  , mkVec2 , mkVec3 , mkVec4+  , project , project' , projectUnsafe , flipNormal+  , householder, householderOrtho+  )+  where+ import Control.Monad import System.Random   import Foreign +-------------------------------------------------------------------------------- -- class declarations  class AbelianGroup g where@@ -20,15 +43,17 @@ vecSum :: AbelianGroup g => [g] -> g vecSum l = foldl (&+) zero l  -class (AbelianGroup r) => -      Ring r where+class MultSemiGroup r where   (.*.) :: r -> r -> r   one   :: r +class (AbelianGroup r, MultSemiGroup r) => Ring r + infixl 7 .*.  -ringProduct :: Ring r => [r] -> r-ringProduct l = foldl (.*.) one l+-- was: ringProduct :: Ring r => [r] -> r+semigroupProduct :: MultSemiGroup r => [r] -> r +semigroupProduct l = foldl (.*.) one l  class LeftModule r m where   lmul :: r -> m -> m@@ -75,7 +100,7 @@   lensqr = normsqr   dotprod :: v -> v -> Flt   normsqr v = (v &. v)  -  norm = sqrt.lensqr+  norm = sqrt . lensqr   dotprod = (&.)  infix 7 &.@@ -112,26 +137,28 @@   fromNormalRadius :: Flt -> u -> v   fromNormalRadius t n = t *& fromNormal n  --- | projects the first vector onto the direction of the second (unit) vector+-- | Projects the first vector onto the direction of the second (unit) vector project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v project' what dir = projectUnsafe what (fromNormal dir) --- | direction (second argument) is assumed to be a /unit/ vector!+-- | Direction (second argument) is assumed to be a /unit/ vector! projectUnsafe :: (Vector v, DotProd v) => v -> v -> v projectUnsafe what dir = what &- dir &* (what &. dir)  project :: (Vector v, DotProd v) => v -> v -> v project what dir = what &- dir &* ((what &. dir) / (dir &. dir)) --- | since unit vectors are not a group, we need a separate function.+-- | Since unit vectors are not a group, we need a separate function. flipNormal :: UnitVector v n => n -> n  flipNormal = toNormalUnsafe . neg . fromNormal  +-- | Cross product class CrossProd v where   crossprod :: v -> v -> v   (&^)      :: v -> v -> v   (&^) = crossprod-  + +-- | Pointwise multiplication  class Pointwise v where   pointwise :: v -> v -> v   (&!)      :: v -> v -> v@@ -166,13 +193,50 @@ "inverse is an involution"    forall m. inverse (inverse m) = m   #-}   +class Matrix m => Orthogonal m o | m->o, o->m where  +  fromOrtho     :: o -> m +  toOrthoUnsafe :: m -> o+  +class (AbelianGroup m, Matrix m) => MatrixNorms m where+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices+  operatorNorm   :: m -> Flt       -- ^ (euclidean) operator norm (not implemented yet)+  matrixDistance m n = frobeniusNorm (n &- m)+  operatorNorm = error "operatorNorm: not implemented yet"+   -- | Outer product (could be unified with Diagonal?) class Tensor t v | t->v where   outer :: v -> v -> t      class Determinant m where   det :: m -> Flt    -    ++class Dimension a where+  dim :: a -> Int+     +-- | Householder matrix, see <http://en.wikipedia.org/wiki/Householder_transformation>.  +-- In plain words, it is the reflection to the hyperplane orthogonal to the input vector.+householder :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m+householder u = idmtx &- (2 *& outer v v) +  where v = fromNormal u++householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o+householderOrtho = toOrthoUnsafe . householder++-- | \"Projective\" matrices have the following form: the top left corner+-- is an any matrix, the bottom right corner is 1, and the top-right+-- column is zero. These describe the affine orthogonal transformation of+-- the space one dimension less.+class (Vector v, Orthogonal n o, Diagonal v n) => Projective v n o m p +    | m->p, p->m, p->o, o->p, p->n, n->p, p->v, v->p, n->o, n->v, v->n where+  fromProjective     :: p -> m+  toProjectiveUnsafe :: m -> p+  orthogonal         :: o -> p+  linear             :: n -> p+  translation        :: v -> p+  scaling            :: v -> p++-------------------------------------------------------------------------------- -- Vec / Mat datatypes   data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt @@ -182,16 +246,16 @@ data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt    deriving (Read,Show) --- | these are /row/ vectors +-- | The components are /row/ vectors  data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show) data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show) data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)  -- | The assumption when dealing with these is always that they are of unit length. -- Also, interpolation works differently.-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) -newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) -newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension)   mkVec2 :: (Flt,Flt) -> Vec2 mkVec3 :: (Flt,Flt,Flt) -> Vec3@@ -201,6 +265,19 @@ mkVec3 (x,y,z)   = Vec3 x y z mkVec4 (x,y,z,w) = Vec4 x y z w +-- | Orthogonal matrices.+--+-- Note: the "Random" instances generates orthogonal matrices with determinant 1+-- (that is, orientation-preserving orthogonal transformations)!+newtype Ortho2 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)++-- | Projective matrices, encoding affine transformations in dimension one less.+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)++-------------------------------------------------------------------------------- -- Unit vectors    instance UnitVector Vec2 Normal2 where@@ -218,47 +295,163 @@   fromNormal (Normal4 v) = v    toNormalUnsafe = Normal4 -rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)-rndUnit g = +_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+_rndUnit g =    if d > 0.01     then ( v &* (1.0/d) , h )-    else rndUnit h+    else _rndUnit h   where     (v,h) = random g     d = norm v      instance Random Normal2 where-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal2 v, h)     randomR _ = random  instance Random Normal3 where-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal3 v, h)     randomR _ = random  instance Random Normal4 where-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal4 v, h)     randomR _ = random -{--instance Storable Normal2 where-  alignment _ = alignment (undefined::Vec2)-  sizeOf    _ = sizeOf    (undefined::Vec2)-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)-  poke p (Normal2 v) = poke (castPtr p) v  - -instance Storable Normal3 where-  alignment _ = alignment (undefined::Vec3)-  sizeOf    _ = sizeOf    (undefined::Vec3)-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)-  poke p (Normal3 v) = poke (castPtr p) v  +instance CrossProd Normal3 where+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w) -instance Storable Normal4 where-  alignment _ = alignment (undefined::Vec4)-  sizeOf    _ = sizeOf    (undefined::Vec4)-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)-  poke p (Normal4 v) = poke (castPtr p) v  --}+--------------------------------------------------------------------------------+-- Orthogonal matrices +instance Orthogonal Mat2 Ortho2 where+  fromOrtho (Ortho2 o) = o+  toOrthoUnsafe = Ortho2++instance Orthogonal Mat3 Ortho3 where+  fromOrtho (Ortho3 o) = o+  toOrthoUnsafe = Ortho3 ++instance Orthogonal Mat4 Ortho4 where+  fromOrtho (Ortho4 o) = o+  toOrthoUnsafe = Ortho4++------++instance Matrix Ortho2 where+  transpose (Ortho2 o) = Ortho2 (transpose o)+  idmtx = Ortho2 idmtx+  inverse = transpose++instance Matrix Ortho3 where+  transpose (Ortho3 o) = Ortho3 (transpose o)+  idmtx = Ortho3 idmtx+  inverse = transpose++instance Matrix Ortho4 where+  transpose (Ortho4 o) = Ortho4 (transpose o)+  idmtx = Ortho4 idmtx+  inverse = transpose++------++instance Random Ortho2 where+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)+  randomR _ = random++instance Random Ortho3 where+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)+  randomR _ = random++instance Random Ortho4 where+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)+  randomR _ = random++------++-- determinant will be -1+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)+_rndOrtho2 g = (h2, g1) where+  h2 = householder u2 :: Mat2 +  (u2,g1) = random g   ++-- generates a uniformly random orthogonal 3x3 matrix +-- /with determinant +1/, with respect to the Haar measure of SO3.+--+-- see Theorem 4 in:+-- Francesco Mezzadri: How to Generate Random Matrices from the Classical Compact Groups +-- Notices of the AMS, May 2007 issue+-- <http://www.ams.org/notices/200705/fea-mezzadri-web.ps>+_rndOrtho3 :: RandomGen g => g -> (Mat3, g) +_rndOrtho3 g = ( (h3 .*. m3), g2) where+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 +  h3 = householder u3 :: Mat3+  (u3,g1) = random g+  (o2,g2) = _rndOrtho2 g1++-- determinant will be -1+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) +_rndOrtho4 g = ( (h4 .*. m4), g2) where+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 +  h4 = householder u4 :: Mat4+  (u4,g1) = random g+  (o3,g2) = _rndOrtho3 g1++------++_flip1stRow2 :: Mat2 -> Mat2+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b++_flip1stRow3 :: Mat3 -> Mat3+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c++_flip1stRow4 :: Mat4 -> Mat4+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d++--------------------------------------------------------------------------------+-- projective matrices+  +instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where+  fromProjective (Proj3 m) = m+  toProjectiveUnsafe = Proj3+  orthogonal = Proj3 . extendWith 1 . fromOrtho+  linear     = Proj3 . extendWith 1+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)+  scaling     v = Proj3 $ diag (extendWith 1 v)+  +instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where+  fromProjective (Proj4 m) = m+  toProjectiveUnsafe = Proj4+  orthogonal = Proj4 . extendWith 1 . fromOrtho +  linear     = Proj4 . extendWith 1+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)+  scaling     v = Proj4 $ diag (extendWith 1 v)++instance Matrix Proj3 where+  idmtx = Proj3 idmtx+  transpose (Proj3 m) = Proj3 (transpose m)+  inverse = _invertProj3++instance Matrix Proj4 where+  idmtx = Proj4 idmtx+  transpose (Proj4 m) = Proj4 (transpose m)+  inverse = _invertProj4++_invertProj3 :: Proj3 -> Proj3+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = +  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec2) .* invm2 +    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)++-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.+-- We assume that the bottom-right corner is 1.+_invertProj4 :: Proj4 -> Proj4+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = +  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec3) .* invm3 +    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)++-------------------------------------------------------------------------------- -- Vec2 instances  instance HasCoordinates Vec2 Flt where@@ -314,7 +507,10 @@         k = sizeOf (undefined::Flt)     poke        p   x     pokeByteOff p k y-               ++instance Dimension Vec2 where dim _ = 2++--------------------------------------------------------------------------------                     -- Mat2 instances  instance HasCoordinates Mat2 Vec2 where@@ -336,19 +532,21 @@   (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)   (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)   neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  -  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)-+  zero = Mat2 zero zero  +   instance Vector Mat2 where   scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s   mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f -instance Ring Mat2 where+instance MultSemiGroup Mat2 where   (.*.) (Mat2 r1 r2) n =      let (Mat2 c1 c2) = transpose n     in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))             (Vec2 (r2 &. c1) (r2 &. c2))   one = idmtx  +instance Ring Mat2+ instance LeftModule Mat2 Vec2 where   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v)    @@ -362,11 +560,6 @@   outer (Vec2 a b) (Vec2 x y) = Mat2     (Vec2 (a*x) (a*y))     (Vec2 (b*x) (b*y))-{--  outer v w = -    let full = Mat2 (Vec2 1 1) (Vec2 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat2 where   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c @@ -393,6 +586,24 @@     poke        p   r1     pokeByteOff p k r2 +instance Random Mat2 where+  random = randomR (Mat2 v1 v1 , Mat2 v2 v2) where +    v1 = Vec2 (-1) (-1) +    v2 = Vec2   1    1+  randomR (Mat2 a b, Mat2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Mat2 x y, gen2)+          +instance Dimension Mat2 where dim _ = 2+     +instance MatrixNorms Mat2 where +  frobeniusNorm (Mat2 r1 r2) =  +    sqrt $+      normsqr r1 + +      normsqr r2+     +--------------------------------------------------------------------------------      -- Vec3 instances  instance HasCoordinates Vec3 Flt where@@ -454,7 +665,10 @@     poke        p       x     pokeByteOff p (k  ) y     pokeByteOff p (k+k) z-   ++instance Dimension Vec3 where dim _ = 3++--------------------------------------------------------------------------------    -- Mat3 instances  instance HasCoordinates Mat3 Vec3 where@@ -495,13 +709,13 @@   (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)   (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)   neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) -  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)+  zero = Mat3 zero zero zero   instance Vector Mat3 where   scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s   mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f -instance Ring Mat3 where+instance MultSemiGroup Mat3 where   (.*.) (Mat3 r1 r2 r3) n =      let (Mat3 c1 c2 c3) = transpose n     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))@@ -509,6 +723,8 @@             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))   one = idmtx  +instance Ring Mat3+ instance LeftModule Mat3 Vec3 where   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)   @@ -523,11 +739,6 @@     (Vec3 (a*x) (a*y) (a*z))     (Vec3 (b*x) (b*y) (b*z))     (Vec3 (c*x) (c*y) (c*z))-{--  outer v w = -    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat3 where   det (Mat3 r1 r2 r3) = det (r1,r2,r3)@@ -556,6 +767,26 @@     pokeByteOff p (k  ) r2     pokeByteOff p (k+k) r3 +instance Random Mat3 where+  random = randomR (Mat3 v1 v1 v1 , Mat3 v2 v2 v2) where+    v1 = Vec3 (-1) (-1) (-1)+    v2 = Vec3   1    1    1+  randomR (Mat3 a b c, Mat3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Mat3 x y z, gen3)+   +instance Dimension Mat3 where dim _ = 3+  +instance MatrixNorms Mat3 where +  frobeniusNorm (Mat3 r1 r2 r3)  = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 +    +-------------------------------------------------------------------------------- -- Vec4 instances  instance HasCoordinates Vec4 Flt where@@ -615,6 +846,9 @@     pokeByteOff p (k+k) z     pokeByteOff p (3*k) w +instance Dimension Vec4 where dim _ = 4++-------------------------------------------------------------------------------- -- Mat4 instances  instance HasCoordinates Mat4 Vec4 where@@ -642,7 +876,7 @@   scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s   mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f -instance Ring Mat4 where+instance MultSemiGroup Mat4 where   (.*.) (Mat4 r1 r2 r3 r4) n =      let (Mat4 c1 c2 c3 c4) = transpose n     in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))@@ -651,6 +885,8 @@             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))   one = idmtx  +instance Ring Mat4+ instance LeftModule Mat4 Vec4 where   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)   @@ -666,14 +902,10 @@     (Vec4 (b*x) (b*y) (b*z) (b*w))     (Vec4 (c*x) (c*y) (c*z) (c*w))     (Vec4 (d*x) (d*y) (d*z) (d*w))-{--  outer v w = -    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)-    in  (diag v) .*. full .*. (diag w)--} ---instance Determinant Mat4 where---  det (Mat4 r1 r2 r3 r4)+instance Determinant Mat4 where+  det = error "det/Mat4: not implemented yet" +  -- det (Mat4 r1 r2 r3 r4) =   {- instance Show Mat4 where@@ -701,6 +933,28 @@     pokeByteOff p (k+k) r3     pokeByteOff p (3*k) r4 +instance Random Mat4 where+  random = randomR (Mat4 v1 v1 v1 v1, Mat4 v2 v2 v2 v2) where+    v1 = Vec4 (-1) (-1) (-1) (-1)+    v2 = Vec4   1    1    1    1+  randomR (Mat4 a b c d, Mat4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Mat4 x y z w, gen4)+    +instance Dimension Mat4 where dim _ = 4+   +instance MatrixNorms Mat4 where +  frobeniusNorm (Mat4 r1 r2 r3 r4) = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 + +      normsqr r4  +    +-------------------------------------------------------------------------------- -- Extend instances  instance Extend Vec2 Vec3 where@@ -719,17 +973,20 @@   trim (Vec4 x y z _)       = Vec3 x y z  instance Extend Mat2 Mat3 where-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith w (Mat2 p q) = Mat3 (extendZero p) (extendZero q) (Vec3 0 0 w)   trim (Mat3 p q _) = Mat2 (trim p) (trim q)  instance Extend Mat2 Mat4 where-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith w (Mat2 p q) = Mat4 (extendZero p) (extendZero q) (Vec4 0 0 w 0) (Vec4 0 0 0 w)   trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)  instance Extend Mat3 Mat4 where-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith w (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) (Vec4 0 0 0 w)   trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+  +--------------------------------------------------------------------------------+
Data/Vect/Float/GramSchmidt.hs view
@@ -10,7 +10,7 @@  import Data.Vect.Flt.Base --------------------------------------------------------+--------------------------------------------------------------------------------  liftPair :: (a -> b) -> (a,a) -> (b,b) liftPair f (x,y) = (f x, f y)@@ -21,7 +21,7 @@ liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b) liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w) --------------------------------------------------------+--------------------------------------------------------------------------------      -- | produces orthogonal\/orthonormal vectors from a set of vectors     class GramSchmidt a where@@ -33,7 +33,7 @@ "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a    #-} --------------------------------------------------------+--------------------------------------------------------------------------------  instance GramSchmidt (Vec2,Vec2) where   gramSchmidt = gramSchmidtPair
Data/Vect/Float/OpenGL.hs view
@@ -3,21 +3,92 @@ -- TODO: the pointer versions of these functions should be really implemented  -- via the pointer versions of the original opengl functions... --- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.   module Data.Vect.Flt.OpenGL where  import Control.Monad import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Projective import qualified Graphics.Rendering.OpenGL as GL  import Foreign -import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)+import Graphics.Rendering.OpenGL hiding +  ( Normal3 , rotate , translate , scale+  , matrix , currentMatrix , withMatrix , multMatrix +  ) --------------------------------------------------------+-------------------------------------------------------------------------------- -{-# SPECIALISE radianToDegrees :: Float -> Float #-}+-- | There should be a big warning here about the different conventions, +-- hidden transpositions, and all the confusion this will inevitably cause...+--+-- As it stands, +--+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3+-- +-- has the same result as+--+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)+--+-- because at the interface of OpenGL and this library there is a transposition+-- to compensate for the different conventions. (This transposition is implicit+-- in the code, because the way the matrices are stored in the memory is also+-- different: OpenGL stores them column-major, and we store them row-major).++class ToOpenGLMatrix m where+  makeGLMatrix :: m -> IO (GLmatrix Flt)++class FromOpenGLMatrix m where+  peekGLMatrix :: GLmatrix Flt -> IO m+  +setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x+ +getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix++matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)++currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m+currentMatrix = matrix Nothing++multMatrix :: ToOpenGLMatrix m => m -> IO ()+multMatrix m = makeGLMatrix m >>= GL.multMatrix++instance ToOpenGLMatrix Mat4 where+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) + +instance FromOpenGLMatrix Mat4 where+  -- huh? GL.withMatrix is strange+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)+  +instance ToOpenGLMatrix Mat3 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)+ +instance ToOpenGLMatrix Mat2 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)++instance ToOpenGLMatrix Ortho4 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)++instance ToOpenGLMatrix Ortho3 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)++instance ToOpenGLMatrix Ortho2 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)++instance ToOpenGLMatrix Proj4 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)++instance ToOpenGLMatrix Proj3 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)+  +--------------------------------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float  -> Float  #-} {-# SPECIALISE radianToDegrees :: Double -> Double #-} radianToDegrees :: RealFrac a => a -> a radianToDegrees x = x * 57.295779513082322@@ -28,20 +99,61 @@ degreesToRadian x = x * 1.7453292519943295e-2  -- | The angle is in radians. (WARNING: OpenGL uses degrees!)-rotate :: Flt -> Vec3 -> IO ()-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)+glRotate :: Flt -> Vec3 -> IO ()+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z) -translate :: Vec3 -> IO ()-translate (Vec3 x y z) = GL.translate (Vector3 x y z)+glTranslate :: Vec3 -> IO ()+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z) -scale3 :: Vec3 -> IO ()-scale3 (Vec3 x y z) = GL.scale x y z+glScale3 :: Vec3 -> IO ()+glScale3 (Vec3 x y z) = GL.scale x y z -scale :: Flt -> IO ()-scale x = GL.scale x x x+glScale :: Flt -> IO ()+glScale x = GL.scale x x x --------------------------------------------------------+--------------------------------------------------------------------------------+ +-- | \"Orthogonal projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1)+orthoMatrix +  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +orthoMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2/(r-l)) 0 0 0)+  (Vec4 0 (2/(t-b)) 0 0)+  (Vec4 0 0 (-2/(f-n)) 0)+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)+  +-- | The same as "orthoMatrix", but with a different parametrization.+orthoMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f) +-- | \"Perspective projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1).+frustumMatrix+  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +frustumMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2*n/(r-l)) 0 0 0)+  (Vec4 0 (2*n/(t-b)) 0 0)+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))+  (Vec4 0 0 (-2*f*n*(f-n)) 0)+  +-- | The same as "frustumMatrix", but with a different parametrization.+frustumMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)++-------------------------------------------------------------------------------- -- Vertex instances  instance GL.Vertex Vec2 where@@ -56,17 +168,20 @@   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)   vertexv p = peek p >>= vertex    ---------------------------------------------------------+-------------------------------------------------------------------------------- -- the Normal instance -- note that there is no Normal2\/Normal4 in the OpenGL binding  instance GL.Normal Normal3 where-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normal u = GL.normal (GL.Normal3 x y z) +    where Vec3 x y z = fromNormal u    normalv p = peek p >>= normal  --------------------------------------------------------+instance GL.Normal Vec3 where+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) +  normalv p = peek p >>= normal  +-------------------------------------------------------------------------------- -- Color instances    instance GL.Color Vec3 where@@ -88,8 +203,7 @@   secondaryColorv p = peek p >>= secondaryColor -} ---------------------------------------------------------+-------------------------------------------------------------------------------- -- TexCoord instances  instance GL.TexCoord Vec2 where@@ -110,8 +224,7 @@   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)   multiTexCoordv unit p = peek p >>= multiTexCoord unit ---------------------------------------------------------    +-------------------------------------------------------------------------------- -- Vertex Attributes (experimental)  class VertexAttrib' a where@@ -130,16 +243,18 @@   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w   instance VertexAttrib' Normal2 where-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y+  vertexAttrib loc u = GL.vertexAttrib2 loc x y+    where Vec2 x y = fromNormal u   instance VertexAttrib' Normal3 where-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z+    where Vec3 x y z = fromNormal u   instance VertexAttrib' Normal4 where-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w----------------------------------------------------------+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w+    where Vec4 x y z w = fromNormal u +   +-------------------------------------------------------------------------------- -- Uniform (again, experimental)  -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; @@ -181,3 +296,5 @@   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)      #endif++
Data/Vect/Float/Util/Dim2.hs view
@@ -4,12 +4,13 @@  import Data.Vect.Flt.Base --- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@. structVec2 :: [Flt] -> [Vec2] structVec2 [] = [] structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls)  structVec2 _ = error "structVec2" +-- | The opposite of "structVec2". destructVec2 :: [Vec2] -> [Flt] destructVec2 [] = [] destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  @@ -49,10 +50,13 @@ angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt angle2' = angle2 . fromNormal --- |Rotation matrix by a given angle (in radians), counterclockwise.+-- | Rotation matrix by a given angle (in radians), counterclockwise. rotMatrix2 :: Flt -> Mat2 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a +rotMatrixOrtho2 :: Flt -> Ortho2+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2+ rotate2 :: Flt -> Vec2 -> Vec2 rotate2 a v = v .* (rotMatrix2 a)  @@ -63,3 +67,5 @@ -- |Rotates clockwise by 90 degrees. rotateCW :: Vec2 -> Vec2 rotateCW (Vec2 x y) = Vec2 y (-x)++
Data/Vect/Float/Util/Dim3.hs view
@@ -4,18 +4,26 @@  import Data.Vect.Flt.Base +--------------------------------------------------------------------------------++-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@. structVec3 :: [Flt] -> [Vec3] structVec3 [] = [] structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls)  structVec3 _ = error "structVec3" +-- | The opposite of "structVec3". destructVec3 :: [Vec3] -> [Flt] destructVec3 [] = [] destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)   +--------------------------------------------------------------------------------+ det3 :: Vec3 -> Vec3 -> Vec3 -> Flt det3 u v w = det (u,v,w) +--------------------------------------------------------------------------------+ translate3X :: Flt -> Vec3 -> Vec3 translate3Y :: Flt -> Vec3 -> Vec3 translate3Z :: Flt -> Vec3 -> Vec3@@ -43,26 +51,34 @@ rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a -rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)-         -> Normal3   -- ^ axis (should be a /unit/ vector!) -         -> Vec3      -- ^ vector-         -> Vec3      -- ^ result+--------------------------------------------------------------------------------++rotate3' {- ' CPP is sensitive to primes -} +  :: Flt       -- ^ angle (in radians)+  -> Normal3   -- ^ axis (should be a /unit/ vector!) +  -> Vec3      -- ^ vector+  -> Vec3      -- ^ result rotate3' angle axis v = v .* (rotMatrix3' axis angle) -rotate3 :: Flt    -- ^ angle (in radians)-        -> Vec3   -- ^ axis (arbitrary nonzero vector)-        -> Vec3   -- ^ vector-        -> Vec3   -- ^ result+rotate3 +  :: Flt    -- ^ angle (in radians)+  -> Vec3   -- ^ axis (arbitrary nonzero vector)+  -> Vec3   -- ^ vector+  -> Vec3   -- ^ result rotate3 angle axis v = v .* (rotMatrix3 axis angle)       --- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3 :: Vec3 -> Flt -> Mat3 rotMatrix3 v a = rotMatrix3' (mkNormal v) a --- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3+rotMatrixOrtho3 v a = toOrthoUnsafe $ rotMatrix3 v a++-- | Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3-rotMatrix3' (Normal3 v) a = -  let c = cos a+rotMatrix3' u a = +  let v = fromNormal u+      c = cos a       s = sin a       m1 = scalarMul (1-c) (outer v v)       x = _1 v@@ -73,4 +89,65 @@                 (Vec3 ( s*y) (-s*x)   c   )   in (m1 &+ m2) +rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3+rotMatrixOrtho3' u a = toOrthoUnsafe $ rotMatrix3' u a +--------------------------------------------------------------------------------++-- | Reflects a vector to an axis: that is, the result of @reflect n v@ is+-- 2\<n,v\>n - v+reflect :: Normal3 -> Vec3 -> Vec3+reflect u v = (s *& n) &- v where +  n = fromNormal u+  s = 2 * (n &. v)++reflect' :: Normal3 -> Normal3 -> Normal3+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)+  +refract :: Flt -> Normal3 -> Vec3 -> Vec3+refract eta u v = s *& fromNormal w where+  s = norm v +  w = refract' eta u (toNormalUnsafe $ v &* (1.0/s))+  +-- | Refraction.+-- First parameter (@eta@) is the relative refraction index +--+-- >        refl_inside+-- > eta = --------------+-- >        refl_outside+--+-- where \"inside\" is the direction of the second argument +-- (to vector normal to plane which models the boundary +-- between the two materials). That is, total internal reflection+-- can occur when @eta>1@.+--+-- The convention is that the origin is the point of intersection+-- of the ray and the surface, and all the vectors \"point away\"+-- from here (unlike, say, GLSL's @refract@, where the incident+-- vector \"points towards\" the material)+refract' {- ' CPP is sensitive to primes -} +  :: Flt -> Normal3 -> Normal3 -> Normal3+refract' eta u i = +  if k<0+    then reflect' u i +    else toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++-- | When total internal reflection would occur, we return "Nothing".+refractOnly' {- ' CPP is sensitive to primes -}  +  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3+refractOnly' eta u i = +  if k<0+    then Nothing +    else Just $ toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++--------------------------------------------------------------------------------
Data/Vect/Float/Util/Dim4.hs view
@@ -40,7 +40,7 @@ vec4Z = Vec4 0 0 1 0 vec4W = Vec4 0 0 0 1 ----------------------------------------------------------------------------+--------------------------------------------------------------------------------  -- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code) -- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.@@ -76,8 +76,9 @@ -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors. -- Intended for multiplication on the /right/! rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +rotMatrix4' angle (u1,u2) = m1 &+ (s *& m2) &+ m3    where+    v = fromNormal u1 ; w = fromNormal u2     c = cos angle ; s = sin angle     m1 = scalarMul (1-c) ( outer v v  &+  outer w w )     m2 = biVector4AsTensor v w
Data/Vect/Float/Util/Projective.hs view
@@ -1,6 +1,7 @@ {-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} --- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- | Classic 4x4 projective matrices, encoding the affine transformations of R^3.+-- Our convention is that they are intended for multiplication on -- the /right/, that is, they are of the form -- -- >     _____@@ -13,6 +14,8 @@ -- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices  -- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory, -- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.+--+-- Warning: The naming conventions will probably change in the future.  module Data.Vect.Flt.Util.Projective where @@ -21,68 +24,73 @@  import qualified Data.Vect.Flt.Util.Dim4 as Dim4 -class ExtendProjective v e | v->e where-  extendProj     :: v -> e-  extendProjWith :: Flt -> v -> e-  extendProj = extendProjWith 1-  -instance ExtendProjective Vec2 Vec4 where-  extendProj       (Vec2 x y) = Vec4 x y 0 1-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w-  -instance ExtendProjective Vec3 Vec4 where-  extendProj       (Vec3 x y z) = Vec4 x y z 1-  extendProjWith w (Vec3 x y z) = Vec4 x y z w--instance ExtendProjective Vec4 Vec4 where-  extendProj = id-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w+-------------------------------------------------------------------------------- -instance ExtendProjective Mat2 Mat4 where-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle -instance ExtendProjective Mat3 Mat4 where-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)+rotMatrixProj4 :: Flt -> Vec3 -> Proj4+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle -rotMatrixProj :: Flt -> Normal3 -> Mat4-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle+-- | synonym for "rotateAfterProj4"+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateProj4 = rotateAfterProj4 -rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis)  -translMatrixProj :: Vec3 -> Mat4-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m  --- | we assume that the bottom-right corner is 1.-translWithProj :: Vec3 -> Mat4 -> Mat4-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)+--------------- -scaleMatrixProj :: Vec3 -> Mat4-scaleMatrixProj v = diag $ extendProj v+--scalingUniformProj3 :: Flt -> Proj3+--scalingUniformProj3 x = scaling (Vec2 x x) -scaleMatrixUniformProj :: Flt -> Mat4-scaleMatrixUniformProj s = diag (Vec4 s s s 1)+scalingUniformProj4 :: Flt -> Proj4+scalingUniformProj4 x = scaling (Vec3 x x x) -class ProjectiveAction v where-  actProj :: v -> Mat4 -> v- -instance ProjectiveAction Vec3 where-  actProj v m = trim $ (extendProj v) .* m +-- | Equivalent to @\m -> scaling v .*. m@.+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4+scaleBeforeProj4 (Vec3 u v w) p4 =  +  toProjectiveUnsafe $ +    Mat4 (u*&a) (v*&b) (w*&c) t+  where+    Mat4 a b c t = fromProjective p4 -instance ProjectiveAction Vec4 where-  actProj v m = v .* m +-- | Equivalent to @\m -> m .*. scaling v@.+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4+scaleAfterProj4 v p4 =+  toProjectiveUnsafe $ +    Mat4 (a&!w) (b&!w) (c&!w) (t&!w)+  where+    w = extendWith 1 v+    Mat4 a b c t = fromProjective p4+    +--------------- --- | When acting on unit vectors, we ignore the translation part.-instance ProjectiveAction Normal3 where-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))+-- | Synonym for "translateAfter4"+translate4 :: Vec3 -> Proj4 -> Proj4+translate4 = translateAfter4 --- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,--- and the bottom-right corner is 1.-invertProj :: Mat4 -> Mat4-invertProj mat@(Mat4 u v w t) = -  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+-- | Equivalent to @\m -> m .*. translation v@.+translateAfter4 :: Vec3 -> Proj4 -> Proj4+translateAfter4 v p4 = +  toProjectiveUnsafe $+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)   where-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)-    +    Mat4 r1 r2 r3 r4 = fromProjective p4 ++-- | Equivalent to @\m -> translation v .*. m@.+translateBefore4 :: Vec3 -> Proj4 -> Proj4+translateBefore4 v p4 = +  toProjectiveUnsafe $ +    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) +  where +   u = v .* (trim mat :: Mat3) +   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4+   +---------------+
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2008, Balazs Komuves+Copyright (c) 2008-2009, Balazs Komuves All rights reserved.  Redistribution and use in source and binary forms, with or without@@ -15,15 +15,15 @@ may be used to endorse or promote products derived from this software without specific prior written permission.  -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER -OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
src/flt/Base.hs view
@@ -1,10 +1,33 @@ -module Data.Vect.Flt.Base where+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, GeneralizedNewtypeDeriving #-} +module Data.Vect.Flt.Base+  ( AbelianGroup(..) , vecSum+  , MultSemiGroup(..) , Ring , semigroupProduct+  , LeftModule(..) , RightModule(..)+  , Vector(..) , DotProd(..) , CrossProd(..)+  , normalize , distance , angle , angle'+  , UnitVector(..)+  , Pointwise(..)+  , Extend(..) , HasCoordinates(..) , Dimension(..)+  , Matrix(..) , Tensor(..) , Diagonal (..) , Determinant(..)+  , Orthogonal(..) , Projective(..) , MatrixNorms(..)+  , Vec2(..) , Vec3(..) , Vec4(..)+  , Mat2(..) , Mat3(..) , Mat4(..)+  , Ortho2 , Ortho3 , Ortho4+  , Normal2 , Normal3 , Normal4+  , Proj3 , Proj4+  , mkVec2 , mkVec3 , mkVec4+  , project , project' , projectUnsafe , flipNormal+  , householder, householderOrtho+  )+  where+ import Control.Monad import System.Random   import Foreign +-------------------------------------------------------------------------------- -- class declarations  class AbelianGroup g where@@ -19,15 +42,17 @@ vecSum :: AbelianGroup g => [g] -> g vecSum l = foldl (&+) zero l  -class (AbelianGroup r) => -      Ring r where+class MultSemiGroup r where   (.*.) :: r -> r -> r   one   :: r +class (AbelianGroup r, MultSemiGroup r) => Ring r + infixl 7 .*.  -ringProduct :: Ring r => [r] -> r-ringProduct l = foldl (.*.) one l+-- was: ringProduct :: Ring r => [r] -> r+semigroupProduct :: MultSemiGroup r => [r] -> r +semigroupProduct l = foldl (.*.) one l  class LeftModule r m where   lmul :: r -> m -> m@@ -74,7 +99,7 @@   lensqr = normsqr   dotprod :: v -> v -> Flt   normsqr v = (v &. v)  -  norm = sqrt.lensqr+  norm = sqrt . lensqr   dotprod = (&.)  infix 7 &.@@ -111,26 +136,28 @@   fromNormalRadius :: Flt -> u -> v   fromNormalRadius t n = t *& fromNormal n  --- | projects the first vector onto the direction of the second (unit) vector+-- | Projects the first vector onto the direction of the second (unit) vector project' :: (Vector v, UnitVector v u, DotProd v) => v -> u -> v project' what dir = projectUnsafe what (fromNormal dir) --- | direction (second argument) is assumed to be a /unit/ vector!+-- | Direction (second argument) is assumed to be a /unit/ vector! projectUnsafe :: (Vector v, DotProd v) => v -> v -> v projectUnsafe what dir = what &- dir &* (what &. dir)  project :: (Vector v, DotProd v) => v -> v -> v project what dir = what &- dir &* ((what &. dir) / (dir &. dir)) --- | since unit vectors are not a group, we need a separate function.+-- | Since unit vectors are not a group, we need a separate function. flipNormal :: UnitVector v n => n -> n  flipNormal = toNormalUnsafe . neg . fromNormal  +-- | Cross product class CrossProd v where   crossprod :: v -> v -> v   (&^)      :: v -> v -> v   (&^) = crossprod-  + +-- | Pointwise multiplication  class Pointwise v where   pointwise :: v -> v -> v   (&!)      :: v -> v -> v@@ -165,13 +192,50 @@ "inverse is an involution"    forall m. inverse (inverse m) = m   #-}   +class Matrix m => Orthogonal m o | m->o, o->m where  +  fromOrtho     :: o -> m +  toOrthoUnsafe :: m -> o+  +class (AbelianGroup m, Matrix m) => MatrixNorms m where+  frobeniusNorm  :: m -> Flt       -- ^ the frobenius norm (= euclidean norm in the space of matrices)+  matrixDistance :: m -> m -> Flt  -- ^ euclidean distance in the space of matrices+  operatorNorm   :: m -> Flt       -- ^ (euclidean) operator norm (not implemented yet)+  matrixDistance m n = frobeniusNorm (n &- m)+  operatorNorm = error "operatorNorm: not implemented yet"+   -- | Outer product (could be unified with Diagonal?) class Tensor t v | t->v where   outer :: v -> v -> t      class Determinant m where   det :: m -> Flt    -    ++class Dimension a where+  dim :: a -> Int+     +-- | Householder matrix, see <http://en.wikipedia.org/wiki/Householder_transformation>.  +-- In plain words, it is the reflection to the hyperplane orthogonal to the input vector.+householder :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v) => u -> m+householder u = idmtx &- (2 *& outer v v) +  where v = fromNormal u++householderOrtho :: (Vector v, UnitVector v u, Matrix m, Vector m, Tensor m v, Orthogonal m o) => u -> o+householderOrtho = toOrthoUnsafe . householder++-- | \"Projective\" matrices have the following form: the top left corner+-- is an any matrix, the bottom right corner is 1, and the top-right+-- column is zero. These describe the affine orthogonal transformation of+-- the space one dimension less.+class (Vector v, Orthogonal n o, Diagonal v n) => Projective v n o m p +    | m->p, p->m, p->o, o->p, p->n, n->p, p->v, v->p, n->o, n->v, v->n where+  fromProjective     :: p -> m+  toProjectiveUnsafe :: m -> p+  orthogonal         :: o -> p+  linear             :: n -> p+  translation        :: v -> p+  scaling            :: v -> p++-------------------------------------------------------------------------------- -- Vec / Mat datatypes   data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt @@ -181,16 +245,16 @@ data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt    deriving (Read,Show) --- | these are /row/ vectors +-- | The components are /row/ vectors  data Mat2 = Mat2 !Vec2 !Vec2              deriving (Read,Show) data Mat3 = Mat3 !Vec3 !Vec3 !Vec3        deriving (Read,Show) data Mat4 = Mat4 !Vec4 !Vec4 !Vec4 !Vec4  deriving (Read,Show)  -- | The assumption when dealing with these is always that they are of unit length. -- Also, interpolation works differently.-newtype Normal2 = Normal2 Vec2 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) -newtype Normal3 = Normal3 Vec3 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable,CrossProd) -newtype Normal4 = Normal4 Vec4 deriving ({-AbelianGroup,Vector,-}Read,Show,DotProd,Storable) +newtype Normal2 = Normal2 Vec2 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal3 = Normal3 Vec3 deriving (Read,Show,Storable,DotProd,Dimension) +newtype Normal4 = Normal4 Vec4 deriving (Read,Show,Storable,DotProd,Dimension)   mkVec2 :: (Flt,Flt) -> Vec2 mkVec3 :: (Flt,Flt,Flt) -> Vec3@@ -200,6 +264,19 @@ mkVec3 (x,y,z)   = Vec3 x y z mkVec4 (x,y,z,w) = Vec4 x y z w +-- | Orthogonal matrices.+--+-- Note: the "Random" instances generates orthogonal matrices with determinant 1+-- (that is, orientation-preserving orthogonal transformations)!+newtype Ortho2 = Ortho2 Mat2 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho3 = Ortho3 Mat3 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)+newtype Ortho4 = Ortho4 Mat4 deriving (Read,Show,Storable,MultSemiGroup,Determinant,Dimension)++-- | Projective matrices, encoding affine transformations in dimension one less.+newtype Proj3 = Proj3 Mat3 deriving (Read,Show,Storable,MultSemiGroup)+newtype Proj4 = Proj4 Mat4 deriving (Read,Show,Storable,MultSemiGroup)++-------------------------------------------------------------------------------- -- Unit vectors    instance UnitVector Vec2 Normal2 where@@ -217,47 +294,163 @@   fromNormal (Normal4 v) = v    toNormalUnsafe = Normal4 -rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)-rndUnit g = +_rndUnit :: (RandomGen g, Random v, Vector v, DotProd v) => g -> (v,g)+_rndUnit g =    if d > 0.01     then ( v &* (1.0/d) , h )-    else rndUnit h+    else _rndUnit h   where     (v,h) = random g     d = norm v      instance Random Normal2 where-  random g = let (v,h) = rndUnit g in (Normal2 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal2 v, h)     randomR _ = random  instance Random Normal3 where-  random g = let (v,h) = rndUnit g in (Normal3 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal3 v, h)     randomR _ = random  instance Random Normal4 where-  random g = let (v,h) = rndUnit g in (Normal4 v, h)  +  random g = let (v,h) = _rndUnit g in (Normal4 v, h)     randomR _ = random -{--instance Storable Normal2 where-  alignment _ = alignment (undefined::Vec2)-  sizeOf    _ = sizeOf    (undefined::Vec2)-  peek p = liftM (\v -> Normal2 v) (peek $ castPtr p)-  poke p (Normal2 v) = poke (castPtr p) v  - -instance Storable Normal3 where-  alignment _ = alignment (undefined::Vec3)-  sizeOf    _ = sizeOf    (undefined::Vec3)-  peek p = liftM (\v -> Normal3 v) (peek $ castPtr p)-  poke p (Normal3 v) = poke (castPtr p) v  +instance CrossProd Normal3 where+  crossprod (Normal3 v) (Normal3 w) = mkNormal (crossprod v w) -instance Storable Normal4 where-  alignment _ = alignment (undefined::Vec4)-  sizeOf    _ = sizeOf    (undefined::Vec4)-  peek p = liftM (\v -> Normal4 v) (peek $ castPtr p)-  poke p (Normal4 v) = poke (castPtr p) v  --}+--------------------------------------------------------------------------------+-- Orthogonal matrices +instance Orthogonal Mat2 Ortho2 where+  fromOrtho (Ortho2 o) = o+  toOrthoUnsafe = Ortho2++instance Orthogonal Mat3 Ortho3 where+  fromOrtho (Ortho3 o) = o+  toOrthoUnsafe = Ortho3 ++instance Orthogonal Mat4 Ortho4 where+  fromOrtho (Ortho4 o) = o+  toOrthoUnsafe = Ortho4++------++instance Matrix Ortho2 where+  transpose (Ortho2 o) = Ortho2 (transpose o)+  idmtx = Ortho2 idmtx+  inverse = transpose++instance Matrix Ortho3 where+  transpose (Ortho3 o) = Ortho3 (transpose o)+  idmtx = Ortho3 idmtx+  inverse = transpose++instance Matrix Ortho4 where+  transpose (Ortho4 o) = Ortho4 (transpose o)+  idmtx = Ortho4 idmtx+  inverse = transpose++------++instance Random Ortho2 where+  random g = let (o,h) = _rndOrtho2 g in (toOrthoUnsafe (_flip1stRow2 o), h)+  randomR _ = random++instance Random Ortho3 where+  random g = let (o,h) = _rndOrtho3 g in (toOrthoUnsafe (             o), h)+  randomR _ = random++instance Random Ortho4 where+  random g = let (o,h) = _rndOrtho4 g in (toOrthoUnsafe (_flip1stRow4 o), h)+  randomR _ = random++------++-- determinant will be -1+_rndOrtho2 :: RandomGen g => g -> (Mat2, g)+_rndOrtho2 g = (h2, g1) where+  h2 = householder u2 :: Mat2 +  (u2,g1) = random g   ++-- generates a uniformly random orthogonal 3x3 matrix +-- /with determinant +1/, with respect to the Haar measure of SO3.+--+-- see Theorem 4 in:+-- Francesco Mezzadri: How to Generate Random Matrices from the Classical Compact Groups +-- Notices of the AMS, May 2007 issue+-- <http://www.ams.org/notices/200705/fea-mezzadri-web.ps>+_rndOrtho3 :: RandomGen g => g -> (Mat3, g) +_rndOrtho3 g = ( (h3 .*. m3), g2) where+  m3 = (extendWith :: Flt -> Mat2 -> Mat3) 1 o2 +  h3 = householder u3 :: Mat3+  (u3,g1) = random g+  (o2,g2) = _rndOrtho2 g1++-- determinant will be -1+_rndOrtho4 :: RandomGen g => g -> (Mat4, g) +_rndOrtho4 g = ( (h4 .*. m4), g2) where+  m4 = (extendWith :: Flt -> Mat3 -> Mat4) 1 o3 +  h4 = householder u4 :: Mat4+  (u4,g1) = random g+  (o3,g2) = _rndOrtho3 g1++------++_flip1stRow2 :: Mat2 -> Mat2+_flip1stRow2 (Mat2 a b) = Mat2 (neg a) b++_flip1stRow3 :: Mat3 -> Mat3+_flip1stRow3 (Mat3 a b c) = Mat3 (neg a) b c++_flip1stRow4 :: Mat4 -> Mat4+_flip1stRow4 (Mat4 a b c d) = Mat4 (neg a) b c d++--------------------------------------------------------------------------------+-- projective matrices+  +instance Projective Vec2 Mat2 Ortho2 Mat3 Proj3 where+  fromProjective (Proj3 m) = m+  toProjectiveUnsafe = Proj3+  orthogonal = Proj3 . extendWith 1 . fromOrtho+  linear     = Proj3 . extendWith 1+  translation v = Proj3 $ Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (extendWith 1 v)+  scaling     v = Proj3 $ diag (extendWith 1 v)+  +instance Projective Vec3 Mat3 Ortho3 Mat4 Proj4 where+  fromProjective (Proj4 m) = m+  toProjectiveUnsafe = Proj4+  orthogonal = Proj4 . extendWith 1 . fromOrtho +  linear     = Proj4 . extendWith 1+  translation v = Proj4 $ Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (extendWith 1 v)+  scaling     v = Proj4 $ diag (extendWith 1 v)++instance Matrix Proj3 where+  idmtx = Proj3 idmtx+  transpose (Proj3 m) = Proj3 (transpose m)+  inverse = _invertProj3++instance Matrix Proj4 where+  idmtx = Proj4 idmtx+  transpose (Proj4 m) = Proj4 (transpose m)+  inverse = _invertProj4++_invertProj3 :: Proj3 -> Proj3+_invertProj3 (Proj3 mat@(Mat3 _ _ t)) = +  Proj3 $ Mat3 (extendZero a) (extendZero b) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec2) .* invm2 +    invm2@(Mat2 a b) = inverse $ (trim mat :: Mat2)++-- Inverts a projective 4x4 matrix. But you can simply use "inverse" instead.+-- We assume that the bottom-right corner is 1.+_invertProj4 :: Proj4 -> Proj4+_invertProj4 (Proj4 mat@(Mat4 _ _ _ t)) = +  Proj4 $ Mat4 (extendZero a) (extendZero b) (extendZero c) (extendWith 1 t') +  where+    t' = neg $ (trim t :: Vec3) .* invm3 +    invm3@(Mat3 a b c) = inverse $ (trim mat :: Mat3)++-------------------------------------------------------------------------------- -- Vec2 instances  instance HasCoordinates Vec2 Flt where@@ -313,7 +506,10 @@         k = sizeOf (undefined::Flt)     poke        p   x     pokeByteOff p k y-               ++instance Dimension Vec2 where dim _ = 2++--------------------------------------------------------------------------------                     -- Mat2 instances  instance HasCoordinates Mat2 Vec2 where@@ -335,19 +531,21 @@   (&+) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &+ s1) (r2 &+ s2)   (&-) (Mat2 r1 r2) (Mat2 s1 s2) = Mat2 (r1 &- s1) (r2 &- s2)   neg  (Mat2 r1 r2)              = Mat2 (neg r1) (neg r2)  -  zero = Mat2 zero zero   -- (zero::Vec2) (zero::Vec2)-+  zero = Mat2 zero zero  +   instance Vector Mat2 where   scalarMul s (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = scalarMul s   mapVec    f (Mat2 r1 r2) = Mat2 (g r1) (g r2) where g = mapVec f -instance Ring Mat2 where+instance MultSemiGroup Mat2 where   (.*.) (Mat2 r1 r2) n =      let (Mat2 c1 c2) = transpose n     in Mat2 (Vec2 (r1 &. c1) (r1 &. c2))             (Vec2 (r2 &. c1) (r2 &. c2))   one = idmtx  +instance Ring Mat2+ instance LeftModule Mat2 Vec2 where   lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v)    @@ -361,11 +559,6 @@   outer (Vec2 a b) (Vec2 x y) = Mat2     (Vec2 (a*x) (a*y))     (Vec2 (b*x) (b*y))-{--  outer v w = -    let full = Mat2 (Vec2 1 1) (Vec2 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat2 where   det (Mat2 (Vec2 a b) (Vec2 c d)) = a*d - b*c @@ -392,6 +585,24 @@     poke        p   r1     pokeByteOff p k r2 +instance Random Mat2 where+  random = randomR (Mat2 v1 v1 , Mat2 v2 v2) where +    v1 = Vec2 (-1) (-1) +    v2 = Vec2   1    1+  randomR (Mat2 a b, Mat2 c d) gen = +    let (x,gen1) = randomR (a,c) gen+        (y,gen2) = randomR (b,d) gen1+    in (Mat2 x y, gen2)+          +instance Dimension Mat2 where dim _ = 2+     +instance MatrixNorms Mat2 where +  frobeniusNorm (Mat2 r1 r2) =  +    sqrt $+      normsqr r1 + +      normsqr r2+     +--------------------------------------------------------------------------------      -- Vec3 instances  instance HasCoordinates Vec3 Flt where@@ -453,7 +664,10 @@     poke        p       x     pokeByteOff p (k  ) y     pokeByteOff p (k+k) z-   ++instance Dimension Vec3 where dim _ = 3++--------------------------------------------------------------------------------    -- Mat3 instances  instance HasCoordinates Mat3 Vec3 where@@ -494,13 +708,13 @@   (&+) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3)   (&-) (Mat3 r1 r2 r3) (Mat3 s1 s2 s3) = Mat3 (r1 &- s1) (r2 &- s2) (r3 &- s3)   neg  (Mat3 r1 r2 r3)                 = Mat3 (neg r1) (neg r2) (neg r3) -  zero = Mat3 zero zero zero   -- (zero::Vec3) (zero::Vec3) (zero::Vec3)+  zero = Mat3 zero zero zero   instance Vector Mat3 where   scalarMul s (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = scalarMul s   mapVec    f (Mat3 r1 r2 r3) = Mat3 (g r1) (g r2) (g r3) where g = mapVec f -instance Ring Mat3 where+instance MultSemiGroup Mat3 where   (.*.) (Mat3 r1 r2 r3) n =      let (Mat3 c1 c2 c3) = transpose n     in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3))@@ -508,6 +722,8 @@             (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3))   one = idmtx  +instance Ring Mat3+ instance LeftModule Mat3 Vec3 where   lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v)   @@ -522,11 +738,6 @@     (Vec3 (a*x) (a*y) (a*z))     (Vec3 (b*x) (b*y) (b*z))     (Vec3 (c*x) (c*y) (c*z))-{--  outer v w = -    let full = Mat3 (Vec3 1 1 1) (Vec3 1 1 1) (Vec3 1 1 1)-    in  (diag v) .*. full .*. (diag w)--}  instance Determinant Mat3 where   det (Mat3 r1 r2 r3) = det (r1,r2,r3)@@ -555,6 +766,26 @@     pokeByteOff p (k  ) r2     pokeByteOff p (k+k) r3 +instance Random Mat3 where+  random = randomR (Mat3 v1 v1 v1 , Mat3 v2 v2 v2) where+    v1 = Vec3 (-1) (-1) (-1)+    v2 = Vec3   1    1    1+  randomR (Mat3 a b c, Mat3 d e f) gen = +    let (x,gen1) = randomR (a,d) gen+        (y,gen2) = randomR (b,e) gen1+        (z,gen3) = randomR (c,f) gen2  +    in (Mat3 x y z, gen3)+   +instance Dimension Mat3 where dim _ = 3+  +instance MatrixNorms Mat3 where +  frobeniusNorm (Mat3 r1 r2 r3)  = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 +    +-------------------------------------------------------------------------------- -- Vec4 instances  instance HasCoordinates Vec4 Flt where@@ -614,6 +845,9 @@     pokeByteOff p (k+k) z     pokeByteOff p (3*k) w +instance Dimension Vec4 where dim _ = 4++-------------------------------------------------------------------------------- -- Mat4 instances  instance HasCoordinates Mat4 Vec4 where@@ -641,7 +875,7 @@   scalarMul s (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = scalarMul s   mapVec    f (Mat4 r1 r2 r3 r4) = Mat4 (g r1) (g r2) (g r3) (g r4) where g = mapVec f -instance Ring Mat4 where+instance MultSemiGroup Mat4 where   (.*.) (Mat4 r1 r2 r3 r4) n =      let (Mat4 c1 c2 c3 c4) = transpose n     in Mat4 (Vec4 (r1 &. c1) (r1 &. c2) (r1 &. c3) (r1 &. c4))@@ -650,6 +884,8 @@             (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4))   one = idmtx  +instance Ring Mat4+ instance LeftModule Mat4 Vec4 where   lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v)   @@ -665,14 +901,10 @@     (Vec4 (b*x) (b*y) (b*z) (b*w))     (Vec4 (c*x) (c*y) (c*z) (c*w))     (Vec4 (d*x) (d*y) (d*z) (d*w))-{--  outer v w = -    let full = Mat4 (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1) (Vec4 1 1 1 1)-    in  (diag v) .*. full .*. (diag w)--} ---instance Determinant Mat4 where---  det (Mat4 r1 r2 r3 r4)+instance Determinant Mat4 where+  det = error "det/Mat4: not implemented yet" +  -- det (Mat4 r1 r2 r3 r4) =   {- instance Show Mat4 where@@ -700,6 +932,28 @@     pokeByteOff p (k+k) r3     pokeByteOff p (3*k) r4 +instance Random Mat4 where+  random = randomR (Mat4 v1 v1 v1 v1, Mat4 v2 v2 v2 v2) where+    v1 = Vec4 (-1) (-1) (-1) (-1)+    v2 = Vec4   1    1    1    1+  randomR (Mat4 a b c d, Mat4 e f g h) gen = +    let (x,gen1) = randomR (a,e) gen+        (y,gen2) = randomR (b,f) gen1+        (z,gen3) = randomR (c,g) gen2  +        (w,gen4) = randomR (d,h) gen3  +    in (Mat4 x y z w, gen4)+    +instance Dimension Mat4 where dim _ = 4+   +instance MatrixNorms Mat4 where +  frobeniusNorm (Mat4 r1 r2 r3 r4) = +    sqrt $+      normsqr r1 + +      normsqr r2 + +      normsqr r3 + +      normsqr r4  +    +-------------------------------------------------------------------------------- -- Extend instances  instance Extend Vec2 Vec3 where@@ -718,17 +972,20 @@   trim (Vec4 x y z _)       = Vec3 x y z  instance Extend Mat2 Mat3 where-  extendZero (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat3 (extendZero p) (extendZero q) zero+  extendWith w (Mat2 p q) = Mat3 (extendZero p) (extendZero q) (Vec3 0 0 w)   trim (Mat3 p q _) = Mat2 (trim p) (trim q)  instance Extend Mat2 Mat4 where-  extendZero (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat2 p q) = Mat4 (extendZero p) (extendZero q) zero zero+  extendWith w (Mat2 p q) = Mat4 (extendZero p) (extendZero q) (Vec4 0 0 w 0) (Vec4 0 0 0 w)   trim (Mat4 p q _ _) = Mat2 (trim p) (trim q)  instance Extend Mat3 Mat4 where-  extendZero (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero-  extendWith _ _ = error "extendWith is meaningless for matrices"+  extendZero   (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) zero+  extendWith w (Mat3 p q r) = Mat4 (extendZero p) (extendZero q) (extendZero r) (Vec4 0 0 0 w)   trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r)+  +--------------------------------------------------------------------------------+
src/flt/GramSchmidt.hs view
@@ -9,7 +9,7 @@  import Data.Vect.Flt.Base --------------------------------------------------------+--------------------------------------------------------------------------------  liftPair :: (a -> b) -> (a,a) -> (b,b) liftPair f (x,y) = (f x, f y)@@ -20,7 +20,7 @@ liftQuadruple :: (a -> b) -> (a,a,a,a) -> (b,b,b,b) liftQuadruple f (x,y,z,w) = (f x, f y, f z, f w) --------------------------------------------------------+--------------------------------------------------------------------------------      -- | produces orthogonal\/orthonormal vectors from a set of vectors     class GramSchmidt a where@@ -32,7 +32,7 @@ "gramSchmidtNormalize is idempotent"  forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a    #-} --------------------------------------------------------+--------------------------------------------------------------------------------  instance GramSchmidt (Vec2,Vec2) where   gramSchmidt = gramSchmidtPair
src/flt/OpenGL.hs view
@@ -2,21 +2,92 @@ -- TODO: the pointer versions of these functions should be really implemented  -- via the pointer versions of the original opengl functions... --- | OpenGL support, inclduing 'vertex', 'texCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.+-- | OpenGL support, including 'Vertex', 'TexCoord', etc instances for 'Vec2', 'Vec3' and 'Vec4'.   module Data.Vect.Flt.OpenGL where  import Control.Monad import Data.Vect.Flt.Base+import Data.Vect.Flt.Util.Projective import qualified Graphics.Rendering.OpenGL as GL  import Foreign -import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale)+import Graphics.Rendering.OpenGL hiding +  ( Normal3 , rotate , translate , scale+  , matrix , currentMatrix , withMatrix , multMatrix +  ) --------------------------------------------------------+-------------------------------------------------------------------------------- -{-# SPECIALISE radianToDegrees :: Float -> Float #-}+-- | There should be a big warning here about the different conventions, +-- hidden transpositions, and all the confusion this will inevitably cause...+--+-- As it stands, +--+-- > glRotate t1 axis1 >> glRotate t2 axis2 >> glRotate t3 axis3+-- +-- has the same result as+--+-- > multMatrix (rotMatrixProj4 t3 axis3 .*. rotMatrixProj4 t2 axis2 .*. rotMatrixProj4 t1 axis1)+--+-- because at the interface of OpenGL and this library there is a transposition+-- to compensate for the different conventions. (This transposition is implicit+-- in the code, because the way the matrices are stored in the memory is also+-- different: OpenGL stores them column-major, and we store them row-major).++class ToOpenGLMatrix m where+  makeGLMatrix :: m -> IO (GLmatrix Flt)++class FromOpenGLMatrix m where+  peekGLMatrix :: GLmatrix Flt -> IO m+  +setMatrix :: ToOpenGLMatrix m => Maybe MatrixMode -> m -> IO ()+setMatrix mode m = makeGLMatrix m >>= \x -> GL.matrix mode $= x+ +getMatrix :: FromOpenGLMatrix m => Maybe MatrixMode -> IO m+getMatrix mode = get (GL.matrix mode) >>= peekGLMatrix++matrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => Maybe MatrixMode -> StateVar m+matrix mode = makeStateVar (getMatrix mode) (setMatrix mode)++currentMatrix :: (ToOpenGLMatrix m, FromOpenGLMatrix m) => StateVar m+currentMatrix = matrix Nothing++multMatrix :: ToOpenGLMatrix m => m -> IO ()+multMatrix m = makeGLMatrix m >>= GL.multMatrix++instance ToOpenGLMatrix Mat4 where+  makeGLMatrix m = GL.withNewMatrix GL.ColumnMajor (flip poke m . castPtr) + +instance FromOpenGLMatrix Mat4 where+  -- huh? GL.withMatrix is strange+  peekGLMatrix x = GL.withMatrix x $ \_ p -> peek (castPtr p)+  +instance ToOpenGLMatrix Mat3 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)+ +instance ToOpenGLMatrix Mat2 where+  makeGLMatrix m = makeGLMatrix (extendWith 1 m :: Mat4)++instance ToOpenGLMatrix Ortho4 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat4)++instance ToOpenGLMatrix Ortho3 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat3)++instance ToOpenGLMatrix Ortho2 where+  makeGLMatrix m = makeGLMatrix (fromOrtho m :: Mat2)++instance ToOpenGLMatrix Proj4 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat4)++instance ToOpenGLMatrix Proj3 where+  makeGLMatrix m = makeGLMatrix (fromProjective m :: Mat3)+  +--------------------------------------------------------------------------------++{-# SPECIALISE radianToDegrees :: Float  -> Float  #-} {-# SPECIALISE radianToDegrees :: Double -> Double #-} radianToDegrees :: RealFrac a => a -> a radianToDegrees x = x * 57.295779513082322@@ -27,20 +98,61 @@ degreesToRadian x = x * 1.7453292519943295e-2  -- | The angle is in radians. (WARNING: OpenGL uses degrees!)-rotate :: Flt -> Vec3 -> IO ()-rotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z)+glRotate :: Flt -> Vec3 -> IO ()+glRotate angle (Vec3 x y z) = GL.rotate (radianToDegrees angle) (Vector3 x y z) -translate :: Vec3 -> IO ()-translate (Vec3 x y z) = GL.translate (Vector3 x y z)+glTranslate :: Vec3 -> IO ()+glTranslate (Vec3 x y z) = GL.translate (Vector3 x y z) -scale3 :: Vec3 -> IO ()-scale3 (Vec3 x y z) = GL.scale x y z+glScale3 :: Vec3 -> IO ()+glScale3 (Vec3 x y z) = GL.scale x y z -scale :: Flt -> IO ()-scale x = GL.scale x x x+glScale :: Flt -> IO ()+glScale x = GL.scale x x x --------------------------------------------------------+--------------------------------------------------------------------------------+ +-- | \"Orthogonal projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1)+orthoMatrix +  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +orthoMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2/(r-l)) 0 0 0)+  (Vec4 0 (2/(t-b)) 0 0)+  (Vec4 0 0 (-2/(f-n)) 0)+  (Vec4 (-(r+l)/(r-l)) (-(t+b)/(t-b)) (-(f+n)/(f-n)) 1)+  +-- | The same as "orthoMatrix", but with a different parametrization.+orthoMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +orthoMatrix2 (Vec3 l t n) (Vec3 r b f) = orthoMatrix (l,r) (b,t) (n,f) +-- | \"Perspective projecton\" matrix, a la OpenGL +-- (the corresponding functionality is removed in OpenGL 3.1).+frustumMatrix+  :: (Flt,Flt)   -- ^ (left,right)+  -> (Flt,Flt)   -- ^ (bottom,top)+  -> (Flt,Flt)   -- ^ (near,far)+  -> Mat4 +frustumMatrix (l,r) (b,t) (n,f) = Mat4+  (Vec4 (2*n/(r-l)) 0 0 0)+  (Vec4 0 (2*n/(t-b)) 0 0)+  (Vec4 ((r+l)/(r-l)) ((t+b)/(t-b)) (-(f+n)/(f-n)) (-1))+  (Vec4 0 0 (-2*f*n*(f-n)) 0)+  +-- | The same as "frustumMatrix", but with a different parametrization.+frustumMatrix2 {- ' CPP is sensitive to primes -}+  :: Vec3     -- ^ (left,top,near)+  -> Vec3     -- ^ (right,bottom,far)+  -> Mat4 +frustumMatrix2 (Vec3 l t n) (Vec3 r b f) = frustumMatrix (l,r) (b,t) (n,f)++-------------------------------------------------------------------------------- -- Vertex instances  instance GL.Vertex Vec2 where@@ -55,17 +167,20 @@   vertex (Vec4 x y z w) = GL.vertex (GL.Vertex4 x y z w)   vertexv p = peek p >>= vertex    ---------------------------------------------------------+-------------------------------------------------------------------------------- -- the Normal instance -- note that there is no Normal2\/Normal4 in the OpenGL binding  instance GL.Normal Normal3 where-  normal (Normal3 (Vec3 x y z)) = GL.normal (GL.Normal3 x y z)+  normal u = GL.normal (GL.Normal3 x y z) +    where Vec3 x y z = fromNormal u    normalv p = peek p >>= normal  --------------------------------------------------------+instance GL.Normal Vec3 where+  normal (Vec3 x y z) = GL.normal (GL.Normal3 x y z) +  normalv p = peek p >>= normal  +-------------------------------------------------------------------------------- -- Color instances    instance GL.Color Vec3 where@@ -87,8 +202,7 @@   secondaryColorv p = peek p >>= secondaryColor -} ---------------------------------------------------------+-------------------------------------------------------------------------------- -- TexCoord instances  instance GL.TexCoord Vec2 where@@ -109,8 +223,7 @@   multiTexCoord unit (Vec4 u v w z) = GL.multiTexCoord unit (GL.TexCoord4 u v w z)   multiTexCoordv unit p = peek p >>= multiTexCoord unit ---------------------------------------------------------    +-------------------------------------------------------------------------------- -- Vertex Attributes (experimental)  class VertexAttrib' a where@@ -129,16 +242,18 @@   vertexAttrib loc (Vec4 x y z w) = GL.vertexAttrib4 loc x y z w   instance VertexAttrib' Normal2 where-  vertexAttrib loc (Normal2 (Vec2 x y)) = GL.vertexAttrib2 loc x y+  vertexAttrib loc u = GL.vertexAttrib2 loc x y+    where Vec2 x y = fromNormal u   instance VertexAttrib' Normal3 where-  vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z+  vertexAttrib loc u = GL.vertexAttrib3 loc x y z+    where Vec3 x y z = fromNormal u   instance VertexAttrib' Normal4 where-  vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w----------------------------------------------------------+  vertexAttrib loc u = GL.vertexAttrib4 loc x y z w+    where Vec4 x y z w = fromNormal u +   +-------------------------------------------------------------------------------- -- Uniform (again, experimental)  -- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; @@ -180,3 +295,5 @@   uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt)      #endif++
src/flt/Util/Dim2.hs view
@@ -3,12 +3,13 @@  import Data.Vect.Flt.Base --- |example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@.+-- | Example: @structVec2 [1,2,3,4] = [ Vec2 1 2 , Vec2 3 4 ]@. structVec2 :: [Flt] -> [Vec2] structVec2 [] = [] structVec2 (x:y:ls) = (Vec2 x y):(structVec2 ls)  structVec2 _ = error "structVec2" +-- | The opposite of "structVec2". destructVec2 :: [Vec2] -> [Flt] destructVec2 [] = [] destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls)  @@ -48,10 +49,13 @@ angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt angle2' = angle2 . fromNormal --- |Rotation matrix by a given angle (in radians), counterclockwise.+-- | Rotation matrix by a given angle (in radians), counterclockwise. rotMatrix2 :: Flt -> Mat2 rotMatrix2 a = Mat2 (Vec2 c s) (Vec2 (-s) c) where c = cos a; s = sin a +rotMatrixOrtho2 :: Flt -> Ortho2+rotMatrixOrtho2 = toOrthoUnsafe . rotMatrix2+ rotate2 :: Flt -> Vec2 -> Vec2 rotate2 a v = v .* (rotMatrix2 a)  @@ -62,3 +66,5 @@ -- |Rotates clockwise by 90 degrees. rotateCW :: Vec2 -> Vec2 rotateCW (Vec2 x y) = Vec2 y (-x)++
src/flt/Util/Dim3.hs view
@@ -3,18 +3,26 @@  import Data.Vect.Flt.Base +--------------------------------------------------------------------------------++-- | Example: @structVec3 [1,2,3,4,5,6] = [ Vec3 1 2 3 , Vec3 4 5 6]@. structVec3 :: [Flt] -> [Vec3] structVec3 [] = [] structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls)  structVec3 _ = error "structVec3" +-- | The opposite of "structVec3". destructVec3 :: [Vec3] -> [Flt] destructVec3 [] = [] destructVec3 ((Vec3 x y z):ls) = x:y:z:(destructVec3 ls)   +--------------------------------------------------------------------------------+ det3 :: Vec3 -> Vec3 -> Vec3 -> Flt det3 u v w = det (u,v,w) +--------------------------------------------------------------------------------+ translate3X :: Flt -> Vec3 -> Vec3 translate3Y :: Flt -> Vec3 -> Vec3 translate3Z :: Flt -> Vec3 -> Vec3@@ -42,26 +50,34 @@ rotMatrixY a = Mat3 (Vec3 c 0 (-s)) (Vec3 0 1 0) (Vec3 s 0 c) where c = cos a; s = sin a rotMatrixX a = Mat3 (Vec3 1 0 0) (Vec3 0 c s) (Vec3 0 (-s) c) where c = cos a; s = sin a -rotate3' :: {- ' CPP is sensitive to primes -} Flt       -- ^ angle (in radians)-         -> Normal3   -- ^ axis (should be a /unit/ vector!) -         -> Vec3      -- ^ vector-         -> Vec3      -- ^ result+--------------------------------------------------------------------------------++rotate3' {- ' CPP is sensitive to primes -} +  :: Flt       -- ^ angle (in radians)+  -> Normal3   -- ^ axis (should be a /unit/ vector!) +  -> Vec3      -- ^ vector+  -> Vec3      -- ^ result rotate3' angle axis v = v .* (rotMatrix3' axis angle) -rotate3 :: Flt    -- ^ angle (in radians)-        -> Vec3   -- ^ axis (arbitrary nonzero vector)-        -> Vec3   -- ^ vector-        -> Vec3   -- ^ result+rotate3 +  :: Flt    -- ^ angle (in radians)+  -> Vec3   -- ^ axis (arbitrary nonzero vector)+  -> Vec3   -- ^ vector+  -> Vec3   -- ^ result rotate3 angle axis v = v .* (rotMatrix3 axis angle)       --- |Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +-- | Rotation around an arbitrary 3D vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3 :: Vec3 -> Flt -> Mat3 rotMatrix3 v a = rotMatrix3' (mkNormal v) a --- |Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/. +rotMatrixOrtho3 :: Vec3 -> Flt -> Ortho3+rotMatrixOrtho3 v a = toOrthoUnsafe $ rotMatrix3 v a++-- | Rotation around an arbitrary 3D /unit/ vector. The resulting 3x3 matrix is intended for multiplication on the /right/.  rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3-rotMatrix3' (Normal3 v) a = -  let c = cos a+rotMatrix3' u a = +  let v = fromNormal u+      c = cos a       s = sin a       m1 = scalarMul (1-c) (outer v v)       x = _1 v@@ -72,4 +88,65 @@                 (Vec3 ( s*y) (-s*x)   c   )   in (m1 &+ m2) +rotMatrixOrtho3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Ortho3+rotMatrixOrtho3' u a = toOrthoUnsafe $ rotMatrix3' u a +--------------------------------------------------------------------------------++-- | Reflects a vector to an axis: that is, the result of @reflect n v@ is+-- 2\<n,v\>n - v+reflect :: Normal3 -> Vec3 -> Vec3+reflect u v = (s *& n) &- v where +  n = fromNormal u+  s = 2 * (n &. v)++reflect' :: Normal3 -> Normal3 -> Normal3+reflect' u x = toNormalUnsafe $ reflect u (fromNormal x)+  +refract :: Flt -> Normal3 -> Vec3 -> Vec3+refract eta u v = s *& fromNormal w where+  s = norm v +  w = refract' eta u (toNormalUnsafe $ v &* (1.0/s))+  +-- | Refraction.+-- First parameter (@eta@) is the relative refraction index +--+-- >        refl_inside+-- > eta = --------------+-- >        refl_outside+--+-- where \"inside\" is the direction of the second argument +-- (to vector normal to plane which models the boundary +-- between the two materials). That is, total internal reflection+-- can occur when @eta>1@.+--+-- The convention is that the origin is the point of intersection+-- of the ray and the surface, and all the vectors \"point away\"+-- from here (unlike, say, GLSL's @refract@, where the incident+-- vector \"points towards\" the material)+refract' {- ' CPP is sensitive to primes -} +  :: Flt -> Normal3 -> Normal3 -> Normal3+refract' eta u i = +  if k<0+    then reflect' u i +    else toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++-- | When total internal reflection would occur, we return "Nothing".+refractOnly' {- ' CPP is sensitive to primes -}  +  :: Flt -> Normal3 -> Normal3 -> Maybe Normal3+refractOnly' eta u i = +  if k<0+    then Nothing +    else Just $ toNormalUnsafe $ ((-eta) *& v) &- (- eta*c + sqrt k) *& n +  where+    n = fromNormal u+    v = fromNormal i+    c = n &. v+    k = 1 - eta*eta*(1-c*c)++--------------------------------------------------------------------------------
src/flt/Util/Dim4.hs view
@@ -39,7 +39,7 @@ vec4Z = Vec4 0 0 1 0 vec4W = Vec4 0 0 0 1 ----------------------------------------------------------------------------+--------------------------------------------------------------------------------  -- |If @(x,y,u,v)@ is an orthonormal system, then (written in pseudo-code) -- @biVector4 (x,y) = plusMinus (reverse $ biVector4 (u,v))@.@@ -75,8 +75,9 @@ -- | Rotation matrix around a plane specified by two normalized and /orthogonal/ vectors. -- Intended for multiplication on the /right/! rotMatrix4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Mat4-rotMatrix4' angle (Normal4 v, Normal4 w) = m1 &+ (s *& m2) &+ m3 +rotMatrix4' angle (u1,u2) = m1 &+ (s *& m2) &+ m3    where+    v = fromNormal u1 ; w = fromNormal u2     c = cos angle ; s = sin angle     m1 = scalarMul (1-c) ( outer v v  &+  outer w w )     m2 = biVector4AsTensor v w
src/flt/Util/Projective.hs view
@@ -1,5 +1,6 @@ --- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on+-- | Classic 4x4 projective matrices, encoding the affine transformations of R^3.+-- Our convention is that they are intended for multiplication on -- the /right/, that is, they are of the form -- -- >     _____@@ -12,6 +13,8 @@ -- store them by rows; but OpenGL also use the opposite convention (so the OpenGL projective matrices  -- are intended for multiplication on the /left/). So in effect, they are the same when stored in the memory, -- say with @poke :: Ptr Mat4 -> Mat4 -> IO ()@.+--+-- Warning: The naming conventions will probably change in the future.  module Data.Vect.Flt.Util.Projective where @@ -20,68 +23,73 @@  import qualified Data.Vect.Flt.Util.Dim4 as Dim4 -class ExtendProjective v e | v->e where-  extendProj     :: v -> e-  extendProjWith :: Flt -> v -> e-  extendProj = extendProjWith 1-  -instance ExtendProjective Vec2 Vec4 where-  extendProj       (Vec2 x y) = Vec4 x y 0 1-  extendProjWith w (Vec2 x y) = Vec4 x y 0 w-  -instance ExtendProjective Vec3 Vec4 where-  extendProj       (Vec3 x y z) = Vec4 x y z 1-  extendProjWith w (Vec3 x y z) = Vec4 x y z w--instance ExtendProjective Vec4 Vec4 where-  extendProj = id-  extendProjWith w (Vec4 x y z w') = let s = w/w' in Vec4 (s*x) (s*y) (s*z) w+-------------------------------------------------------------------------------- -instance ExtendProjective Mat2 Mat4 where-  extendProj       (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 1)-  extendProjWith w (Mat2 r1 r2) = Mat4 (extendZero r1) (extendZero r2) (Dim4.vec4Z) (Vec4 0 0 0 w)+rotMatrixProj4' :: {- ' CPP is sensitive to primes -}  Flt -> Normal3 -> Proj4+rotMatrixProj4' angle axis = linear $ rotMatrix3' axis angle -instance ExtendProjective Mat3 Mat4 where-  extendProj       (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 1)-  extendProjWith w (Mat3 r1 r2 r3) = Mat4 (extendZero r1) (extendZero r2) (extendZero r3) (Vec4 0 0 0 w)+rotMatrixProj4 :: Flt -> Vec3 -> Proj4+rotMatrixProj4 angle axis = linear $ rotMatrix3 axis angle -rotMatrixProj :: Flt -> Normal3 -> Mat4-rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle+-- | synonym for "rotateAfterProj4"+rotateProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateProj4 = rotateAfterProj4 -rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4-rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle+-- | Synonym for @\m -> m .*. rotMatrixProj4 angle axis@.+rotateAfterProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateAfterProj4 angle axis m = m .*. (rotMatrixProj4' angle axis)  -translMatrixProj :: Vec3 -> Mat4-translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v)+-- | Synonym for @\m -> rotMatrixProj4 angle axis .*. m@.+rotateBeforeProj4 :: Flt -> Normal3 -> Proj4 -> Proj4+rotateBeforeProj4 angle axis m = (rotMatrixProj4' angle axis) .*. m  --- | we assume that the bottom-right corner is 1.-translWithProj :: Vec3 -> Mat4 -> Mat4-translWithProj v mat@(Mat4 r1 r2 r3 r4) = Mat4 r1 r2 r3 (extendProjWith 0 v &+ r4)+--------------- -scaleMatrixProj :: Vec3 -> Mat4-scaleMatrixProj v = diag $ extendProj v+--scalingUniformProj3 :: Flt -> Proj3+--scalingUniformProj3 x = scaling (Vec2 x x) -scaleMatrixUniformProj :: Flt -> Mat4-scaleMatrixUniformProj s = diag (Vec4 s s s 1)+scalingUniformProj4 :: Flt -> Proj4+scalingUniformProj4 x = scaling (Vec3 x x x) -class ProjectiveAction v where-  actProj :: v -> Mat4 -> v- -instance ProjectiveAction Vec3 where-  actProj v m = trim $ (extendProj v) .* m +-- | Equivalent to @\m -> scaling v .*. m@.+scaleBeforeProj4 :: Vec3 -> Proj4 -> Proj4+scaleBeforeProj4 (Vec3 u v w) p4 =  +  toProjectiveUnsafe $ +    Mat4 (u*&a) (v*&b) (w*&c) t+  where+    Mat4 a b c t = fromProjective p4 -instance ProjectiveAction Vec4 where-  actProj v m = v .* m +-- | Equivalent to @\m -> m .*. scaling v@.+scaleAfterProj4 :: Vec3 -> Proj4 -> Proj4+scaleAfterProj4 v p4 =+  toProjectiveUnsafe $ +    Mat4 (a&!w) (b&!w) (c&!w) (t&!w)+  where+    w = extendWith 1 v+    Mat4 a b c t = fromProjective p4+    +--------------- --- | When acting on unit vectors, we ignore the translation part.-instance ProjectiveAction Normal3 where-  actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3))+-- | Synonym for "translateAfter4"+translate4 :: Vec3 -> Proj4 -> Proj4+translate4 = translateAfter4 --- | Inverts a projective 4x4 matrix, assuming that the top-left 3x3 part is /orthogonal/,--- and the bottom-right corner is 1.-invertProj :: Mat4 -> Mat4-invertProj mat@(Mat4 u v w t) = -  translWithProj t' $ extendProj $ transpose $ (trim mat :: Mat3)+-- | Equivalent to @\m -> m .*. translation v@.+translateAfter4 :: Vec3 -> Proj4 -> Proj4+translateAfter4 v p4 = +  toProjectiveUnsafe $+    Mat4 r1 r2 r3 (extendWith 0 v &+ r4)   where-    t' = Vec3 (- u &. t) (- v &. t) (- w &. t)-    +    Mat4 r1 r2 r3 r4 = fromProjective p4 ++-- | Equivalent to @\m -> translation v .*. m@.+translateBefore4 :: Vec3 -> Proj4 -> Proj4+translateBefore4 v p4 = +  toProjectiveUnsafe $ +    Mat4 r1 r2 r3 (extendWith 0 u &+ r4) +  where +   u = v .* (trim mat :: Mat3) +   mat@(Mat4 r1 r2 r3 r4) = fromProjective p4+   +---------------+
vect.cabal view
@@ -1,5 +1,5 @@ Name:                vect-Version:             0.4.0+Version:             0.4.5 Synopsis:            A low-dimensional linear algebra library, tailored to computer graphics. Description:         A low-dimensional (2, 3 and 4) linear algebra library,                       with lots of useful functions. Intended usage is primarily