Enable higher-precision rotation.

1 files changed, 81 insertions, 60 deletions

diff --git a/MeshSketch.hs b/MeshSketch.hs
index 8f38863..273c3dc 100644
--- a/MeshSketch.hs
+++ b/MeshSketch.hs
@@ -1,7 +1,8 @@
-{-# LANGUAGE OverloadedLabels #-}
+{-# LANGUAGE CPP               #-}
+{-# LANGUAGE LambdaCase        #-}
+{-# LANGUAGE OverloadedLabels  #-}
 {-# LANGUAGE OverloadedStrings #-}
-{-# LANGUAGE LambdaCase #-}
-{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE RecordWildCards   #-}
 module MeshSketch where
 
 import Codec.Picture             as Juicy
@@ -16,6 +17,7 @@ import Data.Text (Text)
 import Data.Map.Strict (Map)
 import qualified Data.Map.Strict as Map
 import qualified Data.Vector     as V
+import qualified Data.Vector.Generic as G (init)
 import Foreign.Marshal.Array
 import Foreign.Storable
 import GI.Gdk
@@ -24,7 +26,7 @@ import GI.GLib.Constants
 import GI.Gtk hiding (IsWindow,windowFullscreen,windowUnfullscreen)
 import LambdaCube.GL             as LC
 import LambdaCube.GL.Mesh        as LC
-import Numeric.LinearAlgebra hiding ((<>))
+import Numeric.LinearAlgebra as Math hiding ((<>))
 import System.Environment
 import System.IO
 import System.IO.Error
@@ -86,16 +88,17 @@ data State = State
 
 data Camera = Camera
     { camHeightAngle   :: Float
-    , camTarget        :: Vector Float
-    , camDirection     :: Vector Float
+    , camTarget        :: Vector Float -- 3-vector
+    , camDirection     :: Vector Float -- 3-vector
     , camDistance      :: Float
     , camWidth         :: Float
     , camHeight        :: Float
-    , camUp            :: Vector Float
-    , camWorldToScreen :: Maybe (Matrix Float)
-    , camScreenToWorld :: Maybe (Matrix Float)
+    , camUp            :: Vector Float -- 3-vector
+    , camWorldToScreen :: Maybe (Matrix Float) -- 4×4
+    , camScreenToWorld :: Maybe (Matrix Float) -- 4×4
     }
 
+camPos :: Camera -> Vector Float
 camPos c = camTarget c - scale (camDistance c) (camDirection c)
 
 initCamera :: Camera
@@ -122,6 +125,13 @@ viewProjection c
     pos = camTarget c - scale (camDistance c) (camDirection c)
     proj = perspective 0.1 100 (camHeightAngle c) (camWidth c / camHeight c)
 
+realToFracVector :: ( Real a
+                    , Fractional b
+                    , Storable a
+                    , Storable b
+                    ) => Vector a -> Vector b
+realToFracVector v = Math.fromList $ map realToFrac $ Math.toList v
+
 realToFracMatrix :: (Real a, Fractional t, Element t, Element a) => Matrix a -> Matrix t
 realToFracMatrix m = fromLists $ map realToFrac <$> toLists m
 
@@ -337,12 +347,12 @@ onRender w realized gl = do
 
 onResize :: GLArea -> Realized -> Int32 -> Int32 -> IO ()
 onResize glarea realized w h = do
-    let storage = stStorage realized
     -- Plenty of options here.  I went with the last one.
     -- 1. gLContextGetWindow :: HasCallStack => GLContext -> IO (Maybe Window)
     -- 2. getGLContextWindow ::                 GLContext -> IO (Maybe Window)
     -- 3. widgetGetWindow    :: HasCallStack => GLArea    -> IO (Maybe Window)
-    widgetGetWindow glarea >>= mapM_ (\win -> do
+    mwin <- widgetGetWindow glarea
+    forM_ mwin $ \win -> do
         (wd,ht) <- do wd <- windowGetWidth win
                       ht <- windowGetHeight win
                       return (fromIntegral wd,fromIntegral ht)
@@ -352,11 +362,19 @@ onResize glarea realized w h = do
                             , camWorldToScreen = Nothing
                             , camScreenToWorld = Nothing
                             }
-        LC.setScreenSize (stStorage realized) wd ht)
+        LC.setScreenSize (stStorage realized) wd ht
 
 -- This computes a point in world coordinates on the view screen if
 -- we assume the camera is located at the origin.
 computeDirection :: Camera -> Double -> Double -> Vector Float
+computeDirection cam h k | Just pv <- camScreenToWorld cam  =
+    let d0 = fromList [ 2 * realToFrac h/camWidth cam - 1
+                      , 1 - 2 * realToFrac k/camHeight cam
+                      , 1
+                      , 1
+                      ] :: Vector Float
+        d1 = pv #> d0
+    in scale (1 /(d1!3) ) $ G.init d1
 computeDirection cam h k =
     let d̂ = camDirection cam -- forward
         û = camUp cam        -- upward
@@ -366,7 +384,11 @@ computeDirection cam h k =
         xd = (camHeight cam / 2) / tan (camHeightAngle cam / 2)
     in scale xr r̂ + scale xu û + scale xd d̂
 
-rotate :: Float -> Vector Float -> Matrix Float
+rotate :: ( Floating a
+          , Math.Container Vector a
+          , Indexable (Vector a) a
+          , Normed (Vector a)
+          ) => a -> Vector a -> Matrix a
 rotate cosθ u = (3><3)
   [ cosθ + ux² mcosθ       , (uy.uy)mcosθ - uz sinθ , (ux.uz)mcosθ + uy sinθ
   , (uy.ux)mcosθ + uz sinθ , cosθ + uy² mcosθ       , (uy.uz)mcosθ - ux sinθ
@@ -383,19 +405,27 @@ rotate cosθ u = (3><3)
     uy² = uy . uy
     uz² = uz . uz
 
+updateCameraRotation :: IsWidget a => a -> State -> Double -> Double -> IO ()
 updateCameraRotation w st h k = do
     m <- readIORef (stDragFrom st)
-    forM_ m $ \(df,cam) -> do
+    forM_ m $ \(df0,cam) -> do
         let d̂ = camDirection cam -- forward
             û = camUp cam        -- upward
-            r̂ = d̂ `cross` û      -- rightward
-            -- fr = df `dot` r̂
-            -- fu = df `dot` û
-            -- fd = df `dot` d̂
-            dt = computeDirection cam h k
-            -- tr = dt `dot` r̂
-            -- tu = dt `dot` û
-            -- td = dt `dot` d̂
+            -- r̂ = d̂ `cross` û      -- rightward
+#if 0
+            -- This turned out to be pointless.
+            promote :: Vector Float -> Vector Double
+            promote = realToFracVector
+            demote :: Vector Double -> Vector Float
+            demote = realToFracVector
+#else
+            promote = id
+            demote = id
+            {-# INLINE promote #-}
+            {-# INLINE demote #-}
+#endif
+            df = promote df0
+            dt = promote $ computeDirection cam h k
             cosθ = dot df dt / realToFrac (norm_2 df) / realToFrac (norm_2 dt)
             axis0 = df `cross` dt
             small x = abs x < 0.00001
@@ -404,8 +434,8 @@ updateCameraRotation w st h k = do
                         then fromList [0,1,0]
                         else axis0
             cam' = cam
-                    { camDirection     = d̂ <# rotate cosθ axis
-                    , camUp            = û <# rotate cosθ axis
+                    { camDirection     = demote $ promote d̂ <# rotate cosθ axis
+                    , camUp            = demote $ promote û <# rotate cosθ axis
                     , camWorldToScreen = Nothing
                     , camScreenToWorld = Nothing
                     }
@@ -414,6 +444,7 @@ updateCameraRotation w st h k = do
         forM_ mwin $ \win ->
             windowInvalidateRect win Nothing False
 
+sanitizeCamera :: State -> IO ()
 sanitizeCamera st = do
     modifyIORef (stCamera st) $ \cam ->
         let d = camDirection cam
@@ -435,10 +466,29 @@ sanitizeCamera st = do
             , camScreenToWorld = Nothing
             }
 
+worldCoordinates :: State -> Double -> Double -> IO (Vector Float)
+worldCoordinates st h k = do
+    pv <- atomicModifyIORef' (stCamera st) projectionView
+    cam <- readIORef (stCamera st)
+    let d0 = fromList [ 2 * realToFrac h/camWidth cam - 1
+                      , 1 - 2 * realToFrac k/camHeight cam
+                      , 1
+                      , 1
+                      ] :: Vector Float
+        d1 = pv #> d0
+    return $ scale (1 /(d1!3) ) d1
+
+pushRing :: IsWidget w => w -> State -> Double -> Double -> IO (Vector Float)
+pushRing w st h k = do
+    d <- worldCoordinates st h k
+    Just win <- getWidgetWindow w
+    pushBack (stRingBuffer st) (d!0) (d!1) (d!2)
+    windowInvalidateRect win Nothing False
+    return d
+
 
 onEvent :: IsWidget w => w -> Realized -> Event -> IO Bool
 onEvent w realized ev = do
-    let st = stState realized
     msrc <- eventGetSourceDevice ev
     inputSource <- forM msrc $ \src -> do
         src <- get src #inputSource
@@ -446,7 +496,7 @@ onEvent w realized ev = do
     etype <- get ev #type
     -- putStrLn $ "onEvent! " ++ show (etype,inputSource)
     let put x = putStrLn (show inputSource ++ " " ++ show x)
-    let st = stState realized
+        st = stState realized
     case etype of
 
         EventTypeMotionNotify -> do
@@ -457,19 +507,7 @@ onEvent w realized ev = do
                         mev <- get ev #motion
                         h <- get mev #x
                         k <- get mev #y
-                        -- let d = camPos cam + computeDirection cam h k
-                        pv <- atomicModifyIORef' (stCamera st) projectionView
-                        cam <- readIORef (stCamera st)
-                        let d0 = fromList [ 2 * realToFrac h/camWidth cam - 1
-                                          , 1 - 2 * realToFrac k/camHeight cam
-                                          , 1
-                                          , 1
-                                          ] :: Vector Float
-                            d1 = pv #> d0
-                            d = scale (1 /(d1!3) ) d1
-                        pushBack (stRingBuffer st) (d!0) (d!1) (d!2)
-                        Just win <- getWidgetWindow w
-                        windowInvalidateRect win Nothing False
+                        d <- pushRing w st h k
                         put (etype,(h,k),d)
                 _ -> do
                     mev <- get ev #motion
@@ -486,16 +524,7 @@ onEvent w realized ev = do
                     bev <- get ev #button
                     h <- get bev #x
                     k <- get bev #y
-                    pv <- atomicModifyIORef' (stCamera st) projectionView
-                    cam <- readIORef (stCamera st)
-                    let d0 = fromList [ 2 * realToFrac h/camWidth cam - 1
-                                      , 1 - 2 * realToFrac k/camHeight cam
-                                      , 1
-                                      , 1
-                                      ] :: Vector Float
-                        d1 = pv #> d0
-                        d = scale (1 /(d1!3) ) d1
-                    pushBack (stRingBuffer st) (d!0) (d!1) (d!2)
+                    d <- pushRing w st h k
                     Just win <- getWidgetWindow w
                     windowInvalidateRect win Nothing False
                     put (etype,(h,k),d)
@@ -503,8 +532,9 @@ onEvent w realized ev = do
                     bev <- get ev #button
                     h <- get bev #x
                     k <- get bev #y
+                    d <- G.init <$> worldCoordinates st h k
                     cam <- readIORef (stCamera st)
-                    let d = computeDirection cam h k
+                    -- let d = computeDirection cam h k
                     writeIORef (stDragFrom st) $ Just (d,cam)
                     put (etype,(h,k),d)
                     return ()
@@ -516,16 +546,7 @@ onEvent w realized ev = do
                     bev <- get ev #button
                     h <- get bev #x
                     k <- get bev #y
-                    pv <- atomicModifyIORef' (stCamera st) projectionView
-                    cam <- readIORef (stCamera st)
-                    let d0 = fromList [ 2 * realToFrac h/camWidth cam - 1
-                                      , 1 - 2 * realToFrac k/camHeight cam
-                                      , 1
-                                      , 1
-                                      ] :: Vector Float
-                        d1 = pv #> d0
-                        d = scale (1 /(d1!3) ) d1
-                    pushBack (stRingBuffer st) (d!0) (d!1) (d!2)
+                    d <- pushRing w st h k
                     Just win <- getWidgetWindow w
                     windowInvalidateRect win Nothing False
                 _ -> do