{-
   Standard geometric entities and their associated operations
   Sven Panne 2000.   mailto:Sven.Panne@informatik.uni-muenchen.de
-}
module Geometry where

-- To be mathematically correct and make things clearer, we distinguish
-- between vectors and points. This is in contrast to a lot of other
-- graphic libs. All representations are relative to a fixed world frame,
-- see worldFrame below. All components are marked as strict for similar
-- reasons as for the standard Complex library.

---------------------------------------------------------------------------
-- scalars

class Real s => Scalar s

instance Scalar Int
instance Scalar Integer
instance Scalar Float
instance Scalar Double

---------------------------------------------------------------------------
-- Vector space 

data Scalar s => Vector s = Vector !s !s !s
   deriving (Eq, Show)

scaledBy :: Scalar s => Vector s -> s -> Vector s
(Vector x y z) `scaledBy` a = Vector (x*a) (y*a) (z*a)

vAdd :: Scalar s => Vector s -> Vector s -> Vector s
(Vector x1 y1 z1) `vAdd` (Vector x2 y2 z2) = Vector (x1+x2) (y1+y2) (z1+z2)

vSub :: Scalar s => Vector s -> Vector s -> Vector s
(Vector x1 y1 z1) `vSub` (Vector x2 y2 z2) = Vector (x1-x2) (y1-y2) (z1-z2)

---------------------------------------------------------------------------
-- Affine space

data Scalar s => Point s = Point !s !s !s
   deriving (Eq, Show)

pAdd :: Scalar s => Point s -> Vector s -> Point s
(Point x1 y1 z1) `pAdd` (Vector x2 y2 z2) =  Point (x1+x2) (y1+y2) (z1+z2)

pSub :: Scalar s => Point s -> Point s -> Vector s
(Point x1 y1 z1) `pSub` (Point  x2 y2 z2) = Vector (x1-x2) (y1-y2) (z1-z2)

---------------------------------------------------------------------------
-- Euclidean space

dot :: Scalar s => Vector s -> Vector s -> s
(Vector x1 y1 z1) `dot` (Vector x2 y2 z2) = x1*x2 + y1*y2 + z1*z2

cross :: Scalar s => Vector s -> Vector s -> Vector s
(Vector x1 y1 z1) `cross` (Vector x2 y2 z2) =
   Vector (y1*z2 - z1*y2) (z1*x2 - x1*z2) (x1*y2 - y1*x2)

---------------------------------------------------------------------------
-- Angles. NOTE: We use degrees, not radians.

type Angle = Float

addAngle :: Angle -> Angle -> Angle
addAngle a b = until (<360) (+ (-360)) . until (>0) (+360) $ a+b

degreeToRadians :: Angle -> Angle
degreeToRadians = ((pi / 180) *)

radiansToDegree :: Angle -> Angle
radiansToDegree = ((180 / pi) *)

sinD :: Fractional a => Angle -> a
sinD = realToFrac . sin . degreeToRadians

cosD :: Fractional a => Angle -> a
cosD = realToFrac . cos . degreeToRadians

tanD :: Fractional a => Angle -> a
tanD = realToFrac . tan . degreeToRadians

---------------------------------------------------------------------------
-- Frames of reference

data Scalar s => Basis s = Basis {
   basisRight, basisUp, basisBack :: !(Vector s)
 } deriving (Eq, Show)

data Scalar s => Frame s = Frame {
   origin :: !(Point s),
   basis  :: !(Basis s)
 } deriving (Eq, Show)

worldOrigin :: Scalar s => Point s
worldOrigin = Point 0 0 0

worldBasis :: Scalar s => Basis s
worldBasis = Basis {
   basisRight = Vector 1 0 0,
   basisUp    = Vector 0 1 0,
   basisBack  = Vector 0 0 1 }

worldFrame :: Scalar s => Frame s
worldFrame = Frame { origin = worldOrigin, basis = worldBasis }

vWorldToFrame :: Scalar s => Frame s -> Vector s -> Vector s
vWorldToFrame (Frame { basis = b }) vector =
   Vector (vector `dot` basisRight b)
          (vector `dot` basisUp    b)
          (vector `dot` basisBack  b)

-- NOTE: The basis of the frame should always be orthonormal for this
-- function. Under this assumption the inverse of the transformation matrix
-- is its transposition and there is no need for real matrix inversion here.
vFrameToWorld :: Scalar s => Frame s -> Vector s -> Vector s
vFrameToWorld (Frame { basis = b }) vector =
   Vector (vector `dot` Vector rX uX bX)
          (vector `dot` Vector rY uY bY)
          (vector `dot` Vector rZ uZ bZ)
   where Vector rX rY rZ = basisRight b
         Vector uX uY uZ = basisUp    b
         Vector bX bY bZ = basisBack  b

---------------------------------------------------------------------------
-- A plane is represented by a normalized normal and the smallest distance
-- to the origin of the world. This is actually the point/normal
-- representation in disguise.

data Scalar s => Plane s = Plane {
   planeNormal   :: !(Vector s),   -- always normalized
   planeDistance :: !s
 } deriving (Eq, Show)

pointPlaneDistance :: Scalar s => Point s -> Plane s -> s
pointPlaneDistance point plane =
   (point `pSub` worldOrigin) `dot` planeNormal plane

isInFrontOf :: Scalar s => Point s -> Plane s -> Bool
point `isInFrontOf` plane =
   pointPlaneDistance point plane >= planeDistance plane

---------------------------------------------------------------------------
-- An axis aligned bounding box is represented by 2 points in 3D space,
-- the minimum and maxium x,y,z values of all points within the AABB.

data Scalar s => AABBox s = AABBox {
   aABBoxMin, aABBoxMax :: !(Point s)
 } deriving (Eq, Show)

data BBoxPlaneStatus = InFrontOfPlane | IntersectsPlane | OnBackOfPlane
   deriving (Eq, Show)

-- Cunning trick by Greene/Hoff
aABBoxPlaneStatus :: Scalar s => AABBox s -> Plane s -> BBoxPlaneStatus
aABBoxPlaneStatus (AABBox { aABBoxMin = Point minX minY minZ,
                            aABBoxMax = Point maxX maxY maxZ })
                  plane@(Plane { planeNormal = Vector pX pY pZ }) =
   if posFarPt `isInFrontOf` plane
      then if negFarPt `isInFrontOf` plane
              then InFrontOfPlane
              else IntersectsPlane
      else OnBackOfPlane
   where posFarPt = Point (if pX >= 0 then maxX else minX)
                          (if pY >= 0 then maxY else minY)
                          (if pZ >= 0 then maxZ else minZ)
         negFarPt = Point (if pX >= 0 then minX else maxX)
                          (if pY >= 0 then minY else maxY)
                          (if pZ >= 0 then minZ else maxZ)

---------------------------------------------------------------------------
-- A view frustum, represented as a list of planes with inward pointing
-- normals. This can represent volumes different from a real frustum, e.g.
-- we use only right/left/top/bottom planes.

newtype Scalar s => Frustum s = Frustum [Plane s]
   deriving (Eq, Show)

-- Normalized view frustum vectors in view space
frustumNormals :: (Scalar s, Fractional s) => Angle -> Angle -> [Vector s]
frustumNormals fovX fovY = [
   Vector (-cX)   0  (-sX),   -- right
   Vector   cX    0  (-sX),   -- left
   Vector    0 (-cY) (-sY),   -- top
   Vector    0   cY  (-sY) ]  -- bottom
   where cX = cosD (fovX / 2)
         sX = sinD (fovX / 2)
         cY = cosD (fovY / 2)
         sY = sinD (fovY / 2)

calcViewFrustum :: (Scalar s, Fractional s) => Frame s -> Angle -> Angle -> Frustum s
calcViewFrustum frame@(Frame { origin = pos }) fovX fovY = Frustum $
   zipWith (\n d -> Plane { planeNormal = n, planeDistance = d })
           frustumNormalsWS
           distances
   where frustumNormalsWS = map (vFrameToWorld frame) (frustumNormals fovX fovY)
         distances        = map ((pos `pSub` worldOrigin) `dot`) frustumNormalsWS

data BBoxFrustumStatus = InsideFrustum | IntersectsFrustum | OutsideFrustum
   deriving (Eq, Show)

aABBoxFrustumStatus :: Scalar s => AABBox s -> Frustum s -> BBoxFrustumStatus
aABBoxFrustumStatus bbox (Frustum planes) = checkPlanes planes
   where checkPlanes []           = InsideFrustum
         checkPlanes (plane:rest) = case aABBoxPlaneStatus bbox plane of
                                       InFrontOfPlane  -> checkPlanes rest
                                       IntersectsPlane -> IntersectsFrustum
                                       OnBackOfPlane   -> OutsideFrustum
