1 files changed, 79 insertions, 0 deletions

diff --git a/Camera.hs b/Camera.hs
new file mode 100644
index 0000000..6bf4dd8
--- /dev/null
+++ b/Camera.hs
@@ -0,0 +1,79 @@
+module Camera where
+
+import qualified Data.Vector.Generic as G
+import Numeric.LinearAlgebra as Math hiding ((<>))
+
+import Matrix
+import LambdaCube.GL.HMatrix ()
+
+data Camera = Camera
+    { camHeightAngle   :: Float
+    , camTarget        :: Vector Float -- 3-vector
+    , camDirection     :: Vector Float -- 3-vector
+    , camDistance      :: Float
+    , camWidth         :: Float
+    , camHeight        :: Float
+    , camUp            :: Vector Float -- 3-vector
+    , camWorldToScreen :: Maybe (Matrix Float) -- 4×4
+    , camScreenToWorld :: Maybe (Matrix Float) -- 4×4
+    }
+
+-- | Compute the height of a pixel at the given 3d point.
+pixelDelta :: Camera -> Vector Float -> Float
+pixelDelta cam x = realToFrac $ frustumHeight eyeToPoint / realToFrac (camHeight cam)
+ where
+    eyeToPoint = norm_2 (x - camPos cam)
+    frustumHeight d = 2 * d * tan (realToFrac $ camHeightAngle cam / 2)
+
+camPos :: Camera -> Vector Float
+camPos c = camTarget c - scale (camDistance c) (camDirection c)
+
+camWorldCoordinates :: Camera -> Double -> Double -> Maybe (Vector Float) -> Vector Float
+camWorldCoordinates cam h k mplane = case mplane of
+        -- Write on the plane.
+        Just plane -> let n̂ = G.init plane
+                          c = plane!3
+                          a = (c - dot p n̂) / dot d̂ n̂
+                       in p + scale a d̂
+
+        -- Write on the camDistance sphere.
+        Nothing    -> p + scale (camDistance cam) d̂
+ where
+    q0 = fromList [ 2 * realToFrac h/camWidth cam - 1
+                  , 1 - 2 * realToFrac k/camHeight cam
+                  , 1
+                  , 1
+                  ] :: Vector Float
+    pv = snd $ projectionView cam
+    q1 = pv #> q0
+    q2 = scale (1 /(q1!3)) $ G.init q1
+    p = camPos cam
+    d = q2 - p
+    d̂ = unit d
+
+projectionView :: Camera -> (Camera,Matrix Float)
+projectionView c
+    | Just m <- camScreenToWorld c = (c,m)
+    | Just w <- camWorldToScreen c = projectionView c{ camScreenToWorld = Just $ invFloat w }
+    | otherwise                    = projectionView $ fst $ viewProjection c
+
+viewProjection :: Camera -> (Camera,(Matrix Float,Vector Float))
+viewProjection c
+    | Just m <- camWorldToScreen c = (c,(m,pos))
+    | otherwise                    = (c { camWorldToScreen = Just m' }, (m',pos))
+ where
+    m' = proj <> cam
+    cam = lookat pos (camTarget c) (camUp c)
+    pos = camPos c
+    proj = perspective 0.1 100 (camHeightAngle c) (camWidth c / camHeight c)
+
+invFloat :: Matrix Float -> Matrix Float
+invFloat m = realToFracMatrix $ inv (realToFracMatrix m :: Matrix Double)
+
+realToFracMatrix :: (Real a, Fractional t, Element t, Element a) => Matrix a -> Matrix t
+realToFracMatrix m = fromLists $ map realToFrac <$> toLists m
+
+unit :: (Linear t c, Fractional t, Normed (c t)) => c t -> c t
+unit v = scale (1/realToFrac (norm_2 v)) v
+
+