geomancy 0.2.6.0 → 0.3.0.0

raw patch · 17 files changed

+1174/−590 lines, 17 filesdep +webcolor-labelsPVP ok

version bump matches the API change (PVP)

Dependencies added: webcolor-labels

API changes (from Hackage documentation)

Files

• ChangeLog.md +25/−0
• README.md +42/−6
• bench/Main.hs +11/−0
• cbits/mat4.c +24/−24
• geomancy.cabal +7/−5
• src/Geomancy/Gl/Funs.hs +1/−0
• src/Geomancy/Interpolate.hs +3/−2
• src/Geomancy/Mat4.hs +163/−355
• src/Geomancy/Mat4/Internal.hs +286/−0
• src/Geomancy/Quaternion.hs +64/−8
• src/Geomancy/Transform.hs +62/−40
• src/Geomancy/UVec3.hs +12/−0
• src/Geomancy/Vec3.hs +16/−2
• src/Geomancy/Vec4.hs +21/−5
• src/Geomancy/Vulkan/Projection.hs +63/−56
• src/Geomancy/Vulkan/View.hs +37/−29
• test/Spec.hs +337/−58

@@ -1,5 +1,30 @@ # Changelog for geomancy +## 0.3.0.0++> `[0]` releases without row/col-major, handedness, and counter/clockwise confusion.++This release brings sanity and consistency around the transform stack.+- `Transform`s have proper *col-major* definitions.+- Operation ordering follows the GLSL convention.+  * The composition order is `(p <> v <> m <> ...)` (from global to local).+    Reverse your transforms to adjust.+  * The application order is `m !* v` (with the operators becoming infixr 5, just under the `<>`).+    It will do `Mᵀ * v` inside for SIMD reasons and to match what GLSL does for its `M * v` ops.+- `Geomancy.Vulkan.Projection` is right-handed, BUT produces *reverse-depth* range ([1; 0], In 3D infinite-Z converges to 0).+  * Replace your depth tests with `OP_GREATER` and clear to `0.0` - get better precision for your migration troubles.+- `Geomancy.Vulkan.View` is *right-handed*, with +Z being forward.+  * The intended up vector is still `vec3 0 (-1) 0` -- +Y down.+    Silly as it sounds, this matches the XY plane of the window with XY plane in front of a "first person" camera.+- Axis rotations (using `rotateQ`) will appear clockwise when looking along the axis.+- Angle rotations follow Tait-Bryan angles (heading/elevation/bank or yaw/pitch/roll) in the y-x-z frame.+  * `rotateZ (time * rate)` will follow the clock hands in 2D scenes and roll in 3D.+  * `rotateX` will follow the sun from sunrise to sunset, pitching UP / increasing elevation.+  * `rotateY` will turn you right, increasing yaw / heading eastwards.+- You're of course free to define your own transforms, just copy the modules and tune to your liking.+  Just make sure that you use matching row/column constructors and the math layer will do the rest, fast.+- Added `webcolor-labels` instances for UVec3/Vec3/Vec4.+ ## 0.2.6.0  * `Geomancy.Gl.Block` extract to `gl-block` package as `Graphics.Gl.Block`.

@@ -6,12 +6,6 @@ * `Mat4` and `Vec4` are `ByteArray#`. * `Mat4`x`Mat4` and `Mat4`x`Vec4` is done with SIMD. -## Matrix layout--CPU-side matrices compose in MVP order, optimized for `mconcat (local1 : local2 : ... : root)` operation.--GPU-side, in GLSL, it is `PVM * v`.- ## The Numbers  Storing a list of 1000 transformations (e.g. rendering instance data):@@ -64,3 +58,45 @@ Add that time to the poking that'll follow.  Sure, it is in the lower microseconds range, but this budget can be used elsewhere.++## Conventions++### Matrix layout++Transforms produced, composed, and applied to mimic the GLSL order (col-major):++- `vec4 vPosOut = P * V * M * vPosIn;`+- `vPosOut = (p <> v <> m) !* vPosIn`++This way you don't have to transpose your transforms or fiddle with layout annotations.++### Projections / Views++`Geomancy.Vulkan.Projection` is using the "reverse-depth" trick that remaps the vulkan default `[0; 1]` range to `[1; 0]`.+This grants extra precision with one less parameter to specify (you only need "near" now), but makes handedness reasoning tricky.+The default depth range the coordinate is left-handed (+X right, +Y down, +Z forward).+But after reversing the depth it has to be paired with a right-handed view function like `Geomancy.Vulkan.View.lookAtRH`.++The intended up vector is still `vec3 0 (-1) 0` -- +Y down.+Silly as it sounds, this matches the XY plane of the window with XY plane in front of a "first person" camera.++### Rotations++Axis rotations (using `rotateQ`) will appear clockwise when looking along the axis.++Angle rotations follow Tait-Bryan angles (heading/elevation/bank or yaw/pitch/roll) in the y-x-z frame.+- `rotateZ (time * rate)` will follow the clock hands in 2D scenes and roll in 3D.+- `rotateX` will follow the sun from sunrise to sunset, increasing elevation / pitching UP.+- `rotateY` will turn you right, increasing yaw/heading eastwards.++Using `Geomancy.Quaternion.intrinsic roll pitch yaw` will make a rotation from the 3 angles in one go.+You can use it to `rotate` a point directly (e.g `vec3 0 0 1` to get a direction vector from Quaternion) or commit to a matrix using `Transform.rotateQ`.++![yaw-pitch-roll](./yaw-pitch-roll.svg)++You're of course free to define your own transforms, just copy the modules and tune to your liking.+Just make sure that you use matching row/column constructors and the math layer will do the rest, fast.++## GLSL-like functions++To further facilitate conversion between the host and shader code `Geomancy.Gl.Funs` provides common functions like `glFract` and `smoothstep`.

@@ -10,6 +10,7 @@  import qualified Geomancy import qualified Geomancy.Mat4+import qualified Geomancy.Quaternion  import qualified Linear import qualified Linear.Matrix@@ -26,6 +27,10 @@   , bgroup "4x4 poke" $ flip map [10, 100, 1000, 10000] \size ->       env (clones size mempty Linear.identity) $         mat4poke size++  , bgroup "quaternion"+      [ bgroup "euler" quatFromEuler+      ]   ]  mat4transpose :: [Benchmark]@@ -73,3 +78,9 @@     pokeLinearTranspose =       Foreign.withArray (map Linear.transpose m44s) \_ptr ->         pure ()++quatFromEuler :: [Benchmark]+quatFromEuler =+  [ bench "fused" $ whnf (\a -> Geomancy.Quaternion.intrinsic a a a) (pi/3)+  , bench "composed" $ whnf (\a -> Geomancy.Quaternion.extrinsic a a a) (pi/3)+  ]

@@ -3,15 +3,15 @@ #include <xmmintrin.h>  void Mat4xMat4_SIMD(float *A, float *B, float *O) {-  __m128 row0 = _mm_load_ps(&B[0]);-  __m128 row1 = _mm_load_ps(&B[4]);-  __m128 row2 = _mm_load_ps(&B[8]);-  __m128 row3 = _mm_load_ps(&B[12]);+  __m128 row0 = _mm_load_ps(&A[0]);+  __m128 row1 = _mm_load_ps(&A[4]);+  __m128 row2 = _mm_load_ps(&A[8]);+  __m128 row3 = _mm_load_ps(&A[12]);   for(int i=0; i<4; i++) {-    __m128 brod0 = _mm_set1_ps(A[4*i + 0]);-    __m128 brod1 = _mm_set1_ps(A[4*i + 1]);-    __m128 brod2 = _mm_set1_ps(A[4*i + 2]);-    __m128 brod3 = _mm_set1_ps(A[4*i + 3]);+    __m128 brod0 = _mm_set1_ps(B[4*i + 0]);+    __m128 brod1 = _mm_set1_ps(B[4*i + 1]);+    __m128 brod2 = _mm_set1_ps(B[4*i + 2]);+    __m128 brod3 = _mm_set1_ps(B[4*i + 3]);     __m128 row = _mm_add_ps(       _mm_add_ps(         _mm_mul_ps(brod0, row0),@@ -86,28 +86,28 @@     float32x4_t C3 = vmovq_n_f32(0);      // Multiply accumulate in 4x1 blocks to output row-    C0 = vfmaq_laneq_f32(C0, B0, A0, 0);-    C0 = vfmaq_laneq_f32(C0, B1, A0, 1);-    C0 = vfmaq_laneq_f32(C0, B2, A0, 2);-    C0 = vfmaq_laneq_f32(C0, B3, A0, 3);+    C0 = vfmaq_laneq_f32(C0, A0, B0, 0);+    C0 = vfmaq_laneq_f32(C0, A1, B0, 1);+    C0 = vfmaq_laneq_f32(C0, A2, B0, 2);+    C0 = vfmaq_laneq_f32(C0, A3, B0, 3);     vst1q_f32(O, C0); -    C1 = vfmaq_laneq_f32(C1, B0, A1, 0);-    C1 = vfmaq_laneq_f32(C1, B1, A1, 1);-    C1 = vfmaq_laneq_f32(C1, B2, A1, 2);-    C1 = vfmaq_laneq_f32(C1, B3, A1, 3);+    C1 = vfmaq_laneq_f32(C1, A0, B1, 0);+    C1 = vfmaq_laneq_f32(C1, A1, B1, 1);+    C1 = vfmaq_laneq_f32(C1, A2, B1, 2);+    C1 = vfmaq_laneq_f32(C1, A3, B1, 3);     vst1q_f32(&O[4], C1); -    C2 = vfmaq_laneq_f32(C2, B0, A2, 0);-    C2 = vfmaq_laneq_f32(C2, B1, A2, 1);-    C2 = vfmaq_laneq_f32(C2, B2, A2, 2);-    C2 = vfmaq_laneq_f32(C2, B3, A2, 3);+    C2 = vfmaq_laneq_f32(C2, A0, B2, 0);+    C2 = vfmaq_laneq_f32(C2, A1, B2, 1);+    C2 = vfmaq_laneq_f32(C2, A2, B2, 2);+    C2 = vfmaq_laneq_f32(C2, A3, B2, 3);     vst1q_f32(&O[8], C2); -    C3 = vfmaq_laneq_f32(C3, B0, A3, 0);-    C3 = vfmaq_laneq_f32(C3, B1, A3, 1);-    C3 = vfmaq_laneq_f32(C3, B2, A3, 2);-    C3 = vfmaq_laneq_f32(C3, B3, A3, 3);+    C3 = vfmaq_laneq_f32(C3, A0, B3, 0);+    C3 = vfmaq_laneq_f32(C3, A1, B3, 1);+    C3 = vfmaq_laneq_f32(C3, A2, B3, 2);+    C3 = vfmaq_laneq_f32(C3, A3, B3, 3);     vst1q_f32(&O[12], C3);  }

@@ -1,18 +1,18 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.35.2.+-- This file has been generated from package.yaml by hpack version 0.38.1. -- -- see: https://github.com/sol/hpack  name:           geomancy-version:        0.2.6.0-synopsis:       Geometry and matrix manipulation+version:        0.3.0.0+synopsis:       Vectors and matrix manipulation description:    Sometimes it is unavoidable you have to do stuff on CPU.                 Let's at least do it faster. category:       Graphics-author:         Alexander Bondarenko+author:         IC Rainbow maintainer:     aenor.realm@gmail.com-copyright:      2021 Alexander Bondarenko+copyright:      2021-2025 IC Rainbow license:        BSD3 license-file:   LICENSE build-type:     Simple@@ -34,6 +34,7 @@       Geomancy.IVec3       Geomancy.IVec4       Geomancy.Mat4+      Geomancy.Mat4.Internal       Geomancy.Point       Geomancy.Quaternion       Geomancy.Swizzle@@ -64,6 +65,7 @@     , mono-traversable     , ptrdiff     , simple-affine-space+    , webcolor-labels   default-language: Haskell2010  test-suite geomancy-test

@@ -6,6 +6,7 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}  module Geomancy.Gl.Funs   ( GlClamp(..)

@@ -1,8 +1,9 @@ {-# LANGUAGE BlockArguments #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE PolyKinds #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}  module Geomancy.Interpolate   ( linear

@@ -1,58 +1,65 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE BlockArguments #-}-{-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeApplications #-} {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UnliftedFFITypes #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -Wno-orphans #-} -- XXX: until there are order-specific types+ -- | General matrix storage and operations.  module Geomancy.Mat4   ( Mat4--  , rowMajor-  , withRowMajor-  , toListRowMajor-  , toListRowMajor2d-  , fromRowMajor2d--  , colMajor-  , withColMajor-  , toListColMajor-  , toListColMajor2d-+    -- * Operations   , identity   , transpose   , inverse+  , det   , pointwise   , zipWith   , matrixProduct   , scalarMultiply   , (!*)++    -- * Construction++    -- ** Col-major+  , colMajor+  , withColMajor+  , toListColMajor+  , toListColMajor2d+  , showColumns++    -- ** Row-major+  , rowMajor+  , withRowMajor+  , toListRowMajor+  , toListRowMajor2d+  , fromRowMajor2d+  , showRows   ) where  import Prelude hiding (zipWith) import GHC.Exts hiding (VecCount(..), toList)+import Geomancy.Mat4.Internal -import Control.DeepSeq (NFData(rnf))-import Foreign (Storable(..))-import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)-import GHC.IO (IO(..)) import System.IO.Unsafe (unsafePerformIO) import Text.Printf (printf)  import qualified Data.Foldable as Foldable-import qualified Data.List as List -import Graphics.Gl.Block (Block(..))-import Geomancy.Vec4 (Vec4(..), unsafeNewVec4)+import Geomancy.Vec4 (Vec4(..), newVec4) -data Mat4 = Mat4 ByteArray#+toListColMajor :: Coercible a Mat4 => a -> [Float]+toListColMajor = toListTrans . coerce +toListColMajor2d :: Coercible a Mat4 => a -> [[Float]]+toListColMajor2d = toList2dTrans . coerce+ {- | Construct 'Mat4' from @row@ notation. -} rowMajor@@ -62,28 +69,17 @@   -> Float -> Float -> Float -> Float   -> Float -> Float -> Float -> Float   -> a-rowMajor = coerce mat4--{- | Reduce 'Mat4' with a function with @row@ notation of arguments.--}-withRowMajor-  :: Coercible a Mat4-  => a-  ->-    ( Float -> Float -> Float -> Float ->-      Float -> Float -> Float -> Float ->-      Float -> Float -> Float -> Float ->-      Float -> Float -> Float -> Float ->-      r-    )-  -> r-withRowMajor m = withMat4 (coerce m)--toListRowMajor :: Coercible a Mat4 => a -> [Float]-toListRowMajor = toList . coerce--toListRowMajor2d :: Coercible a Mat4 => a -> [[Float]]-toListRowMajor2d = toList2d . coerce+rowMajor+  e0 e1 e2 e3+  e4 e5 e6 e7+  e8 e9 eA eB+  eC eD eE eF =+    -- just store+    coerce $ fromMemory+      e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF  {- |   Build a Mat4 from a list-of-lists kind of container@@ -103,15 +99,15 @@ fromRowMajor2d rows =   case Foldable.toList rows of     [r0, r1, r2, r3] ->-      withRow r0 \m00 m01 m02 m03 ->-      withRow r1 \m10 m11 m12 m13 ->-      withRow r2 \m20 m21 m22 m23 ->-      withRow r3 \m30 m31 m32 m33 ->-        Just . coerce $ mat4-          m00 m01 m02 m03-          m10 m11 m12 m13-          m20 m21 m22 m23-          m30 m31 m32 m33+      withRow r0 \e0 e1 e2 e3 ->+      withRow r1 \e4 e5 e6 e7 ->+      withRow r2 \e8 e9 eA eB ->+      withRow r3 \eC eD eE eF ->+        Just . coerce @a $ rowMajor+          e0 e1 e2 e3+          e4 e5 e6 e7+          e8 e9 eA eB+          eC eD eE eF     _ ->       Nothing   where@@ -122,31 +118,10 @@         _ ->           Nothing -{- | Construct a 'Mat4' from @column@ notation.--}-{-# INLINE colMajor #-}-colMajor-  :: Coercible Mat4 a-  => Float -> Float -> Float -> Float-  -> Float -> Float -> Float -> Float-  -> Float -> Float -> Float -> Float-  -> Float -> Float -> Float -> Float-  -> a-colMajor-  m00 m01 m02 m03-  m10 m11 m12 m13-  m20 m21 m22 m23-  m30 m31 m32 m33 =-    coerce $ mat4-      m00 m10 m20 m30-      m01 m11 m21 m31-      m02 m12 m22 m32-      m03 m13 m23 m33--{- | Reduce 'Mat4' with a function with @column@ notation for arguments.+{- | Reduce 'Mat4' with a function with @row@ notation for arguments. -}-{-# INLINE withColMajor #-}-withColMajor+{-# INLINE withRowMajor #-}+withRowMajor   :: Coercible a Mat4   => a   ->@@ -157,70 +132,44 @@       r     )   -> r-withColMajor m f = withMat4 (coerce m)-  \ m00 m01 m02 m03-    m10 m11 m12 m13-    m20 m21 m22 m23-    m30 m31 m32 m33 ->-  f-    m00 m10 m20 m30-    m01 m11 m21 m31-    m02 m12 m22 m32-    m03 m13 m23 m33--toListColMajor :: Coercible a Mat4 => a -> [Float]-toListColMajor = toListTrans . coerce--toListColMajor2d :: Coercible a Mat4 => a -> [[Float]]-toListColMajor2d = toList2dTrans . coerce+withRowMajor m f = withMemory (coerce m)+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  f -- preserved element order+    e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF -{- | Construct 'Mat4' from elements in memory order.+{- | Construct a 'Mat4' from @column@ notation. -}-{-# INLINE mat4 #-}-mat4-  :: Float -> Float -> Float -> Float+{-# INLINE colMajor #-}+colMajor+  :: forall a . Coercible Mat4 a+  => Float -> Float -> Float -> Float   -> Float -> Float -> Float -> Float   -> Float -> Float -> Float -> Float   -> Float -> Float -> Float -> Float-  -> Mat4-mat4-  (F# m00) (F# m01) (F# m02) (F# m03)-  (F# m10) (F# m11) (F# m12) (F# m13)-  (F# m20) (F# m21) (F# m22) (F# m23)-  (F# m30) (F# m31) (F# m32) (F# m33) =-  runRW# \world ->-    let-      !(# world_, arr #) = newAlignedPinnedByteArray# 64# 16# world--      world00 = writeFloatArray# arr 0x0# m00 world_-      world01 = writeFloatArray# arr 0x1# m01 world00-      world02 = writeFloatArray# arr 0x2# m02 world01-      world03 = writeFloatArray# arr 0x3# m03 world02--      world10 = writeFloatArray# arr 0x4# m10 world03-      world11 = writeFloatArray# arr 0x5# m11 world10-      world12 = writeFloatArray# arr 0x6# m12 world11-      world13 = writeFloatArray# arr 0x7# m13 world12--      world20 = writeFloatArray# arr 0x8# m20 world13-      world21 = writeFloatArray# arr 0x9# m21 world20-      world22 = writeFloatArray# arr 0xA# m22 world21-      world23 = writeFloatArray# arr 0xB# m23 world22--      world30 = writeFloatArray# arr 0xC# m30 world23-      world31 = writeFloatArray# arr 0xD# m31 world30-      world32 = writeFloatArray# arr 0xE# m32 world31-      world33 = writeFloatArray# arr 0xF# m33 world32--      !(# _world', arr' #) = unsafeFreezeByteArray# arr world33-    in-      Mat4 arr'+  -> a+colMajor+  e0 e4 e8 eC+  e1 e5 e9 eD+  e2 e6 eA eE+  e3 e7 eB eF =+    coerce $ fromMemory+      e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF -{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.+{- | Reduce 'Mat4' with a function with @column@ notation for arguments. -}-{-# INLINE withMat4 #-}-withMat4-  :: Mat4+{-# INLINE withColMajor #-}+withColMajor+  :: Coercible a Mat4+  => a   ->     ( Float -> Float -> Float -> Float ->       Float -> Float -> Float -> Float ->@@ -229,68 +178,64 @@       r     )   -> r-withMat4 (Mat4 arr) f =-  f-    (F# (indexFloatArray# arr 0x0#))-    (F# (indexFloatArray# arr 0x1#))-    (F# (indexFloatArray# arr 0x2#))-    (F# (indexFloatArray# arr 0x3#))--    (F# (indexFloatArray# arr 0x4#))-    (F# (indexFloatArray# arr 0x5#))-    (F# (indexFloatArray# arr 0x6#))-    (F# (indexFloatArray# arr 0x7#))--    (F# (indexFloatArray# arr 0x8#))-    (F# (indexFloatArray# arr 0x9#))-    (F# (indexFloatArray# arr 0xA#))-    (F# (indexFloatArray# arr 0xB#))--    (F# (indexFloatArray# arr 0xC#))-    (F# (indexFloatArray# arr 0xD#))-    (F# (indexFloatArray# arr 0xE#))-    (F# (indexFloatArray# arr 0xF#))+withColMajor m f = withMemory (coerce m)+  \ e0 e4 e8 eC+    e1 e5 e9 eD+    e2 e6 eA eE+    e3 e7 eB eF ->+  f -- transposed element order+    e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF -{- | @I@, the identity matrix.+toListRowMajor :: Coercible a Mat4 => a -> [Float]+toListRowMajor = toListMemory . coerce -Neutral element of its monoid, so you can use 'mempty'.--}-{-# INLINE identity #-}-identity :: Mat4-identity = mat4-  1 0 0 0-  0 1 0 0-  0 0 1 0-  0 0 0 1+toListRowMajor2d :: Coercible a Mat4 => a -> [[Float]]+toListRowMajor2d = toList2dMemory . coerce -{-# INLINE transpose #-}-transpose :: Mat4 -> Mat4-transpose =-  flip withMat4+det :: forall a . (Coercible Mat4 a, Coercible Mat4 a) => a -> Float+det m =+  withRowMajor (coerce @_ @a m)     \ m00 m01 m02 m03       m10 m11 m12 m13       m20 m21 m22 m23       m30 m31 m32 m33 ->-    mat4-      m00 m10 m20 m30-      m01 m11 m21 m31-      m02 m12 m22 m32-      m03 m13 m23 m33+        let+          s0 = m00 * m11 - m10 * m01+          s1 = m00 * m12 - m10 * m02+          s2 = m00 * m13 - m10 * m03+          s3 = m01 * m12 - m11 * m02+          s4 = m01 * m13 - m11 * m03+          s5 = m02 * m13 - m12 * m03 +          c5 = m22 * m33 - m32 * m23+          c4 = m21 * m33 - m31 * m23+          c3 = m21 * m32 - m31 * m22+          c2 = m20 * m33 - m30 * m23+          c1 = m20 * m32 - m30 * m22+          c0 = m20 * m31 - m30 * m21+        in+            s0 * c5+          - s1 * c4+          + s2 * c3+          + s3 * c2+          - s4 * c1+          + s5 * c0+ {- | Compute an inverse matrix, slowly. -}-inverse :: (Coercible Mat4 a, Coercible Mat4 a) => a -> a+inverse :: forall a . (Coercible Mat4 a, Coercible Mat4 a) => a -> a inverse m =-  coerce $ withMat4 (coerce m)+  coerce @a $ withRowMajor (coerce @_ @a m)     \ m00 m01 m02 m03       m10 m11 m12 m13       m20 m21 m22 m23       m30 m31 m32 m33 ->         let-          invDet = recip det--          det-            = s0 * c5+          invDet = recip $+              s0 * c5             - s1 * c4             + s2 * c3             + s3 * c2@@ -331,143 +276,36 @@           i32 = (-m30 * s3 + m31 * s1 - m32 * s0) * invDet           i33 = ( m20 * s3 - m21 * s1 + m22 * s0) * invDet         in-          mat4+          rowMajor             i00 i01 i02 i03             i10 i11 i12 i13             i20 i21 i22 i23             i30 i31 i32 i33 -pointwise :: Mat4 -> Mat4 -> (Float -> Float -> Float) -> Mat4-pointwise a b f =-  withMat4 a-    \ a00 a01 a02 a03-      a10 a11 a12 a13-      a20 a21 a22 a23-      a30 a31 a32 a33 ->-  withMat4 b-    \ b00 b01 b02 b03-      b10 b11 b12 b13-      b20 b21 b22 b23-      b30 b31 b32 b33 ->-  mat4-    (f a00 b00) (f a01 b01) (f a02 b02) (f a03 b03)-    (f a10 b10) (f a11 b11) (f a12 b12) (f a13 b13)-    (f a20 b20) (f a21 b21) (f a22 b22) (f a23 b23)-    (f a30 b30) (f a31 b31) (f a32 b32) (f a33 b33)--toList :: Mat4 -> [Float]-toList = flip withMat4-    \ a00 a01 a02 a03-      a10 a11 a12 a13-      a20 a21 a22 a23-      a30 a31 a32 a33 ->-    [ a00, a01, a02, a03-    , a10, a11, a12, a13-    , a20, a21, a22, a23-    , a30, a31, a32, a33-    ]--toList2d :: Mat4 -> [[Float]]-toList2d = flip withMat4-    \ a00 a01 a02 a03-      a10 a11 a12 a13-      a20 a21 a22 a23-      a30 a31 a32 a33 ->-    [ [a00, a01, a02, a03]-    , [a10, a11, a12, a13]-    , [a20, a21, a22, a23]-    , [a30, a31, a32, a33]-    ]--toListTrans :: Mat4 -> [Float]-toListTrans = flip withMat4-    \ a00 a01 a02 a03-      a10 a11 a12 a13-      a20 a21 a22 a23-      a30 a31 a32 a33 ->-    [ a00, a10, a20, a30-    , a01, a11, a21, a31-    , a02, a12, a22, a32-    , a03, a13, a23, a33-    ]--toList2dTrans :: Mat4 -> [[Float]]-toList2dTrans = flip withMat4-    \ a00 a01 a02 a03-      a10 a11 a12 a13-      a20 a21 a22 a23-      a30 a31 a32 a33 ->-    [ [a00, a10, a20, a30]-    , [a01, a11, a21, a31]-    , [a02, a12, a22, a32]-    , [a03, a13, a23, a33]-    ]--zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]-zipWith f a b = List.zipWith f (toList a) (toList b)--foreign import ccall unsafe "Mat4xMat4_SIMD" m4m4simd :: Addr# -> Addr# -> Addr# -> IO ()+foreign import ccall unsafe "Mat4xMat4_SIMD" m4m4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO ()  {-# INLINE matrixProduct #-} matrixProduct :: Mat4 -> Mat4 -> Mat4 matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do-  result@(Mat4 m) <- unsafeNewMat4-  m4m4simd-    (byteArrayContents# l)-    (byteArrayContents# r)-    (byteArrayContents# m)+  result@(Mat4 m) <- newMat4+  m4m4simd l r m   pure result -{-# INLINE unsafeNewMat4 #-}-unsafeNewMat4 :: IO Mat4-unsafeNewMat4 =-  IO \world ->-    let-      !(# world_, arr_ #) = newAlignedPinnedByteArray# 64# 16# world-      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_-    in-      (# world, Mat4 arr #)--{-# INLINE scalarMultiply #-}-scalarMultiply :: Float -> Mat4 -> Mat4-scalarMultiply x m =-  withMat4 m-    \ m00 m01 m02 m03-      m10 m11 m12 m13-      m20 m21 m22 m23-      m30 m31 m32 m33 ->-      mat4-        (m00 * x) (m10 * x) (m20 * x) (m30 * x)-        (m01 * x) (m11 * x) (m21 * x) (m31 * x)-        (m02 * x) (m12 * x) (m22 * x) (m32 * x)-        (m03 * x) (m13 * x) (m23 * x) (m33 * x)+foreign import ccall unsafe "Mat4xVec4_SIMD" m4v4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO () -foreign import ccall unsafe "Mat4xVec4_SIMD" m4v4simd :: Addr# -> Addr# -> Addr# -> IO ()+{- | Matrix - column vector multiplication (post) --- | Matrix - column vector multiplication+@+vOut = p <> v <> m !* vIn+@+-} (!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4 (!*) (coerce -> Mat4 m) (Vec4 v) = unsafePerformIO do-  result@(Vec4 o) <- unsafeNewVec4-  m4v4simd-    (byteArrayContents# m)-    (byteArrayContents# v)-    (byteArrayContents# o)+  result@(Vec4 o) <- newVec4+  !() <- m4v4simd m v o   pure result -  -- withVec4 vec \v1 v2 v3 v4 ->-  --   withColMajor mat-  --     \ m11 m12 m13 m14-  --       m21 m22 m23 m24-  --       m31 m32 m33 m34-  --       m41 m42 m43 m44 ->-  --         vec4-  --           (m11 * v1 + m12 * v2 + m13 * v3 + m14 * v4)-  --           (m21 * v1 + m22 * v2 + m23 * v3 + m24 * v4)-  --           (m31 * v1 + m32 * v2 + m33 * v3 + m34 * v4)-  --           (m41 * v1 + m42 * v2 + m43 * v3 + m44 * v4)--instance NFData Mat4 where-  rnf Mat4{} = ()+infixr 5 !*  instance Semigroup Mat4 where   {-# INLINE (<>) #-}@@ -478,60 +316,30 @@   mempty = identity  instance Show Mat4 where-  show = flip withMat4-    \ m00 m01 m02 m03-      m10 m11 m12 m13-      m20 m21 m22 m23-      m30 m31 m32 m33 ->-    unlines-      [ printf "| %.4f %.4f %.4f %.4f |" m00 m01 m02 m03-      , printf "| %.4f %.4f %.4f %.4f |" m10 m11 m12 m13-      , printf "| %.4f %.4f %.4f %.4f |" m20 m21 m22 m23-      , printf "| %.4f %.4f %.4f %.4f |" m30 m31 m32 m33-      ]--instance Storable Mat4 where-  sizeOf _mat4 = 64--  alignment _mat4 = 16--  {-# INLINE poke #-}-  poke (Ptr addr) (Mat4 arr) = IO \world ->-    let-      world' = copyByteArrayToAddr# arr 0# addr 64# world-    in-      (# world', () #)+  show = showColumns -  {-# INLINE peek #-}-  peek (Ptr addr) = IO \world ->-    let-      !(# world0, arr #)  = newAlignedPinnedByteArray# 64# 16# world-      world1              = copyAddrToByteArray# addr arr 0# 64# world0-      !(# world', arr' #) = unsafeFreezeByteArray# arr world1-    in-      (# world', Mat4 arr' #)+showColumns :: Mat4 -> String+showColumns cm = withColMajor cm+  \ e0 e4 e8 eC+    e1 e5 e9 eD+    e2 e6 eA eE+    e3 e7 eB eF ->+  unlines+    [ printf  "/ %.4f %.4f %.4f %.4f \\" e0 e4 e8 eC+    , printf  "| %.4f %.4f %.4f %.4f |"  e1 e5 e9 eD+    , printf  "| %.4f %.4f %.4f %.4f |"  e2 e6 eA eE+    , printf "\\ %.4f %.4f %.4f %.4f /"  e3 e7 eB eF+    ] -instance Block Mat4 where-  type PackedSize Mat4 = 64-  alignment140 _  = 16-  sizeOf140       = sizeOfPacked-  alignment430    = alignment140-  sizeOf430       = sizeOf140-  isStruct _      = False-  read140     = peekDiffOff-  write140    = pokeDiffOff-  read430     = read140-  write430    = write140-  readPacked  = read140-  writePacked = write140-  {-# INLINE alignment140 #-}-  {-# INLINE sizeOf140 #-}-  {-# INLINE alignment430 #-}-  {-# INLINE sizeOf430 #-}-  {-# INLINE isStruct #-}-  {-# INLINE read140 #-}-  {-# INLINE write140 #-}-  {-# INLINE read430 #-}-  {-# INLINE write430 #-}-  {-# INLINE readPacked #-}-  {-# INLINE writePacked #-}+showRows :: Mat4 -> String+showRows rm = withRowMajor rm+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  unlines+    [ printf "[ %.4f %.4f %.4f %.4f |" e0 e1 e2 e3+    , printf "| %.4f %.4f %.4f %.4f |" e4 e5 e6 e7+    , printf "| %.4f %.4f %.4f %.4f |" e8 e9 eA eB+    , printf "| %.4f %.4f %.4f %.4f ]" eC eD eE eF+    ]

@@ -0,0 +1,286 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UnboxedTuples #-}++-- | General matrix storage and operations.++module Geomancy.Mat4.Internal+  ( Mat4(..)+  , fromMemory+  , withMemory+  , newMat4+    -- * Order-independent operations+  , identity+  , transpose+  , pointwise+  , scalarMultiply+    -- * List operations in memory order+  , toListMemory+  , toList2dMemory+  , toListTrans+  , toList2dTrans+  , zipWith+  ) where++import Prelude hiding (zipWith)+import GHC.Exts hiding (VecCount(..), toList)++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..))+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.IO (IO(..))++import qualified Data.List as List++import Graphics.Gl.Block (Block(..))++data Mat4 = Mat4 ByteArray#++instance NFData Mat4 where+  rnf Mat4{} = ()++{- | Construct 'Mat4' from elements in memory order.+-}+{-# INLINE fromMemory #-}+fromMemory+  :: Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Mat4+fromMemory+  (F# e0) (F# e1) (F# e2) (F# e3)+  (F# e4) (F# e5) (F# e6) (F# e7)+  (F# e8) (F# e9) (F# ea) (F# eb)+  (F# ec) (F# ed) (F# ee) (F# ef) =+  runRW# \world ->+    let+      !(# world_, arr #) = newAlignedPinnedByteArray# 64# 16# world++      world0 = writeFloatArray# arr 0x0# e0 world_+      world1 = writeFloatArray# arr 0x1# e1 world0+      world2 = writeFloatArray# arr 0x2# e2 world1+      world3 = writeFloatArray# arr 0x3# e3 world2++      world4 = writeFloatArray# arr 0x4# e4 world3+      world5 = writeFloatArray# arr 0x5# e5 world4+      world6 = writeFloatArray# arr 0x6# e6 world5+      world7 = writeFloatArray# arr 0x7# e7 world6++      world8 = writeFloatArray# arr 0x8# e8 world7+      world9 = writeFloatArray# arr 0x9# e9 world8+      worldA = writeFloatArray# arr 0xA# ea world9+      worldB = writeFloatArray# arr 0xB# eb worldA++      worldC = writeFloatArray# arr 0xC# ec worldB+      worldD = writeFloatArray# arr 0xD# ed worldC+      worldE = writeFloatArray# arr 0xE# ee worldD+      worldF = writeFloatArray# arr 0xF# ef worldE++      !(# _world', arr' #) = unsafeFreezeByteArray# arr worldF+    in+      Mat4 arr'++{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.+-}+{-# INLINE withMemory #-}+withMemory+  :: Mat4+  ->+    ( Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      r+    )+  -> r+withMemory (Mat4 arr) f =+  f+    (F# (indexFloatArray# arr 0x0#))+    (F# (indexFloatArray# arr 0x1#))+    (F# (indexFloatArray# arr 0x2#))+    (F# (indexFloatArray# arr 0x3#))++    (F# (indexFloatArray# arr 0x4#))+    (F# (indexFloatArray# arr 0x5#))+    (F# (indexFloatArray# arr 0x6#))+    (F# (indexFloatArray# arr 0x7#))++    (F# (indexFloatArray# arr 0x8#))+    (F# (indexFloatArray# arr 0x9#))+    (F# (indexFloatArray# arr 0xA#))+    (F# (indexFloatArray# arr 0xB#))++    (F# (indexFloatArray# arr 0xC#))+    (F# (indexFloatArray# arr 0xD#))+    (F# (indexFloatArray# arr 0xE#))+    (F# (indexFloatArray# arr 0xF#))++{- | @I@, the identity matrix.++Neutral element of its monoid, so you can use 'mempty'.+-}+{-# INLINE identity #-}+identity :: Mat4+identity = fromMemory+  1 0 0 0+  0 1 0 0+  0 0 1 0+  0 0 0 1++-- TODO: simdify+{-# INLINE transpose #-}+transpose :: Mat4 -> Mat4+transpose =+  flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    fromMemory+      e0 e4 e8 eC+      e1 e5 e9 eD+      e2 e6 eA eE+      e3 e7 eB eF++pointwise :: Mat4 -> Mat4 -> (Float -> Float -> Float) -> Mat4+pointwise a b f =+  withMemory a+    \ a0 a1 a2 a3+      a4 a5 a6 a7+      a8 a9 aA aB+      aC aD aE aF ->+  withMemory b+    \ b0 b1 b2 b3+      b4 b5 b6 b7+      b8 b9 bA bB+      bC bD bE bF ->+  fromMemory+    (f a0 b0) (f a1 b1) (f a2 b2) (f a3 b3)+    (f a4 b4) (f a5 b5) (f a6 b6) (f a7 b7)+    (f a8 b8) (f a9 b9) (f aA bA) (f aB bB)+    (f aC bC) (f aD bD) (f aE bE) (f aF bF)++{-# INLINE newMat4 #-}+newMat4 :: IO Mat4+newMat4 =+  IO \world ->+    let+      !(# world_, arr_ #) = newAlignedPinnedByteArray# 64# 16# world+      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_+    in+      (# world, Mat4 arr #)++{-# INLINE scalarMultiply #-}+scalarMultiply :: Float -> Mat4 -> Mat4+scalarMultiply x m =+  withMemory m+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+      fromMemory+      (e0 * x) (e1 * x) (e2 * x) (e3 * x)+      (e4 * x) (e5 * x) (e6 * x) (e7 * x)+      (e8 * x) (e9 * x) (eA * x) (eB * x)+      (eC * x) (eD * x) (eE * x) (eF * x)++instance Storable Mat4 where+  sizeOf _mat4 = 64++  alignment _mat4 = 16++  {-# INLINE poke #-}+  poke (Ptr addr) (Mat4 arr) = IO \world ->+    let+      world' = copyByteArrayToAddr# arr 0# addr 64# world+    in+      (# world', () #)++  {-# INLINE peek #-}+  peek (Ptr addr) = IO \world ->+    let+      !(# world0, arr #)  = newAlignedPinnedByteArray# 64# 16# world+      world1              = copyAddrToByteArray# addr arr 0# 64# world0+      !(# world', arr' #) = unsafeFreezeByteArray# arr world1+    in+      (# world', Mat4 arr' #)++instance Block Mat4 where+  type PackedSize Mat4 = 64+  alignment140 _  = 16+  sizeOf140       = sizeOfPacked+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++toListMemory :: Mat4 -> [Float]+toListMemory = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ e0, e1, e2, e3+    , e4, e5, e6, e7+    , e8, e9, eA, eB+    , eC, eD, eE, eF+    ]++toList2dMemory :: Mat4 -> [[Float]]+toList2dMemory = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ [e0, e1, e2, e3]+    , [e4, e5, e6, e7]+    , [e8, e9, eA, eB]+    , [eC, eD, eE, eF]+    ]++toListTrans :: Mat4 -> [Float]+toListTrans = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ e0, e4, e8, eC+    , e1, e5, e9, eD+    , e2, e6, eA, eE+    , e3, e7, eB, eF+    ]++toList2dTrans :: Mat4 -> [[Float]]+toList2dTrans = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ [e0, e4, e8, eC]+    , [e1, e5, e9, eD]+    , [e2, e6, eA, eE]+    , [e3, e7, eB, eF]+    ]++zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]+zipWith f a b = List.zipWith f (toListMemory a) (toListMemory b)

@@ -10,6 +10,8 @@   , withQuaternion    , axisAngle+  , intrinsic+  , extrinsic   , rotate   , rotatePoint   , rotationBetween@@ -149,14 +151,14 @@       (d - h)    {-# INLINE (*) #-}-  Quaternion a b c d * Quaternion e f g h =-    withVec3 v \y z w ->-      Quaternion x y z w+  Quaternion as ax ay az * Quaternion bs bx by bz =+    withVec3 v \x y z ->+      Quaternion s x y z     where-      x = a * e - Vec3.dot v1 v2-      v = Vec3.cross v1 v2 + v2 Vec3.^* a + v1 Vec3.^* e-      v1 = vec3 b c d-      v2 = vec3 f g h+      s = as * bs - Vec3.dot v1 v2+      v = Vec3.cross v1 v2 + v2 Vec3.^* as + v1 Vec3.^* bs+      v1 = vec3 ax ay az+      v2 = vec3 bx by bz    {-# INLINE fromInteger #-}   fromInteger x = Quaternion (fromInteger x) 0 0 0@@ -217,7 +219,10 @@     where       ptr' = castPtr ptr --- | Quaternion construction from axis and angle.+{- | Quaternion construction from axis and angle.++In a right-handed system, the rotation would appear clockwise when looking in the axis direction.+-} {-# INLINE axisAngle #-} axisAngle :: Vec3 -> Float -> Quaternion axisAngle axis rads =@@ -226,6 +231,55 @@   where     half = rads / 2 +{- A composition of roll-pitch-yaw (@Z-X'-Y''@) rotations using local/object's own frame.++Useful for airplane-like controls and first-person camera directions.++The order of application is:+- Roll is applied first, turning around forward and adjusting the "up" and the "right".+- Pitch is applied second, turning around the *new* "right" and adjusting the "up" again.+- Yaw is applied last, turning around the twice-adjusted "up".++Apply the resulting `Quaternion` to the previous local basis to get an updated local basis.+-}+intrinsic :: Float -> Float -> Float -> Quaternion+intrinsic roll pitch yaw =+  -- XXX: fused construction is a bit faster than chaining axisAngles, but less legible.+  -- see Spec.hs for reference+  quaternion+    -- scalar+    (cx * cy * cz + sx * sy * sz)+    -- vector+    (sx * cy * cz + cx * sy * sz)+    (cx * sy * cz - sx * cy * sz)+    (cx * cy * sz - sx * sy * cz)+  where+    cx = cos (pitch * 0.5)+    sx = sin (pitch * 0.5)+    cy = cos (yaw * 0.5)+    sy = sin (yaw * 0.5)+    cz = cos (roll * 0.5)+    sz = sin (roll * 0.5)++{- A composition of rotations using global/parent/world frame.++Useful for kinematics problems. Use intrinsic if unsure.++The order of application is:+- Heading is applied first, turning around the "gravity" axis (west-east rotation).+- Elevation is applied second, adding some down/up rotation wrt. the original right direction.+- Tilt is applied last, adding some down/up rotation wrt. the original forward direction.++Apply the resulting `Quaternion` to the canonical/parent basis to get a local basis.+-}+extrinsic :: Float -> Float -> Float -> Quaternion+extrinsic heading elevation tilt =+  axisAngle (vec3 0 0 1) tilt *+  axisAngle (vec3 1 0 0) elevation *+  axisAngle (vec3 0 1 0) heading++{- | Rotate point with a unit quaternion.+-} {-# INLINE rotate #-} rotate :: Quaternion -> Vec3 -> Vec3 rotate q v = withQuaternion q' \_a b c d -> vec3 b c d@@ -233,6 +287,8 @@     q' = withVec3 v \x y z ->       q * quaternion 0 x y z * conjugate q +{- | Rotate point around another point with a unit quaternion.+-} {-# INLINE rotatePoint #-} rotatePoint :: Quaternion -> Vec3 -> Vec3 -> Vec3 rotatePoint q origin point =

@@ -2,11 +2,12 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE GeneralisedNewtypeDeriving #-}+{-# LANGUAGE ImportQualifiedPost #-} {-# LANGUAGE TypeFamilies #-}  module Geomancy.Transform   ( Transform(..)-  , inverse+  , Mat4.inverse    , apply   , (!.)@@ -25,19 +26,21 @@   , scaleZ   , scaleXY   , scale3+  , scaleV    , dirPos+  , node   ) where  import Foreign (Storable(..)) import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff) -import Geomancy.Mat4 (Mat4, colMajor, inverse)+import Geomancy.Mat4 (Mat4, colMajor) import Geomancy.Quaternion (Quaternion, withQuaternion) import Geomancy.Vec3 (Vec3, vec3, withVec3) import Geomancy.Vec4 (fromVec3, withVec4)+import Geomancy.Mat4 qualified as Mat4 -import qualified Geomancy.Mat4 as Mat4 import Graphics.Gl.Block (Block(..))  newtype Transform = Transform { unTransform :: Mat4 }@@ -72,7 +75,12 @@ apply :: Vec3 -> Transform -> Vec3 apply = flip (!.) --- | Matrix - column vector multiplication with perspective division+{- | Matrix - row vector multiplication with perspective division++@+vOut = pv <> translate !. vIn+@+-} (!.) :: Transform -> Vec3 -> Vec3 (!.) mat vec =   withVec4 res \x y z w ->@@ -80,6 +88,8 @@   where     res = mat Mat4.!* fromVec3 vec 1.0 +infixr 5 !.+ -- ** Translation  {-# INLINE translate #-}@@ -124,55 +134,62 @@ scaleXY :: Float -> Float -> Transform scaleXY x y = scale3 x y 1 --- ** Euler angle rotations+{-# INLINE scaleV #-}+scaleV :: Vec3 -> Transform+scaleV s = withVec3 s scale3 +-- ** Rotation++{- | Clockwise rotation around positive X axis.++Matches @\a -> rotateQ (axisAngle (vec3 1 0 0) a)@.+-} {-# INLINE rotateX #-} rotateX :: Float -> Transform-rotateX rads = colMajor-  1 0   0   0-  0 t11 t21 0-  0 t12 t22 0-  0 0   0   1+rotateX rads =+  colMajor+    1 0   0  0+    0 c (-s) 0+    0 s   c  0+    0 0   0  1   where-    t11 = cost-    t12 = -sint-    t21 = sint-    t22 = cost+    c = cos rads+    s = sin rads -    cost = cos rads-    sint = sin rads+{- | Clockwise rotation around positive Y axis. +Matches @\a -> rotateQ (axisAngle (vec3 0 1 0) a)@.+-} {-# INLINE rotateY #-} rotateY :: Float -> Transform-rotateY rads = colMajor-  t00 0 t20 0-  0   1 0   0-  t02 0 t22 0-  0   0 0   1+rotateY rads =+  colMajor+    c  0 s 0+    0  1 0 0+  (-s) 0 c 0+    0  0 0 1   where-    cost = cos rads-    sint = sin rads+    c = cos rads+    s = sin rads -    t00 = cost-    t02 = sint-    t20 = -sint-    t22 = cost+{- | Clockwise rotation around positive Z axis. +Can be used for 2D rotation in the XY plane.+Matches @\a -> rotateQ (axisAngle (vec3 0 0 1) a)@.++In the right-handed "window coordinates" (e.g. top-left corner is 0,0) "right" becomes "down" after 90deg turn.+-} {-# INLINE rotateZ #-} rotateZ :: Float -> Transform-rotateZ rads = colMajor-  t00 t10 0 0-  t01 t11 0 0-  0   0   1 0-  0   0   0 1+rotateZ rads =+  colMajor+    c (-s) 0 0+    s   c  0 0+    0   0  1 0+    0   0  0 1   where-   t00 = cost-   t01 = -sint-   t10 = sint-   t11 = cost--   cost = cos rads-   sint = sin rads+   c = cos rads+   s = sin rads  {-# INLINE rotateQ #-} rotateQ :: Quaternion -> Transform@@ -198,4 +215,9 @@         (1 - 2 * (y2 + z2)) (    2 * (xy - zw)) (    2 * (xz + yw)) tx         (    2 * (xy + zw)) (1 - 2 * (x2 + z2)) (    2 * (yz - xw)) ty         (    2 * (xz - yw)) (    2 * (yz + xw)) (1 - 2 * (x2 + y2)) tz-         0                   0                   0                  1+        0                   0                   0                    1++node :: Vec3 -> Quaternion -> Vec3 -> Transform+node t r s+ | s == 1.0 = dirPos r t+ | otherwise = dirPos r t <> scaleV s -- TODO: finish derivation

@@ -1,8 +1,12 @@ {-# LANGUAGE BlockArguments #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ViewPatterns #-} @@ -28,6 +32,8 @@ import Foreign (Storable(..)) import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff) import GHC.Ix (Ix(..))+import GHC.OverloadedLabels (IsLabel(..))+import WebColor.Labels (IsWebColor(..))  import Geomancy.Elementwise (Elementwise(..)) import Graphics.Gl.Block (Block(..))@@ -67,6 +73,12 @@ dot :: UVec3 -> UVec3 -> Word32 dot (UVec3 l1 l2 l3) (UVec3 r1 r2 r3) =   l1 * r1 + l2 * r2 + l3 * r3++instance IsWebColor s => IsLabel s UVec3 where+  {-# INLINE fromLabel #-}+  fromLabel =+    webColor @s \r g b ->+      uvec3 (fromIntegral r) (fromIntegral g) (fromIntegral b)  instance NFData UVec3 where   rnf UVec3{} = ()

@@ -1,13 +1,16 @@ {-# LANGUAGE BlockArguments #-} {-# LANGUAGE DataKinds #-}+{-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ImportQualifiedPost #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ViewPatterns #-}-{-# LANGUAGE DerivingStrategies #-}  -- | Specialized and inlined @V3 Float@. @@ -40,9 +43,11 @@ import Data.Coerce (Coercible, coerce) import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..)) import Data.VectorSpace (VectorSpace)+import Data.VectorSpace qualified as VectorSpace import Foreign (Storable(..), castPtr) import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)-import qualified Data.VectorSpace as VectorSpace+import GHC.OverloadedLabels (IsLabel(..))+import WebColor.Labels (IsWebColor(..))  import Geomancy.Elementwise (Elementwise(..)) import Graphics.Gl.Block (Block(..))@@ -85,6 +90,15 @@ {-# INLINE fromTuple #-} fromTuple :: Coercible Vec3 a => (Float, Float, Float) -> a fromTuple (x, y, z) = coerce (vec3 x y z)++instance IsWebColor s => IsLabel s Vec3 where+  {-# INLINE fromLabel #-}+  fromLabel =+    webColor @s \r g b ->+      vec3+        (fromIntegral r / 255)+        (fromIntegral g / 255)+        (fromIntegral b / 255)  instance NFData Vec3 where   rnf Vec3{} = ()

@@ -2,9 +2,13 @@ {-# LANGUAGE BlockArguments #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ImportQualifiedPost #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UnliftedFFITypes #-}@@ -30,7 +34,7 @@   , dot   , normalize -  , unsafeNewVec4+  , newVec4   ) where  import GHC.Exts hiding (VecCount(..), toList)@@ -38,11 +42,13 @@ import Control.DeepSeq (NFData(rnf)) import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..)) import Data.VectorSpace (VectorSpace)+import Data.VectorSpace qualified as VectorSpace import Foreign (Storable(..)) import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff) import GHC.IO (IO(..))+import GHC.OverloadedLabels (IsLabel(..)) import Text.Printf (printf)-import qualified Data.VectorSpace as VectorSpace+import WebColor.Labels (IsWebColorAlpha(..))  import Geomancy.Elementwise (Elementwise(..)) import Graphics.Gl.Block (Block(..))@@ -136,6 +142,16 @@ fromTuple :: (Float, Float, Float, Float) -> Vec4 fromTuple (x, y, z, w) = vec4 x y z w +instance IsWebColorAlpha s => IsLabel s Vec4 where+  {-# INLINE fromLabel #-}+  fromLabel =+    webColorAlpha @s \r g b a ->+      vec4+        (fromIntegral r / 255)+        (fromIntegral g / 255)+        (fromIntegral b / 255)+        (fromIntegral a / 255)+ instance NFData Vec4 where   rnf Vec4{} = () @@ -402,9 +418,9 @@      nearZero a = abs a <= 1e-6 -{-# INLINE unsafeNewVec4 #-}-unsafeNewVec4 :: IO Vec4-unsafeNewVec4 =+{-# INLINE newVec4 #-}+newVec4 :: IO Vec4+newVec4 =   IO \world ->     let       !(# world_, arr_ #) = newAlignedPinnedByteArray# 16# 16# world

@@ -1,72 +1,79 @@ module Geomancy.Vulkan.Projection-  ( perspective-  , infinitePerspective-  , orthoOffCenter+  ( reverseDepthRH+  , reverseDepthOrthoRH+  , orthoRH   ) where  import Geomancy.Mat4 (colMajor) import Geomancy.Transform (Transform(..)) -perspective-  :: Integral side-  => Float+-- | Construct a view-to-NDC transformation.+--+-- This will shove a camera frustum (an expanding pyramid) into a Vulkan "clip space" box+-- with the dimensions [-1; 1] left-to-right, [-1; 1] top-to-bottom, and [1; 0] into the screen.+--+-- That is, things further away in the view will be pulled into [0; 0; 0] point on the *back* of the box.+-- And things that on the near plane set by the argument to this function will be scaled to match+-- the aspect ratio and the field of view.+--+-- When using FoV @pi/2@ and @width@=@height the points with Z=near will keep their positions in the XY plane.+--+-- NOTE: To update your code using a vanilla perspective:+--+-- Change your depth buffer's clear value.+--    Instead of clearing to 1.0 (farthest), you now clear to 0.0.+-- Change your depth comparison function.+--    Instead of VK_COMPARE_OP_LESS, you must use VK_COMPARE_OP_GREATER.+reverseDepthRH+  :: Float+  -> Float   -> Float -> Float-  -> side -> side   -> Transform-perspective fovRads near far width height = colMajor-  x 0   0   0-  0 y   0   0-  0 0   z w23-  0 0 w32   1-+reverseDepthRH fovRads zn width height =+  colMajor+    sx  0 0  0+    0  sy 0  0+    0   0 0 zn+    0   0 1  0   where-    x = cotFoV * aspectX-    y = -cotFoV-    z = far / (near - far)-    w23 = near * far / (near - far)-    w32 = -1--    cotFoV = recip . tan $ 0.5 * fovRads+    sx = sy * recipAspect+      where+        recipAspect = height / width+    sy = 1 / tan (fovRads / 2) -    aspectX = fromIntegral height / fromIntegral width+{- | Vanilla orthographic projection centered on @0,0@ -infinitePerspective-  :: Integral side-  => Float-  -> side-  -> side-  -> Transform-infinitePerspective fovRads width height = colMajor-  x   0   0  0-  0 (-y)  0  0-  0   0 (-1) w-  0   0 (-1) 0+@orthoRH 0 1 2 2@ gives identity transform.+@orthoRH 0 1 800 600@ will map @vec3 400 300 0@ to @vec3 0 0 0@.+-}+orthoRH :: Float -> Float -> Float -> Float -> Transform+orthoRH near far width height =+  colMajor+    sx 0 0 0+    0 sy 0 0+    0  0 z w+    0  0 0 1   where-    (x, y) =-      if width > height then-        ( cotFoV / camAspect-        , cotFoV-        )-      else-        ( cotFoV-        , cotFoV * camAspect-        )-    camAspect = fromIntegral width / fromIntegral height-    cotFoV = recip . tan $ 0.5 * fovRads+    sx = 2 / width+    sy = 2 / height+    z = 1 / (far - near)+    w = near / (near - far) -    w = -2 * near-    near = 1/128 -- 2048+{- | Reverse-depth orthographic projection centered on @0,0@ -orthoOffCenter :: Integral side => Float -> Float -> side -> side -> Transform-orthoOffCenter near far width height = colMajor-  x 0 0 0-  0 y 0 0-  0 0 z w-  0 0 0 1+Can be used in the same render pass with reverseDepthRH/VK_COMPARE_OP_GREATER pipelines. +@reverseDepthOrthoRH 0 1 800 600@ will map @vec3 400 300 0@ to @vec3 0 0 1@.+-}+reverseDepthOrthoRH :: Float -> Float -> Float -> Float -> Transform+reverseDepthOrthoRH near far width height =+  colMajor+    sx 0 0 0+    0 sy 0 0+    0  0 z w+    0  0 0 1   where-    x = 2 / fromIntegral width-    y = 2 / fromIntegral height-    z = 1 / (far - near)--    w = near * (near - far)+    sx = 2 / width+    sy = 2 / height+    z = 1 / (near - far)+    w = far / (far - near)

@@ -1,41 +1,49 @@ {-# LANGUAGE BlockArguments #-}  module Geomancy.Vulkan.View-  ( orthoFitScreen-  , lookAt+  ( lookAtRH+  , lookAtRH_   ) where -import Geomancy.Mat4 (Mat4, rowMajor)-import Geomancy.Vec3 (Vec3, withVec3)+import Geomancy.Mat4 (colMajor) import Geomancy.Transform (Transform(..))+import Geomancy.Vec3 (Vec3, cross, dot, normalize, vec3, withVec3) -import qualified Geomancy.Vec3 as Vec3+{- | Construct a right-handed world-to-view transformation. -lookAt :: Vec3 -> Vec3 -> Vec3 -> Transform-lookAt eye center up =-  withVec3 xa \xaX xaY xaZ ->-  withVec3 ya \yaX yaY yaZ ->-  withVec3 za \zaX zaY zaZ ->-  rowMajor-    xaX yaX (-zaX) 0-    xaY yaY (-zaY) 0-    xaZ yaZ (-zaZ) 0-    xd  yd    zd   1+This will rotate and translate the world around a static camera.+This would produce an identity when looking from @vec3 0 0 0@ to @vec3 0 0 1@ using (-Y up).++The fallback "up" should be used when you're looking close to your UP in either direction.+-}+lookAtRH+  :: Vec3 -- ^ Eye position+  -> Vec3 -- ^ Target+  -> Vec3 -- ^ Up direction+  -> Transform+lookAtRH eye target up =+  withVec3 rgt \rx ry rz ->+  withVec3 up' \ux uy uz ->+  withVec3 fwd \fx fy fz ->+  colMajor+    rx ry rz (-er)+    ux uy uz (-eu)+    fx fy fz (-ef)+     0  0  0    1   where-    xa = Vec3.normalize $ Vec3.cross za up-    ya = Vec3.cross xa za-    za = Vec3.normalize $ center - eye+    fwd = normalize (target - eye)+    rgt = normalize (cross fwd up)+    up' = normalize (cross fwd rgt)+    er = dot eye rgt+    eu = dot eye up'+    ef = dot eye fwd -    xd = - Vec3.dot xa eye-    yd = - Vec3.dot ya eye-    zd =   Vec3.dot za eye+{- | A shortcut for using the -Y as the up axis. -orthoFitScreen :: Float -> Float -> Float -> Float -> Mat4-orthoFitScreen screenWidth screenHeight targetWidth targetHeight =-  rowMajor-    s 0 0 0-    0 s 0 0-    0 0 1 0-    0 0 0 1+Looking forward from the origin does nothing.+-}+{-# INLINE lookAtRH_ #-}+lookAtRH_ :: Vec3 -> Vec3 -> Transform+lookAtRH_ eye target = lookAtRH eye target yNeg   where-    s = min (screenWidth / targetWidth) (screenHeight / targetHeight)+    yNeg = vec3 0 (-1) 0

@@ -1,27 +1,31 @@ {-# LANGUAGE BlockArguments #-}+{-# LANGUAGE ImportQualifiedPost #-} {-# LANGUAGE NumericUnderscores #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE RecordWildCards #-} {-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeApplications #-} +import Geomancy import Hedgehog-import qualified Hedgehog.Gen as Gen-import qualified Hedgehog.Range as Range -import Control.Monad (unless)-import Control.Monad.IO.Class (MonadIO(..))+import Control.Monad (when) import Data.Foldable (toList)-import Data.Maybe (catMaybes) import GHC.Stack (withFrozenCallStack)+import Hedgehog.Gen qualified as Gen+import Hedgehog.Range qualified as Range+import Linear qualified+import Linear.Quaternion qualified import System.Exit (exitFailure, exitSuccess) import Text.Printf (printf) -import Geomancy.Transform (Transform(..))--import qualified Foreign-import qualified Linear-import qualified Geomancy-import qualified Geomancy.Mat4-import qualified Geomancy.Transform as Transform+import Geomancy.Mat4 (withRowMajor, (!*))+import Geomancy.Mat4 qualified+import Geomancy.Quaternion qualified+import Geomancy.Transform ((!.))+import Geomancy.Transform qualified as Transform+import Geomancy.Vulkan.Projection qualified as Projection+import Geomancy.Vulkan.View qualified as View  import Linear ((!*!)) @@ -41,43 +45,286 @@ pattern PROP_TESTS_BRUTAL :: TestLimit pattern PROP_TESTS_BRUTAL = 10_000_000 -prop_assoc_multiply :: Property-prop_assoc_multiply = withTests PROP_TESTS $ property do-  (a, a_) <- forAllTransform-  (b, b_) <- forAllTransform-  (c, c_) <- forAllTransform+prop_mat4_assoc :: Property+prop_mat4_assoc = withTests PROP_TESTS $ property do+  (p, p_) <- forAllTransform+  (v, v_) <- forAllTransform+  (m, m_) <- forAllTransform    let-    ab'c = (a <> b) <> c-    a'bc = a <> (b <> c)-    delta' = nearlyEqualMat4 ab'c a'bc-  annotateShow delta'+    pv'm = (p <> v) <> m+    p'vm = p <> (v <> m)+    delta' = nearlyEqualMat4 pv'm p'vm+  replicate 16 Nothing === delta'    let-    ab_c = (a_ !*! b_) !*! c_-    a_bc = a_ !*! (b_ !*! c_)-    delta_ = nearlyEqualM44 ab_c a_bc-  annotateShow delta_--  -- Intra-library transitivity-  unless (null $ catMaybes delta') do-    -- XXX: check only if there is some outstanding error-    delta' === delta_+    mv_p = (m_ !*! v_) !*! p_+    m_vp = m_ !*! (v_ !*! p_)+    delta_ = nearlyEqualM44 mv_p m_vp+  replicate 16 Nothing === delta_    -- Inter-library calculated values nearlyEqual-  ab'c_ <- toM44 ab'c-  catMaybes (nearlyEqualM44 ab_c ab'c_) === []+  replicate 16 Nothing === nearlyEqualM44 m_vp (toM44 p'vm)+  replicate 16 Nothing === nearlyEqualM44 mv_p (toM44 pv'm) -  a'bc_ <- toM44 a'bc-  catMaybes (nearlyEqualM44 a_bc a'bc_) === []+prop_mat4_order :: Property+prop_mat4_order = withTests 1 $ property do+  let+    cm = Geomancy.Mat4.colMajor @Geomancy.Mat4.Mat4+      0 1 2 3+      4 5 6 7+      8 9 0 1+      2 3 4 5+    rm = Geomancy.Mat4.rowMajor @Geomancy.Mat4.Mat4+      0 4 8 2+      1 5 9 3+      2 6 0 4+      3 7 1 5+  Geomancy.Mat4.toListColMajor2d cm === Geomancy.Mat4.toListColMajor2d rm +  let+    cmT = Geomancy.Mat4.transpose cm+    rmT = Geomancy.Mat4.transpose rm+  Geomancy.Mat4.toListColMajor2d cmT === Geomancy.Mat4.toListColMajor2d rmT++prop_transform_det :: Property+prop_transform_det = withTests 1 $ property do+  1     === Geomancy.Mat4.det Geomancy.Mat4.identity+  1     === Geomancy.Mat4.det (Transform.rotateY (pi/6))  -- area-preserving+  1     === Geomancy.Mat4.det (Transform.translate 1 2 3) -- area-preserving too+  0     === Geomancy.Mat4.det (Transform.scaleX 0)        -- space collapse, no inverse+  -1    === Geomancy.Mat4.det (Transform.scaleZ (-1))     -- flips the space+  8     === Geomancy.Mat4.det (Transform.scale 2)         -- stretches a little+  (1/8) === Geomancy.Mat4.det (Transform.scale (1/2))     -- stretches backwards++prop_projection_reverseDepthRH :: Property+prop_projection_reverseDepthRH = withTests 1 $ property do+  let fov90 = pi / 2+  let zNear = 0.1+  let Transform p = Projection.reverseDepthRH fov90 zNear 800 600+  annotateShow p+  let+    -- stating the expected in transposed order for testing+    Transform expected = Geomancy.Mat4.rowMajor+      0.7500 0.0000 0.0000 0.0000+      0.0000 1.0000 0.0000 0.0000+      0.0000 0.0000 0.0000 1.0000+      0.0000 0.0000  zNear 0.0000+  annotateShow expected+  let delta_ = nearlyEqualMat4 p expected+  replicate 16 Nothing === delta_++  let+    centerNearIn = vec4 0 0 zNear 1+    centerNearOut = p !* centerNearIn+    centerNearNDC = toNDC centerNearOut+  annotateShow (centerNearIn, centerNearOut)+  vec3 0 0 1 === centerNearNDC++prop_projection_orthoRH_identity :: Property+prop_projection_orthoRH_identity = withTests 1 $ property do+  let Transform pNative = Projection.orthoRH 0 1 2 2+  annotateShow pNative+  let delta_ = nearlyEqualMat4 pNative mempty+  replicate 16 Nothing === delta_++prop_projection_reverseDepthOrthoRH :: Property+prop_projection_reverseDepthOrthoRH = withTests 1 $ property do+  let t = Projection.reverseDepthOrthoRH 1 10 800 600+  annotateShow t++  vec4 0 0 1 1 === t !* vec4 0 0 1 1 -- center/near+  vec4 0 0 0 1 === t !* vec4 0 0 10 1 -- center/far+  vec4 1 1 0 1 === t !* vec4 400 300 10 1 -- edge/far+  vec4 (-1/400) (-1/300) 1 1 === t !* vec4 (-1) (-1) 1 1 -- step/near+  vec4 (1/400) (1/300) 0 1 === t !* vec4 1 1 10 1 -- step/far++toNDC :: Vec4 -> Vec3+toNDC v = withVec4 v \x y z w -> vec3 (x/w) (y/w) (z/w)++prop_view_lookAtRH_identity :: Property+prop_view_lookAtRH_identity = property do+  let alreadyThere = View.lookAtRH_ (vec3 0 0 0) (vec3 0 0 1)+  annotateShow alreadyThere++  p <- forAll genVec3++  p === Transform.apply p alreadyThere++xPos :: Vec3+xPos = vec3 1 0 0++xNeg :: Vec3+xNeg = vec3 (-1) 0 0++yPos :: Vec3+yPos = vec3 0 1 0++yNeg :: Vec3+yNeg = vec3 0 (-1) 0++zPos :: Vec3+zPos = vec3 0 0 1++zNeg :: Vec3+zNeg = vec3 0 0 (-1)++prop_transform_rotateX :: Property+prop_transform_rotateX = property do+  let angle = pi / 2+  let t = Transform.rotateX angle+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 1 0 0) angle)+  let yNegT = t !. zPos+  annotateShow yNegT+  let yNegQ = q !. zPos+  annotateShow yNegQ+  replicate 3 Nothing === nearlyEqualVec3 yNegT yNeg+  replicate 3 Nothing === nearlyEqualVec3 yNegQ yNeg++prop_transform_rotateY :: Property+prop_transform_rotateY = property do+  let angle = pi / 2+  let t = Transform.rotateY angle+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 0 1 0) angle)+  let xPosT = t !. zPos+  annotateShow xPosT+  let xPosQ = q !. zPos+  annotateShow xPosQ+  replicate 3 Nothing === nearlyEqualVec3 xPosT xPos+  replicate 3 Nothing === nearlyEqualVec3 xPosQ xPos++prop_transform_rotateZ :: Property+prop_transform_rotateZ = property do+  let angle = pi / 2+  let t = Transform.rotateZ angle+  let q = Transform.rotateQ (Geomancy.Quaternion.axisAngle (vec3 0 0 1) angle)+  let yPosT = t !. xPos+  annotateShow yPosT+  let yPosQ = q !. xPos+  annotateShow yPosQ+  replicate 3 Nothing === nearlyEqualVec3 yPosT yPos+  replicate 3 Nothing === nearlyEqualVec3 yPosQ yPos++prop_quaternion_ref :: Property+prop_quaternion_ref = property do+  annotate "Construction"+  (axis, angle) <- forAll $ (,) <$> genVec3 <*> genAngle+  let qg = Geomancy.Quaternion.axisAngle axis angle+  let ql = Linear.Quaternion.axisAngle (withVec3 axis Linear.V3) angle+  show qg === show (fromLQ ql)++  annotate "Rotation"+  p <- forAll genVec3+  let rpg = Geomancy.Quaternion.rotate qg p+  let rpl = Linear.Quaternion.rotate ql (withVec3 p Linear.V3)+  show rpg === show (fromV3 rpl)++  annotate "Multiplication"+  (axisB, angleB) <- forAll $ (,) <$> genVec3 <*> genAngle+  let qgB = Geomancy.Quaternion.axisAngle axisB angleB+  let qlB = Linear.Quaternion.axisAngle (withVec3 axisB Linear.V3) angleB+  let prodg = qg * qgB+  let prodl = ql * qlB+  show prodg === show (fromLQ prodl)++prop_quaternion_associativity :: Property+prop_quaternion_associativity = property do+  (a, b, c) <- forAll $ (,,)+    <$> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)+    <*> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)+    <*> (Geomancy.Quaternion.axisAngle <$> genVec3 <*> genAngle)+  let ab'c = (a * b) * c+  annotateShow ab'c+  let a'bc = a * (b * c)+  annotateShow a'bc+  replicate 4 Nothing === nearlyEqualQ ab'c a'bc++prop_quaternion_intrinsic :: Property+prop_quaternion_intrinsic = withTests 1 $ property do+  let yaw = pi / 2+  let pitch = pi / 2+  let roll = 0+  let q = Geomancy.Quaternion.intrinsic roll pitch yaw+  let forward' = Geomancy.Quaternion.rotate q zPos+  let right' = Geomancy.Quaternion.rotate q xPos+  let up' = Geomancy.Quaternion.rotate q yNeg+  replicate 3 Nothing === nearlyEqualVec3 yNeg forward'+  replicate 3 Nothing === nearlyEqualVec3 zNeg right'+  replicate 3 Nothing === nearlyEqualVec3 xNeg up'++prop_quaternion_extrinsic :: Property+prop_quaternion_extrinsic = withTests 1 $ property do+  let heading = pi / 2+  let elevation = pi / 2+  let tilt = 0+  let q = Geomancy.Quaternion.extrinsic heading elevation tilt+  let forward' = Geomancy.Quaternion.rotate q zPos+  let right' = Geomancy.Quaternion.rotate q xPos+  let up' = Geomancy.Quaternion.rotate q yNeg+  replicate 3 Nothing === nearlyEqualVec3 xPos forward'+  replicate 3 Nothing === nearlyEqualVec3 yPos right'+  replicate 3 Nothing === nearlyEqualVec3 zNeg up'++prop_quaternion_intrinsic_ref :: Property+prop_quaternion_intrinsic_ref = property do+  (roll, pitch, yaw) <- forAll $ (,,) <$> genAngle <*> genAngle <*> genAngle+  let fused = Geomancy.Quaternion.intrinsic roll pitch yaw+  annotateShow fused+  let composed = intrinsicComposed roll pitch yaw+  annotateShow composed+  let viaFused = Geomancy.Quaternion.rotate fused zPos+  annotateShow viaFused+  let viaComposed = Geomancy.Quaternion.rotate composed zPos+  replicate 3 Nothing === nearlyEqualVec3 viaFused viaComposed++intrinsicComposed :: Float -> Float -> Float -> Quaternion+intrinsicComposed roll pitch yaw =+  Geomancy.Quaternion.axisAngle (vec3 0 1 0) yaw *+  Geomancy.Quaternion.axisAngle (vec3 1 0 0) pitch *+  Geomancy.Quaternion.axisAngle (vec3 0 0 1) roll++prop_transform_node :: Property+prop_transform_node = withTests 1 $ property do+  let t = vec3 2 3 5+  let r = intrinsicComposed 0 0 (pi/7)+  let s = 11+  let dps = Transform.node t r s+  let node = trs t r s+  let v = vec3 0 0 1+  replicate 3 Nothing === nearlyEqualVec3 (dps !. v) (node !. v)++trs :: Vec3 -> Quaternion -> Vec3 -> Transform+trs t r s = mconcat+  [ Transform.translateV t+  , Transform.rotateQ r+  , Transform.scaleV s+  ]++prop_transform_order :: Property+prop_transform_order = withTests 1 $ property do+  let parent = Transform.node (vec3 4 5 6) (intrinsicComposed 0 0 (pi/6)) 2+  let local = Transform.node (vec3 7 8 9) (intrinsicComposed 0 (pi/6) 0) (1/3)+  let precomp = parent <> local !. vec3 0 0 1+  let nocomp = parent !. local !. vec3 0 0 1++  replicate 3 Nothing === nearlyEqualVec3 precomp nocomp+ forAllTransform :: PropertyT IO (Geomancy.Mat4.Mat4, Linear.M44 Float) forAllTransform = withFrozenCallStack do   (_name, Transform g) <- forAllWith fst genTransform-  l <- toM44 g-  pure (g, l)+  pure (g, toM44 g) -genTransform :: Gen ([Char], Transform)+genVec3 :: Gen Vec3+genVec3 =+  vec3+    <$> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+    <*> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+    <*> Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)++genAngle :: Gen Float+genAngle = Gen.float (Range.linearFracFrom 0.0 (-8*pi) (8*pi))++genTransform :: Gen ([Char], Transform) -- TODO: cross-check with Linear genTransform = Gen.choice   [ genIdentity   , genTranslate@@ -91,6 +338,7 @@       x <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)       y <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)       z <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+      when (abs x == 0 && abs y == 0 && abs z == 0) Gen.discard       pure         ( printf "translate %0.4f %0.4f %0.4f" x y z         , Transform.translate x y z@@ -102,7 +350,8 @@         , ("rotate/y", Transform.rotateY)         , ("rotate/z", Transform.rotateZ)         ]-      angle <- Gen.float (Range.linearFracFrom 0.0 (-4 * pi) (4 * pi))+      angle <- genAngle+      when (abs angle == 0) Gen.discard       pure         ( printf "%s %0.4f" name angle         , axis angle@@ -112,40 +361,70 @@       x <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)       y <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)       z <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)+      when (abs x == 1 && abs y == 1 && abs z == 1) Gen.discard       pure         ( printf "scale %0.4f %0.4f %0.4f" x y z         , Transform.scale3 x y z         ) --- toVulkan :: Linear.M44 Float -> Linear.M44 Float--- toVulkan = {- Linear.transpose . -} (correction Linear.!*!)---   where---     correction = Linear.V4---       (Linear.V4 1   0  0   0)---       (Linear.V4 0 (-1) 0   0)---       (Linear.V4 0   0  0.5 0.5)---       (Linear.V4 0   0  0   1)+fromV3 :: Linear.V3 Float -> Vec3+fromV3 (Linear.V3 x y z) = Geomancy.vec3 x y z -toM44 :: MonadIO io => Geomancy.Mat4 -> io (Linear.M44 Float)-toM44 mat4 =-  liftIO $-    Foreign.with mat4 $-      Foreign.peek . Foreign.castPtr+fromLQ :: Linear.Quaternion.Quaternion Float -> Geomancy.Quaternion+fromLQ (Linear.Quaternion.Quaternion s (Linear.V3 x y z)) = Geomancy.quaternion s x y z -type NearlyEqual = Maybe (Float, Float, Float)+toM44 :: Geomancy.Mat4 -> Linear.M44 Float+toM44 mat4 =+  withRowMajor mat4+    \ m00 m10 m20 m30+      m01 m11 m21 m31+      m02 m12 m22 m32+      m03 m13 m23 m33 ->+        Linear.V4+          (Linear.V4 m00 m10 m20 m30)+          (Linear.V4 m01 m11 m21 m31)+          (Linear.V4 m02 m12 m22 m32)+          (Linear.V4 m03 m13 m23 m33) -nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [NearlyEqual]+nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [Maybe Oops] nearlyEqualMat4 a b = Geomancy.Mat4.zipWith nearlyEqual a b -nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [NearlyEqual]+nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [Maybe Oops] nearlyEqualM44 a b = zipWith nearlyEqual (concatMap toList a) (concatMap toList b) -nearlyEqual :: Float -> Float -> NearlyEqual-nearlyEqual x y =-  if x == y || abs (1 - x / y) < 0.001 then+nearlyEqualVec3 :: Vec3 -> Vec3 -> [Maybe Oops]+nearlyEqualVec3 a b =+  withVec3 a \ax ay az ->+  withVec3 b \bx by bz ->+    [ nearlyEqual ax bx+    , nearlyEqual ay by+    , nearlyEqual az bz+    ]++nearlyEqualQ :: Geomancy.Quaternion -> Geomancy.Quaternion -> [Maybe Oops]+nearlyEqualQ a b =+  Geomancy.withQuaternion a \as ax ay az ->+    Geomancy.withQuaternion b \bs bx by bz ->+      [ nearlyEqual as bs+      , nearlyEqual ax bx+      , nearlyEqual ay by+      , nearlyEqual az bz+      ]++nearlyEqual :: Float -> Float -> Maybe Oops+nearlyEqual lhs rhs =+  if lhs == rhs || absDiff < 1e-4 || relDiff < 1e-4 then     Nothing   else-    Just (x, y, abs (1 - x / y))+    Just Oops{..}+  where+    absDiff = abs $ lhs - rhs+    relDiff = absDiff / (abs lhs + abs rhs)+    -- fltEps = 1.19209290e-07+    -- fltMin = 1.175494e-38++data Oops = Oops {lhs :: Float, rhs :: Float, absDiff :: Float, relDiff :: Float}+  deriving (Eq, Show)  discovered :: Group discovered = $$(discover)