diff --git a/Data/Vect.hs b/Data/Vect.hs new file mode 100644 --- /dev/null +++ b/Data/Vect.hs @@ -0,0 +1,5 @@ + +-- | Importing this module is equivalent to importing "Data.Vect.Float". +module Data.Vect ( module Data.Vect.Float ) where +import Data.Vect.Float + diff --git a/Data/Vect/Double.hs b/Data/Vect/Double.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double.hs @@ -0,0 +1,23 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +module Data.Vect.Flt + ( module Data.Vect.Flt.Base + , module Data.Vect.Flt.Interpolate + , module Data.Vect.Flt.Util.Dim2 + , module Data.Vect.Flt.Util.Dim3 + , module Data.Vect.Flt.Util.Projective +#ifdef VECT_OPENGL + , module Data.Vect.Flt.OpenGL +#endif + ) where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Interpolate + +import Data.Vect.Flt.Util.Dim2 +import Data.Vect.Flt.Util.Dim3 +import Data.Vect.Flt.Util.Projective + +#ifdef VECT_OPENGL +import Data.Vect.Flt.OpenGL +#endif diff --git a/Data/Vect/Double/Base.hs b/Data/Vect/Double/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Base.hs @@ -0,0 +1,735 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +module Data.Vect.Flt.Base where + +import Control.Monad +import System.Random +import Foreign + +-- class declarations + +class AbelianGroup g where + (&+) :: g -> g -> g + (&-) :: g -> g -> g + neg :: g -> g + zero :: g + +infixl 6 &+ +infixl 6 &- + +vecSum :: AbelianGroup g => [g] -> g +vecSum l = foldl (&+) zero l + +class (AbelianGroup r) => + Ring r where + (.*.) :: r -> r -> r + one :: r + +infixl 7 .*. + +ringProduct :: Ring r => [r] -> r +ringProduct l = foldl (.*.) one l + +class LeftModule r m where + lmul :: r -> m -> m + (*.) :: r -> m -> m + (*.) = lmul + +class RightModule m r where + rmul :: m -> r -> m + (.*) :: m -> r -> m + (.*) = rmul + +-- I'm not really sure about this.. may actually degrade the performance in some cases? +{- RULES +"matrix multiplication left" forall m n x. (n .*. m) *. x = n *. (m *. x) +"matrix multiplication right" forall m n x. x .* (m .*. n) = (x .* m) .* n + -} + +infixr 7 *. +infixl 7 .* + +class AbelianGroup v => Vector v where + mapVec :: (Flt -> Flt) -> v -> v + scalarMul :: Flt -> v -> v + (*&) :: Flt -> v -> v + (&*) :: v -> Flt -> v + (*&) s v = scalarMul s v + (&*) v s = scalarMul s v + +infixr 7 *& +infixl 7 &* + +{-# RULES +"scalar multiplication left" forall s t x. t *& (s *& x) = (t*s) *& x +"scalar multiplication right" forall s t x. (x &* s) &* t = x &* (s*t) + #-} + +class DotProd v where + (&.) :: v -> v -> Flt + norm :: v -> Flt + normsqr :: v -> Flt + len :: v -> Flt + lensqr :: v -> Flt + len = norm + lensqr = normsqr + dotprod :: v -> v -> Flt + normsqr v = (v &. v) + norm = sqrt.lensqr + dotprod = (&.) + +infix 7 &. + +{-# RULES +"len/square 1" forall x. (len x)*(len x) = lensqr x +"len/square 2" forall x. (len x)^2 = lensqr x +"norm/square 1" forall x. (norm x)*(norm x) = normsqr x +"norm/square 2" forall x. (norm x)^2 = normsqr x + #-} + +normalize :: (Vector v, DotProd v) => v -> v +normalize v = scalarMul (1.0/(len v)) v + +distance :: (Vector v, DotProd v) => v -> v -> Flt +distance x y = norm (x &- y) + +-- | the angle between two vectors +angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y) + +-- | the angle between two unit vectors +angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y) + +{-# RULES +"normalize is idempotent" forall x. normalize (normalize x) = normalize x + #-} + +class (Vector v, DotProd v) => UnitVector v u | v->u, u->v where + mkNormal :: v -> u -- ^ normalizes the input + toNormalUnsafe :: v -> u -- ^ does not normalize the input! + fromNormal :: u -> v + fromNormalRadius :: Flt -> u -> v + fromNormalRadius t n = t *& fromNormal n + +-- | 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! +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. +flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal + +class CrossProd v where + crossprod :: v -> v -> v + (&^) :: v -> v -> v + (&^) = crossprod + +class Pointwise v where + pointwise :: v -> v -> v + (&!) :: v -> v -> v + (&!) = pointwise + +infix 7 &^ +infix 7 &! + +class HasCoordinates v x | v->x where + _1 :: v -> x + _2 :: v -> x + _3 :: v -> x + _4 :: v -> x + +-- | conversion between vectors (and matrices) of different dimensions +class Extend u v where + extendZero :: u -> v -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@ + extendWith :: Flt -> u -> v -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@ + trim :: v -> u -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@ + +-- | makes a diagonal matrix from a vector +class Diagonal s t | t->s where + diag :: s -> t + +class Matrix m where + transpose :: m -> m + inverse :: m -> m + idmtx :: m + +{-# RULES +"transpose is an involution" forall m. transpose (transpose m) = m +"inverse is an involution" forall m. inverse (inverse m) = m + #-} + +-- | 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 + +-- Vec / Mat datatypes + +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) + +-- | these 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) + +mkVec2 :: (Flt,Flt) -> Vec2 +mkVec3 :: (Flt,Flt,Flt) -> Vec3 +mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4 + +mkVec2 (x,y) = Vec2 x y +mkVec3 (x,y,z) = Vec3 x y z +mkVec4 (x,y,z,w) = Vec4 x y z w + +-- Unit vectors + +instance UnitVector Vec2 Normal2 where + mkNormal v = Normal2 (normalize v) + fromNormal (Normal2 v) = v + toNormalUnsafe = Normal2 + +instance UnitVector Vec3 Normal3 where + mkNormal v = Normal3 (normalize v) + fromNormal (Normal3 v) = v + toNormalUnsafe = Normal3 + +instance UnitVector Vec4 Normal4 where + mkNormal v = Normal4 (normalize v) + fromNormal (Normal4 v) = v + toNormalUnsafe = Normal4 + +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 + where + (v,h) = random g + d = norm v + +instance Random Normal2 where + 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) + randomR _ = random + +instance Random Normal4 where + 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 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 +-} + +-- Vec2 instances + +instance HasCoordinates Vec2 Flt where + _1 (Vec2 x _) = x + _2 (Vec2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance AbelianGroup Vec2 where + (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) + (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2) + neg (Vec2 x y) = Vec2 (-x) (-y) + zero = Vec2 0 0 + +instance Vector Vec2 where + scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y) + mapVec f (Vec2 x y) = Vec2 (f x) (f y) + +instance DotProd Vec2 where + (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2 + +instance Pointwise Vec2 where + pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2) + +instance Determinant (Vec2,Vec2) where + det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1 + +{- +instance Show Vec2 where + show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )" +-} + +instance Random Vec2 where + random = randomR (Vec2 (-1) (-1),Vec2 1 1) + randomR (Vec2 a b, Vec2 c d) gen = + let (x,gen1) = randomR (a,c) gen + (y,gen2) = randomR (b,d) gen1 + in (Vec2 x y, gen2) + +instance Storable Vec2 where + sizeOf _ = 2 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p k + return (Vec2 x y) + + poke q (Vec2 x y) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p k y + +-- Mat2 instances + +instance HasCoordinates Mat2 Vec2 where + _1 (Mat2 x _) = x + _2 (Mat2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance Matrix Mat2 where + transpose (Mat2 row1 row2) = + Mat2 (Vec2 (_1 row1) (_1 row2)) + (Vec2 (_2 row1) (_2 row2)) + idmtx = Mat2 (Vec2 1 0) (Vec2 0 1) + inverse (Mat2 (Vec2 a b) (Vec2 c d)) = + Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) + where r = 1.0 / (a*d - b*c) + +instance AbelianGroup Mat2 where + (&+) (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) + +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 + (.*.) (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 LeftModule Mat2 Vec2 where + lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) + +instance RightModule Vec2 Mat2 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec2 Mat2 where + diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y) + +instance Tensor Mat2 Vec2 where + 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 + +{- +instance Show Mat2 where + show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2 +-} + +instance Storable Mat2 where + sizeOf _ = 2 * sizeOf (undefined::Vec2) + alignment _ = alignment (undefined::Vec2) + + peek q = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + r1 <- peek p + r2 <- peekByteOff p k + return (Mat2 r1 r2) + + poke q (Mat2 r1 r2) = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + poke p r1 + pokeByteOff p k r2 + +-- Vec3 instances + +instance HasCoordinates Vec3 Flt where + _1 (Vec3 x _ _) = x + _2 (Vec3 _ y _) = y + _3 (Vec3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance AbelianGroup Vec3 where + (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) + (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) + neg (Vec3 x y z) = Vec3 (-x) (-y) (-z) + zero = Vec3 0 0 0 + +instance Vector Vec3 where + scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z) + mapVec f (Vec3 x y z) = Vec3 (f x) (f y) (f z) + +instance DotProd Vec3 where + (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2 + +instance Pointwise Vec3 where + pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2) + +{- +instance Show Vec3 where + show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )" +-} + +instance Random Vec3 where + random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1) + randomR (Vec3 a b c, Vec3 d e f) gen = + let (x,gen1) = randomR (a,d) gen + (y,gen2) = randomR (b,e) gen1 + (z,gen3) = randomR (c,f) gen2 + in (Vec3 x y z, gen3) + +instance CrossProd Vec3 where + crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) + +instance Determinant (Vec3,Vec3,Vec3) where + det (u,v,w) = u &. (v &^ w) + +instance Storable Vec3 where + sizeOf _ = 3 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + return (Vec3 x y z) + + poke q (Vec3 x y z) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + +-- Mat3 instances + +instance HasCoordinates Mat3 Vec3 where + _1 (Mat3 x _ _) = x + _2 (Mat3 _ y _) = y + _3 (Mat3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance Matrix Mat3 where + + transpose (Mat3 row1 row2 row3) = + Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) + (Vec3 (_2 row1) (_2 row2) (_2 row3)) + (Vec3 (_3 row1) (_3 row2) (_3 row3)) + + idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1) + + inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = + Mat3 (Vec3 (d11*r) (d21*r) (d31*r)) + (Vec3 (d12*r) (d22*r) (d32*r)) + (Vec3 (d13*r) (d23*r) (d33*r)) + where + r = 1.0 / ( a*d11 + b*d12 + c*d13 ) + + d11 = f*k - g*j + d12 = g*i - e*k + d13 = e*j - f*i + + d31 = b*g - c*f + d32 = c*e - a*g + d33 = a*f - b*e + + d21 = c*j - b*k + d22 = a*k - c*i + d23 = b*i - a*j + +instance AbelianGroup Mat3 where + (&+) (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) + +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 + (.*.) (Mat3 r1 r2 r3) n = + let (Mat3 c1 c2 c3) = transpose n + in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3)) + (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3)) + (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3)) + one = idmtx + +instance LeftModule Mat3 Vec3 where + lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v) + +instance RightModule Vec3 Mat3 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec3 Mat3 where + diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z) + +instance Tensor Mat3 Vec3 where + outer (Vec3 a b c) (Vec3 x y z) = Mat3 + (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) + +{- +instance Show Mat3 where + show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 +-} + +instance Storable Mat3 where + sizeOf _ = 3 * sizeOf (undefined::Vec3) + alignment _ = alignment (undefined::Vec3) + + peek q = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + return (Mat3 r1 r2 r3) + + poke q (Mat3 r1 r2 r3) = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + +-- Vec4 instances + +instance HasCoordinates Vec4 Flt where + _1 (Vec4 x _ _ _) = x + _2 (Vec4 _ y _ _) = y + _3 (Vec4 _ _ z _) = z + _4 (Vec4 _ _ _ w) = w + +instance AbelianGroup Vec4 where + (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2) + (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2) + neg (Vec4 x y z w) = Vec4 (-x) (-y) (-z) (-w) + zero = Vec4 0 0 0 0 + +instance Vector Vec4 where + scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w) + mapVec f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w) + +instance DotProd Vec4 where + (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2 + +instance Pointwise Vec4 where + pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2) + +{- +instance Show Vec4 where + show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )" +-} + +instance Random Vec4 where + random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1) + randomR (Vec4 a b c d, Vec4 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 (Vec4 x y z w, gen4) + +instance Storable Vec4 where + sizeOf _ = 4 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + w <- peekByteOff p (3*k) + return (Vec4 x y z w) + + poke q (Vec4 x y z w) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + pokeByteOff p (3*k) w + +-- Mat4 instances + +instance HasCoordinates Mat4 Vec4 where + _1 (Mat4 x _ _ _) = x + _2 (Mat4 _ y _ _) = y + _3 (Mat4 _ _ z _) = z + _4 (Mat4 _ _ _ w) = w + +instance Matrix Mat4 where + transpose (Mat4 row1 row2 row3 row4) = + Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) + (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) + (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) + (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) + idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1) + inverse = error "inverse/Mat4: not implemented yet" + +instance AbelianGroup Mat4 where + (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4) + (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4) + neg (Mat4 r1 r2 r3 r4) = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) + zero = Mat4 zero zero zero zero + +instance Vector Mat4 where + 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 + (.*.) (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)) + (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4)) + (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4)) + (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4)) + one = idmtx + +instance LeftModule Mat4 Vec4 where + lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v) + +instance RightModule Vec4 Mat4 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec4 Mat4 where + diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w) + +instance Tensor Mat4 Vec4 where + outer (Vec4 a b c d) (Vec4 x y z w) = Mat4 + (Vec4 (a*x) (a*y) (a*z) (a*w)) + (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 Show Mat4 where + show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4 +-} + +instance Storable Mat4 where + sizeOf _ = 4 * sizeOf (undefined::Vec4) + alignment _ = alignment (undefined::Vec4) + + peek q = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + r4 <- peekByteOff p (3*k) + return (Mat4 r1 r2 r3 r4) + + poke q (Mat4 r1 r2 r3 r4) = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + pokeByteOff p (3*k) r4 + +-- Extend instances + +instance Extend Vec2 Vec3 where + extendZero (Vec2 x y) = Vec3 x y 0 + extendWith t (Vec2 x y) = Vec3 x y t + trim (Vec3 x y _) = Vec2 x y + +instance Extend Vec2 Vec4 where + extendZero (Vec2 x y) = Vec4 x y 0 0 + extendWith t (Vec2 x y) = Vec4 x y t t + trim (Vec4 x y _ _) = Vec2 x y + +instance Extend Vec3 Vec4 where + extendZero (Vec3 x y z) = Vec4 x y z 0 + extendWith t (Vec3 x y z) = Vec4 x y z t + 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" + 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" + 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" + trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r) + diff --git a/Data/Vect/Double/GramSchmidt.hs b/Data/Vect/Double/GramSchmidt.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/GramSchmidt.hs @@ -0,0 +1,135 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +-- | Gram-Schmidt orthogonalization. +-- This module is not re-exported by "Data.Vect". + +module Data.Vect.Flt.GramSchmidt + ( GramSchmidt(..) + ) + where + +import Data.Vect.Flt.Base + +------------------------------------------------------- + +liftPair :: (a -> b) -> (a,a) -> (b,b) +liftPair f (x,y) = (f x, f y) + +liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b) +liftTriple f (x,y,z) = (f x, f y, f z) + +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 + gramSchmidt :: a -> a -- ^ does not normalize the vectors! + gramSchmidtNormalize :: a -> a -- ^ normalizes the vectors. + +{-# RULES +"gramSchmidt is idempotent" forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent" forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a + #-} + +------------------------------------------------------- + +instance GramSchmidt (Vec2,Vec2) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec3,Vec3) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec4,Vec4) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +---------- + +instance GramSchmidt (Normal2,Normal2) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal2!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +---------- + +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtPair (u,v) = (u',v') where + u' = u + v' = project v u' + +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtNormalizePair (u,v) = (u',v') where + u' = normalize u + v' = normalize $ projectUnsafe v u' + +---------- + +instance GramSchmidt (Vec3,Vec3,Vec3) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Normal3,Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +instance GramSchmidt (Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +---------- + +gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtTriple (u,v,w) = (u',v',w') where + u' = u + v' = project v u' + w' = project (project w u') v' + +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + +---------- + +instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtQuadruple + gramSchmidtNormalize = gramSchmidtNormalizeQuadruple + +instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal + +---------- + +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where + u' = u + v' = project v u' + w' = project (project w u') v' + z' = project (project (project z u') v') w' + +gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w' + +---------- + diff --git a/Data/Vect/Double/Interpolate.hs b/Data/Vect/Double/Interpolate.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Interpolate.hs @@ -0,0 +1,43 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors. + +module Data.Vect.Flt.Interpolate where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim2 (sinCos',angle2') +import Data.Vect.Flt.Util.Dim3 (rotate3') + +class Interpolate v where + interpolate :: Flt -> v -> v -> v + +instance Interpolate Flt where + interpolate t x y = x + t*(y-x) + +instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x) + +instance Interpolate Normal2 where + interpolate t nx ny = sinCos' $ ax + t*adiff where + ax = angle2' nx + ay = angle2' ny + adiff = helper (ay - ax) + helper d + | d < -pi = d + twopi + | d > pi = d - twopi + | otherwise = d + twopi = 2*pi + +instance Interpolate Normal3 where + interpolate t nx ny = + if maxAngle < 0.001 -- more or less ad-hoc critical angle + then mkNormal $ interpolate t x y + else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where + x = fromNormal nx + y = fromNormal ny + axis = (x &^ y) + maxAngle = acos (x &. y) + + diff --git a/Data/Vect/Double/OpenGL.hs b/Data/Vect/Double/OpenGL.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/OpenGL.hs @@ -0,0 +1,183 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +-- 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'. + +module Data.Vect.Flt.OpenGL where + +import Control.Monad +import Data.Vect.Flt.Base +import qualified Graphics.Rendering.OpenGL as GL + +import Foreign + +import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale) + +------------------------------------------------------- + +{-# SPECIALISE radianToDegrees :: Float -> Float #-} +{-# SPECIALISE radianToDegrees :: Double -> Double #-} +radianToDegrees :: RealFrac a => a -> a +radianToDegrees x = x * 57.295779513082322 + +{-# SPECIALIZE degreesToRadian :: Float -> Float #-} +{-# SPECIALIZE degreesToRadian :: Double -> Double #-} +degreesToRadian :: Floating a => a -> a +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) + +translate :: Vec3 -> IO () +translate (Vec3 x y z) = GL.translate (Vector3 x y z) + +scale3 :: Vec3 -> IO () +scale3 (Vec3 x y z) = GL.scale x y z + +scale :: Flt -> IO () +scale x = GL.scale x x x + +------------------------------------------------------- + +-- Vertex instances + +instance GL.Vertex Vec2 where + vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec3 where + vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec4 where + 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) + normalv p = peek p >>= normal + +------------------------------------------------------- + +-- Color instances + +instance GL.Color Vec3 where + color (Vec3 r g b) = GL.color (GL.Color3 r g b) + colorv p = peek p >>= color + +instance GL.Color Vec4 where + color (Vec4 r g b a) = GL.color (GL.Color4 r g b a) + colorv p = peek p >>= color + +instance GL.SecondaryColor Vec3 where + secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b) + secondaryColorv p = peek p >>= secondaryColor + +{- +-- there is no such thing? +instance GL.SecondaryColor Vec4 where + secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a) + secondaryColorv p = peek p >>= secondaryColor +-} + +------------------------------------------------------- + +-- TexCoord instances + +instance GL.TexCoord Vec2 where + texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec3 where + texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec4 where + texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z) + texCoordv p = peek p >>= texCoord + 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 + vertexAttrib :: GL.AttribLocation -> a -> IO () + +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where + vertexAttrib loc x = GL.vertexAttrib1 loc x + +instance VertexAttrib' Vec2 where + vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y + +instance VertexAttrib' Vec3 where + vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Vec4 where + 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 + +instance VertexAttrib' Normal3 where + vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Normal4 where + vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w + +------------------------------------------------------- + +-- Uniform (again, experimental) + +-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names) + +{- +toFloat :: Flt -> Float +toFloat = realToFrac + +fromFloat :: Float -> Flt +fromFloat = realToFrac +-} + +-- Uniforms are always floats... +#ifdef VECT_Float + +instance GL.Uniform Flt where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc) + setter x = ($=) (uniform loc) (Index1 x) + uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt)) + +instance GL.Uniform Vec2 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc) + setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) + uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec3 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc) + setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) + uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec4 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc) + setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) + uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt) + +#endif diff --git a/Data/Vect/Double/Util/Dim2.hs b/Data/Vect/Double/Util/Dim2.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Util/Dim2.hs @@ -0,0 +1,65 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +module Data.Vect.Flt.Util.Dim2 where + +import Data.Vect.Flt.Base + +-- |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" + +destructVec2 :: [Vec2] -> [Flt] +destructVec2 [] = [] +destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls) + +det2 :: Vec2 -> Vec2 -> Flt +det2 u v = det (u,v) + +vec2X :: Vec2 +vec2Y :: Vec2 + +vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 + +translate2X :: Flt -> Vec2 -> Vec2 +translate2Y :: Flt -> Vec2 -> Vec2 + +translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) + +-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW). +-- A more precise name would be @cosSin@, but that sounds bad :) +sinCos :: Flt -> Vec2 +sinCos a = Vec2 (cos a) (sin a) + +sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2 +sinCos' = toNormalUnsafe . sinCos + +sinCosRadius :: Flt -- ^ angle (in radians) + -> Flt -- ^ radius + -> Vec2 +sinCosRadius a r = Vec2 (r * cos a) (r * sin a) + +-- | The angle relative to the positive X axis +angle2 :: Vec2 -> Flt +angle2 (Vec2 x y) = atan2 y x + +angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt +angle2' = angle2 . fromNormal + +-- |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 + +rotate2 :: Flt -> Vec2 -> Vec2 +rotate2 a v = v .* (rotMatrix2 a) + +-- |Rotates counterclockwise by 90 degrees. +rotateCCW :: Vec2 -> Vec2 +rotateCCW (Vec2 x y) = Vec2 (-y) x + +-- |Rotates clockwise by 90 degrees. +rotateCW :: Vec2 -> Vec2 +rotateCW (Vec2 x y) = Vec2 y (-x) diff --git a/Data/Vect/Double/Util/Dim3.hs b/Data/Vect/Double/Util/Dim3.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Util/Dim3.hs @@ -0,0 +1,76 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +module Data.Vect.Flt.Util.Dim3 where + +import Data.Vect.Flt.Base + +structVec3 :: [Flt] -> [Vec3] +structVec3 [] = [] +structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "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 + +translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) + +vec3X :: Vec3 +vec3Y :: Vec3 +vec3Z :: Vec3 + +vec3X = Vec3 1 0 0 +vec3Y = Vec3 0 1 0 +vec3Z = Vec3 0 0 1 + +rotMatrixZ :: Flt -> Mat3 +rotMatrixY :: Flt -> Mat3 +rotMatrixX :: Flt -> Mat3 + +-- These are intended for multiplication on the /right/. +-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@) +rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a +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' angle axis v = v .* (rotMatrix3' axis angle) + +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/. +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/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3 +rotMatrix3' (Normal3 v) a = + let c = cos a + s = sin a + m1 = scalarMul (1-c) (outer v v) + x = _1 v + y = _2 v + z = _3 v + m2 = Mat3 (Vec3 c ( s*z) (-s*y)) + (Vec3 (-s*z) c ( s*x)) + (Vec3 ( s*y) (-s*x) c ) + in (m1 &+ m2) + + diff --git a/Data/Vect/Double/Util/Dim4.hs b/Data/Vect/Double/Util/Dim4.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Util/Dim4.hs @@ -0,0 +1,93 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea. +-- Not very useful for most people, and not re-exported by "Data.Vect". + +module Data.Vect.Flt.Util.Dim4 where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.GramSchmidt + +structVec4 :: [Flt] -> [Vec4] +structVec4 [] = [] +structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4" + +destructVec4 :: [Vec4] -> [Flt] +destructVec4 [] = [] +destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls) + +--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt +--det4 u v w z = det (u,v,w,z) + +translate4X :: Flt -> Vec4 -> Vec4 +translate4Y :: Flt -> Vec4 -> Vec4 +translate4Z :: Flt -> Vec4 -> Vec4 +translate4W :: Flt -> Vec4 -> Vec4 + +translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w +translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) + +vec4X :: Vec4 +vec4Y :: Vec4 +vec4Z :: Vec4 +vec4W :: Vec4 + +vec4X = Vec4 1 0 0 0 +vec4Y = Vec4 0 1 0 0 +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))@. +-- This is a helper function for the 4 dimensional rotation code. +-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is +-- +-- > [ 0 r q p ] +-- > [ -r 0 z -y ] +-- > [ -q -z 0 x ] +-- > [ -p y -x 0 ] +biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt) +biVector4 (Vec4 x y z w) (Vec4 a b c d) = + ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c ) + +-- | the corresponding antisymmetric tensor +biVector4AsTensor :: Vec4 -> Vec4 -> Mat4 +biVector4AsTensor v w = + Mat4 ( Vec4 0 ( r) ( q) ( p) ) + ( Vec4 (-r) 0 ( z) (-y) ) + ( Vec4 (-q) (-z) 0 ( x) ) + ( Vec4 (-p) ( y) (-x) 0 ) + where + (x,y,z,p,q,r) = biVector4 v w + +-- | We assume that the axes are normalized and /orthogonal/ to each other! +rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4 +rotate4' angle axes v = v .* (rotMatrix4' angle axes) + +-- | We assume only that the axes are independent vectors. +rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4 +rotate4 angle axes v = v .* (rotMatrix4 angle axes) + +-- | 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 + where + c = cos angle ; s = sin angle + m1 = scalarMul (1-c) ( outer v v &+ outer w w ) + m2 = biVector4AsTensor v w + m3 = diag (Vec4 c c c c) + +-- | We assume only that the axes are independent vectors. +rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4 +rotMatrix4 angle axes = + rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes + where + liftPair f (x,y) = (f x, f y) + + diff --git a/Data/Vect/Double/Util/Projective.hs b/Data/Vect/Double/Util/Projective.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Double/Util/Projective.hs @@ -0,0 +1,88 @@ +{-# OPTIONS_GHC -DFlt=Double -DVECT_Double #-} + +-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on +-- the /right/, that is, they are of the form +-- +-- > _____ +-- > [ | | 0 ] +-- > [ | 3x3 | 0 ] +-- > [ |_____| 0 ] +-- > [ p q r 1 ] +-- +-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- 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 ()@. + +module Data.Vect.Flt.Util.Projective where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim3 + +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) + +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) + +rotMatrixProj :: Flt -> Normal3 -> Mat4 +rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle + +rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4 +rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle + +translMatrixProj :: Vec3 -> Mat4 +translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v) + +-- | 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 + +scaleMatrixUniformProj :: Flt -> Mat4 +scaleMatrixUniformProj s = diag (Vec4 s s s 1) + +class ProjectiveAction v where + actProj :: v -> Mat4 -> v + +instance ProjectiveAction Vec3 where + actProj v m = trim $ (extendProj v) .* m + +instance ProjectiveAction Vec4 where + actProj v m = v .* m + +-- | When acting on unit vectors, we ignore the translation part. +instance ProjectiveAction Normal3 where + actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3)) + +-- | 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) + where + t' = Vec3 (- u &. t) (- v &. t) (- w &. t) + diff --git a/Data/Vect/Float.hs b/Data/Vect/Float.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float.hs @@ -0,0 +1,23 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +module Data.Vect.Flt + ( module Data.Vect.Flt.Base + , module Data.Vect.Flt.Interpolate + , module Data.Vect.Flt.Util.Dim2 + , module Data.Vect.Flt.Util.Dim3 + , module Data.Vect.Flt.Util.Projective +#ifdef VECT_OPENGL + , module Data.Vect.Flt.OpenGL +#endif + ) where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Interpolate + +import Data.Vect.Flt.Util.Dim2 +import Data.Vect.Flt.Util.Dim3 +import Data.Vect.Flt.Util.Projective + +#ifdef VECT_OPENGL +import Data.Vect.Flt.OpenGL +#endif diff --git a/Data/Vect/Float/Base.hs b/Data/Vect/Float/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Base.hs @@ -0,0 +1,735 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +module Data.Vect.Flt.Base where + +import Control.Monad +import System.Random +import Foreign + +-- class declarations + +class AbelianGroup g where + (&+) :: g -> g -> g + (&-) :: g -> g -> g + neg :: g -> g + zero :: g + +infixl 6 &+ +infixl 6 &- + +vecSum :: AbelianGroup g => [g] -> g +vecSum l = foldl (&+) zero l + +class (AbelianGroup r) => + Ring r where + (.*.) :: r -> r -> r + one :: r + +infixl 7 .*. + +ringProduct :: Ring r => [r] -> r +ringProduct l = foldl (.*.) one l + +class LeftModule r m where + lmul :: r -> m -> m + (*.) :: r -> m -> m + (*.) = lmul + +class RightModule m r where + rmul :: m -> r -> m + (.*) :: m -> r -> m + (.*) = rmul + +-- I'm not really sure about this.. may actually degrade the performance in some cases? +{- RULES +"matrix multiplication left" forall m n x. (n .*. m) *. x = n *. (m *. x) +"matrix multiplication right" forall m n x. x .* (m .*. n) = (x .* m) .* n + -} + +infixr 7 *. +infixl 7 .* + +class AbelianGroup v => Vector v where + mapVec :: (Flt -> Flt) -> v -> v + scalarMul :: Flt -> v -> v + (*&) :: Flt -> v -> v + (&*) :: v -> Flt -> v + (*&) s v = scalarMul s v + (&*) v s = scalarMul s v + +infixr 7 *& +infixl 7 &* + +{-# RULES +"scalar multiplication left" forall s t x. t *& (s *& x) = (t*s) *& x +"scalar multiplication right" forall s t x. (x &* s) &* t = x &* (s*t) + #-} + +class DotProd v where + (&.) :: v -> v -> Flt + norm :: v -> Flt + normsqr :: v -> Flt + len :: v -> Flt + lensqr :: v -> Flt + len = norm + lensqr = normsqr + dotprod :: v -> v -> Flt + normsqr v = (v &. v) + norm = sqrt.lensqr + dotprod = (&.) + +infix 7 &. + +{-# RULES +"len/square 1" forall x. (len x)*(len x) = lensqr x +"len/square 2" forall x. (len x)^2 = lensqr x +"norm/square 1" forall x. (norm x)*(norm x) = normsqr x +"norm/square 2" forall x. (norm x)^2 = normsqr x + #-} + +normalize :: (Vector v, DotProd v) => v -> v +normalize v = scalarMul (1.0/(len v)) v + +distance :: (Vector v, DotProd v) => v -> v -> Flt +distance x y = norm (x &- y) + +-- | the angle between two vectors +angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y) + +-- | the angle between two unit vectors +angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y) + +{-# RULES +"normalize is idempotent" forall x. normalize (normalize x) = normalize x + #-} + +class (Vector v, DotProd v) => UnitVector v u | v->u, u->v where + mkNormal :: v -> u -- ^ normalizes the input + toNormalUnsafe :: v -> u -- ^ does not normalize the input! + fromNormal :: u -> v + fromNormalRadius :: Flt -> u -> v + fromNormalRadius t n = t *& fromNormal n + +-- | 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! +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. +flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal + +class CrossProd v where + crossprod :: v -> v -> v + (&^) :: v -> v -> v + (&^) = crossprod + +class Pointwise v where + pointwise :: v -> v -> v + (&!) :: v -> v -> v + (&!) = pointwise + +infix 7 &^ +infix 7 &! + +class HasCoordinates v x | v->x where + _1 :: v -> x + _2 :: v -> x + _3 :: v -> x + _4 :: v -> x + +-- | conversion between vectors (and matrices) of different dimensions +class Extend u v where + extendZero :: u -> v -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@ + extendWith :: Flt -> u -> v -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@ + trim :: v -> u -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@ + +-- | makes a diagonal matrix from a vector +class Diagonal s t | t->s where + diag :: s -> t + +class Matrix m where + transpose :: m -> m + inverse :: m -> m + idmtx :: m + +{-# RULES +"transpose is an involution" forall m. transpose (transpose m) = m +"inverse is an involution" forall m. inverse (inverse m) = m + #-} + +-- | 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 + +-- Vec / Mat datatypes + +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) + +-- | these 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) + +mkVec2 :: (Flt,Flt) -> Vec2 +mkVec3 :: (Flt,Flt,Flt) -> Vec3 +mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4 + +mkVec2 (x,y) = Vec2 x y +mkVec3 (x,y,z) = Vec3 x y z +mkVec4 (x,y,z,w) = Vec4 x y z w + +-- Unit vectors + +instance UnitVector Vec2 Normal2 where + mkNormal v = Normal2 (normalize v) + fromNormal (Normal2 v) = v + toNormalUnsafe = Normal2 + +instance UnitVector Vec3 Normal3 where + mkNormal v = Normal3 (normalize v) + fromNormal (Normal3 v) = v + toNormalUnsafe = Normal3 + +instance UnitVector Vec4 Normal4 where + mkNormal v = Normal4 (normalize v) + fromNormal (Normal4 v) = v + toNormalUnsafe = Normal4 + +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 + where + (v,h) = random g + d = norm v + +instance Random Normal2 where + 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) + randomR _ = random + +instance Random Normal4 where + 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 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 +-} + +-- Vec2 instances + +instance HasCoordinates Vec2 Flt where + _1 (Vec2 x _) = x + _2 (Vec2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance AbelianGroup Vec2 where + (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) + (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2) + neg (Vec2 x y) = Vec2 (-x) (-y) + zero = Vec2 0 0 + +instance Vector Vec2 where + scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y) + mapVec f (Vec2 x y) = Vec2 (f x) (f y) + +instance DotProd Vec2 where + (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2 + +instance Pointwise Vec2 where + pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2) + +instance Determinant (Vec2,Vec2) where + det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1 + +{- +instance Show Vec2 where + show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )" +-} + +instance Random Vec2 where + random = randomR (Vec2 (-1) (-1),Vec2 1 1) + randomR (Vec2 a b, Vec2 c d) gen = + let (x,gen1) = randomR (a,c) gen + (y,gen2) = randomR (b,d) gen1 + in (Vec2 x y, gen2) + +instance Storable Vec2 where + sizeOf _ = 2 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p k + return (Vec2 x y) + + poke q (Vec2 x y) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p k y + +-- Mat2 instances + +instance HasCoordinates Mat2 Vec2 where + _1 (Mat2 x _) = x + _2 (Mat2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance Matrix Mat2 where + transpose (Mat2 row1 row2) = + Mat2 (Vec2 (_1 row1) (_1 row2)) + (Vec2 (_2 row1) (_2 row2)) + idmtx = Mat2 (Vec2 1 0) (Vec2 0 1) + inverse (Mat2 (Vec2 a b) (Vec2 c d)) = + Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) + where r = 1.0 / (a*d - b*c) + +instance AbelianGroup Mat2 where + (&+) (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) + +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 + (.*.) (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 LeftModule Mat2 Vec2 where + lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) + +instance RightModule Vec2 Mat2 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec2 Mat2 where + diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y) + +instance Tensor Mat2 Vec2 where + 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 + +{- +instance Show Mat2 where + show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2 +-} + +instance Storable Mat2 where + sizeOf _ = 2 * sizeOf (undefined::Vec2) + alignment _ = alignment (undefined::Vec2) + + peek q = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + r1 <- peek p + r2 <- peekByteOff p k + return (Mat2 r1 r2) + + poke q (Mat2 r1 r2) = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + poke p r1 + pokeByteOff p k r2 + +-- Vec3 instances + +instance HasCoordinates Vec3 Flt where + _1 (Vec3 x _ _) = x + _2 (Vec3 _ y _) = y + _3 (Vec3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance AbelianGroup Vec3 where + (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) + (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) + neg (Vec3 x y z) = Vec3 (-x) (-y) (-z) + zero = Vec3 0 0 0 + +instance Vector Vec3 where + scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z) + mapVec f (Vec3 x y z) = Vec3 (f x) (f y) (f z) + +instance DotProd Vec3 where + (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2 + +instance Pointwise Vec3 where + pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2) + +{- +instance Show Vec3 where + show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )" +-} + +instance Random Vec3 where + random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1) + randomR (Vec3 a b c, Vec3 d e f) gen = + let (x,gen1) = randomR (a,d) gen + (y,gen2) = randomR (b,e) gen1 + (z,gen3) = randomR (c,f) gen2 + in (Vec3 x y z, gen3) + +instance CrossProd Vec3 where + crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) + +instance Determinant (Vec3,Vec3,Vec3) where + det (u,v,w) = u &. (v &^ w) + +instance Storable Vec3 where + sizeOf _ = 3 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + return (Vec3 x y z) + + poke q (Vec3 x y z) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + +-- Mat3 instances + +instance HasCoordinates Mat3 Vec3 where + _1 (Mat3 x _ _) = x + _2 (Mat3 _ y _) = y + _3 (Mat3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance Matrix Mat3 where + + transpose (Mat3 row1 row2 row3) = + Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) + (Vec3 (_2 row1) (_2 row2) (_2 row3)) + (Vec3 (_3 row1) (_3 row2) (_3 row3)) + + idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1) + + inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = + Mat3 (Vec3 (d11*r) (d21*r) (d31*r)) + (Vec3 (d12*r) (d22*r) (d32*r)) + (Vec3 (d13*r) (d23*r) (d33*r)) + where + r = 1.0 / ( a*d11 + b*d12 + c*d13 ) + + d11 = f*k - g*j + d12 = g*i - e*k + d13 = e*j - f*i + + d31 = b*g - c*f + d32 = c*e - a*g + d33 = a*f - b*e + + d21 = c*j - b*k + d22 = a*k - c*i + d23 = b*i - a*j + +instance AbelianGroup Mat3 where + (&+) (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) + +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 + (.*.) (Mat3 r1 r2 r3) n = + let (Mat3 c1 c2 c3) = transpose n + in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3)) + (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3)) + (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3)) + one = idmtx + +instance LeftModule Mat3 Vec3 where + lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v) + +instance RightModule Vec3 Mat3 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec3 Mat3 where + diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z) + +instance Tensor Mat3 Vec3 where + outer (Vec3 a b c) (Vec3 x y z) = Mat3 + (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) + +{- +instance Show Mat3 where + show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 +-} + +instance Storable Mat3 where + sizeOf _ = 3 * sizeOf (undefined::Vec3) + alignment _ = alignment (undefined::Vec3) + + peek q = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + return (Mat3 r1 r2 r3) + + poke q (Mat3 r1 r2 r3) = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + +-- Vec4 instances + +instance HasCoordinates Vec4 Flt where + _1 (Vec4 x _ _ _) = x + _2 (Vec4 _ y _ _) = y + _3 (Vec4 _ _ z _) = z + _4 (Vec4 _ _ _ w) = w + +instance AbelianGroup Vec4 where + (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2) + (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2) + neg (Vec4 x y z w) = Vec4 (-x) (-y) (-z) (-w) + zero = Vec4 0 0 0 0 + +instance Vector Vec4 where + scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w) + mapVec f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w) + +instance DotProd Vec4 where + (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2 + +instance Pointwise Vec4 where + pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2) + +{- +instance Show Vec4 where + show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )" +-} + +instance Random Vec4 where + random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1) + randomR (Vec4 a b c d, Vec4 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 (Vec4 x y z w, gen4) + +instance Storable Vec4 where + sizeOf _ = 4 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + w <- peekByteOff p (3*k) + return (Vec4 x y z w) + + poke q (Vec4 x y z w) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + pokeByteOff p (3*k) w + +-- Mat4 instances + +instance HasCoordinates Mat4 Vec4 where + _1 (Mat4 x _ _ _) = x + _2 (Mat4 _ y _ _) = y + _3 (Mat4 _ _ z _) = z + _4 (Mat4 _ _ _ w) = w + +instance Matrix Mat4 where + transpose (Mat4 row1 row2 row3 row4) = + Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) + (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) + (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) + (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) + idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1) + inverse = error "inverse/Mat4: not implemented yet" + +instance AbelianGroup Mat4 where + (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4) + (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4) + neg (Mat4 r1 r2 r3 r4) = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) + zero = Mat4 zero zero zero zero + +instance Vector Mat4 where + 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 + (.*.) (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)) + (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4)) + (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4)) + (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4)) + one = idmtx + +instance LeftModule Mat4 Vec4 where + lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v) + +instance RightModule Vec4 Mat4 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec4 Mat4 where + diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w) + +instance Tensor Mat4 Vec4 where + outer (Vec4 a b c d) (Vec4 x y z w) = Mat4 + (Vec4 (a*x) (a*y) (a*z) (a*w)) + (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 Show Mat4 where + show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4 +-} + +instance Storable Mat4 where + sizeOf _ = 4 * sizeOf (undefined::Vec4) + alignment _ = alignment (undefined::Vec4) + + peek q = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + r4 <- peekByteOff p (3*k) + return (Mat4 r1 r2 r3 r4) + + poke q (Mat4 r1 r2 r3 r4) = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + pokeByteOff p (3*k) r4 + +-- Extend instances + +instance Extend Vec2 Vec3 where + extendZero (Vec2 x y) = Vec3 x y 0 + extendWith t (Vec2 x y) = Vec3 x y t + trim (Vec3 x y _) = Vec2 x y + +instance Extend Vec2 Vec4 where + extendZero (Vec2 x y) = Vec4 x y 0 0 + extendWith t (Vec2 x y) = Vec4 x y t t + trim (Vec4 x y _ _) = Vec2 x y + +instance Extend Vec3 Vec4 where + extendZero (Vec3 x y z) = Vec4 x y z 0 + extendWith t (Vec3 x y z) = Vec4 x y z t + 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" + 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" + 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" + trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r) + diff --git a/Data/Vect/Float/GramSchmidt.hs b/Data/Vect/Float/GramSchmidt.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/GramSchmidt.hs @@ -0,0 +1,135 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +-- | Gram-Schmidt orthogonalization. +-- This module is not re-exported by "Data.Vect". + +module Data.Vect.Flt.GramSchmidt + ( GramSchmidt(..) + ) + where + +import Data.Vect.Flt.Base + +------------------------------------------------------- + +liftPair :: (a -> b) -> (a,a) -> (b,b) +liftPair f (x,y) = (f x, f y) + +liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b) +liftTriple f (x,y,z) = (f x, f y, f z) + +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 + gramSchmidt :: a -> a -- ^ does not normalize the vectors! + gramSchmidtNormalize :: a -> a -- ^ normalizes the vectors. + +{-# RULES +"gramSchmidt is idempotent" forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent" forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a + #-} + +------------------------------------------------------- + +instance GramSchmidt (Vec2,Vec2) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec3,Vec3) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec4,Vec4) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +---------- + +instance GramSchmidt (Normal2,Normal2) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal2!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +---------- + +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtPair (u,v) = (u',v') where + u' = u + v' = project v u' + +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtNormalizePair (u,v) = (u',v') where + u' = normalize u + v' = normalize $ projectUnsafe v u' + +---------- + +instance GramSchmidt (Vec3,Vec3,Vec3) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Normal3,Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +instance GramSchmidt (Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +---------- + +gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtTriple (u,v,w) = (u',v',w') where + u' = u + v' = project v u' + w' = project (project w u') v' + +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + +---------- + +instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtQuadruple + gramSchmidtNormalize = gramSchmidtNormalizeQuadruple + +instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal + +---------- + +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where + u' = u + v' = project v u' + w' = project (project w u') v' + z' = project (project (project z u') v') w' + +gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w' + +---------- + diff --git a/Data/Vect/Float/Interpolate.hs b/Data/Vect/Float/Interpolate.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Interpolate.hs @@ -0,0 +1,43 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors. + +module Data.Vect.Flt.Interpolate where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim2 (sinCos',angle2') +import Data.Vect.Flt.Util.Dim3 (rotate3') + +class Interpolate v where + interpolate :: Flt -> v -> v -> v + +instance Interpolate Flt where + interpolate t x y = x + t*(y-x) + +instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x) + +instance Interpolate Normal2 where + interpolate t nx ny = sinCos' $ ax + t*adiff where + ax = angle2' nx + ay = angle2' ny + adiff = helper (ay - ax) + helper d + | d < -pi = d + twopi + | d > pi = d - twopi + | otherwise = d + twopi = 2*pi + +instance Interpolate Normal3 where + interpolate t nx ny = + if maxAngle < 0.001 -- more or less ad-hoc critical angle + then mkNormal $ interpolate t x y + else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where + x = fromNormal nx + y = fromNormal ny + axis = (x &^ y) + maxAngle = acos (x &. y) + + diff --git a/Data/Vect/Float/OpenGL.hs b/Data/Vect/Float/OpenGL.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/OpenGL.hs @@ -0,0 +1,183 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +-- 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'. + +module Data.Vect.Flt.OpenGL where + +import Control.Monad +import Data.Vect.Flt.Base +import qualified Graphics.Rendering.OpenGL as GL + +import Foreign + +import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale) + +------------------------------------------------------- + +{-# SPECIALISE radianToDegrees :: Float -> Float #-} +{-# SPECIALISE radianToDegrees :: Double -> Double #-} +radianToDegrees :: RealFrac a => a -> a +radianToDegrees x = x * 57.295779513082322 + +{-# SPECIALIZE degreesToRadian :: Float -> Float #-} +{-# SPECIALIZE degreesToRadian :: Double -> Double #-} +degreesToRadian :: Floating a => a -> a +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) + +translate :: Vec3 -> IO () +translate (Vec3 x y z) = GL.translate (Vector3 x y z) + +scale3 :: Vec3 -> IO () +scale3 (Vec3 x y z) = GL.scale x y z + +scale :: Flt -> IO () +scale x = GL.scale x x x + +------------------------------------------------------- + +-- Vertex instances + +instance GL.Vertex Vec2 where + vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec3 where + vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec4 where + 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) + normalv p = peek p >>= normal + +------------------------------------------------------- + +-- Color instances + +instance GL.Color Vec3 where + color (Vec3 r g b) = GL.color (GL.Color3 r g b) + colorv p = peek p >>= color + +instance GL.Color Vec4 where + color (Vec4 r g b a) = GL.color (GL.Color4 r g b a) + colorv p = peek p >>= color + +instance GL.SecondaryColor Vec3 where + secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b) + secondaryColorv p = peek p >>= secondaryColor + +{- +-- there is no such thing? +instance GL.SecondaryColor Vec4 where + secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a) + secondaryColorv p = peek p >>= secondaryColor +-} + +------------------------------------------------------- + +-- TexCoord instances + +instance GL.TexCoord Vec2 where + texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec3 where + texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec4 where + texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z) + texCoordv p = peek p >>= texCoord + 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 + vertexAttrib :: GL.AttribLocation -> a -> IO () + +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where + vertexAttrib loc x = GL.vertexAttrib1 loc x + +instance VertexAttrib' Vec2 where + vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y + +instance VertexAttrib' Vec3 where + vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Vec4 where + 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 + +instance VertexAttrib' Normal3 where + vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Normal4 where + vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w + +------------------------------------------------------- + +-- Uniform (again, experimental) + +-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names) + +{- +toFloat :: Flt -> Float +toFloat = realToFrac + +fromFloat :: Float -> Flt +fromFloat = realToFrac +-} + +-- Uniforms are always floats... +#ifdef VECT_Float + +instance GL.Uniform Flt where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc) + setter x = ($=) (uniform loc) (Index1 x) + uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt)) + +instance GL.Uniform Vec2 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc) + setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) + uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec3 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc) + setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) + uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec4 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc) + setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) + uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt) + +#endif diff --git a/Data/Vect/Float/Util/Dim2.hs b/Data/Vect/Float/Util/Dim2.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Util/Dim2.hs @@ -0,0 +1,65 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +module Data.Vect.Flt.Util.Dim2 where + +import Data.Vect.Flt.Base + +-- |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" + +destructVec2 :: [Vec2] -> [Flt] +destructVec2 [] = [] +destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls) + +det2 :: Vec2 -> Vec2 -> Flt +det2 u v = det (u,v) + +vec2X :: Vec2 +vec2Y :: Vec2 + +vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 + +translate2X :: Flt -> Vec2 -> Vec2 +translate2Y :: Flt -> Vec2 -> Vec2 + +translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) + +-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW). +-- A more precise name would be @cosSin@, but that sounds bad :) +sinCos :: Flt -> Vec2 +sinCos a = Vec2 (cos a) (sin a) + +sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2 +sinCos' = toNormalUnsafe . sinCos + +sinCosRadius :: Flt -- ^ angle (in radians) + -> Flt -- ^ radius + -> Vec2 +sinCosRadius a r = Vec2 (r * cos a) (r * sin a) + +-- | The angle relative to the positive X axis +angle2 :: Vec2 -> Flt +angle2 (Vec2 x y) = atan2 y x + +angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt +angle2' = angle2 . fromNormal + +-- |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 + +rotate2 :: Flt -> Vec2 -> Vec2 +rotate2 a v = v .* (rotMatrix2 a) + +-- |Rotates counterclockwise by 90 degrees. +rotateCCW :: Vec2 -> Vec2 +rotateCCW (Vec2 x y) = Vec2 (-y) x + +-- |Rotates clockwise by 90 degrees. +rotateCW :: Vec2 -> Vec2 +rotateCW (Vec2 x y) = Vec2 y (-x) diff --git a/Data/Vect/Float/Util/Dim3.hs b/Data/Vect/Float/Util/Dim3.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Util/Dim3.hs @@ -0,0 +1,76 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +module Data.Vect.Flt.Util.Dim3 where + +import Data.Vect.Flt.Base + +structVec3 :: [Flt] -> [Vec3] +structVec3 [] = [] +structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "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 + +translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) + +vec3X :: Vec3 +vec3Y :: Vec3 +vec3Z :: Vec3 + +vec3X = Vec3 1 0 0 +vec3Y = Vec3 0 1 0 +vec3Z = Vec3 0 0 1 + +rotMatrixZ :: Flt -> Mat3 +rotMatrixY :: Flt -> Mat3 +rotMatrixX :: Flt -> Mat3 + +-- These are intended for multiplication on the /right/. +-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@) +rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a +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' angle axis v = v .* (rotMatrix3' axis angle) + +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/. +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/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3 +rotMatrix3' (Normal3 v) a = + let c = cos a + s = sin a + m1 = scalarMul (1-c) (outer v v) + x = _1 v + y = _2 v + z = _3 v + m2 = Mat3 (Vec3 c ( s*z) (-s*y)) + (Vec3 (-s*z) c ( s*x)) + (Vec3 ( s*y) (-s*x) c ) + in (m1 &+ m2) + + diff --git a/Data/Vect/Float/Util/Dim4.hs b/Data/Vect/Float/Util/Dim4.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Util/Dim4.hs @@ -0,0 +1,93 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea. +-- Not very useful for most people, and not re-exported by "Data.Vect". + +module Data.Vect.Flt.Util.Dim4 where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.GramSchmidt + +structVec4 :: [Flt] -> [Vec4] +structVec4 [] = [] +structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4" + +destructVec4 :: [Vec4] -> [Flt] +destructVec4 [] = [] +destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls) + +--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt +--det4 u v w z = det (u,v,w,z) + +translate4X :: Flt -> Vec4 -> Vec4 +translate4Y :: Flt -> Vec4 -> Vec4 +translate4Z :: Flt -> Vec4 -> Vec4 +translate4W :: Flt -> Vec4 -> Vec4 + +translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w +translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) + +vec4X :: Vec4 +vec4Y :: Vec4 +vec4Z :: Vec4 +vec4W :: Vec4 + +vec4X = Vec4 1 0 0 0 +vec4Y = Vec4 0 1 0 0 +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))@. +-- This is a helper function for the 4 dimensional rotation code. +-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is +-- +-- > [ 0 r q p ] +-- > [ -r 0 z -y ] +-- > [ -q -z 0 x ] +-- > [ -p y -x 0 ] +biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt) +biVector4 (Vec4 x y z w) (Vec4 a b c d) = + ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c ) + +-- | the corresponding antisymmetric tensor +biVector4AsTensor :: Vec4 -> Vec4 -> Mat4 +biVector4AsTensor v w = + Mat4 ( Vec4 0 ( r) ( q) ( p) ) + ( Vec4 (-r) 0 ( z) (-y) ) + ( Vec4 (-q) (-z) 0 ( x) ) + ( Vec4 (-p) ( y) (-x) 0 ) + where + (x,y,z,p,q,r) = biVector4 v w + +-- | We assume that the axes are normalized and /orthogonal/ to each other! +rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4 +rotate4' angle axes v = v .* (rotMatrix4' angle axes) + +-- | We assume only that the axes are independent vectors. +rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4 +rotate4 angle axes v = v .* (rotMatrix4 angle axes) + +-- | 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 + where + c = cos angle ; s = sin angle + m1 = scalarMul (1-c) ( outer v v &+ outer w w ) + m2 = biVector4AsTensor v w + m3 = diag (Vec4 c c c c) + +-- | We assume only that the axes are independent vectors. +rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4 +rotMatrix4 angle axes = + rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes + where + liftPair f (x,y) = (f x, f y) + + diff --git a/Data/Vect/Float/Util/Projective.hs b/Data/Vect/Float/Util/Projective.hs new file mode 100644 --- /dev/null +++ b/Data/Vect/Float/Util/Projective.hs @@ -0,0 +1,88 @@ +{-# OPTIONS_GHC -DFlt=Float -DVECT_Float #-} + +-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on +-- the /right/, that is, they are of the form +-- +-- > _____ +-- > [ | | 0 ] +-- > [ | 3x3 | 0 ] +-- > [ |_____| 0 ] +-- > [ p q r 1 ] +-- +-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- 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 ()@. + +module Data.Vect.Flt.Util.Projective where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim3 + +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) + +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) + +rotMatrixProj :: Flt -> Normal3 -> Mat4 +rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle + +rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4 +rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle + +translMatrixProj :: Vec3 -> Mat4 +translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v) + +-- | 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 + +scaleMatrixUniformProj :: Flt -> Mat4 +scaleMatrixUniformProj s = diag (Vec4 s s s 1) + +class ProjectiveAction v where + actProj :: v -> Mat4 -> v + +instance ProjectiveAction Vec3 where + actProj v m = trim $ (extendProj v) .* m + +instance ProjectiveAction Vec4 where + actProj v m = v .* m + +-- | When acting on unit vectors, we ignore the translation part. +instance ProjectiveAction Normal3 where + actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3)) + +-- | 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) + where + t' = Vec3 (- u &. t) (- v &. t) (- w &. t) + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,29 @@ +Copyright (c) 2008, Balazs Komuves +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +- Redistributions of source code must retain the above copyright notice, +this list of conditions and the following disclaimer. + +- Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation +and/or other materials provided with the distribution. + +- Neither names of the copyright holders nor the names of the contributors +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 +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 +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,44 @@ +#! /usr/bin/env runhaskell +> +> import Control.Monad +> import Distribution.Simple +> import Distribution.PackageDescription +> import System.IO +> import System.Directory +> +> copyFileWithPrefix src tgt prefix = +> readFile src >>= \txt -> writeFile tgt (prefix ++ txt) +> +> copyFiles srcdir tgtdir prefix = do +> files <- getDirectoryContents srcdir +> forM_ files $ \fname -> do +> let src = srcdir ++ fname +> tgt = tgtdir ++ fname +> doesFileExist src >>= \b -> when b $ copyFileWithPrefix src tgt prefix +> doesDirectoryExist src>>= \b -> when ( b && fname /= "." && fname /= ".." ) $ do +> createDirectoryIfMissing False tgt +> copyFiles (src ++ "/") (tgt ++ "/") prefix +> +> thePrefix flt = "{-# OPTIONS_GHC -DFlt=" ++ flt ++ " -DVECT_" ++ flt ++ " #-}\n" +> +> myPreBuildHook args buildflags = do +> createDirectoryIfMissing False "Data/Vect/Float" +> createDirectoryIfMissing False "Data/Vect/Double" +> copyFileWithPrefix "src/flt.hs" "Data/Vect/Float.hs" (thePrefix "Float") +> copyFileWithPrefix "src/flt.hs" "Data/Vect/Double.hs" (thePrefix "Double") +> copyFiles "src/flt/" "Data/Vect/Float/" (thePrefix "Float") +> copyFiles "src/flt/" "Data/Vect/Double/" (thePrefix "Double") +> return $ emptyHookedBuildInfo +> +> myPostCleanHook args cleanflags pdep mlocalbuildinfo = do +> removeDirectoryRecursive "Data/Vect/Float" +> removeDirectoryRecursive "Data/Vect/Double" +> +> myUserHooks = simpleUserHooks +> { preBuild = myPreBuildHook +> , postClean = myPostCleanHook +> } +> +> main = do +> defaultMainWithHooks myUserHooks +> diff --git a/src/flt.hs b/src/flt.hs new file mode 100644 --- /dev/null +++ b/src/flt.hs @@ -0,0 +1,22 @@ + +module Data.Vect.Flt + ( module Data.Vect.Flt.Base + , module Data.Vect.Flt.Interpolate + , module Data.Vect.Flt.Util.Dim2 + , module Data.Vect.Flt.Util.Dim3 + , module Data.Vect.Flt.Util.Projective +#ifdef VECT_OPENGL + , module Data.Vect.Flt.OpenGL +#endif + ) where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Interpolate + +import Data.Vect.Flt.Util.Dim2 +import Data.Vect.Flt.Util.Dim3 +import Data.Vect.Flt.Util.Projective + +#ifdef VECT_OPENGL +import Data.Vect.Flt.OpenGL +#endif diff --git a/src/flt/Base.hs b/src/flt/Base.hs new file mode 100644 --- /dev/null +++ b/src/flt/Base.hs @@ -0,0 +1,734 @@ + +module Data.Vect.Flt.Base where + +import Control.Monad +import System.Random +import Foreign + +-- class declarations + +class AbelianGroup g where + (&+) :: g -> g -> g + (&-) :: g -> g -> g + neg :: g -> g + zero :: g + +infixl 6 &+ +infixl 6 &- + +vecSum :: AbelianGroup g => [g] -> g +vecSum l = foldl (&+) zero l + +class (AbelianGroup r) => + Ring r where + (.*.) :: r -> r -> r + one :: r + +infixl 7 .*. + +ringProduct :: Ring r => [r] -> r +ringProduct l = foldl (.*.) one l + +class LeftModule r m where + lmul :: r -> m -> m + (*.) :: r -> m -> m + (*.) = lmul + +class RightModule m r where + rmul :: m -> r -> m + (.*) :: m -> r -> m + (.*) = rmul + +-- I'm not really sure about this.. may actually degrade the performance in some cases? +{- RULES +"matrix multiplication left" forall m n x. (n .*. m) *. x = n *. (m *. x) +"matrix multiplication right" forall m n x. x .* (m .*. n) = (x .* m) .* n + -} + +infixr 7 *. +infixl 7 .* + +class AbelianGroup v => Vector v where + mapVec :: (Flt -> Flt) -> v -> v + scalarMul :: Flt -> v -> v + (*&) :: Flt -> v -> v + (&*) :: v -> Flt -> v + (*&) s v = scalarMul s v + (&*) v s = scalarMul s v + +infixr 7 *& +infixl 7 &* + +{-# RULES +"scalar multiplication left" forall s t x. t *& (s *& x) = (t*s) *& x +"scalar multiplication right" forall s t x. (x &* s) &* t = x &* (s*t) + #-} + +class DotProd v where + (&.) :: v -> v -> Flt + norm :: v -> Flt + normsqr :: v -> Flt + len :: v -> Flt + lensqr :: v -> Flt + len = norm + lensqr = normsqr + dotprod :: v -> v -> Flt + normsqr v = (v &. v) + norm = sqrt.lensqr + dotprod = (&.) + +infix 7 &. + +{-# RULES +"len/square 1" forall x. (len x)*(len x) = lensqr x +"len/square 2" forall x. (len x)^2 = lensqr x +"norm/square 1" forall x. (norm x)*(norm x) = normsqr x +"norm/square 2" forall x. (norm x)^2 = normsqr x + #-} + +normalize :: (Vector v, DotProd v) => v -> v +normalize v = scalarMul (1.0/(len v)) v + +distance :: (Vector v, DotProd v) => v -> v -> Flt +distance x y = norm (x &- y) + +-- | the angle between two vectors +angle :: (Vector v, DotProd v) => v -> v -> Flt +angle x y = acos $ (x &. y) / (norm x * norm y) + +-- | the angle between two unit vectors +angle' {- ' CPP is sensitive to primes -} :: (Vector v, UnitVector v u, DotProd v) => u -> u -> Flt +angle' x y = acos (fromNormal x &. fromNormal y) + +{-# RULES +"normalize is idempotent" forall x. normalize (normalize x) = normalize x + #-} + +class (Vector v, DotProd v) => UnitVector v u | v->u, u->v where + mkNormal :: v -> u -- ^ normalizes the input + toNormalUnsafe :: v -> u -- ^ does not normalize the input! + fromNormal :: u -> v + fromNormalRadius :: Flt -> u -> v + fromNormalRadius t n = t *& fromNormal n + +-- | 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! +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. +flipNormal :: UnitVector v n => n -> n +flipNormal = toNormalUnsafe . neg . fromNormal + +class CrossProd v where + crossprod :: v -> v -> v + (&^) :: v -> v -> v + (&^) = crossprod + +class Pointwise v where + pointwise :: v -> v -> v + (&!) :: v -> v -> v + (&!) = pointwise + +infix 7 &^ +infix 7 &! + +class HasCoordinates v x | v->x where + _1 :: v -> x + _2 :: v -> x + _3 :: v -> x + _4 :: v -> x + +-- | conversion between vectors (and matrices) of different dimensions +class Extend u v where + extendZero :: u -> v -- ^ example: @extendZero (Vec2 5 6) = Vec4 5 6 0 0@ + extendWith :: Flt -> u -> v -- ^ example: @extendWith 1 (Vec2 5 6) = Vec4 5 6 1 1@ + trim :: v -> u -- ^ example: @trim (Vec4 5 6 7 8) = Vec2 5 6@ + +-- | makes a diagonal matrix from a vector +class Diagonal s t | t->s where + diag :: s -> t + +class Matrix m where + transpose :: m -> m + inverse :: m -> m + idmtx :: m + +{-# RULES +"transpose is an involution" forall m. transpose (transpose m) = m +"inverse is an involution" forall m. inverse (inverse m) = m + #-} + +-- | 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 + +-- Vec / Mat datatypes + +data Vec2 = Vec2 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec3 = Vec3 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) +data Vec4 = Vec4 {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt {-# UNPACK #-} !Flt + deriving (Read,Show) + +-- | these 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) + +mkVec2 :: (Flt,Flt) -> Vec2 +mkVec3 :: (Flt,Flt,Flt) -> Vec3 +mkVec4 :: (Flt,Flt,Flt,Flt) -> Vec4 + +mkVec2 (x,y) = Vec2 x y +mkVec3 (x,y,z) = Vec3 x y z +mkVec4 (x,y,z,w) = Vec4 x y z w + +-- Unit vectors + +instance UnitVector Vec2 Normal2 where + mkNormal v = Normal2 (normalize v) + fromNormal (Normal2 v) = v + toNormalUnsafe = Normal2 + +instance UnitVector Vec3 Normal3 where + mkNormal v = Normal3 (normalize v) + fromNormal (Normal3 v) = v + toNormalUnsafe = Normal3 + +instance UnitVector Vec4 Normal4 where + mkNormal v = Normal4 (normalize v) + fromNormal (Normal4 v) = v + toNormalUnsafe = Normal4 + +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 + where + (v,h) = random g + d = norm v + +instance Random Normal2 where + 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) + randomR _ = random + +instance Random Normal4 where + 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 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 +-} + +-- Vec2 instances + +instance HasCoordinates Vec2 Flt where + _1 (Vec2 x _) = x + _2 (Vec2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance AbelianGroup Vec2 where + (&+) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1+x2) (y1+y2) + (&-) (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1-x2) (y1-y2) + neg (Vec2 x y) = Vec2 (-x) (-y) + zero = Vec2 0 0 + +instance Vector Vec2 where + scalarMul s (Vec2 x y) = Vec2 (s*x) (s*y) + mapVec f (Vec2 x y) = Vec2 (f x) (f y) + +instance DotProd Vec2 where + (&.) (Vec2 x1 y1) (Vec2 x2 y2) = x1*x2 + y1*y2 + +instance Pointwise Vec2 where + pointwise (Vec2 x1 y1) (Vec2 x2 y2) = Vec2 (x1*x2) (y1*y2) + +instance Determinant (Vec2,Vec2) where + det (Vec2 x1 y1 , Vec2 x2 y2) = x1*y2 - x2*y1 + +{- +instance Show Vec2 where + show (Vec2 x y) = "( " ++ show x ++ " , " ++ show y ++ " )" +-} + +instance Random Vec2 where + random = randomR (Vec2 (-1) (-1),Vec2 1 1) + randomR (Vec2 a b, Vec2 c d) gen = + let (x,gen1) = randomR (a,c) gen + (y,gen2) = randomR (b,d) gen1 + in (Vec2 x y, gen2) + +instance Storable Vec2 where + sizeOf _ = 2 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p k + return (Vec2 x y) + + poke q (Vec2 x y) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p k y + +-- Mat2 instances + +instance HasCoordinates Mat2 Vec2 where + _1 (Mat2 x _) = x + _2 (Mat2 _ y) = y + _3 _ = error "has only 2 coordinates" + _4 _ = error "has only 2 coordinates" + +instance Matrix Mat2 where + transpose (Mat2 row1 row2) = + Mat2 (Vec2 (_1 row1) (_1 row2)) + (Vec2 (_2 row1) (_2 row2)) + idmtx = Mat2 (Vec2 1 0) (Vec2 0 1) + inverse (Mat2 (Vec2 a b) (Vec2 c d)) = + Mat2 (Vec2 (d*r) (-b*r)) (Vec2 (-c*r) (a*r)) + where r = 1.0 / (a*d - b*c) + +instance AbelianGroup Mat2 where + (&+) (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) + +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 + (.*.) (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 LeftModule Mat2 Vec2 where + lmul (Mat2 row1 row2) v = Vec2 (row1 &. v) (row2 &. v) + +instance RightModule Vec2 Mat2 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec2 Mat2 where + diag (Vec2 x y) = Mat2 (Vec2 x 0) (Vec2 0 y) + +instance Tensor Mat2 Vec2 where + 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 + +{- +instance Show Mat2 where + show (Mat2 r1 r2) = show r1 ++ "\n" ++ show r2 +-} + +instance Storable Mat2 where + sizeOf _ = 2 * sizeOf (undefined::Vec2) + alignment _ = alignment (undefined::Vec2) + + peek q = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + r1 <- peek p + r2 <- peekByteOff p k + return (Mat2 r1 r2) + + poke q (Mat2 r1 r2) = do + let p = castPtr q :: Ptr Vec2 + k = sizeOf (undefined::Vec2) + poke p r1 + pokeByteOff p k r2 + +-- Vec3 instances + +instance HasCoordinates Vec3 Flt where + _1 (Vec3 x _ _) = x + _2 (Vec3 _ y _) = y + _3 (Vec3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance AbelianGroup Vec3 where + (&+) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1+x2) (y1+y2) (z1+z2) + (&-) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1-x2) (y1-y2) (z1-z2) + neg (Vec3 x y z) = Vec3 (-x) (-y) (-z) + zero = Vec3 0 0 0 + +instance Vector Vec3 where + scalarMul s (Vec3 x y z) = Vec3 (s*x) (s*y) (s*z) + mapVec f (Vec3 x y z) = Vec3 (f x) (f y) (f z) + +instance DotProd Vec3 where + (&.) (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = x1*x2 + y1*y2 + z1*z2 + +instance Pointwise Vec3 where + pointwise (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1*x2) (y1*y2) (z1*z2) + +{- +instance Show Vec3 where + show (Vec3 x y z) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " )" +-} + +instance Random Vec3 where + random = randomR (Vec3 (-1) (-1) (-1),Vec3 1 1 1) + randomR (Vec3 a b c, Vec3 d e f) gen = + let (x,gen1) = randomR (a,d) gen + (y,gen2) = randomR (b,e) gen1 + (z,gen3) = randomR (c,f) gen2 + in (Vec3 x y z, gen3) + +instance CrossProd Vec3 where + crossprod (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (y1*z2-y2*z1) (z1*x2-z2*x1) (x1*y2-x2*y1) + +instance Determinant (Vec3,Vec3,Vec3) where + det (u,v,w) = u &. (v &^ w) + +instance Storable Vec3 where + sizeOf _ = 3 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + return (Vec3 x y z) + + poke q (Vec3 x y z) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + +-- Mat3 instances + +instance HasCoordinates Mat3 Vec3 where + _1 (Mat3 x _ _) = x + _2 (Mat3 _ y _) = y + _3 (Mat3 _ _ z) = z + _4 _ = error "has only 3 coordinates" + +instance Matrix Mat3 where + + transpose (Mat3 row1 row2 row3) = + Mat3 (Vec3 (_1 row1) (_1 row2) (_1 row3)) + (Vec3 (_2 row1) (_2 row2) (_2 row3)) + (Vec3 (_3 row1) (_3 row2) (_3 row3)) + + idmtx = Mat3 (Vec3 1 0 0) (Vec3 0 1 0) (Vec3 0 0 1) + + inverse (Mat3 (Vec3 a b c) (Vec3 e f g) (Vec3 i j k)) = + Mat3 (Vec3 (d11*r) (d21*r) (d31*r)) + (Vec3 (d12*r) (d22*r) (d32*r)) + (Vec3 (d13*r) (d23*r) (d33*r)) + where + r = 1.0 / ( a*d11 + b*d12 + c*d13 ) + + d11 = f*k - g*j + d12 = g*i - e*k + d13 = e*j - f*i + + d31 = b*g - c*f + d32 = c*e - a*g + d33 = a*f - b*e + + d21 = c*j - b*k + d22 = a*k - c*i + d23 = b*i - a*j + +instance AbelianGroup Mat3 where + (&+) (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) + +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 + (.*.) (Mat3 r1 r2 r3) n = + let (Mat3 c1 c2 c3) = transpose n + in Mat3 (Vec3 (r1 &. c1) (r1 &. c2) (r1 &. c3)) + (Vec3 (r2 &. c1) (r2 &. c2) (r2 &. c3)) + (Vec3 (r3 &. c1) (r3 &. c2) (r3 &. c3)) + one = idmtx + +instance LeftModule Mat3 Vec3 where + lmul (Mat3 row1 row2 row3) v = Vec3 (row1 &. v) (row2 &. v) (row3 &. v) + +instance RightModule Vec3 Mat3 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec3 Mat3 where + diag (Vec3 x y z) = Mat3 (Vec3 x 0 0) (Vec3 0 y 0) (Vec3 0 0 z) + +instance Tensor Mat3 Vec3 where + outer (Vec3 a b c) (Vec3 x y z) = Mat3 + (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) + +{- +instance Show Mat3 where + show (Mat3 r1 r2 r3) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 +-} + +instance Storable Mat3 where + sizeOf _ = 3 * sizeOf (undefined::Vec3) + alignment _ = alignment (undefined::Vec3) + + peek q = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + return (Mat3 r1 r2 r3) + + poke q (Mat3 r1 r2 r3) = do + let p = castPtr q :: Ptr Vec3 + k = sizeOf (undefined::Vec3) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + +-- Vec4 instances + +instance HasCoordinates Vec4 Flt where + _1 (Vec4 x _ _ _) = x + _2 (Vec4 _ y _ _) = y + _3 (Vec4 _ _ z _) = z + _4 (Vec4 _ _ _ w) = w + +instance AbelianGroup Vec4 where + (&+) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1+x2) (y1+y2) (z1+z2) (w1+w2) + (&-) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1-x2) (y1-y2) (z1-z2) (w1-w2) + neg (Vec4 x y z w) = Vec4 (-x) (-y) (-z) (-w) + zero = Vec4 0 0 0 0 + +instance Vector Vec4 where + scalarMul s (Vec4 x y z w) = Vec4 (s*x) (s*y) (s*z) (s*w) + mapVec f (Vec4 x y z w) = Vec4 (f x) (f y) (f z) (f w) + +instance DotProd Vec4 where + (&.) (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = x1*x2 + y1*y2 + z1*z2 + w1*w2 + +instance Pointwise Vec4 where + pointwise (Vec4 x1 y1 z1 w1) (Vec4 x2 y2 z2 w2) = Vec4 (x1*x2) (y1*y2) (z1*z2) (w1*w2) + +{- +instance Show Vec4 where + show (Vec4 x y z w) = "( " ++ show x ++ " , " ++ show y ++ " , " ++ show z ++ " , " ++ show w ++ " )" +-} + +instance Random Vec4 where + random = randomR (Vec4 (-1) (-1) (-1) (-1),Vec4 1 1 1 1) + randomR (Vec4 a b c d, Vec4 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 (Vec4 x y z w, gen4) + +instance Storable Vec4 where + sizeOf _ = 4 * sizeOf (undefined::Flt) + alignment _ = sizeOf (undefined::Flt) + + peek q = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + x <- peek p + y <- peekByteOff p (k ) + z <- peekByteOff p (k+k) + w <- peekByteOff p (3*k) + return (Vec4 x y z w) + + poke q (Vec4 x y z w) = do + let p = castPtr q :: Ptr Flt + k = sizeOf (undefined::Flt) + poke p x + pokeByteOff p (k ) y + pokeByteOff p (k+k) z + pokeByteOff p (3*k) w + +-- Mat4 instances + +instance HasCoordinates Mat4 Vec4 where + _1 (Mat4 x _ _ _) = x + _2 (Mat4 _ y _ _) = y + _3 (Mat4 _ _ z _) = z + _4 (Mat4 _ _ _ w) = w + +instance Matrix Mat4 where + transpose (Mat4 row1 row2 row3 row4) = + Mat4 (Vec4 (_1 row1) (_1 row2) (_1 row3) (_1 row4)) + (Vec4 (_2 row1) (_2 row2) (_2 row3) (_2 row4)) + (Vec4 (_3 row1) (_3 row2) (_3 row3) (_3 row4)) + (Vec4 (_4 row1) (_4 row2) (_4 row3) (_4 row4)) + idmtx = Mat4 (Vec4 1 0 0 0) (Vec4 0 1 0 0) (Vec4 0 0 1 0) (Vec4 0 0 0 1) + inverse = error "inverse/Mat4: not implemented yet" + +instance AbelianGroup Mat4 where + (&+) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &+ s1) (r2 &+ s2) (r3 &+ s3) (r4 &+ s4) + (&-) (Mat4 r1 r2 r3 r4) (Mat4 s1 s2 s3 s4) = Mat4 (r1 &- s1) (r2 &- s2) (r3 &- s3) (r4 &- s4) + neg (Mat4 r1 r2 r3 r4) = Mat4 (neg r1) (neg r2) (neg r3) (neg r4) + zero = Mat4 zero zero zero zero + +instance Vector Mat4 where + 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 + (.*.) (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)) + (Vec4 (r2 &. c1) (r2 &. c2) (r2 &. c3) (r2 &. c4)) + (Vec4 (r3 &. c1) (r3 &. c2) (r3 &. c3) (r3 &. c4)) + (Vec4 (r4 &. c1) (r4 &. c2) (r4 &. c3) (r4 &. c4)) + one = idmtx + +instance LeftModule Mat4 Vec4 where + lmul (Mat4 row1 row2 row3 row4) v = Vec4 (row1 &. v) (row2 &. v) (row3 &. v) (row4 &. v) + +instance RightModule Vec4 Mat4 where + rmul v mt = lmul (transpose mt) v + +instance Diagonal Vec4 Mat4 where + diag (Vec4 x y z w) = Mat4 (Vec4 x 0 0 0) (Vec4 0 y 0 0) (Vec4 0 0 z 0) (Vec4 0 0 0 w) + +instance Tensor Mat4 Vec4 where + outer (Vec4 a b c d) (Vec4 x y z w) = Mat4 + (Vec4 (a*x) (a*y) (a*z) (a*w)) + (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 Show Mat4 where + show (Mat4 r1 r2 r3 r4) = show r1 ++ "\n" ++ show r2 ++ "\n" ++ show r3 ++ "\n" ++ show r4 +-} + +instance Storable Mat4 where + sizeOf _ = 4 * sizeOf (undefined::Vec4) + alignment _ = alignment (undefined::Vec4) + + peek q = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + r1 <- peek p + r2 <- peekByteOff p (k ) + r3 <- peekByteOff p (k+k) + r4 <- peekByteOff p (3*k) + return (Mat4 r1 r2 r3 r4) + + poke q (Mat4 r1 r2 r3 r4) = do + let p = castPtr q :: Ptr Vec4 + k = sizeOf (undefined::Vec4) + poke p r1 + pokeByteOff p (k ) r2 + pokeByteOff p (k+k) r3 + pokeByteOff p (3*k) r4 + +-- Extend instances + +instance Extend Vec2 Vec3 where + extendZero (Vec2 x y) = Vec3 x y 0 + extendWith t (Vec2 x y) = Vec3 x y t + trim (Vec3 x y _) = Vec2 x y + +instance Extend Vec2 Vec4 where + extendZero (Vec2 x y) = Vec4 x y 0 0 + extendWith t (Vec2 x y) = Vec4 x y t t + trim (Vec4 x y _ _) = Vec2 x y + +instance Extend Vec3 Vec4 where + extendZero (Vec3 x y z) = Vec4 x y z 0 + extendWith t (Vec3 x y z) = Vec4 x y z t + 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" + 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" + 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" + trim (Mat4 p q r _) = Mat3 (trim p) (trim q) (trim r) + diff --git a/src/flt/GramSchmidt.hs b/src/flt/GramSchmidt.hs new file mode 100644 --- /dev/null +++ b/src/flt/GramSchmidt.hs @@ -0,0 +1,134 @@ + +-- | Gram-Schmidt orthogonalization. +-- This module is not re-exported by "Data.Vect". + +module Data.Vect.Flt.GramSchmidt + ( GramSchmidt(..) + ) + where + +import Data.Vect.Flt.Base + +------------------------------------------------------- + +liftPair :: (a -> b) -> (a,a) -> (b,b) +liftPair f (x,y) = (f x, f y) + +liftTriple :: (a -> b) -> (a,a,a) -> (b,b,b) +liftTriple f (x,y,z) = (f x, f y, f z) + +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 + gramSchmidt :: a -> a -- ^ does not normalize the vectors! + gramSchmidtNormalize :: a -> a -- ^ normalizes the vectors. + +{-# RULES +"gramSchmidt is idempotent" forall a. gramSchmidt (gramSchmidt a) = gramSchmidt a +"gramSchmidtNormalize is idempotent" forall a. gramSchmidtNormalize (gramSchmidtNormalize a) = gramSchmidtNormalize a + #-} + +------------------------------------------------------- + +instance GramSchmidt (Vec2,Vec2) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec3,Vec3) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +instance GramSchmidt (Vec4,Vec4) where + gramSchmidt = gramSchmidtPair + gramSchmidtNormalize = gramSchmidtNormalizePair + +---------- + +instance GramSchmidt (Normal2,Normal2) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal2!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +instance GramSchmidt (Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftPair toNormalUnsafe . gramSchmidtNormalizePair . liftPair fromNormal + +---------- + +gramSchmidtPair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtPair (u,v) = (u',v') where + u' = u + v' = project v u' + +gramSchmidtNormalizePair :: (Vector v, DotProd v) => (v,v) -> (v,v) +gramSchmidtNormalizePair (u,v) = (u',v') where + u' = normalize u + v' = normalize $ projectUnsafe v u' + +---------- + +instance GramSchmidt (Vec3,Vec3,Vec3) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtTriple + gramSchmidtNormalize = gramSchmidtNormalizeTriple + +instance GramSchmidt (Normal3,Normal3,Normal3) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal3!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +instance GramSchmidt (Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftTriple toNormalUnsafe . gramSchmidtNormalizeTriple . liftTriple fromNormal + +---------- + +gramSchmidtTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtTriple (u,v,w) = (u',v',w') where + u' = u + v' = project v u' + w' = project (project w u') v' + +gramSchmidtNormalizeTriple :: (Vector v, DotProd v) => (v,v,v) -> (v,v,v) +gramSchmidtNormalizeTriple (u,v,w) = (u',v',w') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + +---------- + +instance GramSchmidt (Vec4,Vec4,Vec4,Vec4) where + gramSchmidt = gramSchmidtQuadruple + gramSchmidtNormalize = gramSchmidtNormalizeQuadruple + +instance GramSchmidt (Normal4,Normal4,Normal4,Normal4) where + gramSchmidt = error "use 'gramSchmidtNormalize' for Normal4!" + gramSchmidtNormalize = liftQuadruple toNormalUnsafe . gramSchmidtNormalizeQuadruple . liftQuadruple fromNormal + +---------- + +gramSchmidtQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtQuadruple (u,v,w,z) = (u',v',w',z') where + u' = u + v' = project v u' + w' = project (project w u') v' + z' = project (project (project z u') v') w' + +gramSchmidtNormalizeQuadruple :: (Vector v, DotProd v) => (v,v,v,v) -> (v,v,v,v) +gramSchmidtNormalizeQuadruple (u,v,w,z) = (u',v',w',z') where + u' = normalize $ u + v' = normalize $ projectUnsafe v u' + w' = normalize $ projectUnsafe (projectUnsafe w u') v' + z' = normalize $ projectUnsafe (projectUnsafe (projectUnsafe z u') v') w' + +---------- + diff --git a/src/flt/Interpolate.hs b/src/flt/Interpolate.hs new file mode 100644 --- /dev/null +++ b/src/flt/Interpolate.hs @@ -0,0 +1,42 @@ + +-- | Interpolation of vectors. +-- Note: we interpolate unit vectors differently from ordinary vectors. + +module Data.Vect.Flt.Interpolate where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim2 (sinCos',angle2') +import Data.Vect.Flt.Util.Dim3 (rotate3') + +class Interpolate v where + interpolate :: Flt -> v -> v -> v + +instance Interpolate Flt where + interpolate t x y = x + t*(y-x) + +instance Interpolate Vec2 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec3 where interpolate t x y = x &+ t *& (y &- x) +instance Interpolate Vec4 where interpolate t x y = x &+ t *& (y &- x) + +instance Interpolate Normal2 where + interpolate t nx ny = sinCos' $ ax + t*adiff where + ax = angle2' nx + ay = angle2' ny + adiff = helper (ay - ax) + helper d + | d < -pi = d + twopi + | d > pi = d - twopi + | otherwise = d + twopi = 2*pi + +instance Interpolate Normal3 where + interpolate t nx ny = + if maxAngle < 0.001 -- more or less ad-hoc critical angle + then mkNormal $ interpolate t x y + else toNormalUnsafe $ rotate3' (t*maxAngle) (mkNormal axis) x where + x = fromNormal nx + y = fromNormal ny + axis = (x &^ y) + maxAngle = acos (x &. y) + + diff --git a/src/flt/OpenGL.hs b/src/flt/OpenGL.hs new file mode 100644 --- /dev/null +++ b/src/flt/OpenGL.hs @@ -0,0 +1,182 @@ + +-- 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'. + +module Data.Vect.Flt.OpenGL where + +import Control.Monad +import Data.Vect.Flt.Base +import qualified Graphics.Rendering.OpenGL as GL + +import Foreign + +import Graphics.Rendering.OpenGL hiding (Normal3,rotate,translate,scale) + +------------------------------------------------------- + +{-# SPECIALISE radianToDegrees :: Float -> Float #-} +{-# SPECIALISE radianToDegrees :: Double -> Double #-} +radianToDegrees :: RealFrac a => a -> a +radianToDegrees x = x * 57.295779513082322 + +{-# SPECIALIZE degreesToRadian :: Float -> Float #-} +{-# SPECIALIZE degreesToRadian :: Double -> Double #-} +degreesToRadian :: Floating a => a -> a +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) + +translate :: Vec3 -> IO () +translate (Vec3 x y z) = GL.translate (Vector3 x y z) + +scale3 :: Vec3 -> IO () +scale3 (Vec3 x y z) = GL.scale x y z + +scale :: Flt -> IO () +scale x = GL.scale x x x + +------------------------------------------------------- + +-- Vertex instances + +instance GL.Vertex Vec2 where + vertex (Vec2 x y) = GL.vertex (GL.Vertex2 x y) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec3 where + vertex (Vec3 x y z) = GL.vertex (GL.Vertex3 x y z) + vertexv p = peek p >>= vertex + +instance GL.Vertex Vec4 where + 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) + normalv p = peek p >>= normal + +------------------------------------------------------- + +-- Color instances + +instance GL.Color Vec3 where + color (Vec3 r g b) = GL.color (GL.Color3 r g b) + colorv p = peek p >>= color + +instance GL.Color Vec4 where + color (Vec4 r g b a) = GL.color (GL.Color4 r g b a) + colorv p = peek p >>= color + +instance GL.SecondaryColor Vec3 where + secondaryColor (Vec3 r g b) = GL.secondaryColor (GL.Color3 r g b) + secondaryColorv p = peek p >>= secondaryColor + +{- +-- there is no such thing? +instance GL.SecondaryColor Vec4 where + secondaryColor (Vec4 r g b a) = GL.secondaryColor (GL.Color4 r g b a) + secondaryColorv p = peek p >>= secondaryColor +-} + +------------------------------------------------------- + +-- TexCoord instances + +instance GL.TexCoord Vec2 where + texCoord (Vec2 u v) = GL.texCoord (GL.TexCoord2 u v) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec2 u v) = GL.multiTexCoord unit (GL.TexCoord2 u v) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec3 where + texCoord (Vec3 u v w) = GL.texCoord (GL.TexCoord3 u v w) + texCoordv p = peek p >>= texCoord + multiTexCoord unit (Vec3 u v w) = GL.multiTexCoord unit (GL.TexCoord3 u v w) + multiTexCoordv unit p = peek p >>= multiTexCoord unit + +instance GL.TexCoord Vec4 where + texCoord (Vec4 u v w z) = GL.texCoord (GL.TexCoord4 u v w z) + texCoordv p = peek p >>= texCoord + 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 + vertexAttrib :: GL.AttribLocation -> a -> IO () + +instance VertexAttrib' {- ' CPP is sensitive to primes -} Flt where + vertexAttrib loc x = GL.vertexAttrib1 loc x + +instance VertexAttrib' Vec2 where + vertexAttrib loc (Vec2 x y) = GL.vertexAttrib2 loc x y + +instance VertexAttrib' Vec3 where + vertexAttrib loc (Vec3 x y z) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Vec4 where + 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 + +instance VertexAttrib' Normal3 where + vertexAttrib loc (Normal3 (Vec3 x y z)) = GL.vertexAttrib3 loc x y z + +instance VertexAttrib' Normal4 where + vertexAttrib loc (Normal4 (Vec4 x y z w)) = GL.vertexAttrib4 loc x y z w + +------------------------------------------------------- + +-- Uniform (again, experimental) + +-- (note that the uniform location code in the OpenGL 2.2.1.1 is broken; +-- a work-around is to put a zero character at the end of uniform names) + +{- +toFloat :: Flt -> Float +toFloat = realToFrac + +fromFloat :: Float -> Flt +fromFloat = realToFrac +-} + +-- Uniforms are always floats... +#ifdef VECT_Float + +instance GL.Uniform Flt where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Index1 x) -> x) $ get (uniform loc) + setter x = ($=) (uniform loc) (Index1 x) + uniformv loc cnt ptr = uniformv loc cnt (castPtr ptr :: Ptr (Index1 Flt)) + +instance GL.Uniform Vec2 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex2 x y) -> Vec2 x y) $ get (uniform loc) + setter (Vec2 x y) = ($=) (uniform loc) (Vertex2 x y) + uniformv loc cnt ptr = uniformv loc (2*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec3 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex3 x y z) -> Vec3 x y z) $ get (uniform loc) + setter (Vec3 x y z) = ($=) (uniform loc) (Vertex3 x y z) + uniformv loc cnt ptr = uniformv loc (3*cnt) (castPtr ptr :: Ptr Flt) + +instance GL.Uniform Vec4 where + uniform loc = GL.makeStateVar getter setter where + getter = liftM (\(GL.Vertex4 x y z w) -> Vec4 x y z w) $ get (uniform loc) + setter (Vec4 x y z w) = ($=) (uniform loc) (Vertex4 x y z w) + uniformv loc cnt ptr = uniformv loc (4*cnt) (castPtr ptr :: Ptr Flt) + +#endif diff --git a/src/flt/Util/Dim2.hs b/src/flt/Util/Dim2.hs new file mode 100644 --- /dev/null +++ b/src/flt/Util/Dim2.hs @@ -0,0 +1,64 @@ + +module Data.Vect.Flt.Util.Dim2 where + +import Data.Vect.Flt.Base + +-- |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" + +destructVec2 :: [Vec2] -> [Flt] +destructVec2 [] = [] +destructVec2 ((Vec2 x y):ls) = x:y:(destructVec2 ls) + +det2 :: Vec2 -> Vec2 -> Flt +det2 u v = det (u,v) + +vec2X :: Vec2 +vec2Y :: Vec2 + +vec2X = Vec2 1 0 +vec2Y = Vec2 0 1 + +translate2X :: Flt -> Vec2 -> Vec2 +translate2Y :: Flt -> Vec2 -> Vec2 + +translate2X t (Vec2 x y) = Vec2 (x+t) y +translate2Y t (Vec2 x y) = Vec2 x (y+t) + +-- | unit vector with given angle relative to the positive X axis (in the positive direction, that is, CCW). +-- A more precise name would be @cosSin@, but that sounds bad :) +sinCos :: Flt -> Vec2 +sinCos a = Vec2 (cos a) (sin a) + +sinCos' {- ' CPP is sensitive to primes -} :: Flt -> Normal2 +sinCos' = toNormalUnsafe . sinCos + +sinCosRadius :: Flt -- ^ angle (in radians) + -> Flt -- ^ radius + -> Vec2 +sinCosRadius a r = Vec2 (r * cos a) (r * sin a) + +-- | The angle relative to the positive X axis +angle2 :: Vec2 -> Flt +angle2 (Vec2 x y) = atan2 y x + +angle2' {- ' CPP is sensitive to primes -} :: Normal2 -> Flt +angle2' = angle2 . fromNormal + +-- |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 + +rotate2 :: Flt -> Vec2 -> Vec2 +rotate2 a v = v .* (rotMatrix2 a) + +-- |Rotates counterclockwise by 90 degrees. +rotateCCW :: Vec2 -> Vec2 +rotateCCW (Vec2 x y) = Vec2 (-y) x + +-- |Rotates clockwise by 90 degrees. +rotateCW :: Vec2 -> Vec2 +rotateCW (Vec2 x y) = Vec2 y (-x) diff --git a/src/flt/Util/Dim3.hs b/src/flt/Util/Dim3.hs new file mode 100644 --- /dev/null +++ b/src/flt/Util/Dim3.hs @@ -0,0 +1,75 @@ + +module Data.Vect.Flt.Util.Dim3 where + +import Data.Vect.Flt.Base + +structVec3 :: [Flt] -> [Vec3] +structVec3 [] = [] +structVec3 (x:y:z:ls) = (Vec3 x y z):(structVec3 ls) +structVec3 _ = error "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 + +translate3X t (Vec3 x y z) = Vec3 (x+t) y z +translate3Y t (Vec3 x y z) = Vec3 x (y+t) z +translate3Z t (Vec3 x y z) = Vec3 x y (z+t) + +vec3X :: Vec3 +vec3Y :: Vec3 +vec3Z :: Vec3 + +vec3X = Vec3 1 0 0 +vec3Y = Vec3 0 1 0 +vec3Z = Vec3 0 0 1 + +rotMatrixZ :: Flt -> Mat3 +rotMatrixY :: Flt -> Mat3 +rotMatrixX :: Flt -> Mat3 + +-- These are intended for multiplication on the /right/. +-- Should be consistent with the rotation around an arbitrary axis +-- (eg, @rotMatrixY a == rotate3 a vec3Y@) +rotMatrixZ a = Mat3 (Vec3 c s 0) (Vec3 (-s) c 0) (Vec3 0 0 1) where c = cos a; s = sin a +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' angle axis v = v .* (rotMatrix3' axis angle) + +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/. +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/. +rotMatrix3' :: {- ' CPP is sensitive to primes -} Normal3 -> Flt -> Mat3 +rotMatrix3' (Normal3 v) a = + let c = cos a + s = sin a + m1 = scalarMul (1-c) (outer v v) + x = _1 v + y = _2 v + z = _3 v + m2 = Mat3 (Vec3 c ( s*z) (-s*y)) + (Vec3 (-s*z) c ( s*x)) + (Vec3 ( s*y) (-s*x) c ) + in (m1 &+ m2) + + diff --git a/src/flt/Util/Dim4.hs b/src/flt/Util/Dim4.hs new file mode 100644 --- /dev/null +++ b/src/flt/Util/Dim4.hs @@ -0,0 +1,92 @@ + +-- | Rotation around an arbitrary plane in four dimensions, and other miscellanea. +-- Not very useful for most people, and not re-exported by "Data.Vect". + +module Data.Vect.Flt.Util.Dim4 where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.GramSchmidt + +structVec4 :: [Flt] -> [Vec4] +structVec4 [] = [] +structVec4 (x:y:z:w:ls) = (Vec4 x y z w):(structVec4 ls) +structVec4 _ = error "structVec4" + +destructVec4 :: [Vec4] -> [Flt] +destructVec4 [] = [] +destructVec4 ((Vec4 x y z w):ls) = x:y:z:w:(destructVec4 ls) + +--det4 :: Vec4 -> Vec4 -> Vec4 -> Vec4 -> Flt +--det4 u v w z = det (u,v,w,z) + +translate4X :: Flt -> Vec4 -> Vec4 +translate4Y :: Flt -> Vec4 -> Vec4 +translate4Z :: Flt -> Vec4 -> Vec4 +translate4W :: Flt -> Vec4 -> Vec4 + +translate4X t (Vec4 x y z w) = Vec4 (x+t) y z w +translate4Y t (Vec4 x y z w) = Vec4 x (y+t) z w +translate4Z t (Vec4 x y z w) = Vec4 x y (z+t) w +translate4W t (Vec4 x y z w) = Vec4 x y z (w+t) + +vec4X :: Vec4 +vec4Y :: Vec4 +vec4Z :: Vec4 +vec4W :: Vec4 + +vec4X = Vec4 1 0 0 0 +vec4Y = Vec4 0 1 0 0 +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))@. +-- This is a helper function for the 4 dimensional rotation code. +-- If @(x,y,z,p,q,r) = biVector4 a b@, then the corresponding antisymmetric tensor is +-- +-- > [ 0 r q p ] +-- > [ -r 0 z -y ] +-- > [ -q -z 0 x ] +-- > [ -p y -x 0 ] +biVector4 :: Vec4 -> Vec4 -> (Flt,Flt,Flt,Flt,Flt,Flt) +biVector4 (Vec4 x y z w) (Vec4 a b c d) = + ( x*b-y*a , x*c-z*a , x*d-w*a , y*c-z*b , -y*d+w*b , z*d-w*c ) + +-- | the corresponding antisymmetric tensor +biVector4AsTensor :: Vec4 -> Vec4 -> Mat4 +biVector4AsTensor v w = + Mat4 ( Vec4 0 ( r) ( q) ( p) ) + ( Vec4 (-r) 0 ( z) (-y) ) + ( Vec4 (-q) (-z) 0 ( x) ) + ( Vec4 (-p) ( y) (-x) 0 ) + where + (x,y,z,p,q,r) = biVector4 v w + +-- | We assume that the axes are normalized and /orthogonal/ to each other! +rotate4' :: {- ' CPP is sensitive to primes -} Flt -> (Normal4,Normal4) -> Vec4 -> Vec4 +rotate4' angle axes v = v .* (rotMatrix4' angle axes) + +-- | We assume only that the axes are independent vectors. +rotate4 :: Flt -> (Vec4,Vec4) -> Vec4 -> Vec4 +rotate4 angle axes v = v .* (rotMatrix4 angle axes) + +-- | 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 + where + c = cos angle ; s = sin angle + m1 = scalarMul (1-c) ( outer v v &+ outer w w ) + m2 = biVector4AsTensor v w + m3 = diag (Vec4 c c c c) + +-- | We assume only that the axes are independent vectors. +rotMatrix4 :: Flt -> (Vec4,Vec4) -> Mat4 +rotMatrix4 angle axes = + rotMatrix4' angle $ liftPair toNormalUnsafe $ gramSchmidtNormalize axes + where + liftPair f (x,y) = (f x, f y) + + diff --git a/src/flt/Util/Projective.hs b/src/flt/Util/Projective.hs new file mode 100644 --- /dev/null +++ b/src/flt/Util/Projective.hs @@ -0,0 +1,87 @@ + +-- | Classic 4x4 projective matrices. Our convention is that they are intended for multiplication on +-- the /right/, that is, they are of the form +-- +-- > _____ +-- > [ | | 0 ] +-- > [ | 3x3 | 0 ] +-- > [ |_____| 0 ] +-- > [ p q r 1 ] +-- +-- Please note that by default, OpenGL stores the matrices (in memory) by columns, while we +-- 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 ()@. + +module Data.Vect.Flt.Util.Projective where + +import Data.Vect.Flt.Base +import Data.Vect.Flt.Util.Dim3 + +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) + +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) + +rotMatrixProj :: Flt -> Normal3 -> Mat4 +rotMatrixProj angle axis = extendProj $ rotMatrix3' axis angle + +rotMatrixProj' :: {- ' CPP is sensitive to primes -} Flt -> Vec3 -> Mat4 +rotMatrixProj' angle axis = extendProj $ rotMatrix3 axis angle + +translMatrixProj :: Vec3 -> Mat4 +translMatrixProj v = Mat4 Dim4.vec4X Dim4.vec4Y Dim4.vec4Z (extendProj v) + +-- | 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 + +scaleMatrixUniformProj :: Flt -> Mat4 +scaleMatrixUniformProj s = diag (Vec4 s s s 1) + +class ProjectiveAction v where + actProj :: v -> Mat4 -> v + +instance ProjectiveAction Vec3 where + actProj v m = trim $ (extendProj v) .* m + +instance ProjectiveAction Vec4 where + actProj v m = v .* m + +-- | When acting on unit vectors, we ignore the translation part. +instance ProjectiveAction Normal3 where + actProj (Normal3 v) m = Normal3 (v .* (trim m :: Mat3)) + +-- | 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) + where + t' = Vec3 (- u &. t) (- v &. t) (- w &. t) + diff --git a/vect.cabal b/vect.cabal new file mode 100644 --- /dev/null +++ b/vect.cabal @@ -0,0 +1,78 @@ +Name: vect +Version: 0.4.0 +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 + computer graphics (basic OpenGL support is included). + Projective 4 dimensional operations, as used in eg. + OpenGL, are also supported. + The base field is either Float or Double. +License: BSD3 +License-file: LICENSE +Author: Balazs Komuves +Copyright: (c) 2008-2009 Balazs Komuves +Maintainer: bkomuves (plus) hackage (at) gmail (dot) com +Homepage: http://code.haskell.org/~bkomuves/ +Stability: Experimental +Category: Graphics, Math +Tested-With: GHC == 6.10.1 +Cabal-Version: >= 1.6 +Build-Type: Custom + +extra-source-files: src/flt.hs, + src/flt/Base.hs, + src/flt/Interpolate.hs, + src/flt/OpenGL.hs, + src/flt/GramSchmidt.hs, + src/flt/Util/Dim2.hs, + src/flt/Util/Dim3.hs, + src/flt/Util/Dim4.hs, + src/flt/Util/Projective.hs + +Flag splitBase + Description: Choose the new smaller, split-up base package. + +Flag OpenGL + Description: Compile with OpenGL support + Default: True + +Library + if flag(splitBase) + Build-Depends: base >= 3, random + else + Build-Depends: base < 3 + + if flag(OpenGL) + Build-Depends: OpenGL + cpp-options: -DVECT_OPENGL + Exposed-Modules: Data.Vect.Float.OpenGL, + Data.Vect.Double.OpenGL + + Exposed-Modules: Data.Vect, + + Data.Vect.Float, + Data.Vect.Float.Base, + Data.Vect.Float.Interpolate, + Data.Vect.Float.GramSchmidt, + Data.Vect.Float.Util.Dim2, + Data.Vect.Float.Util.Dim3, + Data.Vect.Float.Util.Dim4, + Data.Vect.Float.Util.Projective, + + Data.Vect.Double, + Data.Vect.Double.Base, + Data.Vect.Double.Interpolate, + Data.Vect.Double.GramSchmidt, + Data.Vect.Double.Util.Dim2, + Data.Vect.Double.Util.Dim3, + Data.Vect.Double.Util.Dim4, + Data.Vect.Double.Util.Projective + + Hs-Source-Dirs: . + Extensions: ForeignFunctionInterface, CPP, + MultiParamTypeClasses, FunctionalDependencies, + FlexibleInstances, TypeSynonymInstances, + GeneralizedNewtypeDeriving + + ghc-options: -Wall +