diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,26 @@ +Copyright (c) 2006, 2007, Christopher Lane Hinson + 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 the name of Christopher Lane Hinson nor the names of its 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/RSAGL/Math.hs b/RSAGL/Math.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math.hs @@ -0,0 +1,25 @@ +module RSAGL.Math ( + module RSAGL.Math.AbstractVector, + module RSAGL.Math.Affine, + module RSAGL.Math.Angle, + module RSAGL.Math.BoundingBox, + module RSAGL.Math.Curve, + module RSAGL.Math.Interpolation, + module RSAGL.Math.Matrix, + module RSAGL.Math.Orthogonal, + module RSAGL.Math.Ray, + module RSAGL.Math.Types, + module RSAGL.Math.Vector) + where + +import RSAGL.Math.AbstractVector +import RSAGL.Math.Affine +import RSAGL.Math.Angle +import RSAGL.Math.BoundingBox +import RSAGL.Math.Curve +import RSAGL.Math.Interpolation +import RSAGL.Math.Matrix +import RSAGL.Math.Orthogonal +import RSAGL.Math.Ray +import RSAGL.Math.Types +import RSAGL.Math.Vector diff --git a/RSAGL/Math/AbstractVector.hs b/RSAGL/Math/AbstractVector.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/AbstractVector.hs @@ -0,0 +1,230 @@ +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, UndecidableInstances #-} + +-- | Provides generic typeclasses for common operations among many types: addition, subtraction, scalar multiplication, magnitude, and zero. +module RSAGL.Math.AbstractVector + (AbstractVector, + AbstractZero(..), + AbstractAdd(..), + AbstractSubtract(..), + AbstractScale(..), + AbstractMagnitude(..), + abstractScaleTo, + abstractSum, + abstractAverage, + abstractDistance) + where + +import Data.Fixed +import Control.Applicative +import RSAGL.Math.Types + +-- | A data type that has an additive identity. +class AbstractZero a where + zero :: a + +-- | A data type that supports addition. +-- +-- * @a `add` zero = a@ +class AbstractAdd p v | p -> v where + add :: p -> v -> p + +-- | A data type that supports subtraction. +-- +-- * @a `sub` a = zero@ +class AbstractSubtract p v | p -> v where + sub :: p -> p -> v + +-- | A data type that supports scalar multiplication. +-- +-- * @scalarMultiply 0 a = zero@ +class AbstractScale v where + scalarMultiply :: RSdouble -> v -> v + +-- | A data type that supports scalar magnitude. +-- +-- * @magnitude (scalarMultiply (recip $ magnitude a) a) = 1@ +class AbstractMagnitude v where + magnitude :: v -> RSdouble + +-- | A convenience class for many vector types. +class (AbstractZero v,AbstractAdd v v,AbstractSubtract v v,AbstractScale v) => AbstractVector v where + +-- Integer + +instance AbstractZero Integer where + zero = 0 + +instance AbstractAdd Integer Integer where + add = (+) + +instance AbstractSubtract Integer Integer where + sub = (-) + +instance AbstractMagnitude Integer where + magnitude = abs . fromInteger + +-- Float + +instance AbstractZero Float where + zero = 0 + +instance AbstractAdd Float Float where + add = (+) + +instance AbstractSubtract Float Float where + sub = (-) + +instance AbstractScale Float where + scalarMultiply d = (f2f d *) + +instance AbstractMagnitude Float where + magnitude = abs . f2f + +instance AbstractVector Float + +-- Double + +instance AbstractAdd Double Double where + add = (+) + +instance AbstractSubtract Double Double where + sub = (-) + +instance AbstractScale Double where + scalarMultiply d = (f2f d *) + +instance AbstractMagnitude Double where + magnitude = abs . f2f + +instance AbstractVector Double + +instance AbstractZero Double where + zero = 0 + +-- Fixed + +instance (HasResolution a) => AbstractZero (Fixed a) where + zero = 0 + +instance (HasResolution a) => AbstractAdd (Fixed a) (Fixed a) where + add = (+) + +instance (HasResolution a) => AbstractSubtract (Fixed a) (Fixed a) where + sub = (-) + +instance (HasResolution a) => AbstractScale (Fixed a) where + scalarMultiply d = (realToFrac d *) + +instance (HasResolution a) => AbstractMagnitude (Fixed a) where + magnitude = abs . realToFrac + +instance (HasResolution a) => AbstractVector (Fixed a) + +-- Tuples + +instance (AbstractZero a,AbstractZero b) => AbstractZero (a,b) where + zero = (zero,zero) + +instance (AbstractAdd a a',AbstractAdd b b') => AbstractAdd (a,b) (a',b') where + add (a,b) (c,d) = (add a c,add b d) + +instance (AbstractSubtract a a',AbstractSubtract b b') => AbstractSubtract (a,b) (a',b') where + sub (a,b) (c,d) = (sub a c,sub b d) + +instance (AbstractScale a,AbstractScale b) => AbstractScale (a,b) where + scalarMultiply d (a,b) = (scalarMultiply d a,scalarMultiply d b) + +instance (AbstractMagnitude a,AbstractMagnitude b) => AbstractMagnitude (a,b) where + magnitude (a,b) = sqrt $ magnitude a ^ 2 + magnitude b ^ 2 + +instance (AbstractVector a,AbstractVector b) => AbstractVector (a,b) + +-- Functions + +instance (AbstractAdd a a') => AbstractAdd ((->) x a) ((->) x a') where + add a b = \x -> a x `add` b x + +instance (AbstractSubtract a a') => AbstractSubtract ((->) x a) ((->) x a') where + sub a b = \x -> a x `sub` b x + +instance (AbstractScale a) => AbstractScale ((->) x a) where + scalarMultiply d f = scalarMultiply d . f + +-- RSfloat + +instance AbstractAdd RSfloat RSfloat where + add = (+) + +instance AbstractSubtract RSfloat RSfloat where + sub = (-) + +instance AbstractScale RSfloat where + scalarMultiply d = (f2f d *) + +instance AbstractMagnitude RSfloat where + magnitude = abs . f2f + +instance AbstractVector RSfloat + +instance AbstractZero RSfloat where + zero = 0 + +-- RSdouble + +instance AbstractAdd RSdouble RSdouble where + add = (+) + +instance AbstractSubtract RSdouble RSdouble where + sub = (-) + +instance AbstractScale RSdouble where + scalarMultiply = (*) + +instance AbstractMagnitude RSdouble where + magnitude = abs . f2f + +instance AbstractVector RSdouble + +instance AbstractZero RSdouble where + zero = 0 + +-- Lists + +instance (AbstractZero a) => AbstractZero [a] where + zero = repeat zero + +instance (AbstractAdd a b) => AbstractAdd [a] [b] where + add = zipWith add + +instance (AbstractSubtract a b) => AbstractSubtract [a] [b] where + sub = zipWith sub + +instance (AbstractScale a) => AbstractScale [a] where + scalarMultiply d = map (scalarMultiply d) + +instance (AbstractMagnitude a) => AbstractMagnitude [a] where + magnitude = sqrt . sum . map ((^2) . magnitude) + +instance (AbstractVector a) => AbstractVector [a] where + +-- Generic functions. + +-- | Force a vector to the specified magnitude. +abstractScaleTo :: (AbstractScale v,AbstractMagnitude v) => RSdouble -> v -> v +abstractScaleTo _ v | magnitude v == 0 = v +abstractScaleTo x v = scalarMultiply (x / magnitude v) v + +-- | Sum of a list. +abstractSum :: (AbstractAdd p v,AbstractZero p) => [v] -> p +abstractSum = foldr (flip add) zero + +-- | Average of a list. +abstractAverage :: (AbstractAdd p v,AbstractSubtract p v,AbstractVector v,AbstractZero p) => [p] -> p +abstractAverage vs = zero `add` scalarMultiply (recip $ fromInteger total_count) total_sum + where f y (i,x) = i `seq` x `seq` (i+1,y `add` x) + (total_count,total_sum) = foldr f (0,zero) $ map (`sub` zero) vs + +-- | Distance between two points, based on the 'magnitude' of the difference. +abstractDistance :: (AbstractMagnitude v,AbstractSubtract p v) => p -> p -> RSdouble +abstractDistance x y = magnitude $ x `sub` y + diff --git a/RSAGL/Math/Affine.hs b/RSAGL/Math/Affine.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Affine.hs @@ -0,0 +1,145 @@ +-- | Affine Transformations of Arbitrary Geometric Objects +module RSAGL.Math.Affine + (AffineTransformable(..), + scale', + inverseTransform, + withTransformation, + transformAbout, + translateToFrom, + rotateToFrom, + scaleAlong) + where + +import Graphics.Rendering.OpenGL.GL as GL hiding (R) +import RSAGL.Math.Vector +import RSAGL.Math.Matrix +import RSAGL.Math.Angle +import RSAGL.Math.Types +import Data.Maybe + +-- | 'AffineTransformable' objects are subject to affine transformations using matrix multiplication. +class AffineTransformable a where + -- | Apply an affine transformation, defined by a 4x4 matrix. (This is the only required method.) + transform :: RSAGL.Math.Matrix.Matrix -> a -> a + -- | Scale an entity along the @x@ @y@ and @z@ axes. For example, @scale (Vector3D 2 3 4)@ will make an object twice as wide, three times as tall, + -- and four times as deep. It may be helpful to think of the vector as a control point on the vertex of a unit cube. + scale :: Vector3D -> a -> a + scale vector = transform $ scaleMatrix vector + -- | Translate an entity along the specified vector. + translate :: Vector3D -> a -> a + translate vector = transform $ translationMatrix vector + -- | Rotate an entity about the origin, using the specified vector as the axis of rotation. + -- See also 'transformAbout' to rotate around an arbitrary point. + rotate :: Vector3D -> Angle -> a -> a + rotate vector angle = transform $ rotationMatrix vector angle + -- | Specific rotation around the x-axis. + rotateX :: Angle -> a -> a + rotateX = RSAGL.Math.Affine.rotate (Vector3D 1 0 0) + -- | Specific rotation around the y-axis. + rotateY :: Angle -> a -> a + rotateY = RSAGL.Math.Affine.rotate (Vector3D 0 1 0) + -- | Specific rotation around the z-axis. + rotateZ :: Angle -> a -> a + rotateZ = RSAGL.Math.Affine.rotate (Vector3D 0 0 1) + +-- | Apply the inverse of an affine transformation, defined by a 4x4 matrix. +{-# INLINE inverseTransform #-} +inverseTransform :: (AffineTransformable a) => RSAGL.Math.Matrix.Matrix -> a -> a +inverseTransform m = transform (matrixInverse m) + +-- | Specific scale preserving proportions. +{-# INLINE scale' #-} +scale' :: (AffineTransformable a) => RSdouble -> a -> a +scale' x = RSAGL.Math.Affine.scale (Vector3D x x x) + +-- | Apply a function under an affine transformation. @withTransformation m id@ is an identity if @m@ is invertable. +{-# INLINE withTransformation #-} +withTransformation :: (AffineTransformable a) => RSAGL.Math.Matrix.Matrix -> (a -> a) -> a -> a +withTransformation m f = inverseTransform m . f . transform m + +-- | Apply a function treating a particular point as the origin. For example, combining 'transformAbout' with 'RSAGL.Math.Affine.rotate' +-- performs a rotation about an arbitrary point rather than the origin. +{-# INLINE transformAbout #-} +transformAbout :: (AffineTransformable a) => Point3D -> (a -> a) -> a -> a +transformAbout center f = withTransformation (translateToFrom origin_point_3d center identity_matrix) f + +-- | Specific translation along the vector between two points. +-- This ordinary use is to set the second point as the center of a model (typically origin_point_3d) +-- and the first point as the desired position of the model. +{-# INLINE translateToFrom #-} +translateToFrom :: (AffineTransformable a) => Point3D -> Point3D -> a -> a +translateToFrom a b = RSAGL.Math.Affine.translate (vectorToFrom a b) + +-- | Specific rotation along the shortest path that brings the second vector in line with the first. +{-# INLINE rotateToFrom #-} +rotateToFrom :: (AffineTransformable a) => Vector3D -> Vector3D -> a -> a +rotateToFrom u v = RSAGL.Math.Affine.rotate c a + where c = vectorNormalize $ vectorScale (-1) $ fromMaybe (fst $ orthos u) $ aNonZeroVector $ crossProduct u v + a = angleBetween u v + +-- | Specific scale along an arbitary axis. +{-# INLINE scaleAlong #-} +scaleAlong :: (AffineTransformable a) => Vector3D -> RSdouble -> a -> a +scaleAlong v u = withTransformation (rotateToFrom (Vector3D 0 1 0) v identity_matrix) (RSAGL.Math.Affine.scale (Vector3D 1 u 1)) + +instance AffineTransformable a => AffineTransformable (Maybe a) where + transform m = fmap (transform m) + +instance AffineTransformable a => AffineTransformable [a] where + transform m = map (transform m) + +instance (AffineTransformable a,AffineTransformable b) => AffineTransformable (a,b) where + transform m (a,b) = (transform m a,transform m b) + +instance (AffineTransformable a,AffineTransformable b,AffineTransformable c) => AffineTransformable (a,b,c) where + transform m (a,b,c) = (transform m a,transform m b,transform m c) + +instance AffineTransformable RSAGL.Math.Matrix.Matrix where + transform mat = matrixMultiply mat + +instance AffineTransformable Vector3D where + transform m (Vector3D x y z) = transformHomogenous x y z 0 Vector3D m + scale (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) = Vector3D (x1*x2) (y1*y2) (z1*z2) + translate _ = id + rotateX a (Vector3D x y z) = Vector3D x (c*y-s*z) (c*z+s*y) + where s = sine a + c = cosine a + rotateY a (Vector3D x y z) = Vector3D (c*x+s*z) y (c*z-s*x) + where s = sine a + c = cosine a + rotateZ a (Vector3D x y z) = Vector3D (c*x-s*y) (c*y+s*x) z + where s = sine a + c = cosine a + +instance AffineTransformable Point3D where + transform m (Point3D x y z) = transformHomogenous x y z 1 Point3D m + scale (Vector3D x1 y1 z1) (Point3D x2 y2 z2) = Point3D (x1*x2) (y1*y2) (z1*z2) + translate (Vector3D x1 y1 z1) (Point3D x2 y2 z2) = Point3D (x1+x2) (y1+y2) (z1+z2) + rotateX a (Point3D x y z) = Point3D x (c*y-s*z) (c*z+s*y) + where s = sine a + c = cosine a + rotateY a (Point3D x y z) = Point3D (c*x+s*z) y (c*z-s*x) + where s = sine a + c = cosine a + rotateZ a (Point3D x y z) = Point3D (c*x-s*y) (c*y+s*x) z + where s = sine a + c = cosine a + +instance AffineTransformable SurfaceVertex3D where + transform m (SurfaceVertex3D p v) = SurfaceVertex3D (RSAGL.Math.Affine.transform m p) (RSAGL.Math.Affine.transform (matrixTranspose $ matrixInverse m) v) + translate vector (SurfaceVertex3D p v) = SurfaceVertex3D (RSAGL.Math.Affine.translate vector p) v + +-- | The IO monad itself is AffineTransformable. This is done by wrapping the IO action in an OpenGL transformation. +instance AffineTransformable (IO a) where + transform mat iofn = preservingMatrix $ do mat' <- newMatrix RowMajor $ map f2f $ concat $ rowMajorForm mat + multMatrix (mat' :: GLmatrix GLdouble) + iofn + translate (Vector3D x y z) iofn = preservingMatrix $ + do GL.translate $ Vector3 (f2f x) (f2f y) (f2f z :: GLdouble) + iofn + scale (Vector3D x y z) iofn = preservingMatrix $ + do GL.scale (f2f x) (f2f y) (f2f z :: GLdouble) + iofn + rotate (Vector3D x y z) angle iofn = preservingMatrix $ + do GL.rotate (f2f $ toDegrees_ angle) (Vector3 (f2f x) (f2f y) (f2f z :: GLdouble)) + iofn diff --git a/RSAGL/Math/Angle.hs b/RSAGL/Math/Angle.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Angle.hs @@ -0,0 +1,193 @@ +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} +module RSAGL.Math.Angle + (Angle, + BoundAngle(..), + fromDegrees, + fromRadians, + fromRotations, + fromTimeOfDayHMS, + fromArcMinutes, + fromArcSeconds, + sine, + arcSine, + cosine, + arcCosine, + tangent, + arcTangent, + cartesianToPolar, + polarToCartesian, + toRadians, + toRadians_, + toDegrees, + toDegrees_, + toRotations, + toRotations_, + scaleAngle, + supplementaryAngle, + zero_angle, + angularIncrements, + angleAdd, + angleSubtract, + angleNegate, + absoluteAngle, + unboundAngle) + where + +import RSAGL.Math.FMod +import RSAGL.Math.AbstractVector +import RSAGL.Math.Types + +-- | An angular value. +newtype Angle = Radians RSdouble deriving (Read,Show) + +-- | An angular value. 'BoundAngle's are always in the range between -180 and 180 degrees, inclusive. +newtype BoundAngle = BoundAngle Angle deriving (Read,Show) + +zero_angle :: Angle +zero_angle = Radians 0 + +instance Eq Angle where + (==) a b = case (toRadians_ a,toRadians_ b) of + (x,y) | abs x == pi && abs y == pi -> True + (x,y) | x == y -> True + _ -> False + +instance Ord Angle where + compare x y = case () of + _ | x == y -> EQ + _ -> compare (toRadians_ x) (toRadians_ y) + +instance AbstractZero Angle where + zero = zero_angle + +instance AbstractZero BoundAngle where + zero = BoundAngle zero_angle + +instance AbstractAdd Angle Angle where + add = angleAdd + +instance AbstractAdd BoundAngle Angle where + add (BoundAngle a) x = BoundAngle $ boundAngle $ a `add` x + +instance AbstractSubtract Angle Angle where + sub = angleSubtract + +instance AbstractSubtract BoundAngle Angle where + sub (BoundAngle a) (BoundAngle b) = boundAngle $ a `sub` b + +instance AbstractScale Angle where + scalarMultiply = scaleAngle + +instance AbstractVector Angle + +instance AbstractMagnitude Angle where + magnitude = toRotations_ . absoluteAngle + +-- | angularIncrements answers n evenly distributed angles from 0 to 2*pi. +angularIncrements :: Integer -> [Angle] +angularIncrements subdivisions = map (fromRadians . (2*pi*) . (/ fromInteger subdivisions) . fromInteger) [0 .. subdivisions - 1] + +-- | There are 2*pi radians in a circle. +fromRadians :: RSdouble -> Angle +fromRadians = Radians + +-- | There are 260 degrees in a circle. +fromDegrees :: RSdouble -> Angle +fromDegrees = Radians . ((*) (pi/180)) + +-- | There is 1 rotation in a circle. +fromRotations :: RSdouble -> Angle +fromRotations = Radians . ((*) (2*pi)) + +-- | Get an angle based on time of day, hours, minutes, seconds, where noon is considered a zero angle. +fromTimeOfDayHMS :: RSdouble -> RSdouble -> RSdouble -> Angle +fromTimeOfDayHMS h m s = fromRotations (((s/60+m)/60+h)/24) + +-- | There are 21600 arc minutes in a circle, 60 arc minutes in a degree. +fromArcMinutes :: RSdouble -> Angle +fromArcMinutes = fromDegrees . (/60) + +-- | There are 1296000 arc seconds in a circle, 60 arc seconds in an arc minutes. +fromArcSeconds :: RSdouble -> Angle +fromArcSeconds = fromArcMinutes . (/60) + +-- | Answers the angle in the range of -180 to 180, inclusive. +toDegrees :: Angle -> RSdouble +toDegrees x = let x' = toRadians x + in x' * 180 / pi + + +-- | 'toDegrees_' answers the angle in degrees with no range limitation. +toDegrees_ :: Angle -> RSdouble +toDegrees_ (Radians x) = x * 180 / pi + +-- | 'toRadians' answers the angle in the range of -pi .. pi, inclusive. +toRadians :: Angle -> RSdouble +toRadians x = let (Radians x') = boundAngle x + in x' + + +-- | toRadians answers the angle in radians with no range limitation. +toRadians_ :: Angle -> RSdouble +toRadians_ (Radians x) = x + +-- | 'toRotations' answers the angle in the range of -0.5 to 0.5, inclusive. +toRotations :: Angle -> RSdouble +toRotations x= let x' = toRadians x + in x' / pi / 2 + + +-- | 'toRotations' answers the angle in rotations with no range limitation. +toRotations_ :: Angle -> RSdouble +toRotations_ (Radians x) = x / pi / 2 + +scaleAngle :: RSdouble -> Angle -> Angle +scaleAngle x = Radians . (*x) . toRadians_ + +supplementaryAngle :: Angle -> Angle +supplementaryAngle (Radians x) = Radians $ pi - x + +angleAdd :: Angle -> Angle -> Angle +angleAdd (Radians x) (Radians y) = Radians $ x + y + +angleSubtract :: Angle -> Angle -> Angle +angleSubtract (Radians x) (Radians y) = Radians $ x - y + +angleNegate :: Angle -> Angle +angleNegate (Radians x) = Radians $ negate x + +-- | Absolute value of an angle. +absoluteAngle :: Angle -> Angle +absoluteAngle (Radians x) = Radians $ abs x + +sine :: Angle -> RSdouble +sine (Radians x) = sin x + +arcSine :: RSdouble -> Angle +arcSine = fromRadians . asin + +cosine :: Angle -> RSdouble +cosine (Radians x) = cos x + +arcCosine :: RSdouble -> Angle +arcCosine = fromRadians . acos + +tangent :: Angle -> RSdouble +tangent (Radians x) = tan x + +arcTangent :: RSdouble -> Angle +arcTangent = fromRadians . atan + +cartesianToPolar :: (RSdouble,RSdouble) -> (Angle,RSdouble) +cartesianToPolar (u,v) = (fromRadians $ atan2 v u,sqrt $ u*u + v*v) + +polarToCartesian :: (Angle,RSdouble) -> (RSdouble,RSdouble) +polarToCartesian (a,d) = (cosine a*d,sine a*d) + +-- | 'boundAngle' forces the angle into the range (-pi..pi). +boundAngle :: Angle -> Angle +boundAngle (Radians x) = Radians $ if bounded > pi then bounded - 2*pi else bounded + where bounded = x `fmod` (2*pi) + +unboundAngle :: BoundAngle -> Angle +unboundAngle (BoundAngle a) = a diff --git a/RSAGL/Math/BoundingBox.hs b/RSAGL/Math/BoundingBox.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/BoundingBox.hs @@ -0,0 +1,87 @@ +module RSAGL.Math.BoundingBox + (BoundingBox, + Bound3D(..), + boundingCenterRadius, + minimalDistanceToBoundingBox) + where + +import RSAGL.Math.Vector +import RSAGL.Math.Interpolation +import RSAGL.Math.Affine +import RSAGL.Math.Types + +-- | A simple bounding box. Operations on bounding boxes +-- are designed to be fast, not accurate. The only +-- guarantee on any bounding box operation is that +-- objects reported to be outside a bounding box, +-- are. +data BoundingBox = BoundingBox { + bbox_bottom, bbox_top, bbox_left, + bbox_right, bbox_far, bbox_near :: !RSdouble } + deriving (Show) + +-- | A convenience class for any finite geometry. +-- In particular, it's easy to concatenate the bounding +-- box of multiple geometries by placing them in a list +-- and taking the bounding box of the entire list. +class Bound3D a where + boundingBox :: a -> BoundingBox + +instance Bound3D Point3D where + boundingBox (Point3D x y z) = BoundingBox { + bbox_bottom = y, + bbox_top = y, + bbox_right = x, + bbox_left = x, + bbox_near = z, + bbox_far = z } + +instance Bound3D SurfaceVertex3D where + boundingBox = boundingBox . sv3d_position + +instance (Bound3D a) => Bound3D [a] where + boundingBox [] = error "instance Bound3D [a], boundingBox: can't construct boundingBox []" + boundingBox (x:xs) = foldr combineBoundingBoxes (boundingBox x) $ xs + +instance Bound3D BoundingBox where + boundingBox = id + +instance AffineTransformable BoundingBox where + transform m = boundingBox . transform m . boundingBoxToPointCloud + +combineBoundingBoxes :: (Bound3D a,Bound3D b) => a -> b -> BoundingBox +combineBoundingBoxes x y = + BoundingBox { + bbox_left = min (bbox_left b1) (bbox_left b2), + bbox_right = max (bbox_right b1) (bbox_right b2), + bbox_bottom = min (bbox_bottom b1) (bbox_bottom b2), + bbox_top = max (bbox_top b1) (bbox_top b2), + bbox_near = min (bbox_near b1) (bbox_near b2), + bbox_far = max (bbox_far b1) (bbox_far b2) } + where b1 = boundingBox x + b2 = boundingBox y + +boundingBoxToPointCloud :: BoundingBox -> [Point3D] +boundingBoxToPointCloud bbox = + [Point3D (bbox_left bbox) (bbox_bottom bbox) (bbox_near bbox), + Point3D (bbox_right bbox) (bbox_bottom bbox) (bbox_near bbox), + Point3D (bbox_left bbox) (bbox_top bbox) (bbox_near bbox), + Point3D (bbox_right bbox) (bbox_top bbox) (bbox_near bbox), + Point3D (bbox_left bbox) (bbox_bottom bbox) (bbox_far bbox), + Point3D (bbox_right bbox) (bbox_bottom bbox) (bbox_far bbox), + Point3D (bbox_left bbox) (bbox_top bbox) (bbox_far bbox), + Point3D (bbox_right bbox) (bbox_top bbox) (bbox_far bbox)] + +-- | View of a bounding box in the form of a bounding spehre. +boundingCenterRadius :: BoundingBox -> (Point3D,RSdouble) +boundingCenterRadius bbox = (lerp 0.5 (nlb,frt),distanceBetween nlb frt / 2) + where nlb = Point3D (bbox_near bbox) (bbox_left bbox) (bbox_bottom bbox) + frt = Point3D (bbox_far bbox) (bbox_right bbox) (bbox_top bbox) + +-- | Estimates distance between a point and the outside surface of a bounding +-- box. If the value is negative, then the point lies inside the bound +-- region. +minimalDistanceToBoundingBox :: Point3D -> BoundingBox -> RSdouble +minimalDistanceToBoundingBox p bbox = distanceBetween p c - r + where (c,r) = boundingCenterRadius bbox + diff --git a/RSAGL/Math/Curve.hs b/RSAGL/Math/Curve.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Curve.hs @@ -0,0 +1,302 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving, ExistentialQuantification, Rank2Types #-} +module RSAGL.Math.Curve + (Curve, + zipCurve, + iterateCurve, + transposeCurve, + curve, + Surface, + surface, + wrapSurface, + unwrapSurface, + transposeSurface, + zipSurface, + iterateSurface, + halfIterateSurface, + flipTransposeSurface, + translateCurve, + scaleCurve, + clampCurve, + loopCurve, + controlCurve, + transformCurve2, + uv_identity, + translateSurface, + scaleSurface, + transformSurface, + transformSurface2, + surfaceDerivative, + curveDerivative, + orientationLoops, + newellCurve, + surfaceNormals3D, + SamplingAlgorithm, + IntervalSample, + intervalRange, + intervalSize, + intervalSample, + intervalSingleIntegral, + linearSamples, + adaptiveMagnitudeSamples, + integrateCurve) + where + +import RSAGL.Math.ListUtils +import RSAGL.Math.Vector +import RSAGL.Math.Angle +import RSAGL.Math.Affine +import Data.List +import Data.Maybe +import Control.Parallel.Strategies +import Control.Applicative +import RSAGL.Math.AbstractVector +import Debug.Trace +import RSAGL.Math.Interpolation +import RSAGL.Math.FMod +import RSAGL.Math.Types +import RSAGL.Math.BoundingBox + +-- | A parametric function that is aware of it's own sampling interval. The first parameter is the sampling interval, while the second is the curve input parameter. +type CurveF a = (RSdouble,RSdouble) -> a +-- | A surface is a curve of curves. +type SurfaceF a = CurveF (CurveF a) +-- | A 'Curve' is a parametric function from a one-dimensional space into a space of an arbitrary datatype. The key feature of a 'Curve' is that it is aware of it's own +-- sampling interval. Using this information and appropriate arithmetic and scalar multiplication functions provided by RSAGL.AbstractVector, a 'Curve' can be differentiated or integrated. +newtype Curve a = Curve { fromCurve :: CurveF a } + +instance Functor Curve where + fmap g (Curve f) = Curve $ g . f + +instance Applicative Curve where + pure a = Curve $ const a + f <*> a = zipCurve ($) f a + +instance (AffineTransformable a) => AffineTransformable (Curve a) where + scale v = fmap (scale v) + translate v = fmap (translate v) + rotate vector angle = fmap (rotate vector angle) + transform m = fmap (transform m) + rotateX angle = fmap (rotateX angle) + rotateY angle = fmap (rotateY angle) + rotateZ angle = fmap (rotateZ angle) + +instance NFData (Curve a) where + rnf (Curve f) = seq f () + +-- | Sample a specific point on a curve given the sampling interval (first parameter) and the point on the curve (second parameter). +sampleCurve :: Curve a -> RSdouble -> RSdouble -> a +sampleCurve (Curve f) = curry f + +-- | Sample a curve at regular intervals in the range 0..1 inclusive. +iterateCurve :: Integer -> Curve x -> [x] +iterateCurve n c = map f $ zeroToOne n + where f = sampleCurve c (0.25/fromInteger n) + +-- | Combine two curves using an arbitrary function. +zipCurve :: (x -> y -> z) -> Curve x -> Curve y -> Curve z +zipCurve f (Curve x) (Curve y) = Curve $ \hu -> f (x hu) (y hu) + +-- | Arbitrarily transform a 'Curve'. +mapCurve :: (CurveF a -> CurveF a) -> Curve a -> Curve a +mapCurve f = Curve . f . fromCurve + +-- | Arbitrarily transform a surface 'Curve'. +mapCurve2 :: (SurfaceF a -> SurfaceF a) -> Curve (Curve a) -> Curve (Curve a) +mapCurve2 f = Curve . (Curve .) . f . (fromCurve .) . fromCurve + +-- | Translate a curve along the axis of the input parameter. +translateCurve :: RSdouble -> Curve a -> Curve a +translateCurve x = mapCurve $ \f (h,u) -> f (h,u-x) + +-- | Scale a curve along the axis of the input parameter. Factors greater than one have a "zoom in" effect, while +-- factors less than one have a "zoom out" effect. +scaleCurve :: RSdouble -> Curve a -> Curve a +scaleCurve s = mapCurve (\f (h,u) -> f (h/s,u/s)) + +-- | Clamp lower and upper bounds of a curve along the axis of the input parameter. +clampCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a +clampCurve (a,b) | b < a = clampCurve (b,a) +clampCurve (a,b) = mapCurve $ \f (h,u) -> f (h,min b $ max a u) + +-- | Loop a curve onto itself at the specified bounds. +loopCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a +loopCurve (a,b) | b < a = loopCurve (b,a) +loopCurve (a,b) = mapCurve $ \f (h,u) -> f (h,(u-a) `fmod` (b-a) + a) + +-- | Transform a curve by manipulating control points. +controlCurve :: (RSdouble,RSdouble) -> (RSdouble,RSdouble) -> Curve a -> Curve a +controlCurve (u,v) (u',v') = translateCurve u' . scaleCurve ((u'-v') / (u-v)) . translateCurve (negate u) + +-- | Transpose the inner and outer components of a curve. +transposeCurve :: Curve (Curve a) -> Curve (Curve a) +transposeCurve = mapCurve2 flip + +-- | Lift two curve transformations onto each axis of a second order curve. +transformCurve2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Curve (Curve a) -> Curve (Curve a) +transformCurve2 fu fv = transposeCurve . fu . transposeCurve . fv + +-- | Define a simple curve. +curve :: (RSdouble -> a) -> Curve a +curve = Curve . uncurry . const + +-- | A 'Surface' is a based on a 'Curve' with an output of another 'Curve'. +newtype Surface a = Surface (Curve (Curve a)) deriving (NFData,AffineTransformable) + +-- | Define a simple surface. +surface :: (RSdouble -> RSdouble -> a) -> Surface a +surface f = Surface $ curve (\u -> curve $ flip f u) + +wrapSurface :: Curve (Curve a) -> Surface a +wrapSurface = Surface + +unwrapSurface :: Surface a -> Curve (Curve a) +unwrapSurface (Surface s) = s + +-- | Transpose the axes of a 'Surface'. +transposeSurface :: Surface a -> Surface a +transposeSurface (Surface s) = Surface $ transposeCurve s + +-- | Sample a surface at regularly spaced lattice points in the range 0..1 inclusive. +iterateSurface :: (Integer,Integer) -> Surface a -> [[a]] +iterateSurface (u,v) (Surface s) = map (iterateCurve u) $ iterateCurve v s + +-- | Sample the outer 'Curve' of a 'Surface' at regularly spaced intervals. +halfIterateSurface :: Integer -> Surface a -> [Curve a] +halfIterateSurface u = iterateCurve u . unwrapSurface + +instance Functor Surface where + fmap f (Surface x) = Surface $ fmap (fmap f) x + +instance Applicative Surface where + pure a = surface (const $ const a) + f <*> a = zipSurface ($) f a + +-- | Combine two surfaces using an arbitrary function. +zipSurface :: (x -> y -> z) -> Surface x -> Surface y -> Surface z +zipSurface f (Surface x) (Surface y) = Surface $ zipCurve (zipCurve f) x y + +-- | Lift a transformation on a second order 'Curve' onto a Surface. +transformSurface :: (Curve (Curve a) -> Curve (Curve a)) -> Surface a -> Surface a +transformSurface f = Surface . f . unwrapSurface + +-- | Lift two curve transformations onto each axis of a Surface. +transformSurface2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Surface a -> Surface a +transformSurface2 fu fv = transformSurface (transformCurve2 fu fv) + +-- | Translate a surface over each of its input parameter axes, as translateCurve. +translateSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a +translateSurface (u,v) = transformSurface2 (translateCurve u) (translateCurve v) + +-- | Scale a surface along each of its input parameter axes, as scaleCurve. +scaleSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a +scaleSurface (u,v) = transformSurface2 (scaleCurve u) (scaleCurve v) + +-- | Transpose a surface while flipping the inner curve, so that that orientable surfaces retain their original orientation. +flipTransposeSurface :: Surface a -> Surface a +flipTransposeSurface = transformSurface (mapCurve2 $ \f (hu,u) (hv,v) -> f (hu,u) (hv,1-v)) . transposeSurface + +-- | An identity 'Surface'. +uv_identity :: Surface (RSdouble,RSdouble) +uv_identity = surface (curry id) + +-- | Take the derivative of a 'Curve'. +curveDerivative :: (AbstractSubtract p v,AbstractScale v) => Curve p -> Curve v +curveDerivative (Curve f) = Curve $ \(h,u) -> scalarMultiply (recip $ 2 * h) $ f (h/2,u+h) `sub` f (h/2,u-h) + +-- | Take the piecewise derivative of a 'Surface' along the inner and outer curves. +surfaceDerivative :: (AbstractSubtract p v,AbstractScale v) => Surface p -> Surface (v,v) +surfaceDerivative s = zipSurface (,) (curvewiseDerivative s) (transposeSurface $ curvewiseDerivative $ transposeSurface s) + where curvewiseDerivative (Surface t) = Surface $ fmap curveDerivative t + +-- | Determine the orientation of a 'Surface' by passing very small circles centered on each sampled point as the parametric input. +-- +-- A gotchya with this operation is that as much as 3/4ths of the orientation loop may lie outside of the 0..1 range that is normally +-- sampled. Depending on how the surface is constructed, this may produce unexpected results. The solution is to clamp the +-- the problematic parametric inputs at 0 and 1 using 'clampSurface'. +-- +-- As a rule, do clamp longitudinal axes that come to a singularity at each end. +-- Do not clamp latitudinal axes that are connected at each end. +-- +orientationLoops :: Surface p -> Surface (Curve p) +orientationLoops (Surface s) = Surface $ Curve $ \(uh,u) -> Curve $ \(vh,v) -> + curve $ \t -> f (uh/2,u + uh*(sine $ fromRotations t)) + (vh/2,v + vh*(cosine $ fromRotations t)) + where f = fromCurve . fromCurve s + +-- | Try to determine the normal vector to a curve. +newellCurve :: Curve Point3D -> Maybe Vector3D +newellCurve c = newell $ iterateCurve 16 c + +-- | Try to determine the normal vectors of a surface using orientation loops. This is usually slower but more successful than 'surfaceNormals3DByPartialDerivatives'. +-- This generate a warning message if it cannot determine the normal vector at a sampled point. +-- See also 'orientationLoops'. +surfaceNormals3DByOrientationLoops :: Surface Point3D -> Surface SurfaceVertex3D +surfaceNormals3DByOrientationLoops s = SurfaceVertex3D <$> s <*> ((\c -> errmsg c (newellCurve c)) <$> orientationLoops s) + where errmsg c = fromMaybe (trace ("surfaceNormals3DByOrientationLoops: zero normal gave up: " ++ show (iterateCurve 16 c)) (Vector3D 0 0 0)) + +-- | Try to determine the normal vectors of a surface using partial derivatives. +surfaceNormals3DByPartialDerivatives :: Surface Point3D -> Surface (Maybe Vector3D) +surfaceNormals3DByPartialDerivatives s = safeCrossProduct <$> surfaceDerivative s + where x = snd $ boundingCenterRadius $ boundingBox $ concat $ iterateSurface (8,8) s + safeCrossProduct (u_,v_) = + do u <- aLargeVector (x/100) u_ + v <- aLargeVector (x/100) v_ + return $ vectorNormalize $ crossProduct u v + +-- | Try to determine the normal vectors of a surface using multiple techniques. +surfaceNormals3D :: Surface Point3D -> Surface SurfaceVertex3D +surfaceNormals3D s = (\p by_pd by_newell -> case by_pd of + Just v -> SurfaceVertex3D p v + Nothing -> by_newell) <$> + s <*> (surfaceNormals3DByPartialDerivatives s) <*> (surfaceNormals3DByOrientationLoops s) + +-- | An interval of a curve, including the curve, lower and upper bounds of the interval, and an instantaneous sample value for that interval. +data IntervalSample a = IntervalSample (Curve a) RSdouble RSdouble a + +intervalSample :: Curve a -> RSdouble -> RSdouble -> IntervalSample a +intervalSample c l h = IntervalSample c l h $ sampleCurve c ((abs $ l - h) / 2) ((l+h) / 2) + +-- | Lower and upper bounds of an 'IntervalSample'. +intervalRange :: IntervalSample a -> (RSdouble,RSdouble) +intervalRange (IntervalSample _ l h _) = (l,h) + +-- | Size of the range of an 'IntervalSample'. +intervalSize :: IntervalSample a -> RSdouble +intervalSize (IntervalSample _ l h _) = abs $ h - l + +-- | Instantaneous sample value of an 'Interval'. +intervalValue :: IntervalSample a -> a +intervalValue (IntervalSample _ _ _ a) = a + +-- | Integral of the sample value over the range of the 'IntervalSample'. +intervalSingleIntegral :: (AbstractScale a) => IntervalSample a -> a +intervalSingleIntegral x = scalarMultiply (intervalSize x) $ intervalValue x + +-- | Split an interval into three equal parts. +splitInterval :: IntervalSample a -> [IntervalSample a] +splitInterval (IntervalSample c l h a) = [intervalSample c l l',IntervalSample c l' h' a,intervalSample c h' h] + where l' = lerp (1/3) (l,h) + h' = lerp (2/3) (l,h) + +type SamplingAlgorithm a = Curve a -> [IntervalSample a] + +-- | Definite integral of a curve. +integrateCurve :: (AbstractAdd p v,AbstractScale v,AbstractZero p) => SamplingAlgorithm v -> Curve v -> p -> p +integrateCurve samplingAlgorithm c initial_value = foldl' add initial_value $ map intervalSingleIntegral $ samplingAlgorithm c + +-- | Sampling algorithm that takes a fixed count of samples. +linearSamples :: Integer -> SamplingAlgorithm a +linearSamples n c = map (\(l,h) -> intervalSample c l h) $ doubles $ zeroToOne (n+1) + +-- | Sampling algorithm that takes increasing numbers of samples over intervals where the magnitude of the sample is large. +adaptiveMagnitudeSamples :: (AbstractMagnitude a) => Integer -> SamplingAlgorithm a +adaptiveMagnitudeSamples n c = resampleLoop (\xs -> if genericLength xs > n then Nothing else Just $ newSamples xs) $ linearSamples (max 1 $ n `div` 10) c + where newSamples xs = let a = abstractAverage $ map intervalMagnitude xs + in flip concatMap xs $ \x -> if intervalMagnitude x >= a then splitInterval x else [x] + intervalMagnitude :: (AbstractMagnitude a) => IntervalSample a -> RSdouble + intervalMagnitude (IntervalSample _ l h a) = magnitude a * (abs $ h-l) + +-- | Loop to keep generating samples until finished. +resampleLoop :: (b -> Maybe b) -> b -> b +resampleLoop nextPass initial_value = f $ initial_value + where f x = maybe x f $ nextPass x diff --git a/RSAGL/Math/CurveExtras.lhs b/RSAGL/Math/CurveExtras.lhs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/CurveExtras.lhs @@ -0,0 +1,110 @@ +\section{Specific and Interpolated Curves} + +\begin{code} +module RSAGL.Math.CurveExtras + (sphericalCoordinates, + cylindricalCoordinates, + toroidalCoordinates, + circularCoordinates, + polarCoordinates, + transformUnitSquareToUnitCircle, + transformUnitCubeToUnitSphere, + circleXY, + regularPolygon, + linearInterpolation, + loopedLinearInterpolation, + smoothCurve, + loopCurve) + where + +import RSAGL.Math.Curve +import RSAGL.Math.Interpolation +import RSAGL.Math.Vector +import RSAGL.Math.Angle +import RSAGL.Math.AbstractVector +import RSAGL.Math.Affine +import RSAGL.Math.ListUtils +import Control.Arrow +import RSAGL.Math.Types +\end{code} + +\subsection{Alternate Coordinate Systems for Models} + +\begin{code} +sphericalCoordinates :: ((Angle,Angle) -> a) -> Surface a +sphericalCoordinates f = transformSurface2 id (clampCurve (0,1)) $ surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,fromRadians $ ((pi/2) - v*pi)))) + +cylindricalCoordinates :: ((Angle,RSdouble) -> a) -> Surface a +cylindricalCoordinates f = transformSurface2 id (clampCurve (0,1)) $ surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,v))) + +toroidalCoordinates :: ((Angle,Angle) -> a) -> Surface a +toroidalCoordinates f = surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,fromRadians $ negate $ v*2*pi))) + +circularCoordinates :: ((RSdouble,RSdouble) -> a) -> Surface a +circularCoordinates f = surface $ curry $ (f . second negate . transformUnitSquareToUnitCircle) + +polarCoordinates :: ((Angle,RSdouble) -> a) -> Surface a +polarCoordinates f = circularCoordinates (f . cartesianToPolar) +\end{code} + +\subsection{Transformations Between Unit Volumes} + +\begin{code} +transformUnitSquareToUnitCircle :: (RSdouble,RSdouble) -> (RSdouble,RSdouble) +transformUnitSquareToUnitCircle (u,v) = (x,z) + where (Point3D x _ z) = transformUnitCubeToUnitSphere (Point3D u 0.5 v) + +transformUnitCubeToUnitSphere :: Point3D -> Point3D +transformUnitCubeToUnitSphere p = + let p_centered@(Point3D x y z) = scale' 2.0 $ translate (Vector3D (-0.5) (-0.5) (-0.5)) p + p_projected = scale' (minimum [recip $ abs x,recip $ abs y,recip $ abs z]) p_centered + k = recip $ distanceBetween origin_point_3d p_projected + w = maximum $ [abs x, abs y, abs z] -- 'w' could be 1, but this gives a smoother tesselation + in if p_centered == origin_point_3d then origin_point_3d else lerp w (p_centered,scale' k p_centered) +\end{code} + +\subsection{Circles} + +\begin{code} +circleXY :: Curve Point3D +circleXY = curve $ \u_ -> let u = fromRotations u_ in Point3D (cosine u) (sine u) 0 +\end{code} + +\subsection{Regular Polygons} + +A regular polygon, centered at the origin, in the XY plane. + +\begin{code} +regularPolygon :: Integer -> Curve Point3D +regularPolygon n = loopedLinearInterpolation $ map (flip rotateZ (Point3D 0 1 0) . fromRotations) $ zeroToOne n +\end{code} + +\subsection{Piecewise Length Normalization} + +\begin{code} +normalizePolyline :: (AbstractSubtract p v,AbstractMagnitude v) => [p] -> [(RSdouble,p)] +normalizePolyline pts = zip (map (/ total_dist) accumulated_dists) pts + where dists = map (uncurry abstractDistance) $ doubles pts + total_dist = last accumulated_dists + accumulated_dists = scanl (+) 0 dists +\end{code} + +\subsection{Interpolated Curves} + +\begin{code} +linearInterpolation :: (AbstractSubtract p v,AbstractAdd p v,AbstractMagnitude v,AbstractScale v) => [p] -> Curve p +linearInterpolation = curve . lerpMap . normalizePolyline + +loopedLinearInterpolation :: (AbstractSubtract p v,AbstractAdd p v,AbstractMagnitude v,AbstractScale v) => [p] -> Curve p +loopedLinearInterpolation = loopCurve (0,1) . linearInterpolation . (\a -> last a:a) +\end{code} + +\subsection{Smoothing Curves} + +\texttt{smoothCurve i h} takes i samples of an h-long piece of a 'Curve' at each point to smooth it. This is not an interpolation function and will tend to shrink shapes +toward their center of gravity. + +\begin{code} +smoothCurve :: (AbstractAdd p v,AbstractSubtract p v,AbstractVector v,AbstractZero p) => Integer -> RSdouble -> Curve p -> Curve p +smoothCurve i h c = curve $ \u -> abstractAverage $ iterateCurve i $ controlCurve (u-h/2,u+h/2) (0,1) c +\end{code} diff --git a/RSAGL/Math/FMod.hs b/RSAGL/Math/FMod.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/FMod.hs @@ -0,0 +1,24 @@ +module RSAGL.Math.FMod + (FMod(..)) + where + +import RSAGL.Math.Types + +class FMod f where + fmod :: f -> f -> f + +instance FMod Float where + fmod n d = n - fromInteger f * d + where f = floor $ n / d + +instance FMod Double where + fmod n d = n - fromInteger f * d + where f = floor $ n / d + +instance FMod RSfloat where + fmod n d = n - fromInteger f * d + where f = floor $ n / d + +instance FMod RSdouble where + fmod n d = n - fromInteger f * d + where f = floor $ n / d diff --git a/RSAGL/Math/Interpolation.lhs b/RSAGL/Math/Interpolation.lhs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Interpolation.lhs @@ -0,0 +1,109 @@ +\section{RSAGL.Interpolation} + +\begin{code} +module RSAGL.Math.Interpolation + (lerp, + lerpClamped, + lerpBetween, + lerpBetweenMutated, + lerpBetweenClamped, + lerpBetweenClampedMutated, + lerp_mutator_continuous_1st, + lerpMap) + where + +import RSAGL.Math.AbstractVector +import Data.Map as Map +import Data.Maybe +import RSAGL.Math.Types +\end{code} + +\subsection{The Lerpable typeclass} + +Implements linear interpolation. + +\begin{code} +{-# INLINE lerp #-} +lerp :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => r -> (p,p) -> p +lerp u (a,b) = a `add` scalarMultiply (f2f u) (b `sub` a) +\end{code} + +\subsection{Non-linear interpolations} + +The mutated versions of the lerp functions takes an arbitrary ``mutator function'' to define non-linear interpolation curves. +The ``between'' versions of the lerp functions allow the u-value to lie between any two numbers, as opposed to between 0 and 1. +The ``clamped'' versions of the lerp functions clamp the u-value to lie between its boundaries. Otherwise, with non-clamped +interpolations, the u-value may lie outside of its boundaries. + +\begin{code} +{-# INLINE lerpClamped #-} + +lerpClamped :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => r -> (p,p) -> p +lerpClamped u = lerpBetweenClamped (0,u,1) + +{-# INLINE lerpBetween #-} +lerpBetween :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => (r,r,r) -> (p,p) -> p +lerpBetween = lerpBetweenMutated id + +{-# INLINE lerpBetweenMutated #-} +lerpBetweenMutated :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => (r -> r) -> (r,r,r) -> (p,p) -> p +lerpBetweenMutated _ (l,_,r) | l == r = lerp 0.5 +lerpBetweenMutated mutator (l,u,r) = lerp $ mutator $ (u-l) / (r-l) + +{-# INLINE lerpBetweenClamped #-} +lerpBetweenClamped :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r,Ord r) => (r,r,r) -> (p,p) -> p +lerpBetweenClamped = lerpBetweenClampedMutated id + +{-# INLINE lerpBetweenClampedMutated #-} +lerpBetweenClampedMutated :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r,Ord r) => (r -> r) -> (r,r,r) -> (p,p) -> p +lerpBetweenClampedMutated mutator (l,u,r) = lerpBetweenMutated (lerp_mutator_clamp . mutator) (l,u,r) +\end{code} + +\subsection{Lerp mutators} + +\texttt{lerp\_mutator\_clamp} implements clamping between 0 and 1. + +\begin{code} +{-# INLINE lerp_mutator_clamp #-} + +lerp_mutator_clamp :: (RealFloat r) => r -> r +lerp_mutator_clamp = min 1 . max 0 +\end{code} + +\texttt{lerp\_mutator\_continuous\_1st} implements clamping between 0 and 1, but such that the 1st derivative of the result is continuous. + +\begin{code} +{-# INLINE lerp_mutator_continuous_1st #-} + +lerp_mutator_continuous_1st :: (RealFloat r) => r -> r +lerp_mutator_continuous_1st x | x < 0 = 0 +lerp_mutator_continuous_1st x | x > 1 = 1 +lerp_mutator_continuous_1st x | x <= 0.5 = 2*x^2 +lerp_mutator_continuous_1st x = 4*x - 2*x^2 - 1 +\end{code} + +\subsection{lerpMap} + +Given many entities, lerp between the two entities closest to the given point +on either side. For example, if we wanted to lerp between colors on a rainbow, +we might use the map [(0,red),(1,orange),(2,yellow),(3,green),(4,blue),(5,indigo),(6,violet)]. +lerpMap 3.5 would result in a blue-green color. + +\begin{code} +{-# INLINE lerpMap #-} + +lerpMap :: (RealFloat r,AbstractScale v,AbstractSubtract p v,AbstractAdd p v) => [(r,p)] -> r -> p +lerpMap pts = lerpMapSorted $ Map.fromList pts + +{-# INLINE lerpMapSorted #-} +lerpMapSorted :: (RealFloat r,AbstractScale v,AbstractSubtract p v,AbstractAdd p v) => Map r p -> r -> p +lerpMapSorted m _ | Map.null m = error "lerpMapSorted: empty map" +lerpMapSorted m u = case () of + () | isJust exactly -> fromJust exactly + () | Map.null lowers -> higher_a + () | Map.null highers -> lower_a + () -> lerpBetween (lower_r,u,higher_r) (lower_a,higher_a) + where (lowers, exactly, highers) = splitLookup u m + (lower_r,lower_a) = findMax lowers + (higher_r,higher_a) = findMin highers +\end{code} diff --git a/RSAGL/Math/ListUtils.hs b/RSAGL/Math/ListUtils.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/ListUtils.hs @@ -0,0 +1,53 @@ +module RSAGL.Math.ListUtils + (doubles, + loopedDoubles, + consecutives, + loopedConsecutives, + zeroToOne) + where + +import RSAGL.Math.Types +import Debug.Trace + +-- | Transforms a list to a list of adjacent elements. +-- +-- @doubles [1,2,3,4,5] = [(1,2),(2,3),(3,4),(4,5)]@ +doubles :: [a] -> [(a,a)] +doubles [] = [] +doubles [_] = [] +doubles (x:y:zs) = (x,y) : doubles (y:zs) + +-- | loopedDoubles transforms a list to a list of adjacent elements, looping +-- back to the beginning of the list. +-- +-- @loopedRSdoubles [1,2,3,4,5] = [(1,2),(2,3),(3,4),(4,5),(5,1)]@ +loopedDoubles :: [a] -> [(a,a)] +loopedDoubles as = loopedDoubles_ (head as) as + where loopedDoubles_ _ [] = [] + loopedDoubles_ a [x] = [(x,a)] + loopedDoubles_ a (x:y:zs) = (x,y) : loopedDoubles_ a (y:zs) + +-- | Answers a list containing every sequence of n consecutive +-- elements in the parameter. +-- +-- @consecutives 3 [1,2,3,4] = [[1,2,3],[2,3,4]]@ +consecutives :: Int -> [a] -> [[a]] +consecutives n xs = let taken = take n xs + in if (length taken == n) + then (taken : (consecutives n $ tail xs)) + else [] + +-- | Answers a list containing every sequence of n consecutive +-- elements in the parameter, looping back to the beginning of the list. +-- +-- @consecutives 3 [1,2,3,4] = [[1,2,3],[2,3,4],[3,4,1],[4,1,2]]@ +loopedConsecutives :: Int -> [a] -> [[a]] +loopedConsecutives n xs = consecutives n $ take (n + length xs - 1) $ cycle xs + +-- | Creates a list of numbers from 0.0 to 1.0, using n steps. +-- This can't be done with the enum-from-to method, due to roundoff errors. +zeroToOne :: Integer -> [RSdouble] +zeroToOne n | n > 100000 = trace ("Warning: zeroToOne was asked for " ++ show n ++ " subdivisions, which seems high. Using 100,000 instead.") zeroToOne 100000 +zeroToOne n = map (*x) [0..(fromInteger $ n-1)] + where x = recip (fromInteger $ n - 1) + diff --git a/RSAGL/Math/Matrix.lhs b/RSAGL/Math/Matrix.lhs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Matrix.lhs @@ -0,0 +1,190 @@ +\section{RSAGL.Matrix} + +\begin{code} + +module RSAGL.Math.Matrix + (Matrix, + matrix, + transformHomogenous, + rowMajorForm, + colMajorForm, + rowAt, + matrixAt, + identity_matrix, + translationMatrix, + rotationMatrix, + scaleMatrix, + xyzMatrix, + matrixAdd, + matrixMultiply, + matrixTranspose, + matrixInverse, + determinant, + matrixInversePrim, + matrixTransposePrim, + matrixInverseTransposePrim, + determinantPrim) + where + +import Data.List as List +import RSAGL.Math.Angle +import RSAGL.Math.Vector +import Data.Vec as Vec +import RSAGL.Math.Types +\end{code} + +A 4-by-4 matrix with cached inverse, transpose, inverse transpose, and determinant. + +\begin{code} +data Matrix = Matrix { matrix_data :: !(Mat44 RSdouble), + matrix_inverse :: Matrix, + matrix_transpose :: Matrix, + matrix_determinant :: RSdouble } + +instance Eq Matrix where + x == y = matrix_data x == matrix_data y + +instance Show Matrix where + show m = show $ rowMajorForm m + +rowMajorForm :: Matrix -> [[RSdouble]] +rowMajorForm = matToLists . matrix_data + +colMajorForm :: Matrix -> [[RSdouble]] +colMajorForm = List.transpose . rowMajorForm +\end{code} + +rowAt answers the nth row of a matrix. + +\begin{code} +rowAt :: Matrix -> Int -> [RSdouble] +rowAt m n = (rowMajorForm m) !! n +\end{code} + +matrixAt answers the (i'th,j'th) element of a matrix. + +\begin{code} +matrixAt :: Matrix -> (Int,Int) -> RSdouble +matrixAt m (i,j) = rowAt m i !! j +\end{code} + +\subsection{Constructing matrices} + +matrix constructs a matrix from row major list form. (Such a list form can be formatted correctly in +monospaced font and haskell syntax, so that it looks like a matrix as it would be normally written.) + +\begin{code} +matrix :: [[RSdouble]] -> Matrix +matrix = uncheckedMatrix . matFromLists + +-- | Generate a column matrix of length 4, perform an affine transformation on it, and produce the resulting value. +{-# INLINE transformHomogenous #-} +transformHomogenous :: RSdouble -> RSdouble -> RSdouble -> RSdouble -> (RSdouble -> RSdouble -> RSdouble -> a) -> Matrix -> a +transformHomogenous x y z w f m = f x' y' z' + where (x':.y':.z':._) = multmv (matrix_data m) (x:.y:.z:.w:.()) + +uncheckedMatrix :: Mat44 RSdouble -> Matrix +uncheckedMatrix dats = m + where m_inverse = (matrixInversePrim m) { matrix_inverse = m, matrix_transpose = m_inverse_transpose, matrix_determinant = recip m_det } + m_transpose = (matrixTransposePrim m) { matrix_inverse = m_inverse_transpose, matrix_transpose = m, matrix_determinant = m_det } + m_inverse_transpose = (matrixTransposePrim m_inverse) { matrix_inverse = m_transpose, matrix_transpose = m_inverse, matrix_determinant = recip m_det } + m_det = determinantPrim m + m = Matrix { matrix_data=dats, + matrix_inverse = m_inverse, + matrix_transpose = m_transpose, + matrix_determinant = m_det } +\end{code} + +identityMatrix constructs the n by n identity matrix + +\begin{code} +identity_matrix :: Matrix +identity_matrix = Matrix { + matrix_data = identity, + matrix_inverse = identity_matrix, + matrix_transpose = identity_matrix, + matrix_determinant = 1 } +\end{code} + +\begin{code} +translationMatrix :: Vector3D -> Matrix +translationMatrix (Vector3D x y z) = uncheckedMatrix $ (1:.0:.0:.x:.()) :. + (0:.1:.0:.y:.()) :. + (0:.0:.1:.z:.()) :. + (0:.0:.0:.1:.()) :. () +\end{code} + +\begin{code} +rotationMatrix :: Vector3D -> Angle -> Matrix +rotationMatrix vector angle = let s = sine angle + c = cosine angle + c' = 1 - c + (Vector3D x y z) = vectorNormalize vector + in uncheckedMatrix $ ((c+c'*x*x) :. (c'*y*x - s*z) :. (c'*z*x + s*y) :. 0 :. ()) :. + ((c'*x*y+s*z) :. (c+c'*y*y) :. (c'*z*y - s*x) :. 0 :. ()) :. + ((c'*x*z-s*y) :. (c'*y*z+s*x) :. (c+c'*z*z) :. 0 :. ()) :. + (0 :. 0 :. 0 :. 1 :. ()) :. () +\end{code} + +\begin{code} +scaleMatrix :: Vector3D -> Matrix +scaleMatrix (Vector3D x y z) = uncheckedMatrix $ (x :. 0 :. 0 :. 0) :. + (0 :. y :. 0 :. 0) :. + (0 :. 0 :. z :. 0) :. + (0 :. 0 :. 0 :. 1) :. () +\end{code} + +\texttt{xyzMatrix} constructs the matrix in which the x y and z axis are transformed to point in the direction of the specified +vectors. + +\begin{code} +xyzMatrix :: Vector3D -> Vector3D -> Vector3D -> Matrix +xyzMatrix (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) (Vector3D x3 y3 z3) = + uncheckedMatrix $ (x1 :. x2 :. x3 :. 0) :. + (y1 :. y2 :. y3 :. 0) :. + (z1 :. z2 :. z3 :. 0) :. + (0 :. 0 :. 0 :. 1.0) :. () +\end{code} + +\subsection{Matrix arithmetic} + +\begin{code} +matrixAdd :: Matrix -> Matrix -> Matrix +matrixAdd m n = uncheckedMatrix $ Vec.zipWith (+) (matrix_data m) (matrix_data n) +\end{code} + +\label{howtoUseMatrixMultiplyImpl} + +\begin{code} +matrixMultiply :: Matrix -> Matrix -> Matrix +matrixMultiply m n = uncheckedMatrix $ multmm (matrix_data m) (matrix_data n) + +matrixTranspose :: Matrix -> Matrix +matrixTranspose = matrix_transpose +\end{code} + +\subsection{Inverse and determinant of a matrix} + +\begin{code} +matrixInverse :: Matrix -> Matrix +matrixInverse = matrix_inverse +\end{code} + +\begin{code} +determinant :: Matrix -> RSdouble +determinant = matrix_determinant +\end{code} + +\begin{code} +matrixInverseTransposePrim :: Matrix -> Matrix +matrixInverseTransposePrim = uncheckedMatrix . Vec.transpose . fst . invertAndDet . matrix_data + +matrixInversePrim :: Matrix -> Matrix +matrixInversePrim = uncheckedMatrix . fst . invertAndDet . matrix_data + +determinantPrim :: Matrix -> RSdouble +determinantPrim = det . matrix_data + +matrixTransposePrim :: Matrix -> Matrix +matrixTransposePrim = uncheckedMatrix . Vec.transpose . matrix_data +\end{code} diff --git a/RSAGL/Math/Orthogonal.hs b/RSAGL/Math/Orthogonal.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Orthogonal.hs @@ -0,0 +1,116 @@ + +-- | It's useful to work with the set of coordinate systems restricted to those +-- that use orthogonal unit-scaled axes, that is, that are subject only to +-- rotation and translation. This is because these coordinate systems are the +-- describe rigid objects. +module RSAGL.Math.Orthogonal + (up,down,left,right,forward,backward, + orthogonalFrame, + modelLookAt, + FUR) + where + +import RSAGL.Math.Affine +import RSAGL.Math.Vector +import RSAGL.Math.Matrix + +-- | FUR stands for Forward Up Right. It's used to specify arbitrary +-- orthogonal coordinate systems given any combination of forward up and right +-- vectors. It also accepts down, left, and backward vectors. + +-- | This modules uses a left-handed coordinate system. Right is positive X, +-- up is positive Y, forward is positive Z. + +-- When specifying FUR coordinate systems, the first vector is fixed +-- while the second vector will be adjusted as little as possible to guarantee +-- that it is orthogonal to the first. The third vector never needs to be +-- specified, it can be deduced. +data FURAxis = ForwardAxis | UpAxis | RightAxis | DownAxis | LeftAxis | BackwardAxis + +data FUR a = FUR FURAxis a + +instance Functor FUR where + fmap f (FUR a x) = FUR a $ f x + +-- | A reference to the +Y axis. +up :: a -> FUR a +up = FUR UpAxis + +-- | A reference to the -Y axis. +down :: a -> FUR a +down = FUR DownAxis + +-- | A reference to the -X axis. +left :: a -> FUR a +left = FUR LeftAxis + +-- | A reference to the +X axis. +right :: a -> FUR a +right = FUR RightAxis + +-- | A reference to the +Z axis. +forward :: a -> FUR a +forward = FUR ForwardAxis + +-- | A reference to the -Z axis. +backward :: a -> FUR a +backward = FUR BackwardAxis + +-- | Combine two axial references to describe a rigid affine transformation. +-- Accepts any combination of non-coaxial references. +-- In the affine transformation, the old axes will be mapped onto the specified +-- freeform axes. +-- +-- The first parameter is absolute, meaning that the source axis will always map +-- perfectly onto the destination axis. The second parameter will be obeyed +-- on a "best effort" basis. +-- +orthogonalFrame :: (AffineTransformable a) => FUR Vector3D -> FUR Vector3D -> a -> a +orthogonalFrame (FUR ForwardAxis f) (FUR RightAxis r) = + let (r',u') = fixOrtho2 f r in transform (xyzMatrix r' u' (vectorNormalize f)) +orthogonalFrame (FUR UpAxis u) (FUR ForwardAxis f) = + let (f',r') = fixOrtho2 u f in transform (xyzMatrix r' (vectorNormalize u) f') +orthogonalFrame (FUR RightAxis r) (FUR UpAxis u) = + let (u',f') = fixOrtho2 r u in transform (xyzMatrix (vectorNormalize r) u' f') +orthogonalFrame (FUR RightAxis r) (FUR ForwardAxis f) = + let (f',u') = fixOrtho2Left r f in transform (xyzMatrix (vectorNormalize r) u' f') +orthogonalFrame (FUR ForwardAxis f) (FUR UpAxis u) = + let (u',r') = fixOrtho2Left f u in transform (xyzMatrix r' u' (vectorNormalize f)) +orthogonalFrame (FUR UpAxis u) (FUR RightAxis r) = + let (r',f') = fixOrtho2Left u r in transform (xyzMatrix r' (vectorNormalize u) f') +orthogonalFrame (FUR ForwardAxis _) (FUR ForwardAxis _) = + error "orthogonalFrame: two forward vectors" +orthogonalFrame (FUR UpAxis _) (FUR UpAxis _) = + error "orthogonalFrame: two up vectors" +orthogonalFrame (FUR RightAxis _) (FUR RightAxis _) = + error "orthogonalFrame: two right vectors" +orthogonalFrame x y = orthogonalFrame (furCorrect x) (furCorrect y) + +-- | Function to transform down, left, or backward axes into +-- up, right, or forward axes. +furCorrect :: FUR Vector3D -> FUR Vector3D +furCorrect (FUR ForwardAxis f) = FUR ForwardAxis f +furCorrect (FUR UpAxis u) = FUR UpAxis u +furCorrect (FUR RightAxis r) = FUR RightAxis r +furCorrect (FUR DownAxis d) = FUR UpAxis $ vectorScale (-1) d +furCorrect (FUR LeftAxis l) = FUR RightAxis $ vectorScale (-1) l +furCorrect (FUR BackwardAxis b) = FUR ForwardAxis $ vectorScale (-1) b + +-- | Translates and rotates a model to aim at a given position or in a +-- given direction from a given vantage point. This is analogous +-- to camera look-at functions, and could be used, for example, to +-- cause a model of an eyeball to track a particular target. +-- The first parameter is the position of the model. Typically the second +-- parameter will be the position of the target, and the third parameter will +-- @(up \$ Vector3D 0 1 0)@. +modelLookAt :: (AffineTransformable a) => + Point3D -> + FUR (Either Point3D Vector3D) -> + FUR (Either Point3D Vector3D) -> + a -> a +modelLookAt pos primaryish secondaryish = + RSAGL.Math.Affine.translate (vectorToFrom pos origin_point_3d) . + orthogonalFrame primary secondary + where primary = fmap (either (`vectorToFrom` pos) id) primaryish + secondary = fmap (either (`vectorToFrom` pos) id) secondaryish + diff --git a/RSAGL/Math/Ray.hs b/RSAGL/Math/Ray.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Ray.hs @@ -0,0 +1,47 @@ +module RSAGL.Math.Ray + (Ray3D(..), + projectRay, + distanceAlong, + angleFrom, + normalizeRay) + where + +import RSAGL.Math.Vector +import RSAGL.Math.Angle +import RSAGL.Math.Affine +import RSAGL.Math.Types + +-- | Rays with endpoints and vectors. +-- +-- Although a ray is isomorphic to a 'SurfaceVertex3D', it does not have the +-- same behavior. +-- +data Ray3D = Ray3D { ray_endpoint :: Point3D, + ray_vector :: Vector3D } + deriving (Read,Show) + +instance AffineTransformable Ray3D where + transform m (Ray3D p v) = Ray3D (transform m p) (transform m v) + +-- | The parametric function of a ray. The parameter is measured as a +-- proportion of the length of the vector. @projectRay 0@ is the +-- endpoint of the ray. @projectRay 1@ is the endpoint offset +-- by the ray's vector. +projectRay :: RSdouble -> Ray3D -> Point3D +projectRay t (Ray3D (Point3D x y z) (Vector3D u v w)) = Point3D (x+u*t) (y+v*t) (z+w*t) + +-- | The inverse operation to 'projectRay'. This could also be +-- understood as the height of the point above the plane defined +-- by the ray. +distanceAlong :: Ray3D -> Point3D -> RSdouble +distanceAlong (Ray3D p v) p' = dotProduct (vectorToFrom p' p) v / vectorLengthSquared v + +-- | The angle between vector of the ray and the vector from the +-- endpoint of the ray to the specified point. +angleFrom :: Ray3D -> Point3D -> Angle +angleFrom (Ray3D p v) p' = angleBetween (vectorToFrom p' p) v + +-- | A ray normalize to a length of 1. +normalizeRay :: Ray3D -> Ray3D +normalizeRay (Ray3D p v) = Ray3D p $ vectorNormalize v + diff --git a/RSAGL/Math/Types.hs b/RSAGL/Math/Types.hs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Types.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +module RSAGL.Math.Types + (RSfloat, + RSdouble, + toGLfloat, + toGLdouble, + fromGLfloat, + fromGLdouble, + f2f) + where + +import Graphics.Rendering.OpenGL.Raw.Core31 +import Data.Vec as Vec +import Data.Vec.OpenGLRaw +import System.Random +import Control.Arrow + +newtype RSfloat = RSfloat GLfloat + deriving (Enum, + Eq, + Floating, + Fractional, + Num, + Ord, + Read, + Real, + RealFloat, + RealFrac, + Show, + NearZero) +newtype RSdouble = RSdouble GLdouble + deriving (Enum, + Eq, + Floating, + Fractional, + Num, + Ord, + Read, + Real, + RealFloat, + RealFrac, + Show, + NearZero) + +toGLfloat :: RSfloat -> GLfloat +toGLfloat (RSfloat x) = x + +toGLdouble :: RSdouble -> GLdouble +toGLdouble (RSdouble x) = x + +fromGLfloat :: GLfloat -> RSfloat +fromGLfloat = RSfloat + +fromGLdouble :: GLdouble -> RSdouble +fromGLdouble = RSdouble + +{-# RULES +"f2f/id" f2f = id #-} +f2f :: (RealFloat a,RealFloat b) => a -> b +f2f = uncurry encodeFloat . decodeFloat +{-# INLINE f2f #-} + +instance Random RSfloat where + randomR (lo,hi) g = first f2f $ randomR (f2f lo, f2f hi :: Float) g + random g = first (f2f :: Float -> RSfloat) $ random g + +instance Random RSdouble where + randomR (lo,hi) g = first f2f $ randomR (f2f lo, f2f hi :: Double) g + random g = first (f2f :: Double -> RSdouble) $ random g + diff --git a/RSAGL/Math/Vector.lhs b/RSAGL/Math/Vector.lhs new file mode 100644 --- /dev/null +++ b/RSAGL/Math/Vector.lhs @@ -0,0 +1,298 @@ +\section{Points and Vectors: RSAGL.Vector} + +\begin{code} +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies, PatternGuards #-} +module RSAGL.Math.Vector + (Point3D(..), + origin_point_3d, + Vector3D(..), + SurfaceVertex3D(..), + zero_vector, + point3d, + points3d, + point2d, + points2d, + vector3d, + dotProduct, + angleBetween, + crossProduct, + distanceBetween, + distanceBetweenSquared, + aNonZeroVector, + aLargeVector, + vectorAdd, + vectorSum, + vectorScale, + vectorScaleTo, + vectorToFrom, + vectorNormalize, + vectorAverage, + vectorLength, + vectorLengthSquared, + newell, + Xyz(..), + XYZ, + vectorString, + randomXYZ, + fixOrtho, + fixOrtho2, + fixOrtho2Left, + orthos) + where + +import Control.Parallel.Strategies +import RSAGL.Math.Angle +import System.Random +import RSAGL.Math.AbstractVector +import RSAGL.Math.Types +import RSAGL.Math.ListUtils +\end{code} + +\subsection{Generic 3-dimensional types and operations} + +\begin{code} +type XYZ = (RSdouble,RSdouble,RSdouble) + +class Xyz a where + toXYZ :: a -> XYZ + fromXYZ :: XYZ -> a + +instance Xyz (RSdouble,RSdouble,RSdouble) where + toXYZ = id + fromXYZ = id + +vectorString :: Xyz a => a -> String +vectorString xyz = let (x,y,z) = toXYZ xyz + in (show x) ++ "," ++ (show y) ++ "," ++ (show z) + +uncurry3d :: (RSdouble -> RSdouble -> RSdouble -> a) -> XYZ -> a +uncurry3d fn (x,y,z) = fn x y z +\end{code} + +\subsection{Points in 3-space} + +\begin{code} +data Point3D = Point3D {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble + deriving (Read,Show,Eq) + +origin_point_3d :: Point3D +origin_point_3d = Point3D 0 0 0 + +point3d :: (RSdouble,RSdouble,RSdouble) -> Point3D +point3d = uncurry3d Point3D + +point2d :: (RSdouble,RSdouble) -> Point3D +point2d (x,y) = point3d (x,y,0) + +points3d :: [(RSdouble,RSdouble,RSdouble)] -> [Point3D] +points3d = map point3d + +points2d :: [(RSdouble,RSdouble)] -> [Point3D] +points2d = map point2d + +instance Xyz Point3D where + toXYZ (Point3D x y z) = (x,y,z) + fromXYZ (x,y,z) = Point3D x y z + +instance AbstractZero Point3D where + zero = origin_point_3d + +instance AbstractAdd Point3D Vector3D where + add = displace + +instance AbstractSubtract Point3D Vector3D where + sub = vectorToFrom + +instance NFData Point3D +\end{code} + +\subsection{Vectors in 3-space} + +\begin{code} +data Vector3D = Vector3D {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble + deriving (Read,Show,Eq) + +zero_vector :: Vector3D +zero_vector = Vector3D 0 0 0 + +vector3d :: (RSdouble,RSdouble,RSdouble) -> Vector3D +vector3d = uncurry3d Vector3D + +instance Xyz Vector3D where + toXYZ (Vector3D x y z) = (x,y,z) + fromXYZ (x,y,z) = Vector3D x y z + +instance NFData Vector3D + +instance AbstractZero Vector3D where + zero = zero_vector + +instance AbstractAdd Vector3D Vector3D where + add = vectorAdd + +instance AbstractSubtract Vector3D Vector3D where + sub = vectorToFrom + +instance AbstractScale Vector3D where + scalarMultiply = vectorScale + +instance AbstractMagnitude Vector3D where + magnitude = vectorLength + +instance AbstractVector Vector3D +\end{code} + +A \texttt{SurfaceVertex3D} is a a point on an orientable surface, having a position and a normal vector. + +\subsection{Surface Vertices in 3-space} + +\begin{code} +data SurfaceVertex3D = SurfaceVertex3D { sv3d_position :: Point3D, + sv3d_normal :: Vector3D } + deriving (Read,Show) + +instance NFData SurfaceVertex3D where + rnf (SurfaceVertex3D p v) = rnf (p,v) +\end{code} + +\subsection{Vector Arithmetic} + +\begin{code} +aNonZeroVector :: Vector3D -> Maybe Vector3D +aNonZeroVector v = case vectorLength v of + x | x <= 0 -> Nothing + x | isDenormalized x -> Nothing + x | isNaN x -> Nothing + x | isInfinite x -> Nothing + _ | otherwise -> Just v + +aLargeVector :: RSdouble -> Vector3D -> Maybe Vector3D +aLargeVector x v_ = + case aNonZeroVector v_ of + Just v | vectorLength v > x -> Just v + _ | otherwise -> Nothing + +dotProduct :: Vector3D -> Vector3D -> RSdouble +dotProduct (Vector3D ax ay az) (Vector3D bx by bz) = + (ax*bx) + (ay*by) + (az*bz) + +angleBetween :: Vector3D -> Vector3D -> Angle +angleBetween a b = arcCosine $ dotProduct (vectorNormalize a) (vectorNormalize b) + +crossProduct :: Vector3D -> Vector3D -> Vector3D +crossProduct (Vector3D ax ay az) (Vector3D bx by bz) = + Vector3D (ay*bz - az*by) (az*bx - ax*bz) (ax*by - ay*bx) + +distanceBetween :: (Xyz xyz) => xyz -> xyz -> RSdouble +distanceBetween a b = vectorLength $ vectorToFrom a b + +distanceBetweenSquared :: (Xyz xyz) => xyz -> xyz -> RSdouble +distanceBetweenSquared a b = vectorLengthSquared $ vectorToFrom a b + +displace :: (Xyz xyz) => xyz -> Vector3D -> xyz +displace xyz (Vector3D x2 y2 z2) = + let (x1,y1,z1) = toXYZ xyz + in fromXYZ (x1+x2,y1+y2,z1+z2) + +vectorAdd :: Vector3D -> Vector3D -> Vector3D +vectorAdd = displace + +vectorSum :: [Vector3D] -> Vector3D +vectorSum vectors = foldr vectorAdd zero_vector vectors + +vectorToFrom :: (Xyz xyz) => xyz -> xyz -> Vector3D +vectorToFrom a b = + let (ax,ay,az) = toXYZ a + (bx,by,bz) = toXYZ b + in Vector3D (ax - bx) (ay - by) (az - bz) + +vectorLength :: Vector3D -> RSdouble +vectorLength = sqrt . vectorLengthSquared + +vectorLengthSquared :: Vector3D -> RSdouble +vectorLengthSquared (Vector3D x y z) = (x*x + y*y + z*z) + +vectorScale :: RSdouble -> Vector3D -> Vector3D +vectorScale s (Vector3D x y z) = Vector3D (x*s) (y*s) (z*s) +\end{code} + +vectorScaleTo forces the length of a vector to a certain value, without changing the vector's direction. + +\begin{code} +vectorScaleTo :: RSdouble -> Vector3D -> Vector3D +vectorScaleTo new_length vector = vectorScale new_length $ vectorNormalize vector +\end{code} + +vectorNormalize forces the length of a vector to 1, without changing the vector's direction. + +\begin{code} +vectorNormalize :: Vector3D -> Vector3D +vectorNormalize v = + let l = vectorLength v + in maybe (Vector3D 0 0 0) (vectorScale (1/l)) $ aNonZeroVector v +\end{code} + +vectorAverage takes the average of a list of vectors. The result is a normalized vector +where the length of the element vectors is not reflected in the result. + +\begin{code} +vectorAverage :: [Vector3D] -> Vector3D +vectorAverage vects = vectorNormalize $ vectorSum $ map vectorNormalize vects +\end{code} + +\subsection{Generating normal vectors} + +\texttt{newell} calculates the normal vector of an arbitrary polygon. If the points specified are non-coplanar, +\texttt{newell} often calculates a reasonable result; if they are colinear or singular \texttt{newell} will fail. + +The result is a normalized vector. + +\begin{code} +newell :: [Point3D] -> Maybe Vector3D +newell points = fmap vectorNormalize $ aNonZeroVector $ vectorSum $ map newell_ $ loopedDoubles points + where newell_ (Point3D x0 y0 z0,Point3D x1 y1 z1) = + (Vector3D + ((y0 - y1)*(z0 + z1)) + ((z0 - z1)*(x0 + x1)) + ((x0 - x1)*(y0 + y1))) +\end{code} + +\subsection{Randomly Generated Coordinates} + +\texttt{randomXYZ} can generate random coordinates within the cube where x, y, and z are each in the range (lo,hi). + +\begin{code} +randomXYZ :: (RandomGen g,Xyz p) => (RSdouble,RSdouble) -> g -> (p,g) +randomXYZ (lo,hi) g = (fromXYZ (f2f x,f2f y,f2f z),g') + where (g_,g') = split g + (x:y:z:_) = randomRs (f2f lo,f2f hi) g_ :: [Double] +\end{code} + +\subsection{Orthagonal Vectors} + +\texttt{fixOrtho a v} finds the vector, orthogonal to a, that has the least angle to v. + +\texttt{fixOrtho2 right up} yields \texttt{(up,forward)}. + +\texttt{fixOrtho2Left right up} yields \texttt{(up,backward)}. + +\texttt{orthos} finds two arbitrary vectors orthogonal to the parameter. + +\begin{code} +fixOrtho :: Vector3D -> Vector3D -> Vector3D +fixOrtho a = fst . fixOrtho2 a + +fixOrtho2 :: Vector3D -> Vector3D -> (Vector3D,Vector3D) +fixOrtho2 a v = (vectorNormalize $ crossProduct a $ vectorScale (-1) b,vectorNormalize b) + where b = crossProduct a v + +fixOrtho2Left :: Vector3D -> Vector3D -> (Vector3D,Vector3D) +fixOrtho2Left a v = (vectorNormalize $ crossProduct a b,vectorNormalize b) + where b = vectorScale (-1) $ crossProduct a v + +orthos :: Vector3D -> (Vector3D,Vector3D) +orthos v@(Vector3D x y z) | abs y >= abs x && abs z >= abs x = fixOrtho2 v (Vector3D (abs x + abs y + abs z) y z) +orthos v@(Vector3D x y z) | abs x >= abs y && abs z >= abs y = fixOrtho2 v (Vector3D x (abs x + abs y + abs z) z) +orthos v@(Vector3D x y z) | abs x >= abs z && abs y >= abs z = fixOrtho2 v (Vector3D x y (abs x + abs y + abs z)) +orthos v = error $ "orthos: (" ++ show v ++ ")" +\end{code} diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,5 @@ +#!/usr/bin/runhaskell + +import Distribution.Simple + +main = defaultMainWithHooks simpleUserHooks diff --git a/rsagl-math.cabal b/rsagl-math.cabal new file mode 100644 --- /dev/null +++ b/rsagl-math.cabal @@ -0,0 +1,52 @@ +name: rsagl-math +version: 0.6.0.0 +license: BSD3 +license-file: LICENSE +author: Christopher Lane Hinson +maintainer: Christopher Lane Hinson + +category: Math +synopsis: The RogueStar Animation and Graphics Library: Mathematics +description: RSAGL, the RogueStar Animation and Graphics Library, + was specifically designed for a computer game called + roguestar, but effort has been made to make it accessable + to other projects that might benefit from it. This + package includes mathematical algorithms to support + the project. +cabal-version: >= 1.2 +homepage: http://roguestar.downstairspeople.org/ + +build-type: Simple +tested-with: GHC==6.12.1 + +Library + exposed-modules: RSAGL.Math + RSAGL.Math.AbstractVector, + RSAGL.Math.Affine, + RSAGL.Math.Angle, + RSAGL.Math.BoundingBox, + RSAGL.Math.Curve, + RSAGL.Math.CurveExtras, + RSAGL.Math.FMod, + RSAGL.Math.Interpolation, + RSAGL.Math.ListUtils, + RSAGL.Math.Matrix, + RSAGL.Math.Orthogonal, + RSAGL.Math.Ray, + RSAGL.Math.Types + RSAGL.Math.Vector + + ghc-options: -fno-warn-type-defaults -fexcess-precision + ghc-prof-options: -prof -auto-all + + build-depends: base>=4 && <5, + random>= 1.0.0.2 && < 1.1, + array>= 0.3.0.0 && < 0.4, + containers>= 0.3.0.0 && < 0.4, + OpenGL>= 2.4.0.1 && < 2.5, + OpenGLRaw>= 1.1.0.1 && < 1.2, + parsec>=3.1.0 && < 3.2, + parallel>=2.2.0.1 && < 2.3, + Vec>=0.9.8 && < 0.10, + Vec-OpenGLRaw>=0.2.0.0 && < 0.3 +