commit a95f0ac66a2e664a9a10cf0344f1e845340dfcb5 from: James Cook date: Mon Jan 17 20:21:37 2022 UTC Initial import. commit - /dev/null commit + a95f0ac66a2e664a9a10cf0344f1e845340dfcb5 blob - /dev/null blob + 7bcb33375e6cdccdd3220b06a00515c9c77d5b62 (mode 644) --- /dev/null +++ Main.html @@ -0,0 +1,9 @@ + + + + + + + + + blob - /dev/null blob + 1074d6df84092d77c3a70417fc29a24690743f57 (mode 644) --- /dev/null +++ README @@ -0,0 +1,17 @@ +To run: +0. You will need idris2 installed with this change that makes List.length + tail-recursive: https://github.com/idris-lang/Idris2/pull/2100 + As of 2022-01-17, it is not merged. +1. Build the precompute executable: + idris2 --build s3d_precompute.ipkg +2. Generate the S3D.Polyhedra.Precomputed module: + build/exec/s3d_precompute > src/S3D/Polyhedra/Precomputed.idr +3. Build by running + idris2 --cg javascript --build s3d.ipkg +4. Open Main.html in a web browser. + +To run the tests: + idris2 --build s3d_test.ipkg + build/exec/s3d_test + +The src directory is shared by all the *.ipkg files. blob - /dev/null blob + d2e9c1d7c5b6a9eb8babef3e4b1a4008871bbaac (mode 644) --- /dev/null +++ s3d.ipkg @@ -0,0 +1,6 @@ +package s3d + +sourcedir = "src" +main = Main +executable = "Main.js" +depends = contrib, dom, sop blob - /dev/null blob + 24234b6b36188ad0b70d5ea615d0fa3504fdc881 (mode 644) --- /dev/null +++ s3d_precompute.ipkg @@ -0,0 +1,6 @@ +package s3d_precompute + +sourcedir = "src" +main = PrecomputeMain +executable = "s3d_precompute" +depends = contrib blob - /dev/null blob + 96f433a161d9e9a32a83a75c52a27407ddf1228d (mode 644) --- /dev/null +++ s3d_test.ipkg @@ -0,0 +1,6 @@ +package s3d_test + +sourcedir = "src" +main = TestMain +executable = "s3d_test" +depends = contrib blob - /dev/null blob + 9396d24407598ea785748a69100529d2db3ce02c (mode 644) --- /dev/null +++ src/Control/Order/Fold.idr @@ -0,0 +1,35 @@ +module Control.Order.Fold + +import Data.Order.Reversed +import Data.Fin +import Data.Vect + +%default total + +minimum' : Ord a => a -> List a -> a +minimum' x [] = x +minimum' x (y :: t) = minimum' (if y < x then y else x) t + +export +minimum : Ord a => List a -> Maybe a +minimum [] = Nothing +minimum (h :: t) = Just $ minimum' h t + +export +maximum : Ord a => List a -> Maybe a +maximum l = + do MkReversed x <- minimum $ map MkReversed l + pure x + +maxAndIndex : Ord a => Vect (S n) a -> (a, Fin (S n)) +maxAndIndex [x] = (x, 0) +maxAndIndex (x :: xs @ (_ :: _)) = + let (maxRest, indexRest) = maxAndIndex xs + in + if x >= maxRest + then (x, 0) + else (maxRest, FS indexRest) + +export +argmax : Ord a => Vect (S n) a -> Fin (S n) +argmax = snd . maxAndIndex blob - /dev/null blob + 1115c327c3c3f909f98fc04b23645690c685bb2c (mode 644) --- /dev/null +++ src/Data/DefaultMap.idr @@ -0,0 +1,98 @@ +{- + +A map with a default value. A key not having a value is the same as a key +having the default value. + +-} + +module Data.DefaultMap + +import Data.List +import Data.SortedMap +import Decidable.Equality + +%default total + +data AnythingBut : {t : Type} -> t -> Type where + MkAnythingBut : (value : t) -> + {0 notForbidden : Not (value = forbidden)} -> + AnythingBut forbidden + +Eq t => (forbidden : t) => Eq (AnythingBut {t} forbidden) where + MkAnythingBut v0 == MkAnythingBut v1 = v0 == v1 + +Show t => Show (AnythingBut {t} x) where + show (MkAnythingBut value) = show value + +maybeAnythingBut : DecEq a => (forbidden : a) -> a -> Maybe (AnythingBut forbidden) +maybeAnythingBut forbidden x = case decEq x forbidden of + Yes _ => Nothing + No p => Just (MkAnythingBut x {notForbidden=p}) + +export +data DefaultMap : (key : Type) -> (value : Type) -> (defaultValue : v) -> Type where + MkDefaultMap : SortedMap key (AnythingBut {t=value} defaultValue) + -> DefaultMap key value defaultValue + +export +Eq key => Eq value => (defaultValue : value) => Eq (DefaultMap key value defaultValue) where + MkDefaultMap m0 == MkDefaultMap m1 = m0 == m1 + +export +Show key => Show value => (defaultValue : value) => Show (DefaultMap key value defaultValue) where + show (MkDefaultMap m) = show m + +export +empty : Ord k => {0 defaultValue : v} -> DefaultMap k v defaultValue +empty = MkDefaultMap empty + +||| Returns the (key, value) pairs where `value` is not `d`. +export +toList : DefaultMap k v d -> List (k, v) +toList (MkDefaultMap m) = + [(k, value) | (k, MkAnythingBut value) <- Data.SortedMap.toList m] + +export +fromList : Ord k => DecEq v => {d: v} -> List (k, v) -> DefaultMap k v d +fromList l = MkDefaultMap $ fromList $ catMaybes $ map (\(x, y) => map (MkPair x) (maybeAnythingBut _ y)) l + +export +lookup : {defaultValue : v} -> k -> DefaultMap k v defaultValue -> v +lookup x (MkDefaultMap m) = + case Data.SortedMap.lookup x m of + Nothing => defaultValue + Just (MkAnythingBut y) => y + +export +change : DecEq v => + {defaultValue : v} -> + k -> v -> DefaultMap k v defaultValue -> + DefaultMap k v defaultValue +change k v (MkDefaultMap m) = + case decEq v defaultValue of + Yes _ => MkDefaultMap (delete k m) + No p => MkDefaultMap (insert k (MkAnythingBut v {notForbidden=p}) m) + +export +addToValue : Semigroup v => DecEq v => + {d : v} -> + k -> v -> DefaultMap k v d -> DefaultMap k v d +addToValue x y m = change x (y <+> lookup x m) m + +||| Returns the least key with a nondefault value, and that value. +export +leftMost : DefaultMap k v d -> Maybe (k, v) +leftMost (MkDefaultMap m) = + map (\(x, MkAnythingBut y) => (x, y)) (leftMost m) + +export +Monoid v => DecEq v => (d : v) => (d = Prelude.Interfaces.neutral) => Semigroup (DefaultMap k v d) where + x <+> (MkDefaultMap y) = + foldl + (\ acc, (k, MkAnythingBut v) => addToValue k v acc) + x + (Data.SortedMap.toList y) + +export +Ord k => Monoid v => DecEq v => (d : v) => (d = Prelude.Interfaces.neutral) => Monoid (DefaultMap k v d) where + neutral = empty blob - /dev/null blob + ff698e6e304920895663fc698f14e9f0f42c6240 (mode 644) --- /dev/null +++ src/Data/Linear/SizedArray/Prim.idr @@ -0,0 +1,36 @@ +-- Backend for Data.Linear.SizedArray. The main goal is to make things fast by +-- avoiding overhead, e.g. converting numeric types. +-- +-- Caveats: +-- * Does no bounds checking. +-- * Doesn't bother converting numeric types, so may break if numbers are too big. +-- * Breaks referential transparency: an ArrayData value might change; results +-- may be different depending on when the compiler chooses to evaluate them. + +module Data.Linear.SizedArray.Prim + +import Data.Fin + +%default total + +export +data ArrayData : Type where + +export +%foreign + "javascript:lambda:(size) => Array(size)" + "scheme:make-vector" +prim__newArray : Nat -> ArrayData + +export +%foreign + "javascript:lambda:(_0, _1, array, i) => array[i]" + "scheme:(lambda (_0 _1 array i) (vector-ref array i))" +prim__read : ArrayData -> Fin n -> t + +-- The input is modified. +export +%foreign + "javascript:lambda:(_0, _1, array, i, x) => { array[i] = x; return array; }" + "scheme:(lambda (_0 _1 array i x) (vector-set! array i x) array)" +prim__write : ArrayData -> Fin n -> t -> ArrayData blob - /dev/null blob + 0c3884c73f88280d67dd12a304fe4d14363c371f (mode 644) --- /dev/null +++ src/Data/Linear/SizedArray.idr @@ -0,0 +1,85 @@ +-- Linear arrays with size determined by the type. +-- +-- For linear types, using an idea from +-- https://github.com/idris-lang/Idris2/issues/613#issuecomment-960392537 + +module Data.Linear.SizedArray + +import Data.Fin +import Data.Linear.SizedArray.Prim +import Data.Vect + +%default total + +export +data SizedIArray : Nat -> Type -> Type where + MkSizedIArray : ArrayData -> SizedIArray n t + +export +read : SizedIArray n t -> Fin n -> t +read (MkSizedIArray arrayData) i = prim__read arrayData i + +export +data SizedLinArray : Nat -> Type -> Type where + MkSizedLinArray : ArrayData -> SizedLinArray n t + +public export +LinearResult : Type -> Type +LinearResult t = {a : Type} -> (1 _ : (1 _ : t) -> Res a (const t)) -> a + +consume : (1 _ : a) -> b -> b +consume x y = assert_linear (const y) x + +resFst : (1 _ : Res a (const (SizedLinArray n t))) -> a +resFst (x # array) = consume array x + +export +newSizedArray : {n : Nat} -> LinearResult (SizedLinArray n t) +newSizedArray f = resFst $ f $ MkSizedLinArray $ prim__newArray n + +export +newSizedArrayFreeze : {n : Nat} -> (1 _ : SizedLinArray n t -> SizedLinArray n t) -> + SizedIArray n t +newSizedArrayFreeze {n} f = + let MkSizedLinArray arrayData = f $ MkSizedLinArray $ prim__newArray n in + MkSizedIArray arrayData + +export +mread : (1 _ : SizedLinArray n t) -> Fin n -> Res t (const (SizedLinArray n t)) +mread = + assert_linear $ + \array@(MkSizedLinArray arrayData), i => + let value : t + value = prim__read arrayData i + in value # array + +export +write : (1 _ : SizedLinArray n t) -> Fin n -> t -> SizedLinArray n t +write (MkSizedLinArray array) i x = MkSizedLinArray $ prim__write array i x + +export +fillFromVect : {n : Nat} -> Vect n t -> (1 _ : SizedLinArray n t) -> SizedLinArray n t +fillFromVect vect array = + let writeElements : (1 array : SizedLinArray n t) -> List (Fin n, t) -> SizedLinArray n t + writeElements array [] = array + writeElements array ((hi, hx) :: t) = + writeElements (write array hi hx) t + in writeElements array $ toList $ zip range vect + +export +fromVect : {n : Nat} -> Vect n t -> LinearResult (SizedLinArray n t) +fromVect vect f = newSizedArray $ \array => f (fillFromVect vect array) + +namespace IArray + export + fromVect : {n : Nat} -> Vect n t -> SizedIArray n t + fromVect vect = newSizedArrayFreeze (\array => fillFromVect vect array) + +||| This should be safe as long as the resulting value of type `a` doesn't keep +||| a reference to the input. +export +unsafeFreeze : (1 _ : SizedLinArray n t) -> (SizedIArray n t -> a) -> Res a (const (SizedLinArray n t)) +unsafeFreeze = + assert_linear $ + \array@(MkSizedLinArray arrayData), f => + f (MkSizedIArray arrayData) # array blob - /dev/null blob + 40142e74adc7f6aebb82c76e921e233503c00710 (mode 644) --- /dev/null +++ src/Data/Order/Reversed.idr @@ -0,0 +1,18 @@ +module Data.Order.Reversed + +%default total + +public export +data Reversed a = MkReversed a + +public export +Eq a => Eq (Reversed a) where + MkReversed x == MkReversed y = x == y + +export +Ord a => Ord (Reversed a) where + compare (MkReversed x) (MkReversed y) = + case compare x y of + LT => GT + EQ => EQ + GT => LT blob - /dev/null blob + f7f80dc46b1d2d5ef549bf751427292b6c2be85f (mode 644) --- /dev/null +++ src/Main.idr @@ -0,0 +1,115 @@ +module Main + +import Data.IORef +import Data.Vect +import JS +import S3D.AnimationFrame +import S3D.ArrayBuffer +import S3D.Draw +import S3D.Figure +import S3D.GLProgram +import S3D.GLProgram.ColourTexture +import S3D.GLProgram.SimpleColour +import S3D.Drawable +import Math.LinearAlgebra +import S3D.Obstacles +import S3D.PhysicalObject +import S3D.PhysicalState +import S3D.Scenes.House +import S3D.Transformable +import S3D.WebUI +import System.Clock +import Web.Dom +import Web.Html +import Web.Raw.Html +import Web.Raw.Webgl + +%default total + +-- Everything the main loop needs to operate. +record Context where + constructor MkContext + ui : UIContext + drawables : List Drawable + obstacles : Obstacles + +-- The compiler seems to have trouble if I just directly say makeProgram. + +makeColourTextureProgram : WebGL2RenderingContext -> JSIO S3D.GLProgram.ColourTexture.Program +makeColourTextureProgram = makeProgram + +makeSimpleColourProgram : WebGL2RenderingContext -> JSIO S3D.GLProgram.SimpleColour.Program +makeSimpleColourProgram = makeProgram + +theTexturePixels : List Bits8 +theTexturePixels = + concat + [ (if mod u 2 == mod v 2 then [255, 255, 255, 255] else [128, 128, 128, 255]) + | u <- [0 .. 7] + , v <- [0 .. 7] + ] + +makeTheTexture : WebGL2RenderingContext -> Bits32 -> JSIO () +makeTheTexture gl textureUnit = + do texture <- unMaybe "createTexture" $ createTexture gl + activeTexture gl textureUnit + bindTexture gl TEXTURE_2D (Just texture) + texParameteri gl TEXTURE_2D TEXTURE_WRAP_S (cast CLAMP_TO_EDGE) + texParameteri gl TEXTURE_2D TEXTURE_WRAP_T (cast CLAMP_TO_EDGE) + texParameteri gl TEXTURE_2D TEXTURE_MIN_FILTER (cast NEAREST) + texParameteri gl TEXTURE_2D TEXTURE_MAG_FILTER (cast NEAREST) + pixels <- makeUInt8Array !(fromListIO theTexturePixels) + let detailLevel = 0 + let internalFormat = RGBA + let sourceFormat = RGBA + let imageDataType = UNSIGNED_BYTE + let border = 0 + texImage2D gl TEXTURE_2D detailLevel (cast internalFormat) + 8 8 border sourceFormat imageDataType + (Just (S (S (S (Z pixels))))) + +swapXZ : Matrix 4 4 Double +swapXZ = + [ [0, 0, 1, 0] + , [0, 1, 0, 0] + , [1, 0, 0, 0] + , [0, 0, 0, 1] + ] + +covering +mainLoop : Context -> PhysicalState -> Clock Monotonic -> JSIO () +mainLoop context physicalState time = + do let gl = context.ui.gl + (physicalState', viewMatrix) <- + stepState context.obstacles time context.ui.state physicalState + drawUniverse gl viewMatrix context.drawables + requestAnimationFrame (runJS . mainLoop context physicalState') + +partial +mainAfterLoad : JSIO () +mainAfterLoad = + do ui <- S3D.WebUI.init + let gl = ui.gl + enable gl DEPTH_TEST + colourTextureProgram <- makeColourTextureProgram gl + simpleColourProgram <- makeSimpleColourProgram gl + scene <- theHouse + drawables <- compileFigure gl colourTextureProgram simpleColourProgram + scene.figure + + makeTheTexture gl TEXTURE0 + useProgram gl (Just colourTextureProgram.program) + uniform1i gl (Just colourTextureProgram.textureLocation) 0 + + let context = + MkContext + { ui = ui + , drawables = drawables + , obstacles = scene.obstacles + } + requestAnimationFrame (runJS . mainLoop context initialState) + +partial +main : IO () +main = runJS $ + do Web.Raw.Html.GlobalEventHandlers.onload !window ?> mainAfterLoad blob - /dev/null blob + d1c3c542b63ed2053d487da3cf1536dd77407c8b (mode 644) --- /dev/null +++ src/Math/LinearAlgebra/Normalized.idr @@ -0,0 +1,22 @@ +module Math.LinearAlgebra.Normalized + +import Math.LinearAlgebra + +%default total + +export +data Normalized vector = MkNormalized vector + +||| Normalize a vector. +export +makeNormalized : DotProduct vector Double => vector -> Normalized vector +makeNormalized v = MkNormalized (normalize v) + +||| Convert to the original vector type. +export +fromNormalized : Normalized vector -> vector +fromNormalized (MkNormalized x) = x + +export +Show vector => Show (Normalized vector) where + show = show . fromNormalized blob - /dev/null blob + e22e65dad2127790b99581c5c31f3b5e551f733d (mode 644) --- /dev/null +++ src/Math/LinearAlgebra/RotateAxes.idr @@ -0,0 +1,25 @@ +module Math.LinearAlgebra.RotateAxes + +import Data.Fin +import Data.Vect +import Math.LinearAlgebra + +%default total + +basisVector : Num coord => {n : Nat} -> (i : Fin n) -> Vector n coord +basisVector {n = S n'} i = vectToVector (insertAt i 1 (replicate n' 0)) + +-- Helper for rotateAxes. +rotateAxesRow : (Neg coord, Num coord) => {n : Nat} -> (a : Fin n) -> (b : Fin n) -> (r : Fin n) -> Vector n coord +rotateAxesRow a b r = + if r == a + then basisVector b + else if r == b + then negateVector (basisVector a) + else basisVector r + +||| Rotate a quarter circle, changing the first axis to the second. +export +rotateAxes : (Neg coord, Num coord) => {n : Nat} -> (a : Fin n) -> (b : Fin n) -> Matrix n n coord +rotateAxes {n = n@(S n')} a b = + vectorsToMatrix (map (rotateAxesRow a b) range) blob - /dev/null blob + 11fa42c9737ae09735607926fbeb1bd4db37e0a7 (mode 644) --- /dev/null +++ src/Math/LinearAlgebra.idr @@ -0,0 +1,460 @@ +{- + +Operator names (!*! etc) are based on the "linear" package for Haskell. + +-} + +module Math.LinearAlgebra + +import Control.Order.Fold +import Data.Vect + +%default total + +infixl 8 ^+^ +infixl 8 ^-^ +infixl 8 !+! + +infixl 9 ^*! +infixl 9 !*^ +infixl 9 !*! +infix 9 *^ +infix 9 *! +infix 9 ^/ +infix 9 ^.^ + +paritySign : (Neg a, Num a) => Fin n -> a +paritySign FZ = 1 +paritySign (FS x) = - paritySign x + +namespace Vector + + ||| Vector spaces. More accurately, modules, since `coord` might be a ring + ||| and not a field: we don't require `Fractional coord`. + public export + interface Neg coord => VectorSpace vector coord | vector where + negateVector : vector -> vector + ||| Vector sum + (^+^) : vector -> vector -> vector + ||| Vector difference + (^-^) : vector -> vector -> vector + ||| Scalar times vector + (*^) : coord -> vector -> vector + + ||| Vector divided by scalar + export + (^/) : Fractional coord => Neg coord => VectorSpace vector coord => vector -> coord -> vector + v ^/ x = (1 / x) *^ v + + public export + interface VectorSpace vector coord => DotProduct vector coord | vector where + (^.^) : vector -> vector -> coord + + export + Neg coord => VectorSpace (Vect n coord) coord where + negateVector = map negate + (^+^) = zipWith (+) + (^-^) = zipWith (-) + (*^) x = map (x *) + + export + Neg coord => DotProduct (Vect n coord) coord where + x ^.^ y = sum (zipWith (*) x y) + + export + norm : DotProduct vector Double => vector -> Double + norm v = sqrt (v ^.^ v) + + export + normalize : DotProduct vector Double => vector -> vector + normalize v = v ^/ norm v + + ||| Deprecated explicit vector type. + export + data Vector : Nat -> Type -> Type where + MkVector : Vect n a -> Vector n a + + export + Functor (Vector n) where + map f (MkVector v) = MkVector (map f v) + + export + vectToVector : Vect n a -> Vector n a + vectToVector = MkVector + + export + vectorToVect : Vector n a -> Vect n a + vectorToVect (MkVector v) = v + + export + Nil : Vector Z a + Nil = MkVector Nil + + export + (::) : a -> Vector n a -> Vector (S n) a + x :: MkVector xs = MkVector (x :: xs) + + export + Show coord => Show (Vector n coord) where + show = show . vectorToVect + + export + Neg coord => VectorSpace (Vector n coord) coord where + negateVector (MkVector v) = MkVector (negateVector v) + MkVector u ^+^ MkVector v = MkVector ((^+^) u v) + MkVector u ^-^ MkVector v = MkVector ((^-^) u v) + x *^ MkVector v = MkVector ((*^) x v) + + export + Neg coord => DotProduct (Vector n coord) coord where + MkVector u ^.^ MkVector v = u ^.^ v + + export + sumVectorsFrom : VectorSpace vector coord => vector -> List vector -> vector + sumVectorsFrom acc Nil = acc + sumVectorsFrom acc (h :: t) = sumVectorsFrom (acc ^+^ h) t + +namespace Matrix + + export + data Matrix : Nat -> Nat -> Type -> Type where + MkMatrix : Vect numRows (Vector numCols a) -> Matrix numRows numCols a + + export + Functor (Matrix a b) where + map f (MkMatrix m) = MkMatrix (map (map f) m) + + export + Show a => Show (Matrix r c a) where + show (MkMatrix m) = show m + + export + vectorsToMatrix : Vect numRows (Vector numCols coord) -> Matrix numRows numCols coord + vectorsToMatrix = MkMatrix + + export + vectsToMatrix : Vect numRows (Vect numCols coord) -> Matrix numRows numCols coord + vectsToMatrix = MkMatrix . map vectToVector + + export + matrixToVectors : Matrix numRows numCols coord -> Vect numRows (Vector numCols coord) + matrixToVectors (MkMatrix v) = v + + export + matrixToVects : Matrix numRows numCols coord -> Vect numRows (Vect numCols coord) + matrixToVects (MkMatrix v) = map vectorToVect v + + export + Nil : Matrix Z n a + Nil = MkMatrix Nil + + export + (::) : Vector n a -> Matrix m n a -> Matrix (S m) n a + x :: MkMatrix xs = MkMatrix (x :: xs) + + export + negateMatrix : Neg coord => Matrix a b coord -> Matrix a b coord + negateMatrix = vectorsToMatrix . map negateVector . matrixToVectors + + diagonalVects : Num coord => {n : Nat} -> Vect n coord -> Vect n (Vect n coord) + diagonalVects Nil = Nil + diagonalVects {n = S n'} (x :: xs) = + (x :: replicate n' 0) + :: + (map (0 ::) (diagonalVects xs)) + + export + diagonalMatrix : Num coord => {n : Nat} -> Vect n coord -> Matrix n n coord + diagonalMatrix = vectsToMatrix . diagonalVects + + export + identityMatrix : Num a => (n : Nat) -> Matrix n n a + identityMatrix n = diagonalMatrix (replicate n 1) + + export + transpose : {b : Nat} -> Matrix a b coord -> Matrix b a coord + transpose = + vectsToMatrix . + Data.Vect.transpose . + matrixToVects + + -- Scalar times matrix + export + (*!) : Num coord => coord -> Matrix a b coord -> Matrix a b coord + x *! m = map (x *) m + + -- Matrix sum + export + (!+!) : Neg coord => Matrix a b coord -> Matrix a b coord -> Matrix a b coord + MkMatrix x !+! MkMatrix y = MkMatrix $ zipWith (^+^) x y + + -- Matrix product + export + (!*!) : Neg coord => {c : Nat} -> Matrix a b coord -> Matrix b c coord -> Matrix a c coord + (!*!) (MkMatrix x) y = + let MkMatrix y' := transpose y in + MkMatrix $ + map + (\xRow => vectToVector (map (xRow ^.^) y')) + x + + -- This implementation is inefficient. More than O(n!) time. + export + determinant : (Neg coord, Num coord) => {n : Nat} -> Matrix n n coord -> coord + determinant (MkMatrix Nil) = 1 + determinant (MkMatrix (v :: vs)) = + sum $ + zipWith + ( \i, vi => + paritySign i * vi * determinant (MkMatrix (map (vectToVector . deleteAt i . vectorToVect) vs)) + ) + range + (vectorToVect v) + + swapRows : Fin n -> Fin n -> Vect n a -> Vect n a + swapRows FZ FZ v = v + swapRows (FS i) (FS j) (h :: t) = h :: swapRows i j t + swapRows FZ (FS j) (h :: t) = index j t :: replaceAt j h t + swapRows (FS i) FZ (h :: t) = index i t :: replaceAt i h t + + subtractRowsBelow : Neg coord => + {m : Nat} -> (n : Nat) -> Vect m coord -> + Vect (n + S m) (Vect k coord) -> + Vect (n + S m) (Vect k coord) + subtractRowsBelow n subtractionCoefficients matrix = + case splitAt {m = S m} n matrix of + (above, rowToSubtract :: below) => + above ++ + ( rowToSubtract :: + ( [ row ^-^ (coefficient *^ rowToSubtract) + | (row, coefficient) <- zip below subtractionCoefficients + ] + ) + ) + + subtractRowsAbove : Neg coord => + (n : Nat) -> Vect n coord -> + Vect (n + S m) (Vect k coord) -> + Vect (n + S m) (Vect k coord) + subtractRowsAbove n subtractionCoefficients matrix = + case splitAt {m = S m} n matrix of + (above, rowToSubtract :: below) => + [ row ^-^ (coefficient *^ rowToSubtract) + | (row, coefficient) <- zip above subtractionCoefficients + ] + ++ + rowToSubtract :: below + + export + matrixInverse : Abs coord => Fractional coord => Neg coord => Ord coord => + {n : Nat} -> Matrix n n coord -> Matrix n n coord + matrixInverse matrix = + let ||| Requires: + ||| * matrix * inverseSoFar = leftToInvert + ||| * leftToInvert is lower-triangular. + ||| * The diagonal of leftToInvert is all 1s. + ||| * The last n - numInverted columns of leftToInvert are the columns of + ||| the identity matrix. + ||| + ||| Returns the inverse of matrix. + finishInvertingLT : + {auto 0 eq : n = n'} -> + (numLeft : Nat) -> (0 lte : numLeft `LTE` n') -> + (inverseSoFar : Vect n (Vect n coord)) -> + (leftToInvert : Vect n (Vect n coord)) -> + Matrix n n coord + finishInvertingLT Z _ inverseSoFar _ = + vectsToMatrix inverseSoFar + finishInvertingLT + numLeft@(S numLeft') lte@(LTESucc lte') inverseSoFar leftToInvert = + let 0 lteN : numLeft `LTE` n + lteN = rewrite eq in lte + workingIndex : Fin n + workingIndex = natToFinLT numLeft' {prf = lteN} + numInverted : Nat + numInverted = minus n numLeft + 0 leftPlusInvertedEqN : (numLeft + numInverted = n) + leftPlusInvertedEqN = + trans + (plusCommutative numLeft numInverted) + (plusMinusLte numLeft n lteN) + toSubtract : Vect numInverted coord + toSubtract = + map (index workingIndex) $ + drop numLeft $ + rewrite leftPlusInvertedEqN in leftToInvert + 0 predLeftPlusSuccInvertedEqN : (numLeft' + S numInverted = n) + predLeftPlusSuccInvertedEqN = + trans + (sym $ plusSuccRightSucc numLeft' numInverted) + leftPlusInvertedEqN + leftToInvert' : Vect n (Vect n coord) + leftToInvert' = + rewrite sym predLeftPlusSuccInvertedEqN in + subtractRowsBelow numLeft' toSubtract $ + rewrite predLeftPlusSuccInvertedEqN in leftToInvert + inverseSoFar' : Vect n (Vect n coord) + inverseSoFar' = + rewrite sym predLeftPlusSuccInvertedEqN in + subtractRowsBelow numLeft' toSubtract $ + rewrite predLeftPlusSuccInvertedEqN in inverseSoFar + in finishInvertingLT numLeft' (lteSuccRight lte') inverseSoFar' leftToInvert' + ||| Requires: + ||| * matrix * inverseSoFar = leftToInvert + ||| * leftToInvert is lower-triangular. + ||| * The diagonal of leftToInvert is all 1s. + ||| + ||| Returns the inverse of matrix. + invertLT : (inverseSoFar : Vect n (Vect n coord)) -> + (leftToInvert : Vect n (Vect n coord)) -> + Matrix n n coord + invertLT = finishInvertingLT n reflexive + ||| Requires: + ||| * matrix * inverseSoFar = leftToInvert + ||| * The last numInverted columns of leftToInvert are lower-triangular. + ||| * The diagonal of leftToInvert ends with numInverted 1s. + ||| + ||| Returns the inverse of matrix. + invertFromPartlyLT : + {auto 0 eq : n = n'} -> + (numLeft : Nat) -> + (0 lte : numLeft `LTE` n') -> + (inverseSoFar : Vect n (Vect n coord)) -> + (leftToInvert : Vect n (Vect n coord)) -> + Matrix n n coord + invertFromPartlyLT Z _ inverseSoFar leftToInvert = + invertLT inverseSoFar leftToInvert + invertFromPartlyLT (S numLeft') (LTESucc lte') inverseSoFar leftToInvert = + let -- We don't just use an @ pattern to define numLeft because of + -- https://github.com/idris-lang/Idris2/issues/2249 + numLeft : Nat + numLeft = S numLeft' + 0 lte : numLeft `LTE` n + lte = rewrite eq in LTESucc lte' + numInverted : Nat + numInverted = minus n numLeft + 0 nEqLeftPlusInverted : (n = numLeft + numInverted) + nEqLeftPlusInverted = + sym $ + trans + (plusCommutative numLeft numInverted) + (plusMinusLte numLeft n lte) + 0 nEqPredLeftPlusSuccInverted : (n = numLeft' + S numInverted) + nEqPredLeftPlusSuccInverted = + trans + nEqLeftPlusInverted + (plusSuccRightSucc numLeft' numInverted) + workingIndex : Fin n + workingIndex = natToFinLT numLeft' {prf = lte} + toSwap : Fin n + toSwap = + rewrite nEqLeftPlusInverted in + weakenN numInverted $ + argmax $ map (abs . Vect.index workingIndex) $ + Vect.take numLeft $ + replace {p = \m => Vect m (Vect n coord)} + nEqLeftPlusInverted leftToInvert + -- Swap rows + inverseSoFar' : Vect n (Vect n coord) + inverseSoFar' = + swapRows workingIndex toSwap inverseSoFar + leftToInvert' : Vect n (Vect n coord) + leftToInvert' = + swapRows workingIndex toSwap leftToInvert + -- Divide the workingIndex-th row by the value at + -- (workingIndex, workingIndex). + workingValue : coord + workingValue = index workingIndex (index workingIndex leftToInvert') + inverseSoFar'' : Vect n (Vect n coord) + inverseSoFar'' = + updateAt workingIndex (map (/ workingValue)) inverseSoFar' + leftToInvert'' : Vect n (Vect n coord) + leftToInvert'' = + updateAt workingIndex (map (/ workingValue)) leftToInvert' + -- Zero out the entries above (workingIndex, workingIndex). + subtractionCoefficients : Vect numLeft' coord + subtractionCoefficients = + map (index workingIndex) $ take numLeft' $ + replace {p = \m => Vect m (Vect n coord)} + nEqPredLeftPlusSuccInverted leftToInvert'' + inverseSoFar''' : Vect n (Vect n coord) + inverseSoFar''' = + replace {p = \m => Vect m (Vect n coord)} + (sym nEqPredLeftPlusSuccInverted) + $ + subtractRowsAbove numLeft' subtractionCoefficients $ + the (Vect (numLeft' + S numInverted) _) $ + replace {p = \m => Vect m (Vect n coord)} + nEqPredLeftPlusSuccInverted + inverseSoFar'' + leftToInvert''' : Vect n (Vect n coord) + leftToInvert''' = + replace {p = \m => Vect m (Vect n coord)} + (sym nEqPredLeftPlusSuccInverted) + $ + subtractRowsAbove numLeft' subtractionCoefficients $ + the (Vect (numLeft' + S numInverted) _) $ + replace {p = \m => Vect m (Vect n coord)} + nEqPredLeftPlusSuccInverted + leftToInvert'' + in invertFromPartlyLT numLeft' (lteSuccRight lte') + inverseSoFar''' leftToInvert''' + in invertFromPartlyLT n reflexive (diagonalVects (replicate n 1)) (matrixToVects matrix) + + ||| This behaves like Euclidean translation near [0, 0, 0, 1]. + export + translate : Vect 3 Double -> Matrix 4 4 Double + translate v = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], vectToVector (v ++ [1])] + +-- Matrix-vector product +export +(!*^) : Neg coord => {a : Nat} -> Matrix a b coord -> Vector b coord -> Vector a coord +(!*^) a x = + let a' = matrixToVectors a + in vectToVector (map (x ^.^) a') + +-- Vector-matrix product +export +(^*!) : Neg coord => {b : Nat} -> Vector a coord -> Matrix a b coord -> Vector b coord +(^*!) x a = transpose a !*^ x + +-- Cross product. +export +cross : (Neg coord, Num coord) => {n : Nat} -> Vect n (Vector (S n) coord) -> Vector (S n) coord +cross vs = + vectToVector $ + tabulate $ + \ i => paritySign i * determinant (vectsToMatrix (map (deleteAt i . vectorToVect) vs)) + +||| Given two nonzero vectors, returns a vector that's not in their span. +notSpanned : Vector 4 Double -> Vector 4 Double -> + Vector 4 Double +notSpanned u v = + let candidates : Vect 3 (Vector 4 Double) + candidates = + [ [1, 0, 0, 0] + , [0, 1, 0, 0] + , [0, 0, 1, 0] + ] + in index (argmax [norm $ cross [u, v, c] | c <- candidates]) candidates + +||| `from` and `to` should not be colinear. +||| +||| This returns a matrix `m` such that `normalize from !*! m = normalize to`. +||| +||| TODO: test +export +rotationMatrix : Vector 4 Double -> Vector 4 Double -> Matrix 4 4 Double +rotationMatrix from to = + let -- Vectors that aren't moved by the rotation. + fixed0 : Vector 4 Double + fixed0 = normalize $ cross [from, to, notSpanned from to] + fixed1 : Vector 4 Double + fixed1 = normalize $ cross [from, to, fixed0] + -- A vector that moves along with `from`. + otherFrom : Vector 4 Double + otherFrom = normalize $ cross [from, fixed0, fixed1] + -- Where otherFrom goes. + otherTo : Vector 4 Double + otherTo = normalize $ cross [to, fixed0, fixed1] + in transpose [normalize from, otherFrom, fixed0, fixed1] !*! [normalize to, otherTo, fixed0, fixed1] blob - /dev/null blob + cb61879cff4b669e76d77a6cd3e608d841e464b5 (mode 644) --- /dev/null +++ src/Math/LinearAlgebraTest.idr @@ -0,0 +1,70 @@ +module Math.LinearAlgebraTest + +import Data.List +import Data.Vect +import Math.LinearAlgebra +import System.Random +import Test.TestTree + +%default total + +almostEqual : Double -> Double -> Bool +almostEqual x y = abs (x - y) < 0.0001 + +vectorsAlmostEqual : Vect n Double -> Vect n Double -> Bool +vectorsAlmostEqual u v = all (uncurry almostEqual) (zip u v) + +isAlmostIdentity : {n : Nat} -> Matrix n n Double -> Bool +isAlmostIdentity matrix = + all (uncurry vectorsAlmostEqual) $ + zip (matrixToVects matrix) (matrixToVects (identityMatrix _)) + +inverseTestCase : {n : Nat} -> String -> IO (Matrix n n Double) -> TestCase +inverseTestCase name makeMatrix = + MkTestCase + { name = name + , run = do matrix <- makeMatrix + let computedInverse = matrixInverse matrix + let product = computedInverse !*! matrix + pure $ + if isAlmostIdentity product + then MkSuccess + else MkFailure $ + fastConcat + [ show computedInverse + , " (inverse) !*! " + , show matrix + , " (original) = " + , show product + , " which is not close to the identity matrix" + ] + } + +randomUpperTriangular : {n : Nat} -> IO (Matrix n n Double) +randomUpperTriangular {n = Z} = pure Nil +randomUpperTriangular {n = S _} = + do firstRow <- map vectToVector $ sequence (replicate _ $ randomRIO (-100, 100)) + rest <- map (vectorsToMatrix . map (0 ::) . matrixToVectors) randomUpperTriangular + pure (firstRow :: rest) + +export +tests : TestTree +tests = + MkTestGroup "LinearAlgebra.inverse" $ map MkTestLeaf $ + [ inverseTestCase "0x0 matrix" $ pure Nil + , inverseTestCase "1x1 identity" $ pure $ identityMatrix 1 + , inverseTestCase "the number 5" $ pure $ [[5]] + , inverseTestCase "5x5 identity" $ pure $ identityMatrix 5 + , inverseTestCase "diagonal" $ pure $ diagonalMatrix [2, 1, 0.1, 4] + , inverseTestCase "010 001 -100" $ pure [[0, 1, 0], [0, 0, 1], [-1, 0, 0]] + , inverseTestCase "another 3x3" $ pure [[-1, -1, -1], [-1, 0, -1], [0, -1, 1]] + ] + ++ + concat + [ replicate 3 $ + inverseTestCase ("random " <+> show n <+> "x" <+> show n) $ + do upper <- randomUpperTriangular {n} + lower <- map transpose randomUpperTriangular + pure (lower !*! upper) + | n <- the (List _) [1, 2, 3, 4, 5] + ] blob - /dev/null blob + dfb5e75f8dbc966d1add1b65bbe4129623499ad8 (mode 644) --- /dev/null +++ src/PrecomputeMain.idr @@ -0,0 +1,16 @@ +module PrecomputeMain + +import S3D.Polyhedra.Precompute +import S3D.Polyhedra.ToPrecompute + +%default total + +partial +forceEither : Either String a -> a +forceEither (Left s) = idris_crash s +forceEither (Right x) = x + +partial +main : IO () +main = putStr $ forceEither $ + precomputePolyhedra toPrecompute "S3D.Polyhedra.Precomputed" blob - /dev/null blob + 48de365bd78ed20623fb3963217dd780bd31d375 (mode 644) --- /dev/null +++ src/S3D/AnimationFrame/Prim.idr @@ -0,0 +1,5 @@ +module S3D.AnimationFrame.Prim + +export +%foreign "javascript:lambda:(callback) => window.requestAnimationFrame(time => callback(time)())" +prim__requestAnimationFrame : (Double -> PrimIO ()) -> PrimIO () blob - /dev/null blob + 4d3508f0ae03f5e67bd053921e01b95e0147f233 (mode 644) --- /dev/null +++ src/S3D/AnimationFrame.idr @@ -0,0 +1,20 @@ +module S3D.AnimationFrame + +import S3D.AnimationFrame.Prim +import System.Clock + +%default total + +fromSeconds : {type : ClockType} -> Double -> Clock type +fromSeconds seconds = + MkClock (cast (floor seconds)) (cast ((seconds - floor seconds) * 1000000000)) + +fromMillis : {type : ClockType} -> Double -> Clock type +fromMillis millis = fromSeconds (millis / 1000) + +export +requestAnimationFrame : HasIO io => (Clock Monotonic -> IO ()) -> io () +requestAnimationFrame callback = + let callback' : Double -> PrimIO () + callback' time = toPrim $ callback $ fromMillis time + in do primIO $ prim__requestAnimationFrame callback' blob - /dev/null blob + f05b8af787e1a429d8e156cdc7f86df7e012d372 (mode 644) --- /dev/null +++ src/S3D/Array.idr @@ -0,0 +1,8 @@ +module S3D.Array + +import JS.Array + +-- A version of JS.Array.fromListIO that doesn't blow the stack. +export +tailRecFromListIO : HasIO io => List a -> io (Array a) +tailRecFromListIO as = liftIO $ fromList' as thaw blob - /dev/null blob + 853ad2d8d7dc02d372ccac3668d0175e17311fe9 (mode 644) --- /dev/null +++ src/S3D/ArrayBuffer/Prim.idr @@ -0,0 +1,15 @@ +module S3D.ArrayBuffer.Prim + +import JS + +export +%foreign "javascript:lambda:(x) => new Float32Array(x)" +prim__newFloat32Array : Array Double -> PrimIO Float32Array + +export +%foreign "javascript:lambda:(x) => new Int16Array(x)" +prim__newInt16Array : Array Int16 -> PrimIO Int16Array + +export +%foreign "javascript:lambda:(x) => new Uint8Array(x)" +prim__newUInt8Array : Array Bits8 -> PrimIO UInt8Array blob - /dev/null blob + a86059b073ec7224122e3b60897c227ce0953152 (mode 644) --- /dev/null +++ src/S3D/ArrayBuffer.idr @@ -0,0 +1,18 @@ +module S3D.ArrayBuffer + +import JS +import S3D.ArrayBuffer.Prim + +%default total + +export +makeFloat32Array : Array Double -> JSIO Float32Array +makeFloat32Array = primIO . prim__newFloat32Array + +export +makeInt16Array : Array Int16 -> JSIO Int16Array +makeInt16Array = primIO . prim__newInt16Array + +export +makeUInt8Array : Array Bits8 -> JSIO UInt8Array +makeUInt8Array = primIO . prim__newUInt8Array blob - /dev/null blob + 53f7951ac3f850aee63142a06b7e10695660dc17 (mode 644) --- /dev/null +++ src/S3D/CSG/DrawablePolyhedron.idr @@ -0,0 +1,40 @@ +||| Constructive Solid Geometry with `DrawablePolyhedron`s. + +module S3D.CSG.DrawablePolyhedron + +import Data.Vect +import S3D.DrawablePolyhedron +import SolidGeometry.Construct + +%default total + +export +complement : DrawablePolyhedron numPlanes attributes False -> + DrawablePolyhedron numPlanes attributes False +complement = { polyhedron $= SolidGeometry.Construct.complement } + +export +intersect : {numPlanes0 : Nat} -> + DrawablePolyhedron numPlanes0 attributes False -> + DrawablePolyhedron numPlanes1 attributes False -> + DrawablePolyhedron (numPlanes0 + numPlanes1) attributes False +intersect drawable0 drawable1 = + MkDrawablePolyhedron + { polyhedron = SolidGeometry.Construct.intersect + drawable0.polyhedron drawable1.polyhedron + , faceAttributes = drawable0.faceAttributes ++ drawable1.faceAttributes + , triangulation = () + } + +export +union : {numPlanes0 : Nat} -> + DrawablePolyhedron numPlanes0 attributes False -> + DrawablePolyhedron numPlanes1 attributes False -> + DrawablePolyhedron (numPlanes0 + numPlanes1) attributes False +union drawable0 drawable1 = + MkDrawablePolyhedron + { polyhedron = SolidGeometry.Construct.union + drawable0.polyhedron drawable1.polyhedron + , faceAttributes =drawable0.faceAttributes ++ drawable1.faceAttributes + , triangulation = () + } blob - /dev/null blob + 691740da9ed6733da7391aeb1fe34598a8fef3de (mode 644) --- /dev/null +++ src/S3D/CSG/PhysicalPolyhedron.idr @@ -0,0 +1,43 @@ +||| Constructive Solid Geometry with `PhysicalPolyhedron`s. + +module S3D.CSG.PhysicalPolyhedron + +import S3D.CSG.DrawablePolyhedron +import S3D.PhysicalPolyhedron +import SolidGeometry.Construct + +%default total + +export +complement : PhysicalPolyhedron n attributes False -> + PhysicalPolyhedron n attributes False +complement = + { drawablePolyhedron $= S3D.CSG.DrawablePolyhedron.complement + , obstaclePolyhedron $= SolidGeometry.Construct.complement + } + +export +intersect : {n0 : Nat} -> + PhysicalPolyhedron n0 attributes False -> + PhysicalPolyhedron n1 attributes False -> + PhysicalPolyhedron (n0 + n1) attributes False +intersect poly0 poly1 = + MkPhysicalPolyhedron + { drawablePolyhedron = + intersect poly0.drawablePolyhedron poly1.drawablePolyhedron + , obstaclePolyhedron = + intersect poly0.obstaclePolyhedron poly1.obstaclePolyhedron + } + +export +union : {n0 : Nat} -> + PhysicalPolyhedron n0 attributes False -> + PhysicalPolyhedron n1 attributes False -> + PhysicalPolyhedron (n0 + n1) attributes False +union poly0 poly1 = + MkPhysicalPolyhedron + { drawablePolyhedron = + union poly0.drawablePolyhedron poly1.drawablePolyhedron + , obstaclePolyhedron = + union poly0.obstaclePolyhedron poly1.obstaclePolyhedron + } blob - /dev/null blob + d623c0c40ad7169342f6bdfc5dea7c1ff365cce9 (mode 644) --- /dev/null +++ src/S3D/Controls.idr @@ -0,0 +1,71 @@ +module S3D.Controls + +import Data.IORef +import Data.SortedSet +import Data.Vect +import Math.LinearAlgebra +import S3D.Obstacles +import S3D.PlayerMovement +import S3D.WebUI + +%default total + +||| `positiveKeyDown` and `negativeKeyDown` should be the state of an opposite +||| pair of keys --- up and down, or left and right. +||| +||| If just the positive key is down, applies `update 1`. If just the negative +||| key: `update (-1)`. Otherwise applies the identity function. +arrowKeyChange : (positiveKeyDown : Bool) -> (negativeKeyDown : Bool) -> + (update : Double -> a -> a) -> a -> a +arrowKeyChange True False f = f 1 +arrowKeyChange False True f = f (-1) +arrowKeyChange _ _ _ = id + +capMagnitude : Double -> Double -> Double +capMagnitude cap x = + if abs x > cap + then x / abs x * cap + else x + +export +stepPlacement : HasIO io => + Obstacles -> PlayerPlacement -> IORef UIState -> io PlayerPlacement +stepPlacement obstacles oldPlacement uiStateRef = + do uiState <- readIORef uiStateRef + modifyIORef uiStateRef { touchDeltas := replicate _ [0, 0] } + pure $ + let [ [touchTurnX, touchTurnY] + , [touchMoveX, touchMoveY] + ] = map vectorToVect uiState.touchDeltas in + let shift = contains "Shift" uiState.keysDown + right = contains "ArrowRight" uiState.keysDown + left = contains "ArrowLeft" uiState.keysDown + up = contains "ArrowUp" uiState.keysDown + down = contains "ArrowDown" uiState.keysDown + -- Limit the movement speed from touch events to stop the player + -- from going through walls. Collision detection only checks whether + -- the player would end up inside something, not whether they jump + -- through something. + touchMoveX' = capMagnitude 0.06 touchMoveX + touchMoveY' = capMagnitude 0.06 touchMoveY + ||| The new placement only taking keyboard controls into account. + placementFromKeyboard : PlayerPlacement + placementFromKeyboard = + if shift + then arrowKeyChange up down oneFrameLookUp $ + arrowKeyChange right left oneFrameStrafe $ + oldPlacement + else arrowKeyChange right left oneFrameTurn $ + arrowKeyChange up down oneFrameMoveForward $ + oldPlacement + ||| Where the player would go if we didn't have collision detection. + attemptedPlacement : PlayerPlacement + attemptedPlacement = + touchMoveForward (-touchMoveY') $ + touchStrafe (-touchMoveX') $ + touchLookUp (-touchTurnY) $ + touchTurn (-touchTurnX) $ + placementFromKeyboard + in if collision (position attemptedPlacement) obstacles + then oldPlacement + else attemptedPlacement blob - /dev/null blob + d65cfc16c88e682652ad108b4963e07180922e2d (mode 644) --- /dev/null +++ src/S3D/DomExtras.idr @@ -0,0 +1,50 @@ +||| Stuff not in idris2-dom. + +module S3D.DomExtras + +import JS +import Web.Dom + +%default total + +export +data Touch : Type where + +export +data TouchList : Type where + +%foreign "javascript:lambda:(x) => x.clientX" +prim__getClientX : Touch -> PrimIO Int32 + +%foreign "javascript:lambda:(x) => x.clientY" +prim__getClientY : Touch -> PrimIO Int32 + +%foreign "javascript:lambda:(x) => x.length" +prim__getNumTouches : TouchList -> PrimIO Int32 + +%foreign "javascript:lambda:(t,i) => t.item(i)" +prim__getTouchItem : TouchList -> Int32 -> PrimIO Touch + +%foreign "javascript:lambda:(x) => x.touches" +prim__getTouchList : Event -> PrimIO TouchList + +export +getClientX : Touch -> JSIO Int32 +getClientX = primIO . prim__getClientX + +export +getClientY : Touch -> JSIO Int32 +getClientY = primIO . prim__getClientY + +export +getNumTouches : TouchList -> JSIO Nat +getNumTouches touchList = + map cast $ primIO $ prim__getNumTouches touchList + +export +getTouchItem : TouchList -> Nat -> JSIO Touch +getTouchItem touchList i = primIO $ prim__getTouchItem touchList $ cast i + +export +getTouchList : Event -> JSIO TouchList +getTouchList = primIO . prim__getTouchList blob - /dev/null blob + 768195839f38f7b9b59668bc6a6482517b7eb3ca (mode 644) --- /dev/null +++ src/S3D/Draw.idr @@ -0,0 +1,25 @@ +module S3D.Draw + +import Data.Bits +import Data.Vect +import JS +import S3D.Drawable +import Math.LinearAlgebra +import Web.Html +import Web.Raw.Webgl + +%default total + +-- Draws everything twice, so things behind the viewer are visible. +-- +-- This function handles shrinking z coordinates so that far-away things aren't clipped. +export +drawUniverse : WebGL2RenderingContext -> Matrix 4 4 Double -> List Drawable -> JSIO () +drawUniverse gl viewMatrix drawables = + do let shrunkViewMatrix = viewMatrix !*! (diagonalMatrix [1, 1, 0.00001, 1]) + let negShrunkViewMatrix = negateMatrix shrunkViewMatrix + clearColor gl 0 0 0 1 + clear gl (COLOR_BUFFER_BIT .|. DEPTH_BUFFER_BIT) + traverse_ (\drawable => draw gl drawable negShrunkViewMatrix) drawables + clear gl DEPTH_BUFFER_BIT + traverse_ (\drawable => draw gl drawable shrunkViewMatrix) drawables blob - /dev/null blob + d692a33f0a663f5797e515f45372bfcf03286b4a (mode 644) --- /dev/null +++ src/S3D/DrawInput.idr @@ -0,0 +1,58 @@ +module S3D.DrawInput + +import Data.List +import Data.Vect +import Math.LinearAlgebra +import S3D.Numbers +import S3D.Transformable +import SolidGeometry.Polygon + +%default total + +public export +record TriangleIndices where + constructor T + indices : Vect 3 Int16 + +export +convertTriangle : Triangle (Fin n) -> TriangleIndices +convertTriangle (MkTriangle a b c) = T (map finToInt16 [a, b, c]) + +-- A sequence of vertex attributes together with a list of triangles to draw +-- given as indices into that sequence. +public export +record DrawInput attributes where + constructor MkDrawInput + triangles : List TriangleIndices + vertexAttributes : List attributes + +export +Transformable attributes => Transformable (DrawInput attributes) where + transform matrix = { vertexAttributes $= map (transform matrix) } + +export +emptyDrawInput : DrawInput attributes +emptyDrawInput = MkDrawInput Nil Nil + +addOffsetToTriangle : Int16 -> TriangleIndices -> TriangleIndices +addOffsetToTriangle offset (T [a, b, c]) = T [a + offset, b + offset, c + offset] + +genericLength : Num a => List b -> a +genericLength = fromInteger . cast . length + +concatDrawInputs' : DrawInput attributes -> Int16 -> List (DrawInput attributes) -> DrawInput attributes +concatDrawInputs' tailDrawInput _ Nil = tailDrawInput +concatDrawInputs' tailDrawInput numVerticesRemaining (nextDrawInput :: reversedDrawInputs) = + let numVerticesRemaining' = numVerticesRemaining - genericLength (vertexAttributes nextDrawInput) + tailDrawInput' = + { triangles $= (map (addOffsetToTriangle numVerticesRemaining') (triangles nextDrawInput) ++) + , vertexAttributes $= (vertexAttributes nextDrawInput ++) + } + tailDrawInput + in concatDrawInputs' tailDrawInput' numVerticesRemaining' reversedDrawInputs + +export +concatDrawInputs : List (DrawInput attributes) -> DrawInput attributes +concatDrawInputs drawInputs = + let numVertices = sum (map (genericLength . vertexAttributes) drawInputs) + in concatDrawInputs' emptyDrawInput numVertices (reverse drawInputs) blob - /dev/null blob + 2e08c1889e83c1f0a8cc5fc93d1afe6bb5a0b39d (mode 644) --- /dev/null +++ src/S3D/Drawable.idr @@ -0,0 +1,34 @@ +module S3D.Drawable + +import Data.Vect +import JS +import S3D.Array +import Math.LinearAlgebra +import Web.Html +import Web.Raw.Webgl + +%default total + +-- A Drawable is a GL vertex array together with everything needed to draw it. +-- We also store the associated buffers so that we can delete them when we're +-- done with the Drawable. +public export +record Drawable where + constructor MkDrawable + program : WebGLProgram + viewMatrixLocation : WebGLUniformLocation + vertexArray : WebGLVertexArrayObject + -- 3 times the number of triangles + numElements : Int32 + buffers : List WebGLBuffer + +export +draw : WebGL2RenderingContext -> Drawable -> Matrix 4 4 Double -> JSIO () +draw gl drawable matrix = + do useProgram gl (Just drawable.program) + bindVertexArray gl (Just drawable.vertexArray) + matrixArray <- tailRecFromListIO (concat (map (toList . vectorToVect) (matrixToVectors matrix))) + uniformMatrix4fv' gl (Just drawable.viewMatrixLocation) False (S (Z matrixArray)) + let count = drawable.numElements + let offset = 0 + drawElements gl TRIANGLES count UNSIGNED_SHORT offset blob - /dev/null blob + 3e9739fdfac47cefcfd0d3a1da8b45b3b2e2d3e9 (mode 644) --- /dev/null +++ src/S3D/DrawablePolyhedron/SimpleColour.idr @@ -0,0 +1,30 @@ +||| DrawablePolyhedrons with SimpleColour attributes. + +module S3D.DrawablePolyhedron.SimpleColour + +import Data.Vect +import S3D.DrawablePolyhedron +import S3D.GLProgram.SimpleColour.Attributes +import Math.LinearAlgebra +import Math.LinearAlgebra.Normalized +import SolidGeometry.Polyhedron + +%default total + +buildAttributes : Vector 4 Double -> VertexCoordinates -> Attributes +buildAttributes colour coordinates = + MkAttributes + { position = fromNormalized coordinates.absolute + , colour = colour + } + +export +makeSimpleColourPolyhedron : + Polyhedron numPlanes -> Vect numPlanes (Vect 2 (Vector 4 Double)) -> + DrawablePolyhedron numPlanes Attributes False +makeSimpleColourPolyhedron polyhedron faceColours = + MkDrawablePolyhedron + { polyhedron = polyhedron + , faceAttributes = map (map buildAttributes) faceColours + , triangulation = () + } blob - /dev/null blob + cf8a941621465c5cc57ae44804643e7ee552bf62 (mode 644) --- /dev/null +++ src/S3D/DrawablePolyhedron.idr @@ -0,0 +1,71 @@ +||| A `DrawablePolyhedron` is a Polyhedron with attribute information +||| associated with each face. These can be converted to `drawInput`s using +||| `drawPolyhedron`. + +module S3D.DrawablePolyhedron + +import Data.Vect +import S3D.DrawInput +import Math.LinearAlgebra +import SolidGeometry.Polyhedron +import SolidGeometry.Sign + +%default total + +||| Attributes that depend on a vertex's absolute location, and its location in +||| the local coordinate system of the face it's on. +public export +AttributeBuilder : Type -> Type +AttributeBuilder attributes = VertexCoordinates -> attributes + +||| If pretriangulated is True, this includes a precomputed triangulation of +||| the polyhedron. +public export +record DrawablePolyhedron (numPlanes : Nat) (attributes : Type) (pretriangulated : Bool) where + constructor MkDrawablePolyhedron + polyhedron : Polyhedron numPlanes + ||| For each face index i, how to compute attributes for (MkFaceId i + ||| Positive) and for (mkFaceId i Negative). + faceAttributes : Vect numPlanes (Vect 2 (AttributeBuilder attributes)) + triangulation : if pretriangulated then List (TriangulatedFace numPlanes) else () + +signToIndex : Sign -> Fin 2 +signToIndex Positive = 0 +signToIndex Negative = 1 + +lookUpAttributes : FaceId numPlanes -> Vect numPlanes (Vect 2 a) -> a +lookUpAttributes face = + index (signToIndex face.outside) . index face.planeIndex + +faceToDrawInput : {0 numPlanes : Nat} -> TriangulatedFace numPlanes -> + AttributeBuilder attributes -> DrawInput attributes +faceToDrawInput face buildAttributes = + MkDrawInput + (map convertTriangle face.triangles) + (map buildAttributes $ toList face.vertices) + +export +addTriangulation : List (TriangulatedFace numPlanes) -> + DrawablePolyhedron numPlanes attributes False -> + DrawablePolyhedron numPlanes attributes True +addTriangulation t = {triangulation := t} + +||| Turn a `DrawablePolyhedron` into a `DrawInput` using the precomputed +||| triangulation. +export +drawPrecomputedPolyhedron : {numPlanes : Nat} -> DrawablePolyhedron numPlanes attributes True -> + DrawInput attributes +drawPrecomputedPolyhedron polyhedron = + concatDrawInputs + [ faceToDrawInput face (lookUpAttributes face.faceId polyhedron.faceAttributes) + | face <- polyhedron.triangulation + ] + +||| Turn a `DrawablePolyhedron` into a `DrawInput` by computing a new +||| triangulation. +export partial +drawPolyhedron : {numPlanes : Nat} -> DrawablePolyhedron numPlanes attributes False -> + Either String (DrawInput attributes) +drawPolyhedron polyhedron = + do faces <- triangulate polyhedron.polyhedron + pure $ drawPrecomputedPolyhedron $ addTriangulation faces polyhedron blob - /dev/null blob + 4b4a1152490d292f162666f6fbb3ec566b7f3c14 (mode 644) --- /dev/null +++ src/S3D/Figure.idr @@ -0,0 +1,58 @@ +{- + +A Figure consists of a DrawInput for each program we are willing to run on +every frame. + +-} + +module S3D.Figure + +import JS +import S3D.DrawInput +import S3D.Drawable +import S3D.GLProgram +import S3D.GLProgram.ColourTexture +import S3D.GLProgram.SimpleColour +import Math.LinearAlgebra +import S3D.Transformable +import Web.Html + +%default total + +public export +record Figure where + constructor MkFigure + colourTextureInput : DrawInput ColourTexture.Attributes + simpleColourInput : DrawInput SimpleColour.Attributes + +export +Transformable Figure where + transform matrix = + { colourTextureInput $= transform matrix + , simpleColourInput $= transform matrix + } + +export +compileFigure : WebGL2RenderingContext + -> S3D.GLProgram.ColourTexture.Program -> S3D.GLProgram.SimpleColour.Program -> Figure + -> JSIO (List Drawable) +compileFigure gl colourTextureProgram simpleColourProgram figure = + sequence + [ compileDrawInput gl colourTextureProgram figure.colourTextureInput + , compileDrawInput gl simpleColourProgram figure.simpleColourInput + ] + +export +concatFigures : List Figure -> Figure +concatFigures fs = + MkFigure + (concatDrawInputs (map colourTextureInput fs)) + (concatDrawInputs (map simpleColourInput fs)) + +export +colourTextureFigure : DrawInput ColourTexture.Attributes -> Figure +colourTextureFigure drawInput = MkFigure drawInput emptyDrawInput + +export +simpleColourFigure : DrawInput SimpleColour.Attributes -> Figure +simpleColourFigure drawInput = MkFigure emptyDrawInput drawInput blob - /dev/null blob + e90a5f4f8d42b0dd39abbb6167ad53a3a20faf8f (mode 644) --- /dev/null +++ src/S3D/GLProgram/ColourTexture.idr @@ -0,0 +1,107 @@ +module S3D.GLProgram.ColourTexture + +import Data.Vect +import JS +import S3D.DrawInput +import S3D.Drawable +import S3D.GLProgram +import S3D.GLProgram.Helpers +import Math.LinearAlgebra +import S3D.Transformable +import Web.Html +import Web.Raw.Webgl + +%default total + +public export +record Program where + constructor MkProgram + program : WebGLProgram + viewMatrixLocation, textureLocation : WebGLUniformLocation + positionLocation, colourLocation, textureCoordsLocation : Bits32 + +-- "ColourTexture": vertices have a position, texture coordinates, and a colour to combine with the texture. +public export +record Attributes where + constructor MkAttributes + position, colour : Vector 4 Double + textureCoords : Vector 2 Double + +export +Transformable Attributes where + transform matrix = { position $= transform matrix } + +vertexShaderSource : String +vertexShaderSource = + "#version 300 es\n" ++ + " " ++ + "in vec4 a_position; " ++ + "uniform mat4 u_matrix; " ++ + "in vec4 a_colour; " ++ + "in vec2 a_texture_coords; " ++ + " " ++ + "out vec4 v_colour; " ++ + "out vec2 v_texture_coords; " ++ + " " ++ + "void main() { " ++ + " gl_Position = u_matrix * a_position; " ++ + " v_colour = a_colour; " ++ + " v_texture_coords = a_texture_coords; " ++ + "}" + +fragmentShaderSource : String +fragmentShaderSource = + "#version 300 es\n" ++ + " " ++ + "precision mediump float; " ++ + " " ++ + "uniform sampler2D u_texture; " ++ + " " ++ + "in vec4 v_colour; " ++ + "in vec2 v_texture_coords; " ++ + " " ++ + "out vec4 colour; " ++ + " " ++ + "void main() { " ++ + " colour = v_colour * texture(u_texture, v_texture_coords); " ++ + "}" + +export +GLProgram Program Attributes where + makeProgram gl = + do program <- compileAndLinkProgram gl vertexShaderSource fragmentShaderSource + matrixLocation <- getAndCheckUniformLocation gl program "u_matrix" + textureLocation <- getAndCheckUniformLocation gl program "u_texture" + positionLocation <- getAndCheckAttribLocation gl program "a_position" + colourLocation <- getAndCheckAttribLocation gl program "a_colour" + textureCoordsLocation <- getAndCheckAttribLocation gl program "a_texture_coords" + pure $ + MkProgram + { program = program + , viewMatrixLocation = matrixLocation + , textureLocation = textureLocation + , positionLocation = positionLocation + , colourLocation = colourLocation + , textureCoordsLocation = textureCoordsLocation + } + + compileDrawInput gl program drawInput = + do vertexArray <- unMaybe "createVertexArray" $ createVertexArray gl + bindVertexArray gl (Just vertexArray) + positionBuffer <- makeAndBindBuffer gl program.positionLocation + (map (vectorToVect . position) drawInput.vertexAttributes) + colourBuffer <- makeAndBindBuffer gl program.colourLocation + (map (vectorToVect . colour) drawInput.vertexAttributes) + textureCoordsBuffer <- + makeAndBindBuffer gl program.textureCoordsLocation + (map (vectorToVect . textureCoords) drawInput.vertexAttributes) + elementBuffer <- makeAndBindElementBuffer gl drawInput.triangles + let numElements = 3 * length drawInput.triangles + pure $ + MkDrawable + { program = program.program + , viewMatrixLocation = program.viewMatrixLocation + , vertexArray = vertexArray + , numElements = cast numElements + , buffers = [positionBuffer, colourBuffer, textureCoordsBuffer, elementBuffer] + } blob - /dev/null blob + 720acfbc86505cff1113fb8a9e4d21f8c1f87929 (mode 644) --- /dev/null +++ src/S3D/GLProgram/Helpers.idr @@ -0,0 +1,84 @@ +module S3D.GLProgram.Helpers + +import Data.Vect +import JS +import S3D.Array +import S3D.ArrayBuffer +import S3D.DrawInput +import Web.Html +import Web.Raw.Webgl + +%default total + +export +makeAndBindBuffer : {componentsPerIteration : Nat} -> WebGL2RenderingContext -> Bits32 -> List (Vect componentsPerIteration Double) -> JSIO WebGLBuffer +makeAndBindBuffer gl location content = + do arrayBuffer <- makeFloat32Array !(tailRecFromListIO (concat (map toList content))) + buffer <- unMaybe "createBuffer" $ createBuffer gl + bindBuffer gl ARRAY_BUFFER (Just buffer) + bufferData1 gl ARRAY_BUFFER (Just (S (S (S (S (S (S (S (Z arrayBuffer))))))))) STATIC_DRAW + enableVertexAttribArray gl location + vertexAttribPointer gl + location + (cast componentsPerIteration) + FLOAT + {- normalize? -} False + {- stride = -} 0 + {- offset = -} 0 + pure buffer + +export +makeAndBindElementBuffer : WebGL2RenderingContext -> List TriangleIndices -> JSIO WebGLBuffer +makeAndBindElementBuffer gl triangleIndices = + do arrayBuffer <- makeInt16Array !(tailRecFromListIO (concat (map (toList . indices) triangleIndices))) + buffer <- unMaybe "createBuffer" $ createBuffer gl + bindBuffer gl ELEMENT_ARRAY_BUFFER (Just buffer) + bufferData1 gl ELEMENT_ARRAY_BUFFER (Just (S (Z arrayBuffer))) STATIC_DRAW + pure buffer + +castVia : (t : Type) -> (SafeCast t, FromFFI b t) => a -> Maybe b +castVia _ x = safeCast x >>= fromFFI + +makeShader : WebGL2RenderingContext -> Bits32 -> String -> JSIO WebGLShader +makeShader gl type source = + do shader <- unMaybe "createShader" $ createShader gl type + shaderSource gl shader source + compileShader gl shader + success <- unMaybe "castVia (getShaderParameter)" $ map (castVia Boolean) $ getShaderParameter gl shader COMPILE_STATUS + if success + then pure shader + else do log <- unMaybe "getShaderInfoLog" $ getShaderInfoLog gl shader + deleteShader gl (Just shader) + throwError $ Caught log + +export +compileAndLinkProgram : WebGL2RenderingContext -> String -> String -> JSIO WebGLProgram +compileAndLinkProgram gl vertexShaderSource fragmentShaderSource = + do vertexShader <- makeShader gl VERTEX_SHADER vertexShaderSource + fragmentShader <- makeShader gl FRAGMENT_SHADER fragmentShaderSource + program <- unMaybe "createProgram" $ createProgram gl + attachShader gl program vertexShader + attachShader gl program fragmentShader + linkProgram gl program + success <- unMaybe "castVia (getProgramParameter)" $ map (castVia Boolean) $ getProgramParameter gl program LINK_STATUS + + if success + then pure program + else do log <- unMaybe "getProgramInfoLog" $ getProgramInfoLog gl program + deleteProgram gl (Just program) + deleteShader gl (Just vertexShader) + deleteShader gl (Just fragmentShader) + throwError $ Caught log + +export +getAndCheckUniformLocation : WebGL2RenderingContext -> WebGLProgram -> String -> JSIO WebGLUniformLocation +getAndCheckUniformLocation gl program attributeName = + unMaybe "getUniformLocation" $ getUniformLocation gl program attributeName + +export +getAndCheckAttribLocation : WebGL2RenderingContext -> WebGLProgram -> String -> JSIO Bits32 +getAndCheckAttribLocation gl program attributeName = + do location <- getAttribLocation gl program attributeName + if location >= 0 + then pure (cast location) + else throwError $ Caught ("bad location for attribute " ++ attributeName) blob - /dev/null blob + e8395223e9db5294f484a194568935f44fc5e321 (mode 644) --- /dev/null +++ src/S3D/GLProgram/SimpleColour/Attributes.idr @@ -0,0 +1,16 @@ +||| We put the SimpleColour.Attributes type in a separate module so it's +||| possible to use it in non-Javascript code. + +module S3D.GLProgram.SimpleColour.Attributes + +import Math.LinearAlgebra + +%default total + +-- Wrap it in a namespace so we can write SimpleColour.Attributes to refer to +-- the Attributes type. +namespace SimpleColour + public export + record Attributes where + constructor MkAttributes + position, colour : Vector 4 Double blob - /dev/null blob + 1c03b5aa5791de846029534c091d11d15728a094 (mode 644) --- /dev/null +++ src/S3D/GLProgram/SimpleColour.idr @@ -0,0 +1,87 @@ +module S3D.GLProgram.SimpleColour + +import public S3D.GLProgram.SimpleColour.Attributes + +import Data.Vect +import JS +import S3D.DrawInput +import S3D.Drawable +import S3D.GLProgram +import S3D.GLProgram.Helpers +import Math.LinearAlgebra +import S3D.Transformable +import Web.Html +import Web.Raw.Webgl + +%default total + +public export +record Program where + constructor MkProgram + program : WebGLProgram + viewMatrixLocation : WebGLUniformLocation + positionLocation, colourLocation : Bits32 + +export +Transformable Attributes where + transform matrix = { position $= transform matrix } + +vertexShaderSource : String +vertexShaderSource = + "#version 300 es\n" ++ + " " ++ + "in vec4 a_position; " ++ + "uniform mat4 u_matrix; " ++ + " " ++ + "in vec4 a_colour; " ++ + "out vec4 v_colour; " ++ + " " ++ + "void main() { " ++ + " gl_Position = u_matrix * a_position; " ++ + " v_colour = a_colour; " ++ + "}" + +fragmentShaderSource : String +fragmentShaderSource = + "#version 300 es\n" ++ + " " ++ + "precision mediump float; " ++ + " " ++ + "in vec4 v_colour; " ++ + "out vec4 colour; " ++ + "void main() { " ++ + " colour = v_colour; " ++ + "}" + +export +GLProgram Program Attributes where + makeProgram gl = + do program <- compileAndLinkProgram gl vertexShaderSource fragmentShaderSource + matrixLocation <- getAndCheckUniformLocation gl program "u_matrix" + positionLocation <- getAndCheckAttribLocation gl program "a_position" + colourLocation <- getAndCheckAttribLocation gl program "a_colour" + pure $ + MkProgram + { program = program + , viewMatrixLocation = matrixLocation + , positionLocation = positionLocation + , colourLocation = colourLocation + } + + compileDrawInput gl program drawInput = + do vertexArray <- unMaybe "createVertexArray" $ createVertexArray gl + bindVertexArray gl (Just vertexArray) + positionBuffer <- makeAndBindBuffer gl program.positionLocation + (map (vectorToVect . position) drawInput.vertexAttributes) + colourBuffer <- makeAndBindBuffer gl program.colourLocation + (map (vectorToVect . colour) drawInput.vertexAttributes) + elementBuffer <- makeAndBindElementBuffer gl drawInput.triangles + let numElements = 3 * length drawInput.triangles + pure $ + MkDrawable + { program = program.program + , viewMatrixLocation = program.viewMatrixLocation + , vertexArray = vertexArray + , numElements = cast numElements + , buffers = [positionBuffer, colourBuffer, elementBuffer] + } blob - /dev/null blob + 2f6ad104f7531f276f6389518f434f79219df001 (mode 644) --- /dev/null +++ src/S3D/GLProgram.idr @@ -0,0 +1,22 @@ +module S3D.GLProgram + +import JS +import S3D.Drawable +import S3D.DrawInput +import Web.Html +import Web.Raw.Webgl + +%default total + +-- An implementation of GLProgram is a program (i.e. collection of shaders) +-- together with types and values related to it. +-- +-- The program type represents an instance of the program, including the +-- locations of attributes. (There's probably no reason to ever make more +-- than one instance of the same program.) +-- +-- The attributes type represents the attributes for a single vertex. +export +interface GLProgram program attributes | program where + makeProgram : WebGL2RenderingContext -> JSIO program + compileDrawInput : WebGL2RenderingContext -> program -> DrawInput attributes -> JSIO Drawable blob - /dev/null blob + 26ab0b4a83ca9b2bc4e673d804bd1c8ef3e00c1d (mode 644) --- /dev/null +++ src/S3D/Lighting/SimpleColour.idr @@ -0,0 +1,15 @@ +module S3D.Lighting.SimpleColour + +import Data.Vect +import Math.LinearAlgebra +import S3D.GLProgram.SimpleColour +import S3D.Lighting + +%default total + +export +UseLightColour Attributes where + useLightColour lightColour attributes = + case (vectorToVect lightColour, vectorToVect attributes.colour) of + ([r, g, b], [r', g', b', a]) => + { colour := [r * r', g * g', b * b', a] } attributes blob - /dev/null blob + 93e90c6a7771c75d69e302964e868b70c69707a6 (mode 644) --- /dev/null +++ src/S3D/Lighting.idr @@ -0,0 +1,97 @@ +||| Primitive lighting: a face on a polyhedron is illuminated depending on its +||| relation to some light sources. +||| +||| Unprincipled and subject to change. + +module S3D.Lighting + +import Data.List +import Data.Vect +import Math.LinearAlgebra +import Math.LinearAlgebra.Normalized +import S3D.DrawablePolyhedron +import S3D.PhysicalPolyhedron +import SolidGeometry.Polyhedron + +%default total + +public export +record LightSource where + constructor MkLightSource + location : Vector 4 Double + colour : Vector 3 Double + +interface UseLightColour attributes where + ||| Modify the input so its total illumination is `lightColour`. + useLightColour : (lightColour : Vector 3 Double) -> attributes -> attributes + +||| `Nothing` means the same as `Just [0, 0, 0]`, but might save some +||| computation. +faceIlluminationFromLight : (faceNormal : Vector 4 Double) -> + (locationOnFace : Normalized (Vector 4 Double)) -> + LightSource -> Maybe (Vector 3 Double) +faceIlluminationFromLight faceNormal locationOnFace lightSource = + let -- An ad-hoc function. Should get dimmer as the light gets further away + -- or if `faceNormal` isn't pointing at the light. + brightness : Double + brightness = + normalize faceNormal ^.^ + normalize (normalize lightSource.location ^-^ fromNormalized locationOnFace) + in if brightness >= 0 + then Just $ brightness *^ lightSource.colour + else Nothing + +||| `faceNormal` should be pointed so that points `p` with +||| `p ^.^ faceNormal > 0` are on the side of the face that will be drawn +||| (should be the outside). +illuminateFace : UseLightColour attributes => + List LightSource -> (faceNormal : Vector 4 Double) -> + AttributeBuilder attributes -> AttributeBuilder attributes +illuminateFace {attributes} + lightSources faceNormal attributeBuilder + vertexCoordinates = + let totalLight = + sumVectorsFrom [0, 0, 0] $ + mapMaybe (faceIlluminationFromLight faceNormal vertexCoordinates.absolute) + lightSources + in useLightColour totalLight $ attributeBuilder vertexCoordinates + +illuminateFacePair : UseLightColour attributes => + List LightSource -> + (plane : Matrix 3 4 Double) -> + (faceAttributes : Vect 2 (AttributeBuilder attributes)) -> + Vect 2 (AttributeBuilder attributes) +illuminateFacePair + lightSources plane [positiveFaceAttributes, negativeFaceAttributes] = + let planeNormal = cross $ matrixToVectors plane + in [ illuminateFace lightSources (negateVector planeNormal) + positiveFaceAttributes + , illuminateFace lightSources planeNormal + negativeFaceAttributes + ] + +export +interface Illuminate t where + illuminate : List LightSource -> t -> t + +export +UseLightColour attributes => + -- This should also work with the `pretriangulated` type argument + -- set to False, but the compiler complains if we replace `True` + -- with a variable here. + Illuminate (DrawablePolyhedron numplanes attributes True) +where + illuminate lightSources polyhedron = + { faceAttributes := + zipWith (illuminateFacePair lightSources) + polyhedron.polyhedron.planes + polyhedron.faceAttributes + } + polyhedron + +export +UseLightColour attributes => + Illuminate (PhysicalPolyhedron numplanes attributes True) +where + illuminate lightSources = + {drawablePolyhedron $= illuminate lightSources} blob - /dev/null blob + 0544d80b4cb390b1d476345bedd1e926a3dcd572 (mode 644) --- /dev/null +++ src/S3D/Numbers.idr @@ -0,0 +1,9 @@ +module S3D.Numbers + +import Data.Fin + +%default total + +export +finToInt16 : Fin n -> Int16 +finToInt16 = cast . finToInteger blob - /dev/null blob + 4b496cf0db5dec9d789e153c0e2a4724924d8bea (mode 644) --- /dev/null +++ src/S3D/Obstacles.idr @@ -0,0 +1,78 @@ +module S3D.Obstacles + +import Data.Fin +import Data.Linear.SizedArray +import Data.Vect +import Math.LinearAlgebra +import S3D.Transformable +import SolidGeometry.HalfSpace +import SolidGeometry.Polyhedron +import SolidGeometry.Sign + +%default total + +||| A `Polyhedron`, except: +||| - `numPlanes` is a data member, not part of the type. +||| - The planes are given with an "implicit" representation --- i.e. a vector +||| you take a dot product with. +record ObstaclePolyhedron where + constructor MkObstaclePolyhedron + {numPlanes : Nat} + planes : Vect numPlanes (Vector 4 Double) + inside : (Fin numPlanes -> Sign) -> InOutSide + +Transformable ObstaclePolyhedron where + transform matrix = {planes $= map (matrixInverse matrix !*^)} + +convertPolyhedron : {numPlanes : Nat} -> Polyhedron numPlanes -> ObstaclePolyhedron +convertPolyhedron polyhedron = + MkObstaclePolyhedron + { numPlanes = _ + , planes = map (negateVector . cross . matrixToVectors) polyhedron.planes + , inside = polyhedron.inside + } + +pointInsidePolyhedron : Vector 4 Double -> ObstaclePolyhedron -> InOutSide +pointInsidePolyhedron position polyhedron = + let sides : SizedIArray polyhedron.numPlanes Sign + sides = fromVect + [ softSign $ position ^.^ plane + | plane <- polyhedron.planes + ] + in polyhedron.inside (read sides) + +export +data Obstacles = MkObstacles (List ObstaclePolyhedron) + +export +Semigroup Obstacles where + MkObstacles x <+> MkObstacles y = MkObstacles (x <+> y) + +export +Monoid Obstacles where + neutral = MkObstacles neutral + +export +Transformable Obstacles where + transform matrix (MkObstacles polyhedra) = + MkObstacles $ map (transform matrix) polyhedra + +export +makeObstacles : List (numPlanes : Nat ** Polyhedron numPlanes) -> Obstacles +makeObstacles polyhedra = + MkObstacles + [ convertPolyhedron polyhedron + | (_ ** polyhedron) <- polyhedra + ] + +export +collision' : Vector 4 Double -> List ObstaclePolyhedron -> Bool +collision' _ Nil = False +collision' position (h :: t) = + case pointInsidePolyhedron position h of + Inside => True + Outside => collision' position t + +export +collision : Vector 4 Double -> Obstacles -> Bool +collision v (MkObstacles o) = collision' v o blob - /dev/null blob + 85e96c5d59d45f6f871d1ac6baea980bd26628e1 (mode 644) --- /dev/null +++ src/S3D/PhysicalObject.idr @@ -0,0 +1,22 @@ +module S3D.PhysicalObject + +import S3D.Figure +import S3D.Obstacles + +%default total + +||| A PhysicalObject is a Figure for drawing together with an Obstacle for +||| collision detection. +public export +record PhysicalObject where + constructor MkPhysicalObject + figure : Figure + obstacles : Obstacles + +export +concatPhysicalObjects : List PhysicalObject -> PhysicalObject +concatPhysicalObjects objects = + MkPhysicalObject + { figure = concatFigures $ map figure objects + , obstacles = concat $ map obstacles objects + } blob - /dev/null blob + 9ecfab44b6795f9373a5bea29b68632052ef25b0 (mode 644) --- /dev/null +++ src/S3D/PhysicalPolyhedron.idr @@ -0,0 +1,32 @@ +module S3D.PhysicalPolyhedron + +import S3D.DrawablePolyhedron +import SolidGeometry.Polyhedron + +%default total + +||| A `PhysicalPolyhedron attributes` is a `DrawablePolyhedron attributes` to +||| draw, together with a `Polyhedron` to use for collision detection. +||| +||| We expose `drawableNumPlanes` as part of the type to make it easier to use +||| `S3D.DrawablePolyhedron.drawPrecomputedPolyhedron` without having an +||| argument with the type checker about whether the number of planes matches. +||| +||| Unlike a `PhysicalObject`, it is possible to perform constructive solid +||| geometry operations on `PhysicalPolyhedron` values. The final result can be +||| converted to a `PhysicalObject`. +public export +record PhysicalPolyhedron (drawableNumPlanes : Nat) (attributes : Type) + (pretriangulated : Bool) +where + constructor MkPhysicalPolyhedron + drawablePolyhedron : DrawablePolyhedron drawableNumPlanes attributes + pretriangulated + {obstacleNumPlanes : Nat} + obstaclePolyhedron : Polyhedron obstacleNumPlanes + +export +addTriangulation : List (TriangulatedFace numPlanes) -> + PhysicalPolyhedron numPlanes attributes False -> + PhysicalPolyhedron numPlanes attributes True +addTriangulation t = {drawablePolyhedron $= addTriangulation t} blob - /dev/null blob + f8d02b74d83c087ff700ffce4432d3334cf3db13 (mode 644) --- /dev/null +++ src/S3D/PhysicalState.idr @@ -0,0 +1,95 @@ +module S3D.PhysicalState + +import Data.IORef +import S3D.Controls +import Math.LinearAlgebra +import S3D.Obstacles +import S3D.PlayerMovement +import S3D.WebUI +import System.Clock + +%default total + +-- How often to update the state. +updatePeriod : Clock Duration +updatePeriod = MkClock 0 100000000 -- 10 fps + +-- A time really far in the past. +reallyEarly : {type : ClockType} -> Clock type +reallyEarly = MkClock (-1000000000) 0 + +toSeconds : Clock Duration -> Double +toSeconds d = fromInteger (seconds d) + fromInteger (nanoseconds d) / 1000000000.0 + +divideDurations : Clock Duration -> Clock Duration -> Double +divideDurations d e = toSeconds d / toSeconds e + +public export +record PhysicalState where + constructor MkPhysicalState + -- The last time physics was computed. + lastUpdateTime : Clock Monotonic + -- The next time it should be computed. Always lastUpdateTime + updatePeriod, except for initialState. + nextUpdateTime : Clock Monotonic + -- The player's state at lastUpdateTime. + lastPlacement : PlayerPlacement + -- The player's state at the next physics frame. Between physics frames, we + -- interpolate. + nextPlacement : PlayerPlacement + -- The corresponding view matrices. + lastViewMatrix, nextViewMatrix : Matrix 4 4 Double + +export +initialState : PhysicalState +initialState = + MkPhysicalState + { lastUpdateTime = reallyEarly + , nextUpdateTime = reallyEarly + , lastPlacement = initialPlacement + , nextPlacement = initialPlacement + , lastViewMatrix = playerViewMatrix initialPlacement + , nextViewMatrix = playerViewMatrix initialPlacement + } + +-- time must be after physicalState.lastUpdateTime. Assuming that, the +-- lastUpdateTime and nextUpdateTime fields of the output are guaranteed to be +-- before (or equal to) and after (or equal to) time. +updatePhysicalState : + HasIO io => + Obstacles -> Clock Monotonic -> IORef UIState -> PhysicalState -> + io PhysicalState +updatePhysicalState obstacles time uiStateRef physicalState = + if time < physicalState.nextUpdateTime + then pure physicalState + else + do nextPlacement' <- stepPlacement obstacles physicalState.nextPlacement + uiStateRef + pure $ + let -- Ideally, we add updatePeriod to both times. But if time would + -- be outside that interval (meaning we're too slow to keep up), + -- we have to add more. + nextUpdateTime' = max time (addDuration physicalState.nextUpdateTime updatePeriod) + lastUpdateTime' = subtractDuration nextUpdateTime' updatePeriod + in MkPhysicalState + lastUpdateTime' + nextUpdateTime' + physicalState.nextPlacement + nextPlacement' + physicalState.nextViewMatrix + (playerViewMatrix nextPlacement') + +-- Returns the new state and a view matrix corresponding to the given +-- timestamp. Between physics updates, the view matrix is interpolated. +export +stepState : + HasIO io => + Obstacles -> Clock Monotonic -> IORef UIState -> PhysicalState -> + io (Pair PhysicalState (Matrix 4 4 Double)) +stepState obstacles time uiStateRef physicalState = + do physicalState' <- updatePhysicalState obstacles time uiStateRef physicalState + pure $ + let progressThroughUpdate = divideDurations (timeDifference time physicalState'.lastUpdateTime) updatePeriod + in ( physicalState' + , (1 - progressThroughUpdate) *! physicalState'.lastViewMatrix + !+! progressThroughUpdate *! physicalState'.nextViewMatrix + ) blob - /dev/null blob + 34c8cde0d11f7f527852eca497bf5cb0ccd096b6 (mode 644) --- /dev/null +++ src/S3D/Player.idr @@ -0,0 +1,8 @@ +module S3D.Player + +%default total + +||| The distance from the floor to the player's eyes, in radians. +export +playerHeight : Double +playerHeight = pi * 3 / 32 blob - /dev/null blob + 0bf2a1fce63a1059a5f1e1d7fec8d4c44e1e897a (mode 644) --- /dev/null +++ src/S3D/PlayerMovement.idr @@ -0,0 +1,191 @@ +{- + +Functions to help compute the player's position and orientation. + +Naming conventions for a movement M (e.g. M = MoveForward), +* oneFrameM means do movement M an appropriate amount for one frame of a + corresponding key being held down. +* touchM means do movement M an appropriate amount given a magnitude of a touch + gesture. + +The -y pole is always down. + +-} + +module S3D.PlayerMovement + +import Data.Vect +import Math.LinearAlgebra +import S3D.Player + +%default total + +{- + +Constants + +-} + +moveStep, turnStepAngle, upDownLevelStep : Double + +-- How far the player moves in one frame, in radians. +moveStep = pi / 60 + +-- How much the player turns in one frame, in radians. +turnStepAngle = pi / 15 + +-- How much the player can turn their head up or down in one frame, in radians. +upDownStep = pi / 20 + + +-- The player's position and orientation. +export +record PlayerPlacement where + constructor MkPlayerPlacement + -- The player's position. Should have unit length. + position : Vector 4 Double + -- The direction the player's facing, not counting looking up or down. Should + -- have unit length and be perpendicular to both position and the y axis. + forward : Vector 4 Double + -- To the player's right. Computed from position and forward in + -- normalizePlacement. + right : Vector 4 Double + -- pi/2 if the player is looking straight up, -pi/2 if they're looking stright + -- down. + upDown : Double + +||| The player's position +export +position : PlayerPlacement -> Vector 4 Double +position player = player.position + +-- makeOrtho u v returns a vector orthogonal to u. u must have norm 1. +makeOrtho : (Neg coord, Num coord) => Vector n coord -> Vector n coord -> Vector n coord +makeOrtho u v = v ^-^ (u ^.^ v) *^ u + +-- Assumes placement.position is already normalized. +normalizePlacement : PlayerPlacement -> PlayerPlacement +normalizePlacement placement = + { forward := normalize (makeOrtho [0, 1, 0, 0] (makeOrtho placement.position placement.forward)) + , right := normalize $ cross [placement.position, placement.forward, [0, -1, 0, 0]] + } + placement + +-- Adjust the player's position so they are playerHeight radians above the ground. +-- +-- Input need not be normalized. Output will be normalized. +adjustHeight : Vector 4 Double -> Vector 4 Double +adjustHeight v = + let [x, y, z, w] = vectorToVect v + nonYNorm = norm $ the (Vector 3 Double) [x, z, w] + nonYMult = cos playerHeight / nonYNorm + in [x * nonYMult, sin playerHeight, z * nonYMult, w * nonYMult] + +export +initialPlacement : PlayerPlacement +initialPlacement = normalizePlacement $ + MkPlayerPlacement + { -- We avoid setting coordinates to zero to avoid an edge case where points + -- end with w=0 in view coordinates. + position = (adjustHeight [0.001, 0, -1, 0.0001]) + , forward = [0, 0, 0, 1] + , right = [1, 0, 0, 0] + , upDown = 0 + } + +-- Move the player the given distance (in radians) in the given direction. +-- Direction should be unit length. If the player was looking up or down, they +-- level out a bit. +move : Double -> Vector 4 Double -> PlayerPlacement -> PlayerPlacement +move distance direction placement = + normalizePlacement $ + { position := adjustHeight $ + cos distance *^ placement.position ^+^ + sin distance *^ direction + } + placement + +-- Move the player forward (sign=1) or backward (sign=-1). +moveForward : Double -> PlayerPlacement -> PlayerPlacement +moveForward amount placement = move amount placement.forward placement + +export +oneFrameMoveForward : Double -> PlayerPlacement -> PlayerPlacement +oneFrameMoveForward sign = moveForward (sign * moveStep) + +export +touchMoveForward : Double -> PlayerPlacement -> PlayerPlacement +touchMoveForward = moveForward + +-- Move the player right or left. +strafe : Double -> PlayerPlacement -> PlayerPlacement +strafe amount placement = move amount placement.right placement + +-- Move the player right (sign=1) or left (sign=-1). +export +oneFrameStrafe : Double -> PlayerPlacement -> PlayerPlacement +oneFrameStrafe sign = strafe (sign * moveStep) + +-- Move the player right or left. +export +touchStrafe : Double -> PlayerPlacement -> PlayerPlacement +touchStrafe = strafe + +-- Turn the player. +turn : Double -> PlayerPlacement -> PlayerPlacement +turn angle placement = + normalizePlacement $ + { forward := cos angle *^ placement.forward ^+^ + sin angle *^ placement.right + } + placement + +-- Turn the player right (sign=1) or left (sign=-1). +export +oneFrameTurn : Double -> PlayerPlacement -> PlayerPlacement +oneFrameTurn sign = turn (sign * turnStepAngle) + +-- Turn by an amount appropriate for a touch gesture. +export +touchTurn : Double -> PlayerPlacement -> PlayerPlacement +touchTurn = turn + +-- The player looks up or down. +lookUp : Double -> PlayerPlacement -> PlayerPlacement +lookUp angle = + { upDown $= \ ud => min (pi/2) (max (-pi/2) (ud + angle)) + } + +-- The player looks up (sign=1) or down (sign=-1). +export +oneFrameLookUp : Double -> PlayerPlacement -> PlayerPlacement +oneFrameLookUp sign = lookUp (sign * upDownStep) + +-- Look up or down by an amount appropriate for a touch gesture. +export +touchLookUp : Double -> PlayerPlacement -> PlayerPlacement +touchLookUp = lookUp + +up : PlayerPlacement -> Vector 4 Double +up placement = + normalize $ cross [placement.right, placement.position, placement.forward] + +export +playerViewMatrix : PlayerPlacement -> Matrix 4 4 Double +playerViewMatrix placement = + -- The inverse of the matrix with columns (right, up, -position, forward). + -- Since everything's orthonormal, the transpose is the inverse. + let placementUp = up placement + cosUpDown = cos placement.upDown + sinUpDown = sin placement.upDown + in + transpose + [ placement.right + , cosUpDown *^ placementUp ^-^ sinUpDown *^ placement.forward + , negateVector placement.position + , cosUpDown *^ placement.forward ^+^ sinUpDown *^ placementUp + ] + +export +Show PlayerPlacement where + show placement = "pos: " ++ show placement.position ++ "; fwd: " ++ show placement.forward ++ "; right: " ++ show placement.right ++ "; up: " ++ show (up placement) ++ "; ud: " ++ show placement.upDown blob - /dev/null blob + c4a6a4bc029b42960e74970423ba0af69a1f1aec (mode 644) --- /dev/null +++ src/S3D/Polyhedra/Drawable.idr @@ -0,0 +1,20 @@ +module S3D.Polyhedra.Drawable + +import Data.Vect +import Math.LinearAlgebra +import S3D.GLProgram.SimpleColour.Attributes +import S3D.DrawablePolyhedron +import S3D.DrawablePolyhedron.SimpleColour +import S3D.Polyhedra + +%default total + +export +monochromeCube : Vector 4 Double -> DrawablePolyhedron ? SimpleColour.Attributes False +monochromeCube colour = + makeSimpleColourPolyhedron cube (replicate _ (replicate 2 colour)) + +export +monochromeTower : Vector 4 Double -> DrawablePolyhedron ? SimpleColour.Attributes False +monochromeTower colour = + makeSimpleColourPolyhedron tower (replicate _ (replicate 2 colour)) blob - /dev/null blob + cd0e8687288f983ae2232890c5be76cd0292de21 (mode 644) --- /dev/null +++ src/S3D/Polyhedra/Precompute.idr @@ -0,0 +1,109 @@ +module S3D.Polyhedra.Precompute + +import Data.Fin +import Data.SortedMap +import Data.Vect +import Math.LinearAlgebra +import Math.LinearAlgebra.Normalized +import SolidGeometry.Polygon +import SolidGeometry.Polyhedron +import SolidGeometry.Sign + +%default total + +encodeList : (a -> String) -> List a -> String +encodeList encodeItem list = + fastConcat [encodeItem item <+> " :: " | item <- list] <+> "Nil" + +encodeFaceId : FaceId numPlanes -> String +encodeFaceId faceId = + fastConcat + [ "MkFaceId { planeIndex = " + , show faceId.planeIndex + , ", outside = " + , show faceId.outside + , "}" + ] + +encodeVertexCoordinates : VertexCoordinates -> String +encodeVertexCoordinates c = + fastConcat + [ "MkVertexCoordinates { absolute = makeNormalized " + , show $ vectorToVect $ fromNormalized $ c.absolute + , ", faceCoordinates = makeNormalized " + , show $ vectorToVect $ fromNormalized $ c.faceCoordinates + , "}" + ] + +encodeTriangle : {0 numVertices : Nat} -> Triangle (Fin numVertices) -> String +encodeTriangle (MkTriangle a b c) = + fastConcat + [ "MkTriangle " + , show a + , " " + , show b + , " " + , show c + ] + +encodeTriangulatedFace : TriangulatedFace numPlanes -> String +encodeTriangulatedFace face = + fastConcat + [ "MkTriangulatedFace { faceId = " + , encodeFaceId face.faceId + , ", vertices = " + , encodeList encodeVertexCoordinates $ toList face.vertices + , ", triangles = " + , encodeList encodeTriangle face.triangles + , "}" + ] + +encodeTriangulation : List (TriangulatedFace numPlanes) -> String +encodeTriangulation = encodeList encodeTriangulatedFace + +partial +precomputeOnePolyhedron : (name : String) -> + (numPlanes : Nat) -> + (polyhedron : Polyhedron numPlanes) -> + Either String String +precomputeOnePolyhedron name numPlanes polyhedron = + pure $ + "export\n" <+> + name <+> " : List (TriangulatedFace " <+> show numPlanes <+> ")\n" <+> + name <+> " = " <+> encodeTriangulation !(triangulate polyhedron) <+> "\n" + +||| Returns Idris source for a module named `moduleName` defining a `List +||| (TriangulatedFace numPlanes)` for each of the polyhedra in the map, using +||| the map keys as identifiers. +||| +||| This allows triangulations to be precomputed and baked into a compiled +||| program. +export partial +precomputePolyhedra : + (polyhedra : SortedMap String (numPlanes : Nat ** Polyhedron numPlanes)) -> + (moduleName : String) -> + Either String String +precomputePolyhedra polyhedra moduleName = + do precomputedPolyhedra <- sequence + [ precomputeOnePolyhedron name numPlanes polyhedron + | (name, (numPlanes ** polyhedron)) <- SortedMap.toList polyhedra + ] + pure $ + fastConcat $ + [ "module " + , moduleName + , "\n\n" + ] + <+> + [ "import " <+> m <+> "\n" + | m <- [ "Data.Fin", "Data.Vect", "Math.LinearAlgebra" + , "Math.LinearAlgebra.Normalized" + , "SolidGeometry.Polygon", "SolidGeometry.Polyhedron" + , "SolidGeometry.Sign" + ] + ] + <+> + [ "\n%default total\n\n" + ] + <+> + precomputedPolyhedra blob - /dev/null blob + ad91fd043c1912ba69b05ce0f38ec9f466e72fa1 (mode 644) --- /dev/null +++ src/S3D/Polyhedra/ToPrecompute.idr @@ -0,0 +1,42 @@ +||| A collection of polyhedra to precompute. + +module S3D.Polyhedra.ToPrecompute + +import Data.SortedMap +import S3D.DrawablePolyhedron +import S3D.PhysicalPolyhedron +import S3D.GLProgram.SimpleColour.Attributes +import S3D.Polyhedra +import S3D.Scenes.House.ToPrecompute.Kitchen +import S3D.Scenes.House.ToPrecompute.LivingRoom +import S3D.Scenes.House.ToPrecompute.Walls +import SolidGeometry.Polyhedron + +%default total + +export partial +toPrecompute : SortedMap String (numPlanes : Nat ** Polyhedron numPlanes) +toPrecompute = + fromList + [ ( "cube" + , (_ ** cube) + ) + , ( "houseCouch" + , (_ ** couch.drawablePolyhedron.polyhedron) + ) + , ( "houseFrontWall" + , (_ ** frontWall.drawablePolyhedron.polyhedron) + ) + , ( "houseKitchenCounter" + , (_ ** counter.drawablePolyhedron.polyhedron) + ) + , ( "wallWithDoor" + , (_ ** wallWithDoor.drawablePolyhedron.polyhedron) + ) + , ( "wallWithOneWindow" + , (_ ** wallWithOneWindow.drawablePolyhedron.polyhedron) + ) + , ( "wallWithWindows" + , (_ ** wallWithWindows.drawablePolyhedron.polyhedron) + ) + ] blob - /dev/null blob + cbc474937f7c29152ec108768d96fc6dfff8e791 (mode 644) --- /dev/null +++ src/S3D/Polyhedra.idr @@ -0,0 +1,108 @@ +-- Some polyhedra. + +module S3D.Polyhedra + +import Data.Fin +import Data.Nat +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.HalfSpace +import SolidGeometry.Polyhedron +import SolidGeometry.Sign + +%default total + +0 lteAddLeft : (n, m : Nat) -> n `LTE` m + n +lteAddLeft n m = plusLteMonotoneRight n Z m LTEZero + +allNegativeFrom : (n : Nat) -> (Fin (m + n) -> Sign) -> Bool +allNegativeFrom Z _ = True +allNegativeFrom (S n') f = + let 0 nLteMPlusN : n `LTE` m + n + nLteMPlusN = lteAddLeft n m + in case f (natToFinLT n' {prf = nLteMPlusN}) of + Positive => False + Negative => allNegativeFrom {m = S m} n' $ + rewrite plusSuccRightSucc m n' in f + +allNegative : {n : Nat} -> (Fin n -> Sign) -> Bool +allNegative = allNegativeFrom {m = 0} n + +||| Intended to be used as the `inside` field for a Polyhedron. Says you're +||| inside iff you're on the `Negative` side of all faces. +planeIntersection : {numPlanes : Nat} -> (Fin numPlanes -> Sign) -> InOutSide +planeIntersection h = + if allNegative h then Inside else Outside + +||| The half-space of points whose y coordinate is positive. [0, 0, 1] in the +||| plane's coordinate system maps to the w pole. +export +positiveY : Polyhedron 1 +positiveY = + MkPolyhedron + { planes = [[[1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]] + , inside = \sides => + case sides FZ of + Positive => Inside + Negative => Outside + } + +||| A tetrahedron centred on [0, 0, 0, 1]. The closer to 0 the input is, the +||| larger the tetrahedron. +export +tetrahedron : Double -> Polyhedron 4 +tetrahedron w = + let a = [ 1, 1, 1, w] in + let b = [ 1, -1, -1, w] in + let c = [-1, 1, -1, w] in + let d = [-1, -1, 1, w] in + MkPolyhedron + [[a, c, b], [a, b, d], [a, d, c], [b, c, d]] + planeIntersection + +namespace Cube + a : Vector 4 Double + a = [1, 0, 0, 0] + b : Vector 4 Double + b = [0, 1, 0, 0] + c : Vector 4 Double + c = [0, 0, 1, 0] + + wVector : Vector 4 Double + wVector = [0, 0, 0, 1] + wMatrix : Matrix 3 4 Double + wMatrix = [wVector, wVector, wVector] + + ||| Given an order [d, e, f] of the three axes, returns two planes + ||| corresponding to d and -d. + sidesForAxes : Vect 3 (Vector 4 Double) -> Vect 2 (Matrix 3 4 Double) + sidesForAxes [d, e, f] = + let -- Three points on one of the sides, without the w part. + g : Vector 4 Double + g = d ^+^ e ^+^ f + h : Vector 4 Double + h = d ^+^ negateVector e ^+^ f + i : Vector 4 Double + i = d ^+^ e ^+^ negateVector f + in [[g, i, h] !+! wMatrix, negateMatrix [g, h, i] !+! wMatrix] + + ||| A cube with corners at [+-1, +-1, +-1, 1]. + export + cube : Polyhedron 6 + cube = + let axisOrders : Vect 3 (Vect 3 (Vector 4 Double)) + axisOrders = [[a, b, c], [b, c, a], [c, a, b]] + sides : Vect 6 (Matrix 3 4 Double) + sides = concat $ map sidesForAxes axisOrders + in MkPolyhedron sides planeIntersection + + ||| Like cube, but ignoring y coordinates, so it extends from the -y pole to + ||| the y pole. + export + tower : Polyhedron 4 + tower = + let axisOrders : Vect 2 (Vect 3 (Vector 4 Double)) + axisOrders = [[a, b, c], [c, a, b]] + sides : Vect 4 (Matrix 3 4 Double) + sides = concat $ map sidesForAxes axisOrders + in MkPolyhedron sides planeIntersection blob - /dev/null blob + a0b48c2c4e0730a3a94052f2505a15d9eca8a304 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/DiningRoom.idr @@ -0,0 +1,65 @@ +module S3D.Scenes.House.DiningRoom + +import Data.Vect +import S3D.DrawablePolyhedron +import S3D.DrawablePolyhedron.SimpleColour +import S3D.GLProgram.SimpleColour +import Math.LinearAlgebra +import S3D.PhysicalPolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Precomputed +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron +import SolidGeometry.Polyhedron + +%default total + +export +monochromeCube : Vector 4 Double -> DrawablePolyhedron ? SimpleColour.Attributes True +monochromeCube colour = + addTriangulation Precomputed.cube $ + makeSimpleColourPolyhedron cube (replicate _ (replicate 2 colour)) + +top : PhysicalPolyhedron ? SimpleColour.Attributes True +top = + transform + ( diagonalMatrix [0.3, 0.01, 0.2, 1] + !*! + translate [0, 0.2, 0.4] + ) + $ + MkPhysicalPolyhedron + { drawablePolyhedron = monochromeCube [0.6, 0.4, 0.3, 1] + , obstaclePolyhedron = tower + } + +leg : PhysicalPolyhedron ? SimpleColour.Attributes True +leg = + transform + ( translate [0, 1, 0] + !*! + diagonalMatrix [0.02, 0.1, 0.02, 1] + ) + $ + MkPhysicalPolyhedron + { drawablePolyhedron = monochromeCube [0.8, 0.7, 0.5, 1] + , obstaclePolyhedron = emptyPolyhedron + } + +export +table : List ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +table = + [ (_ ** top) + , (_ ** transform (translate [ 0.275, 0, 0.575]) leg) + , (_ ** transform (translate [ 0.275, 0, 0.225]) leg) + , (_ ** transform (translate [-0.275, 0, 0.575]) leg) + , (_ ** transform (translate [-0.275, 0, 0.225]) leg) + ] + +export +diningRoom : List ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +diningRoom = table blob - /dev/null blob + bb91f63e172fb3d9aac78202107e5c21ff2260b3 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/FloorLines.idr @@ -0,0 +1,49 @@ +{- + +The floorLines function returns a Figure with bands of colour along the y=0 +plane (the "floor") connecting the +w and -w poles. + +-} + +module S3D.Scenes.House.FloorLines + +import Data.Vect +import S3D.DrawInput +import S3D.Figure +import S3D.GLProgram.SimpleColour +import Math.LinearAlgebra + +%default total + +finToDouble : Fin n -> Double +finToDouble = fromInteger . finToInteger + +floorBand : Double -> Double -> Vector 4 Double -> DrawInput Attributes +floorBand startAngle endAngle colour = + MkDrawInput + [T [0, 2, 3], T [1, 3, 2]] + [ MkAttributes + [0, 0, 0, 1] + colour + , MkAttributes + [0, 0, 0, -1] + colour + , MkAttributes + [cos startAngle, 0, sin startAngle, 0] + colour + , MkAttributes + [cos endAngle, 0, sin endAngle, 0] + colour + ] + +-- Input: the colours of the bands. +export +floorLines : {n : Nat} -> Vect n (Vector 4 Double) -> Figure +floorLines colours = + simpleColourFigure $ + concatDrawInputs $ toList $ + zipWith + ( \i, colour => + floorBand (2 * pi * finToDouble i / cast n) (2 * pi * finToDouble (FS i) / cast n) colour + ) + Data.Vect.Fin.range colours blob - /dev/null blob + 817380a225ac515783e20ad0eb4f2363838eccd4 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/Kitchen.idr @@ -0,0 +1,36 @@ +module S3D.Scenes.House.Kitchen + +import Data.Vect +import Math.LinearAlgebra +import S3D.GLProgram.SimpleColour +import S3D.DrawablePolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Drawable +import S3D.Polyhedra.Precomputed +import S3D.PhysicalPolyhedron +import S3D.Scenes.House.ToPrecompute.Kitchen +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron + +%default total + +fridge : PhysicalPolyhedron ? SimpleColour.Attributes True +fridge = + transform + ( diagonalMatrix [0.15, 0.3, 0.1, 1] !*! + translate [0.65, 0.3, -0.75] + ) + $ + MkPhysicalPolyhedron + { drawablePolyhedron = + addTriangulation Precomputed.cube $ + monochromeCube [0.9, 0.9, 0.9, 1] + , obstaclePolyhedron = tower + } + +export +kitchen : List (numPlanes : Nat ** PhysicalPolyhedron numPlanes SimpleColour.Attributes True) +kitchen = + [ (_ ** addTriangulation houseKitchenCounter counter) + , (_ ** fridge) + ] blob - /dev/null blob + 954bd75eedf924eac6143b3216aae8f3c0c780fa (mode 644) --- /dev/null +++ src/S3D/Scenes/House/LivingRoom.idr @@ -0,0 +1,23 @@ +module S3D.Scenes.House.LivingRoom + +import Data.Vect +import Math.LinearAlgebra +import S3D.GLProgram.SimpleColour +import S3D.PhysicalPolyhedron +import S3D.Polyhedra.Precomputed +import S3D.Scenes.House.ToPrecompute.LivingRoom +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron + +%default total + +export +livingRoom : List ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +livingRoom = + [ ( _ ** + transform (translate [0, 0, -0.68]) + (addTriangulation houseCouch couch) + ) + ] blob - /dev/null blob + 880f9dbfaf44cad10b16edc7f4576e9629149b37 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/Pantry.idr @@ -0,0 +1,54 @@ +module S3D.Scenes.House.Pantry + +import Data.Vect +import Math.LinearAlgebra +import S3D.DrawablePolyhedron +import S3D.DrawablePolyhedron.SimpleColour +import S3D.GLProgram.SimpleColour +import S3D.PhysicalPolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Drawable +import S3D.Polyhedra.Precomputed +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron + +%default total + +box : PhysicalPolyhedron ? SimpleColour.Attributes True +box = + transform + ( diagonalMatrix [0.125, 0.125, 0.125, 1] !*! + translate [0, 0.125, 0] + ) $ + MkPhysicalPolyhedron + { drawablePolyhedron = + addTriangulation Precomputed.cube $ + monochromeCube [0.7, 0.7, 0.85, 1] + , obstaclePolyhedron = tower + } + +export +pantry : List ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +pantry = + [ (_ ** transform (translate [0.7, 0, 0.72]) box) + , ( _ ** + transform + [ [cos (pi/20), 0, sin (pi/20), 0] + , [0, 1, 0, 0] + , [-sin (pi/20), 0, cos (pi/20), 0] + , [0.7, 0.25, 0.7, 1] + ] + box + ) + , ( _ ** + transform + [ [cos (pi/6), 0, sin (pi/6), 0] + , [0, 1, 0, 0] + , [-sin (pi/6), 0, cos (pi/6), 0] + , [-0.6, 0, 0.7, 1] + ] + box + ) + ] blob - /dev/null blob + 249ab345d7c453d20884f58781704c2459c21c3b (mode 644) --- /dev/null +++ src/S3D/Scenes/House/ToPrecompute/Kitchen.idr @@ -0,0 +1,62 @@ +module S3D.Scenes.House.ToPrecompute.Kitchen + +import Data.Vect +import Math.LinearAlgebra +import S3D.CSG.PhysicalPolyhedron +import S3D.GLProgram.SimpleColour.Attributes +import S3D.PhysicalPolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Drawable +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron +import SolidGeometry.Polyhedron + +%default total + +||| Kitchen sink, to be subtracted from the counter +sink : (xOffset : Double) -> PhysicalPolyhedron ? SimpleColour.Attributes False +sink xOffset = + transform + [ [0.125, 0, 0, 0] + , [0, 0.15, 0, 0] + , [0, 0, 0.1, 0] + , [xOffset, 0.3, 0.18, 1] + ] + $ + MkPhysicalPolyhedron + { drawablePolyhedron = monochromeCube [0.9, 0.9, 0.9, 1] + , obstaclePolyhedron = emptyPolyhedron + } + +||| The kitchen counter, placed in the middle of the kitchen. Transformed in +||| S3D.Scenes.House.Kitchen. +export +counter : PhysicalPolyhedron ? SimpleColour.Attributes False +counter = + -- Send [x, y, z, 1] to [x, y, -1+z, 1+z], so the top of the counter is level + -- with the windowsill. Just using LinearAlgebra.translate would make it + -- slanted compared the windowsill. + transform + [ [1, 0, 0, 0] + , [0, 1, 0, 0] + , [0, 0, 1, 1] + , [0, 0, -1, 1] + ] + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sink (-0.2)) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sink (-0.47)) + $ + transform + [ [0.5, 0, 0, 0] + , [0, 0.3, 0, 0] + , [0, 0, 0.15, 0] + , [-0.5, 0, 0.15, 1] + ] + $ + MkPhysicalPolyhedron + { drawablePolyhedron = monochromeCube [0.5, 0.9, 0.5, 1] + , obstaclePolyhedron = tower + } blob - /dev/null blob + f3cb971cee11e17bee7ed918dd327a3af33d39e2 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/ToPrecompute/LivingRoom.idr @@ -0,0 +1,47 @@ +module S3D.Scenes.House.ToPrecompute.LivingRoom + +import Data.Vect +import Math.LinearAlgebra +import S3D.CSG.DrawablePolyhedron +import S3D.GLProgram.SimpleColour.Attributes +import S3D.DrawablePolyhedron +import S3D.PhysicalPolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Drawable +import S3D.Transformable +import S3D.Transformable.DrawablePolyhedron +import S3D.Transformable.Polyhedron +import SolidGeometry.Polyhedron + +%default total + +||| A couch sitting at (0, 0, 0, 1) facing toward (0, 0, 1, 0). +export +couch : PhysicalPolyhedron ? SimpleColour.Attributes False +couch = + MkPhysicalPolyhedron + { drawablePolyhedron = + DrawablePolyhedron.intersect + ( DrawablePolyhedron.complement $ + transform + [ [0.4, 0, 0, 0] + , [0, 0.15, 0, 0] + , [0, 0, 0.2, 0] + , [0, 0.3, 0.1, 1] + ] + $ + monochromeCube [0.8, 0.8, 0.8, 1] + ) + $ + ( transform + [ [0.5, 0, 0, 0] + , [0, 0.15, 0, 0] + , [0, 0, 0.2, 0] + , [0, 0.15, 0, 1] + ] + $ + monochromeCube [0.4, 0.4, 0.9, 1] + ) + , obstaclePolyhedron = + transform (diagonalMatrix [0.5, 1, 0.2, 1]) tower + } blob - /dev/null blob + c53a528067757594398a9656db96366c1391d1ac (mode 644) --- /dev/null +++ src/S3D/Scenes/House/ToPrecompute/Walls.idr @@ -0,0 +1,134 @@ +||| Walls for the house. Lighting is applied in the S3D.Scenes.House.Walls +||| module. + +module S3D.Scenes.House.ToPrecompute.Walls + +import S3D.CSG.PhysicalPolyhedron +import S3D.DrawInput +import S3D.DrawablePolyhedron +import S3D.DrawablePolyhedron.SimpleColour +import S3D.GLProgram.SimpleColour.Attributes +import Math.LinearAlgebra +import S3D.PhysicalPolyhedron +import S3D.Polyhedra +import S3D.Polyhedra.Drawable +import S3D.Transformable +import S3D.Transformable.DrawablePolyhedron +import S3D.Transformable.PhysicalPolyhedron +import S3D.Transformable.Polyhedron +import SolidGeometry.Construct +import SolidGeometry.Polyhedron + +%default total + +doorway : Vector 4 Double -> PhysicalPolyhedron ? SimpleColour.Attributes False +doorway colour = + MkPhysicalPolyhedron + { drawablePolyhedron = + -- Send [+-1, y, +-1, 1] to [+-0.15, {-0.6, 0.6}, 1-y, 1+y]. + transform + [[0.15, 0, 0, 0], [0, 0, -1, 1], [0, 0.6, 0, 0], [0, 0, 1, 1]] + $ + monochromeTower colour + , obstaclePolyhedron = + -- Two planes constraining x coordinate to be between -1/8 and 1/8 at the + -- middle of the wall [0,0,1,1]. + intersect + ( transform + [[0, -1, 0, 0], [1, 0, 0, 0], [-1/16, 0, 1, 0], [-1/16, 0, 0, 1]] + positiveY + ) + ( transform + [[0, 1, 0, 0], [-1, 0, 0, 0], [1/16, 0, 1, 0], [1/16, 0, 0, 1]] + positiveY + ) + } + +||| A window shape to cut out. Ignored for collision detection. +window : (xOffset, yOffset, width, height : Double) -> + PhysicalPolyhedron ? SimpleColour.Attributes False +window xOffset yOffset width height = + MkPhysicalPolyhedron + { drawablePolyhedron = + -- Send [+-1, y, +-1, 1] to [xOffset+-width, yOffset+-height, 1-y, 1+y]. + transform + [[width, 0, 0, 0], [0, 0, -1, 1], [0, height, 0, 0], [xOffset, yOffset, 1, 1]] $ + monochromeTower [1, 1, 1, 1] + , obstaclePolyhedron = emptyPolyhedron + } + +||| A tall window shape to cut out, paramaterized by x coordinate. Ignored for +||| collision detection. +sideWindow : (xOffset : Double) -> PhysicalPolyhedron ? SimpleColour.Attributes False +sideWindow xOffset = window xOffset (1/2) (1/8) (1/4) + +||| A small window shape to cut out, paramaterized by x and y coordinates. +||| Ignored for collision detection. +smallWindow : (xOffset, yOffset : Double) -> PhysicalPolyhedron ? SimpleColour.Attributes False +smallWindow xOffset yOffset = window xOffset yOffset 0.07 0.07 + +||| A wall that counts for drawing and collision detection. +blankWall : PhysicalPolyhedron ? SimpleColour.Attributes False +blankWall = + -- Send [+-1, +-1, +-1, 1] to [+-1, {0,1}, 1, 1] plus or minus [0, 0, 1/32, -1/32]. + transform [[1, 0, 0, 0], [0, 0.5, 0, 0], [0, 0, 1/32, -1/32], [0, 0.5, 1, 1]] $ + MkPhysicalPolyhedron + { drawablePolyhedron = monochromeCube [0.8, 0.8, 0.8, 1] + , obstaclePolyhedron = tower + } + +export +wallWithWindows : PhysicalPolyhedron ? SimpleColour.Attributes False +wallWithWindows = + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ window 0.5 0.6 0.25 0.25) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ window (-0.5) 0.6 0.25 0.25) + $ + blankWall + +export +wallWithOneWindow : PhysicalPolyhedron ? SimpleColour.Attributes False +wallWithOneWindow = + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ window 0.4 0.6 0.25 0.25) + $ + blankWall + +export +wallWithDoor : PhysicalPolyhedron ? SimpleColour.Attributes False +wallWithDoor = + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ doorway [0.8, 0.8, 1, 1]) + $ + blankWall + +export +frontWall : PhysicalPolyhedron ? SimpleColour.Attributes False +frontWall = + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sideWindow (-0.4)) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sideWindow (-0.7)) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sideWindow 0.4) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ sideWindow 0.7) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ smallWindow 0.08 0.72) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ smallWindow (-0.08) 0.72) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ smallWindow 0.08 0.88) + $ + PhysicalPolyhedron.intersect + (PhysicalPolyhedron.complement $ smallWindow (-0.08) 0.88) + $ + wallWithDoor blob - /dev/null blob + 0ace1f546be3879e7ddf09b50201d5cd6912c2c5 (mode 644) --- /dev/null +++ src/S3D/Scenes/House/Walls.idr @@ -0,0 +1,74 @@ +module S3D.Scenes.House.Walls + +import Data.Fin +import Data.Vect +import Math.LinearAlgebra +import Math.LinearAlgebra.RotateAxes +import S3D.DrawablePolyhedron +import S3D.Figure +import S3D.GLProgram.SimpleColour +import S3D.Obstacles +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron +import S3D.PhysicalPolyhedron +import S3D.Polyhedra.Precomputed +import S3D.Scenes.House.ToPrecompute.Walls +import SolidGeometry.Polyhedron + +%default total + +-- Precomputed versions of the walls defined in S3D.Scenes.House.ToPrecompute.Walls. + +frontWall' : PhysicalPolyhedron ? SimpleColour.Attributes True +frontWall' = + addTriangulation Precomputed.houseFrontWall frontWall + +wallWithDoor' : PhysicalPolyhedron ? SimpleColour.Attributes True +wallWithDoor' = + addTriangulation Precomputed.wallWithDoor wallWithDoor + +wallWithOneWindow' : PhysicalPolyhedron ? SimpleColour.Attributes True +wallWithOneWindow' = + addTriangulation Precomputed.wallWithOneWindow wallWithOneWindow + +wallWithWindows' : PhysicalPolyhedron ? SimpleColour.Attributes True +wallWithWindows' = + addTriangulation Precomputed.wallWithWindows wallWithWindows + +-- For each wall we haven't customized, a transformation to apply to the +-- generic wallWithDoor shape. +genericWallTransformations : List (Matrix 4 4 Double) +genericWallTransformations = + [ identityMatrix 4 + , rotateAxes 0 2 + , rotateAxes 0 3 + , rotateAxes 2 0 + , rotateAxes 3 0 + , rotateAxes 3 2 + ] + ++ map (diagonalMatrix [1, 1, -1, -1] !*!) + [ rotateAxes 0 2 + , rotateAxes 2 0 + ] + +export +walls : List + ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +walls = + (_ ** transform (rotateAxes 2 3) frontWall') :: + ( _ ** + transform (diagonalMatrix [1, 1, -1, -1] !*! rotateAxes 3 0) + wallWithOneWindow' + ) :: + [ (_ ** transform t wallWithWindows') + | t <- + map (diagonalMatrix [1, 1, -1, -1] !*!) + [ identityMatrix 4 + , rotateAxes 0 3 + ] + ] ++ + [ (_ ** transform t wallWithDoor') + | t <- genericWallTransformations + ] blob - /dev/null blob + 6f5c761e54fe4c994aa7486315317d0e2f86e42a (mode 644) --- /dev/null +++ src/S3D/Scenes/House/Yard.idr @@ -0,0 +1,101 @@ +module S3D.Scenes.House.Yard + +import Control.Order.Fold +import Data.List +import Data.Vect +import S3D.DrawInput +import S3D.Figure +import S3D.GLProgram.SimpleColour +import Math.LinearAlgebra +import S3D.Transformable + +%default total + +-- System.Random.prim__randomDouble doesn't seem to be implemented (yet?) so +-- roll our own for now. +-- +-- A random Double in (0, 1). +%foreign "javascript:lambda:() => Math.random()" +prim__randomDouble : PrimIO Double +randomDouble : HasIO io => io Double +randomDouble = primIO prim__randomDouble + +moveNearZPole : Vector 3 Double -> Vector 3 Double +moveNearZPole v = + let [x, y, z] = vectorToVect v in + case argmax (map abs (vectorToVect v)) of + 0 => [z, y, abs x] + 1 => [x, z, abs y] + 2 => [x, y, abs z] + +-- Move an object from [0, 0, 0, 1] to the given location. +moveTo : Transformable a => Vector 4 Double -> a -> a +moveTo location = + let newX = normalize (cross [location, [0, 1, 0, 0], [0, 0, 1, 0]]) + newZ = normalize (cross [location, newX, [0, 1, 0, 0]]) + in transform [newX, [0, 1, 0, 0], newZ, location] + +-- Sample from the standard normal distribution using the Box-Muller transform. +sampleStandardNormal : HasIO io => io Double +sampleStandardNormal = + do x <- randomDouble + y <- randomDouble + pure (sqrt (-2 * log x) * cos (2 * pi * y)) + +randomPointOnSphere : HasIO io => {n : Nat} -> io (Vector (S n) Double) +randomPointOnSphere = + do normalDistributionSample <- sequence $ replicate (S n) sampleStandardNormal + pure $ normalize $ vectToVector normalDistributionSample + +-- A random location for a blade of grass. +randomGrassBladeLocation : HasIO io => io (Vector 4 Double) +randomGrassBladeLocation = + do -- Pick a random place on the 2-sphere (fix y=0). Then, apply the + -- moveNearZPole function so that it ends up within the yard. + [x, z, w] <- map (vectorToVect . moveNearZPole) $ randomPointOnSphere {n=2} + pure $ [x, 0, z, w] + +grassColour : Vector 4 Double +grassColour = [0.05, 0.9, 0.1, 1] + +-- A blade of grass +grassBlade : DrawInput Attributes +grassBlade = + MkDrawInput + [T [0, 1, 2]] + [ MkAttributes + (normalize [-0.0015, 0, 0, 1]) + grassColour + , MkAttributes + (normalize [0.0015, 0, 0, 1]) + grassColour + , MkAttributes + (normalize [0, 0.02, 0, 1]) + grassColour + ] + +-- A blade of grass with a random location and orientation. +randomGrassBlade : HasIO io => io (DrawInput Attributes) +randomGrassBlade = + do location <- randomGrassBladeLocation + orientation <- map (2 * pi *) randomDouble + pure $ + moveTo location $ + transform + [ [cos orientation, 0, sin orientation, 0] + , [0, 1, 0, 0] + , [-sin orientation, 0, cos orientation, 0] + , [0, 0, 0, 1] + ] + grassBlade + +grass : HasIO io => io Figure +grass = + let numBlades = 400 + in + map (simpleColourFigure . concatDrawInputs) $ + sequence (replicate numBlades randomGrassBlade) + +export +yard : HasIO io => io Figure +yard = grass blob - /dev/null blob + 93c07be1fac9ace00cd19f6b349273e8591aa721 (mode 644) --- /dev/null +++ src/S3D/Scenes/House.idr @@ -0,0 +1,133 @@ +{- + +The y=0 plane is the floor. + +Each room is initially set up at the w pole, then transformed to its final +location. Each wall is initially set up so its bottom runs along the line from +[-1,0,1,1] to [1,0,1,1], and then is similarly transformed. + +Rooms: +* +x pole: dining room +* -x pole: kitchen +* +z pole: living room +* -z pole: yard +* +w pole: entrance +* -w pole: bedroom + +The typical ceiling for a room at the +w pole is at points of the form +[x, 1, z, 1] where x and z are between -1 and 1. That means the ceiling height +varies from pi/4 radians in the middle of the room down to pi/6 radians at the +corner of the room. + +We'll say 1 metre = pi / 16 radians, so: +- The typical ceiling height ranges from 2.7m to 4m. +- The walls are 6.3m long at the floor and 5.3m at the ceiling. + +-} + +module S3D.Scenes.House + +import Data.Fin +import Data.List +import Data.Vect +import S3D.DrawInput +import S3D.DrawablePolyhedron +import S3D.GLProgram.SimpleColour +import S3D.Figure +import Math.LinearAlgebra +import Math.LinearAlgebra.RotateAxes +import S3D.Lighting +import S3D.Lighting.SimpleColour +import S3D.Obstacles +import S3D.PhysicalObject +import S3D.PhysicalPolyhedron +import S3D.Scenes.House.DiningRoom +import S3D.Scenes.House.FloorLines +import S3D.Scenes.House.Kitchen +import S3D.Scenes.House.LivingRoom +import S3D.Scenes.House.Pantry +import S3D.Scenes.House.Walls +import S3D.Scenes.House.Yard +import S3D.Transformable +import S3D.Transformable.PhysicalPolyhedron +import SolidGeometry.Polyhedron + +%default total + +lights : List LightSource +lights = + [ MkLightSource { location = [0, 1.5, -1, 0] + , colour = [1.2, 1.2, 0.8] + } + , MkLightSource { location = [0, 1.5, 1, 0] + , colour = [0.8, 0.8, 1.2] + } + , MkLightSource { location = [0, 1.5, 0, -1] + , colour = [0.5, 0.3, 0.3] + } + ] + +entrance : List (numPlanes : Nat ** PhysicalPolyhedron numPlanes SimpleColour.Attributes True) +entrance = [] + +floorColours : Vect 32 (Vector 4 Double) +floorColours = + fromList $ concat $ map concat $ the (List (List (List (Vector 4 Double)))) + [ -- dining room + replicate 2 [[1, 0.9, 0.8, 1], [0.8, 0.7, 0.6, 1]] + , -- living room + replicate 4 [[0.3, 0.5, 1, 1], [0.2, 0.4, 0.8, 1]] + , -- kitchen + replicate 4 [[0.9, 0.9, 0.9, 1], [0.7, 0.7, 0.7, 1]] + , -- yard + replicate 4 [[0.3, 0.2, 0.1, 1], [0.1, 0.3, 0.1, 1]] + , -- dining room, again + replicate 2 [[1, 0.9, 0.8, 1], [0.8, 0.7, 0.6, 1]] + ] + +physicalPolyhedra : List + ( numPlanes : Nat ** + PhysicalPolyhedron numPlanes SimpleColour.Attributes True + ) +physicalPolyhedra = + walls ++ + [ (numPlanes ** transform t polyhedron) + | (room, t) <- + the (List (_, Matrix 4 4 Double)) + [ (entrance, identityMatrix 4) + , (livingRoom, rotateAxes 3 0) + , (kitchen, rotateAxes 0 3) + , (pantry, rotateAxes 3 2) + , (diningRoom, diagonalMatrix [1, 1, -1, -1]) + ] + , (numPlanes ** polyhedron) <- room + ] + +physicalPolyhedraFigure : Figure +physicalPolyhedraFigure = + simpleColourFigure $ concatDrawInputs + [ drawPrecomputedPolyhedron $ + illuminate lights p.drawablePolyhedron + | (_ ** p) <- physicalPolyhedra + ] + +physicalPolyhedraObstacles : Obstacles +physicalPolyhedraObstacles = + makeObstacles [(_ ** p.obstaclePolyhedron) | (_ ** p) <- physicalPolyhedra] + +export +theHouse : HasIO io => io PhysicalObject +theHouse = + pure $ + MkPhysicalObject + { figure = + concatFigures + [ physicalPolyhedraFigure + , floorLines floorColours + , transform (rotateAxes 2 3) !yard + ] + , obstacles = physicalPolyhedraObstacles + } + + + blob - /dev/null blob + fa7d97daa8d6a2c315790b746f1dfee1d7860705 (mode 644) --- /dev/null +++ src/S3D/Transformable/DrawablePolyhedron.idr @@ -0,0 +1,22 @@ +module S3D.Transformable.DrawablePolyhedron + +import S3D.DrawablePolyhedron +import Math.LinearAlgebra +import S3D.Transformable +import S3D.Transformable.Polyhedron +import S3D.Transformable.TriangulatedFace +import SolidGeometry.Polyhedron + +%default total + +export +(pretriangulated : Bool) => + Transformable (DrawablePolyhedron numPlanes attributes pretriangulated) +where + transform matrix with (pretriangulated) + transform matrix | False = + {polyhedron $= transform matrix} + transform matrix | True = + { polyhedron $= transform matrix + , triangulation $= map (transform matrix) + } blob - /dev/null blob + 16b12de5e7231d5c7c63ba011895b4d33d08d9b5 (mode 644) --- /dev/null +++ src/S3D/Transformable/PhysicalObject.idr @@ -0,0 +1,13 @@ +module S3D.Transformable.PhysicalObject + +import S3D.Figure +import Math.LinearAlgebra +import S3D.Obstacles +import S3D.PhysicalObject +import S3D.Transformable + +%default total + +export +Transformable PhysicalObject where + transform matrix = {figure $= transform matrix, obstacles $= transform matrix} blob - /dev/null blob + b0efd1b8225f6704c9e75a693d51301033c79f42 (mode 644) --- /dev/null +++ src/S3D/Transformable/PhysicalPolyhedron.idr @@ -0,0 +1,19 @@ +module S3D.Transformable.PhysicalPolyhedron + +import S3D.DrawablePolyhedron +import Math.LinearAlgebra +import S3D.PhysicalPolyhedron +import S3D.Transformable +import S3D.Transformable.DrawablePolyhedron +import S3D.Transformable.Polyhedron +import SolidGeometry.Polyhedron + +%default total + +export +(pretriangulated : Bool) => + Transformable (PhysicalPolyhedron n attributes pretriangulated) where + transform matrix = + { drawablePolyhedron $= transform matrix + , obstaclePolyhedron $= transform matrix + } blob - /dev/null blob + 42e18cd9a887b03057fac2b02e27c596f99cad60 (mode 644) --- /dev/null +++ src/S3D/Transformable/Polyhedron.idr @@ -0,0 +1,14 @@ +-- Transformable instance for `Polyhedron`. + +module S3D.Transformable.Polyhedron + +import Data.Vect +import Math.LinearAlgebra +import S3D.Transformable +import SolidGeometry.Polyhedron + +%default total + +export +Transformable (Polyhedron numPlanes) where + transform matrix = { planes $= map (!*! matrix) } blob - /dev/null blob + 08f7ff362c80450565f51fc03105b71cfbc6ff66 (mode 644) --- /dev/null +++ src/S3D/Transformable/TriangulatedFace.idr @@ -0,0 +1,17 @@ +module S3D.Transformable.TriangulatedFace + +import Data.Vect +import Math.LinearAlgebra +import Math.LinearAlgebra.Normalized +import S3D.Transformable +import SolidGeometry.Polyhedron + +%default total + +export +Transformable VertexCoordinates where + transform matrix = {absolute $= makeNormalized . transform matrix . fromNormalized} + +export +Transformable (TriangulatedFace numPlanes) where + transform matrix = {vertices $= map (transform matrix)} blob - /dev/null blob + a96e056074da55c3d3f0d48e12e51d76cf0a391b (mode 644) --- /dev/null +++ src/S3D/Transformable.idr @@ -0,0 +1,26 @@ +{- + +Class for types that it makes sense to apply a 4D linear transformation to. The +matrix is applied on the right-hand side, since this works better with the GL +point of view. + +Generally, anything position-related is transformed, but other attributes like +colours are unaffected. + +-} + +module S3D.Transformable + +import Math.LinearAlgebra + +%default total + +public export +interface Transformable a where + ||| If positive/negative orientation of the object being transformed matters, + ||| make sure the determinant of the transformation matrix is positive. + transform : Matrix 4 4 Double -> a -> a + +export +Transformable (Vector 4 Double) where + transform matrix vector = vector ^*! matrix blob - /dev/null blob + 789216569320ec6817459648edf4441dda880d6e (mode 644) --- /dev/null +++ src/S3D/WebUI.idr @@ -0,0 +1,142 @@ +module S3D.WebUI + +import Data.IORef +import Data.Vect +import Data.SortedSet +import JS +import Math.LinearAlgebra +import S3D.DomExtras +import Web.Html +import Web.Dom +import Web.Raw.UIEvents + +%default total + +||| We ignore touches if there are more than this many touches active. +public export +maxTouches : Nat +maxTouches = 2 + +-- Mutable state managed by the UI code. The init function sets up event +-- handlers which update the state field of the UIContext record it returns. +public export +record UIState where + constructor MkUIState + keysDown : SortedSet String + ||| The centre of the current touches. + touchCentre : Vector 2 Double + ||| index (numTouches-1) numDeltas is how much touchCentre has changed since + ||| the last physics update while numTouches touches were active. + touchDeltas : Vect WebUI.maxTouches (Vector 2 Double) + +public export +record UIContext where + constructor MkUIContext + canvas : HTMLCanvasElement + gl : WebGL2RenderingContext + state : IORef UIState + +makeCanvas : JSIO HTMLCanvasElement +makeCanvas = + do canvas <- createElement Canvas + ignore (appendChild (!body) canvas) + width canvas .= 300 + height canvas .= 300 + pure canvas + +getGl2Context : HTMLCanvasElement -> JSIO WebGL2RenderingContext +getGl2Context canvas = + do context <- unMaybe "getContext" $ Html.HTMLCanvasElement.getContext canvas "webgl2" Undef + case context of + (S (S (S (Z gl)))) => pure gl + _ => throwError $ Caught "Got wrong context type." + +handleKeyDown : IORef UIState -> Event -> JSIO () +handleKeyDown stateRef event = + do keyboardEvent <- unMaybe "cast keyboard event" $ pure (safeCast event) + theKey <- Web.Raw.UIEvents.KeyboardEvent.key keyboardEvent + modifyIORef stateRef $ { keysDown $= insert theKey } + +handleKeyUp : IORef UIState -> Event -> JSIO () +handleKeyUp stateRef event = + do keyboardEvent <- unMaybe "cast keyboard event" $ pure (safeCast event) + theKey <- Web.Raw.UIEvents.KeyboardEvent.key keyboardEvent + modifyIORef stateRef $ { keysDown $= delete theKey } + +handleFocusout : IORef UIState -> Event -> JSIO () +handleFocusout stateRef event = + modifyIORef stateRef $ { keysDown := empty } + +natsBelow' : List Nat -> Nat -> List Nat +natsBelow' acc Z = acc +natsBelow' acc (S n) = natsBelow' (n :: acc) n + +natsBelow : Nat -> List Nat +natsBelow = natsBelow' [] + +touchListCentre : TouchList -> JSIO (Vector 2 Double) +touchListCentre touchList = + do numTouches <- getNumTouches touchList + locations <- sequence + [ do touch <- getTouchItem touchList i + x <- getClientX touch + y <- getClientY touch + -- TODO: get canvas size instead of hard-coding + pure [cast x / 150, - cast y / 150] + | i <- natsBelow numTouches + ] + pure $ sumVectorsFrom [0, 0] locations ^/ cast numTouches + +handleTouchStartEnd : IORef UIState -> Event -> JSIO () +handleTouchStartEnd stateRef event = + do preventDefault event + touchList <- getTouchList event + centre <- touchListCentre touchList + modifyIORef stateRef {touchCentre := centre} + +handleTouchMove : IORef UIState -> Event -> JSIO () +handleTouchMove stateRef event = + do preventDefault event + touchList <- getTouchList event + newCentre <- touchListCentre touchList + numTouches@(S numTouchesMinusOne) <- getNumTouches touchList + | Z => pure () + case isLTE numTouches maxTouches of + No _ => pure () + Yes lte => + let numTouchesMinusOne' : Fin WebUI.maxTouches + numTouchesMinusOne' = natToFinLT numTouchesMinusOne + in modifyIORef stateRef $ \state => + { touchCentre := newCentre + , touchDeltas $= updateAt numTouchesMinusOne' + (^+^ (newCentre ^-^ state.touchCentre)) + } + state + +export +init : JSIO UIContext +init = + do canvas <- makeCanvas + gl <- getGl2Context canvas + state <- newIORef $ + MkUIState + { keysDown = empty + , touchCentre = [0, 0] + , touchDeltas = replicate _ [0, 0] + } + theDocument <- document + -- Keyboard events + addEventListener' theDocument "keydown" (Just !(callback (handleKeyDown state))) + addEventListener' theDocument "keyup" (Just !(callback (handleKeyUp state))) + addEventListener' theDocument "focusout" (Just !(callback (handleFocusout state))) + -- Touch events + addEventListener' canvas "touchstart" (Just !(callback (handleTouchStartEnd state))) + addEventListener' canvas "touchend" (Just !(callback (handleTouchStartEnd state))) + addEventListener' canvas "touchcancel" (Just !(callback (handleTouchStartEnd state))) + addEventListener' canvas "touchmove" (Just !(callback (handleTouchMove state))) + pure $ + MkUIContext + { canvas = canvas + , gl = gl + , state = state + } blob - /dev/null blob + 455b9621236af4e219aa0253ccfd98a61cfa95dc (mode 644) --- /dev/null +++ src/SolidGeometry/Circle.idr @@ -0,0 +1,213 @@ +module SolidGeometry.Circle + +import Data.Fin +import Data.Linear.SizedArray +import Data.List +import Data.List1 +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.Sign + +%default total + +public export +data CircleEdges : Nat -> Type where + AllOneWay : Sign -> CircleEdges n + Edges : List ((Fin n, Sign), (Fin n, Sign)) -> CircleEdges n + +export +Eq (CircleEdges n) where + AllOneWay s0 == AllOneWay s1 = s0 == s1 + AllOneWay _ == Edges _ = False + Edges _ == AllOneWay _ = False + Edges e0 == Edges e1 = e0 == e1 + +export +Show (CircleEdges n) where + show (AllOneWay s) = "AllOneWay " <+> show s + show (Edges e) = "Edges (" <+> show e <+> ")" + +-- circleToLine and compareLocations define a total order on positions on +-- the circle, starting at [1, 0], then visiting [0, 1], [-1, 0] and [0, +-- -1] in that order. +circleToLine : Vector 2 Double -> Double +circleToLine v = + let [x, y] = vectorToVect v in + if y == 0 && x > 0 + then 0 + else case (softSign x, softSign y) of + (Positive, Positive) => y / (x + y) + (Negative, Positive) => 2 - y / (-x + y) + (Negative, Negative) => 2 - y / (-x - y) + (Positive, Negative) => 4 + y / (x - y) + +compareLocations : (Vector 2 Double, a) -> (Vector 2 Double, a) -> Ordering +compareLocations (v0, _) (v1, _) = compare (circleToLine v0) (circleToLine v1) + +||| Helper for mergeIntervals, for the case where either `snd (head intervals) +||| ==` is Nothing or we've seen the end of the first interval. +||| +||| TODO: Prove it's total. +covering +mergeAfterFirstInterval : + Eq b => + a -> Maybe (a, b) -> List ((a, a), b) -> List1 (a, Maybe b) -> + List ((a, a), b) +mergeAfterFirstInterval x Nothing acc ((hx, Nothing) ::: []) = + reverse acc +mergeAfterFirstInterval x (Just (x', v)) acc ((hx, Nothing) ::: []) = + ((x, x'), v) :: reverse acc +mergeAfterFirstInterval x v acc ((hx, Nothing) ::: h :: t) = + mergeAfterFirstInterval x v acc (h ::: t) +mergeAfterFirstInterval x v acc ((hx, hv@(Just hvv)) ::: t) = + case span (\(_, v') => v' == hv) t of + (matches, []) => + case v of + Nothing => ((hx, x), hvv) :: reverse acc + Just (x', v') => + if v' == hvv + then ((hx, x'), hvv) :: reverse acc + else ((hx, x), hvv) :: ((x, x'), v') :: reverse acc + (matches, h' :: t') => + mergeAfterFirstInterval x v (((hx, fst h'), hvv) :: acc) (h' ::: t') + +||| Helper for mergeIntervals, for the case where `snd (head intervals)` is a +||| Just, and all the `snd`s we've seen so far are that same value. +covering +mergeStillInFirstInterval : Eq b => + a -> b -> List (a, Maybe b) -> + Either b (List ((a, a), b)) +mergeStillInFirstInterval _ v [] = Left v +mergeStillInFirstInterval x v (h@(hx, hv) :: t) = + if hv == Just v + then mergeStillInFirstInterval x v t + else Right $ mergeAfterFirstInterval x (Just (hx, v)) [] (h ::: t) + +covering +mergeIntervals : Eq b => + (intervals : List1 (a, Maybe b)) -> Either b (List ((a, a), b)) +mergeIntervals ((x, Nothing) ::: []) = Right [] +mergeIntervals ((x, Nothing) ::: h :: t) = + Right $ mergeAfterFirstInterval x Nothing [] (h ::: t) +mergeIntervals ((x, Just v) ::: t) = + mergeStillInFirstInterval x v t + +||| Turn an implicit description of edges along a circle into an explicit one. +||| +||| Each element b of `boundaries` describes a half-space: the Positive side of +||| `b` consists of points p where `determinant [b, p]` is positive. +||| +||| `edgeDirection` says which direction edges point within each combination of +||| half-space. A Positive edge is any edge from u to v where +||| `determinant [u, v]` is positive. +||| +||| For example, if +||| +||| boundaries = [[1, 0], [0, 1]] +||| +||| then the circle is divided into quadrants like this: +||| +||| direction of increasing second coordinate +||| ^ +||| | +||| -+- +||| (Negative, Positive) / | \ (Positive, Positive) +||| / | \ +||| | | | +||| -+---+---+---> direction of increasing first coordinate +||| | | | +||| \ | / +||| (Negative, Negative) \ | / (Positive, Negative) +||| -+- +||| | +||| +||| where (sign0, sign1) indicates points in that quadrant have sign0 according +||| to the first boundary, and sign1 according to the second boundary. +||| +||| Then if +||| +||| edgeDirection f = if f 0 == Positive && f 1 == Positive then Just Positive else Nothing +||| +||| we end up with a single edge, from [1, 0] to [0, 1]: +||| +||| | +||| |< +||| | \ +||| | \ +||| | | +||| -----+----- +||| | +||| | +||| | +||| | +||| +||| The return value `AllOneWay s` means the entire circle points in the same +||| direction. E.g. this can happen if `edgeDirection` is `const (Just s)`. +||| +||| Currently can't handle edges that go halfway around the circle or more. +||| (TODO: fix that.) +export covering +boundariesToEdges : {n : Nat} -> + (edgeDirection : (Fin n -> Sign) -> Maybe Sign) -> + (boundaries : Vect n (Vector 2 Double)) -> + CircleEdges n +boundariesToEdges edgeDirection boundaries = + let verticesInOrder : List (Fin n, Sign) + verticesInOrder = + map snd $ + sortBy compareLocations + [ (signAsDouble *^ loc, (i, sign)) + | (i, loc) <- toList $ zip range boundaries + , (sign, signAsDouble) <- [(Positive, 1.0), (Negative, -1.0)] + ] + -- Which side of the given half-space [1, 0] is on. As a special case, + -- it's on the Negative side of [1, 0] and the Positive side of [-1, 0]. + locationToSide : Vector 2 Double -> Sign + locationToSide v = + let [x, y] = vectorToVect v in + if y == 0 + then softSign (-x) + else softSign (-y) + 0 SideAssignment : Type + SideAssignment = SizedLinArray n Sign + startSideLin : LinearResult SideAssignment + startSideLin = fromVect $ map locationToSide boundaries + flipSignAt : Fin n -> (1 _ : SideAssignment) -> SideAssignment + flipSignAt i assignment = + case mread assignment i of + sign # assignment' => + write assignment' i $ flipSign sign + -- Compute a list of (vertex index, direction in the following segment). + computeIntervalsFrom : List ((Fin n, Sign), Maybe Sign) -> (1 assignment : SideAssignment) -> + List (Fin n, Sign) -> Res (List ((Fin n, Sign), Maybe Sign)) (const SideAssignment) + computeIntervalsFrom acc assignment [] = reverse acc # assignment + computeIntervalsFrom acc assignment (h@(hi, _) :: t) = + case unsafeFreeze (flipSignAt hi assignment) (\frozenAssignment => edgeDirection (read frozenAssignment)) of + theEdgeDirection # assignment' => + computeIntervalsFrom + ((h, theEdgeDirection) :: acc) + assignment' + t + in + startSideLin $ + \startSide => + case unsafeFreeze startSide (\frozen => edgeDirection (read frozen)) of + startEdgeDirection # startSide' => + case computeIntervalsFrom [] startSide' verticesInOrder of + intervals # leftoverArray => + let result = case List1.fromList intervals of + Nothing => + case startEdgeDirection of + Nothing => Edges neutral + Just sign => AllOneWay sign + Just intervals' => + case mergeIntervals intervals' of + Left sign => AllOneWay sign + Right merged => + Edges $ + [ case sign of + Positive => (u, v) + Negative => (v, u) + | ((u, v), sign) <- merged + ] + in result # leftoverArray blob - /dev/null blob + 823a9b11b8e4a19f45d5239c17ec902b21a8b342 (mode 644) --- /dev/null +++ src/SolidGeometry/CircleTest.idr @@ -0,0 +1,192 @@ +module SolidGeometry.CircleTest + +import Data.Fin +import Data.List +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.Circle +import SolidGeometry.Sign +import Test.TestTree + +%default total + +record BoundariesToEdgesTestCase where + constructor MkBoundariesToEdgesTestCase + name : String + {numHalfSpaces : Nat} + edgeDirection : (Fin numHalfSpaces -> Sign) -> Maybe Sign + boundaries : Vect numHalfSpaces (Vector 2 Double) + expectedResult : CircleEdges numHalfSpaces + +covering +boundariesToEdgesTestCase : BoundariesToEdgesTestCase -> TestCase +boundariesToEdgesTestCase testCase = + let actualResult : CircleEdges testCase.numHalfSpaces + actualResult = + boundariesToEdges testCase.edgeDirection testCase.boundaries + actualEqualsExpected : Bool + actualEqualsExpected = + case (actualResult, testCase.expectedResult) of + (AllOneWay s0, AllOneWay s1) => s0 == s1 + (Edges e0, Edges e1) => sort e0 == e1 + _ => False + in + MkTestCase testCase.name $ pure $ + if actualEqualsExpected + then MkSuccess + else MkFailure $ + "Expected: " <+> show testCase.expectedResult <+> + "; actual: " <+> show actualResult + +empty : BoundariesToEdgesTestCase +empty = + MkBoundariesToEdgesTestCase + { name = "empty" + , edgeDirection = const Nothing + , boundaries = [] + , expectedResult = Edges neutral + } + +emptyOneWay : BoundariesToEdgesTestCase +emptyOneWay = + MkBoundariesToEdgesTestCase + { name = "empty, one way" + , edgeDirection = const (Just Positive) + , boundaries = [] + , expectedResult = AllOneWay Positive + } + +nonEmptyOneWay : BoundariesToEdgesTestCase +nonEmptyOneWay = + MkBoundariesToEdgesTestCase + { name = "non-empty, one way" + , edgeDirection = const (Just Negative) + , boundaries = [[1, 0], [0, -1], [1, -1]] + , expectedResult = AllOneWay Negative + } + +quarterEdge : Double -> Sign -> BoundariesToEdgesTestCase +quarterEdge startingY sign = + let edgeDirection : (Fin 2 -> Sign) -> Maybe Sign + edgeDirection f = if f 0 == Positive && f 1 == Negative + then Just sign + else Nothing + expectedEdges : List ((Fin 2, Sign), (Fin 2, Sign)) + expectedEdges = + case sign of + Positive => [((0, Positive), (1, Positive))] + Negative => [((1, Positive), (0, Positive))] + in + MkBoundariesToEdgesTestCase + { name = "one quarter-circle edge; start at [1, " <+> show startingY <+> + "]; sign: " <+> show sign + , edgeDirection = edgeDirection + , boundaries = [[1, startingY], [0, 1]] + , expectedResult = Edges expectedEdges + } + +edgeTo10 : BoundariesToEdgesTestCase +edgeTo10 = + let edgeDirection : (Fin 2 -> Sign) -> Maybe Sign + edgeDirection f = if f 0 == Negative && f 1 == Negative + then Just Positive + else Nothing + expectedEdges : List ((Fin 2, Sign), (Fin 2, Sign)) + expectedEdges = [((1, Negative), (0, Positive))] + in + MkBoundariesToEdgesTestCase + { name = "edge ending at [1, 0]" + , edgeDirection = edgeDirection + , boundaries = [[1, 0], [0, 1]] + , expectedResult = Edges expectedEdges + } + +edgeThrough10 : BoundariesToEdgesTestCase +edgeThrough10 = + let edgeDirection : (Fin 2 -> Sign) -> Maybe Sign + edgeDirection f = if f 0 == Negative && f 1 == Positive + then Just Positive + else Nothing + expectedEdges : List ((Fin 2, Sign), (Fin 2, Sign)) + expectedEdges = [((1, Positive), (0, Positive))] + in + MkBoundariesToEdgesTestCase + { name = "edge through [1, 0]" + , edgeDirection = edgeDirection + , boundaries = [[1, 1], [1, -1]] + , expectedResult = Edges expectedEdges + } + +||| [1, 2] +||| [-1, 1] | [1, 1] +||| \<-|<-/ +||| \ | /^ +||| \|/ | +||| ------*---- [2, 1] +||| ^ /|\ +||| |/ | \ +||| /->|<-\ +manyEdges : BoundariesToEdgesTestCase +manyEdges = + let edgeDirection : (Fin 4 -> Sign) -> Maybe Sign + edgeDirection f = + case f 0 of + Positive => + case f 1 of + Positive => Nothing + Negative => Just Positive + Negative => + case f 1 of + Negative => Nothing + Positive => + case f 3 of + Negative => Just Negative + Positive => + case f 2 of + Positive => Just Negative + Negative => Just Positive + expectedEdges : List ((Fin 4, Sign), (Fin 4, Sign)) + expectedEdges = + [ ((0, Positive), (1, Positive)) + , ((1, Negative), (3, Negative)) + , ((2, Negative), (0, Negative)) + , ((2, Negative), (3, Negative)) + ] + in + MkBoundariesToEdgesTestCase + { name = "many edges" + , edgeDirection =edgeDirection + , boundaries = [[2, 1], [-1, 1], [1, 1], [1, 2]] + , expectedResult = Edges expectedEdges + } + +||| Two copies of [1, 0] and two copies of [0, 1]. edgeDirection is only +||| non-Nothing in an imposible case. +degenerateCubeEdges : BoundariesToEdgesTestCase +degenerateCubeEdges = + let edgeDirection : (Fin 4 -> Sign) -> Maybe Sign + edgeDirection f = + if map f range == the (Vect 4 Sign) [Positive, Negative, Positive, Negative] + then Just Positive + else Nothing + in + MkBoundariesToEdgesTestCase + { name = "degenerate case from a cube" + , edgeDirection = edgeDirection + , boundaries = [[1, 0], [1, 0], [0, 1], [0, 1]] + , expectedResult = Edges neutral + } + +export covering +tests : TestTree +tests = + MkTestGroup "Circle.boundariesToEdges" $ map (MkTestLeaf . boundariesToEdgesTestCase) $ + [ empty, emptyOneWay, nonEmptyOneWay + , edgeTo10, edgeThrough10, manyEdges + , degenerateCubeEdges + ] + ++ + [ quarterEdge startAtY sign + | startAtY <- [-1, 0, 1] + , sign <- [Positive, Negative] + ] blob - /dev/null blob + 183544b467093f741b546d59ed691e15ec110d04 (mode 644) --- /dev/null +++ src/SolidGeometry/Construct.idr @@ -0,0 +1,82 @@ +||| Tools for constructing polyhedra out of other polyhedra. + +module SolidGeometry.Construct + +import Data.Fin +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.HalfSpace +import SolidGeometry.Polyhedron +import SolidGeometry.Sign + +%default total + +export +complement : Polyhedron numPlanes -> Polyhedron numPlanes +complement = { inside $= (flipSide .) } + +0 finLtLimit : (n : Nat) -> (x : Fin n) -> finToNat x `LT` n +finLtLimit (S _) FZ = LTESucc LTEZero +finLtLimit (S n') (FS x') = LTESucc $ finLtLimit n' x' + +||| A faster version of Data.Fin.shift. +fastShift : (m : Nat) -> Fin n -> Fin (m + n) +fastShift m x = + let xAsNat : Nat + xAsNat = finToNat x + resultAsNat : Nat + resultAsNat = m + xAsNat + 0 mPlusSxLteMPlusN : m + S xAsNat `LTE` m + n + mPlusSxLteMPlusN = plusLteMonotoneLeft m (S xAsNat) n (finLtLimit n x) + 0 resultLtMPlusN : resultAsNat `LT` m + n + resultLtMPlusN = rewrite plusSuccRightSucc m xAsNat in mPlusSxLteMPlusN + in natToFinLT resultAsNat {prf = resultLtMPlusN} + +0 finToNatToFin : (x : Fin n) -> {prf : finToNat x `LT` n} -> natToFinLT (finToNat x) {prf} = x +finToNatToFin {prf = LTESucc _} FZ = Refl +finToNatToFin {n = S n'} {prf = LTESucc prf'} (FS x') = + cong FS $ finToNatToFin x' {prf = prf'} + +0 fastShiftByZero : (x : Fin n) -> fastShift 0 x = x +fastShiftByZero x = finToNatToFin x + +0 lteEqual : (p, p' : x `LTE` y) -> p = p' +lteEqual {x = Z} LTEZero LTEZero = Refl +lteEqual {x = S x'} (LTESucc p) (LTESucc p') = cong LTESucc $ lteEqual p p' + +0 natToFinLTSucc : {prf : S x `LT` S n} -> {prf' : x `LT` n} -> natToFinLT (S x) {prf} = FS (natToFinLT x {prf = prf'}) +natToFinLTSucc {prf = LTESucc _} {prf'} = rewrite lteEqual prf' _ in Refl + +0 fastShiftSucc : (m : Nat) -> (x : Fin n) -> fastShift (S m) x = FS (fastShift m x) +fastShiftSucc m x = natToFinLTSucc + +0 fastShiftIsShift : (m : Nat) -> (x : Fin n) -> fastShift m x = shift m x +fastShiftIsShift Z x = fastShiftByZero x +fastShiftIsShift (S m') x = + trans (fastShiftSucc m' x) $ cong FS $ fastShiftIsShift m' x + +export +intersect : {numPlanes0 : Nat} -> + Polyhedron numPlanes0 -> Polyhedron numPlanes1 -> + Polyhedron (numPlanes0 + numPlanes1) +intersect polyhedron0 polyhedron1 = + let inside : (Fin (numPlanes0 + numPlanes1) -> Sign) -> InOutSide + inside planeSides = + let planeSides0 : Fin numPlanes0 -> Sign + planeSides0 i = planeSides (weakenN numPlanes1 i) + planeSides1 : Fin numPlanes1 -> Sign + planeSides1 i = planeSides (fastShift numPlanes0 i) + in intersectSides + (polyhedron0.inside planeSides0) + (polyhedron1.inside planeSides1) + in + MkPolyhedron + (polyhedron0.planes ++ polyhedron1.planes) + inside + +export +union : {numPlanes0 : Nat} -> + Polyhedron numPlanes0 -> Polyhedron numPlanes1 -> + Polyhedron (numPlanes0 + numPlanes1) +union polyhedron0 polyhedron1 = + complement (intersect (complement polyhedron0) (complement polyhedron1)) blob - /dev/null blob + f6443f6a2956d2b26daf79e6851fae583a3b0a12 (mode 644) --- /dev/null +++ src/SolidGeometry/EdgeNetwork.idr @@ -0,0 +1,107 @@ +module SolidGeometry.EdgeNetwork + +import Data.DefaultMap +import Data.List +import Data.So +import Decidable.Equality -- DecEq Int is needed for Monoid IntMap instance. + +%default total + +record SortedEdge (vertex : Type) {auto vertexOrder : Ord vertex} where + constructor MkSortedEdge + first, second : vertex + {auto 0 sorted : So (first < second)} + +Ord vertex => Show vertex => Show (SortedEdge vertex) where + show e = "(" <+> show e.first <+> ", " <+> show e.second <+> ")" + +Ord vertex => Eq (SortedEdge vertex) where + MkSortedEdge u0 v0 == MkSortedEdge u1 v1 = u0 == u1 && v0 == v1 + +Ord vertex => Ord (SortedEdge vertex) where + compare (MkSortedEdge u0 v0) (MkSortedEdge u1 v1) = + case compare u0 u1 of + LT => LT + GT => GT + EQ => compare v0 v1 + +public export +IntMap : Type -> Type +IntMap key = DefaultMap key Int (the Int 0) + +export +record EdgeNetwork (vertex : Type) {auto vertexOrder : Ord vertex} where + constructor MkEdgeNetwork + -- The negative of an edge is the edge in the opposite direction. + edges : IntMap (SortedEdge vertex) + +export +Ord vertex => Eq (EdgeNetwork vertex) where + MkEdgeNetwork n0 == MkEdgeNetwork n1 = n0 == n1 + +export +Ord vertex => Show vertex => Show (EdgeNetwork vertex) where + show (MkEdgeNetwork edges) = + let %hint + defaultValue : Int + defaultValue = 0 + in show edges + +export +Ord vertex => Semigroup (EdgeNetwork vertex) where + MkEdgeNetwork x <+> MkEdgeNetwork y = + let %hint + monoidInt : Monoid Int + monoidInt = Additive + in MkEdgeNetwork (x <+> y) + +export +Ord vertex => Monoid (EdgeNetwork vertex) where + neutral = + let %hint + monoidInt : Monoid Int + monoidInt = Additive + in MkEdgeNetwork neutral + +edgeIfSorted : Ord vertex => (u, v : vertex) -> Int -> EdgeNetwork vertex +edgeIfSorted u v sign = + case choose (u < v) of + Left _ => MkEdgeNetwork (fromList [(MkSortedEdge u v, sign)]) + Right _ => neutral + +export +makeSignedEdge : Ord vertex => (u, v : vertex) -> EdgeNetwork vertex +makeSignedEdge u v = + edgeIfSorted u v 1 <+> edgeIfSorted v u (-1) + +||| Get all edges that are present at least once, in the forward direction. +export +getEdgesForward : Ord vertex => EdgeNetwork vertex -> List (vertex, vertex) +getEdgesForward (MkEdgeNetwork edges) = + map (\(MkSortedEdge u v, value) => if value > 0 then (u, v) else (v, u)) + (toList edges) + +||| Returns an edge with positive weight, or Nothing if the input is empty. +export +getAForwardEdge : Ord vertex => EdgeNetwork vertex -> Maybe (vertex, vertex) +getAForwardEdge = + -- TODO: Implement in a way that doesn't require reading the whole map. + head' . getEdgesForward + +||| getNeighbours net gives +w to v when (u, v) has weight w in `net`, or equivalently (v, u) has weight -w. +||| +||| TODO: Store an index so this function doesn't have to scan the whole edge network. +export +getNeighbours : Ord vertex => EdgeNetwork vertex -> vertex -> IntMap vertex +getNeighbours (MkEdgeNetwork edges) u = + let %hint + monoidInt : Monoid Int + monoidInt = Additive + in fromList [(v, w) | (MkSortedEdge u' v, w) <- toList edges, u' == u] + <+> + fromList [(v, -w) | (MkSortedEdge v u', w) <- toList edges, u' == u] + +export +getForwardNeighbours : Ord vertex => EdgeNetwork vertex -> vertex -> List vertex +getForwardNeighbours net u = + [v | (v, w) <- toList $ getNeighbours net u, w > 0] blob - /dev/null blob + 2602fcdee7c4cd9a33c6418a8efee2ccecd9fc83 (mode 644) --- /dev/null +++ src/SolidGeometry/HalfSpace.idr @@ -0,0 +1,40 @@ +||| An (affine) half-spaces is a plane where one side is "inside" and the other +||| is "outside". +||| +||| For now, this is written just for 3-D geometry, but it could probably be +||| generalized to higher dimensions. + +module SolidGeometry.HalfSpace + +%default total + +public export +data InOutSide = Inside | Outside + +export +Eq InOutSide where + Inside == Inside = True + Outside == Outside = True + _ == _ = False + +export +Ord InOutSide where + compare Inside Inside = EQ + compare Inside Outside = LT + compare Outside Inside = GT + compare Outside Outside = EQ + +export +Show InOutSide where + show Inside = "Inside" + show Outside = "Outside" + +export +flipSide : InOutSide -> InOutSide +flipSide Inside = Outside +flipSide Outside = Inside + +export +intersectSides : InOutSide -> InOutSide -> InOutSide +intersectSides Inside Inside = Inside +intersectSides _ _ = Outside blob - /dev/null blob + 37f0fca30fd9adbf1ac88838a40870225147dd4a (mode 644) --- /dev/null +++ src/SolidGeometry/Polygon.idr @@ -0,0 +1,131 @@ +module SolidGeometry.Polygon + +import Control.Monad.Maybe +import Control.Order.Fold +import Data.List +import Data.SortedSet +import Data.Vect +import Decidable.Equality +import Math.LinearAlgebra +import SolidGeometry.EdgeNetwork +import SolidGeometry.PolygonAsPaths + +%default total + +public export +data Triangle vertex = MkTriangle vertex vertex vertex + +export +Eq vertex => Eq (Triangle vertex) where + (MkTriangle a0 b0 c0) == (MkTriangle a1 b1 c1) = + (a0, b0, c0) == (a1, b1, c1) + +export +Ord vertex => Ord (Triangle vertex) where + compare (MkTriangle a0 b0 c0) (MkTriangle a1 b1 c1) = + compare (a0, b0, c0) (a1, b1, c1) + +export +Show vertex => Show (Triangle vertex) where + show (MkTriangle a b c) = show (a, b, c) + +-- Permute the vertices so that the triangle's orientation is changed. +negateTriangle : Triangle vertex -> Triangle vertex +negateTriangle (MkTriangle a b c) = MkTriangle b a c + +triangleBoundary : Ord vertex => Triangle vertex -> EdgeNetwork vertex +triangleBoundary (MkTriangle a b c) = + concat + [makeSignedEdge u v | (u, v) <- the (List _) [(a, b), (b, c), (c, a)]] + +subtractTriangle : Ord vertex => + EdgeNetwork vertex -> Triangle vertex -> EdgeNetwork vertex +subtractTriangle net triangle = + let %hint + monoidInt : Monoid Int + monoidInt = Additive + in net <+> triangleBoundary (negateTriangle triangle) + +maybeToList : Maybe a -> List a +maybeToList Nothing = [] +maybeToList (Just x) = [x] + +covering +unfoldrM' : Monad m => List a -> (b -> m (Maybe (a, b))) -> b -> m (List a) +unfoldrM' acc f x = + do Just (y, x') <- f x + | Nothing => pure (reverse acc) + unfoldrM' (y :: acc) f x' + +covering +unfoldrM : Monad m => (b -> m (Maybe (a, b))) -> b -> m (List a) +unfoldrM = unfoldrM' [] + +maybeToMaybeT : Applicative m => Maybe a -> MaybeT m a +maybeToMaybeT Nothing = nothing +maybeToMaybeT (Just x) = just x + +parameters (pos : vertex -> Vector 3 Double) + isLeftTurn : (vertex, vertex, vertex) -> Bool + isLeftTurn (a, b, c) = determinant (vectorsToMatrix $ map pos [a, b, c]) > 0 + + ||| Finds a convex angle in the polygon, i.e. a pair of edges that turns left rather than right. + findConvexAngle : PolygonAsPaths vertex -> Maybe (vertex, vertex, vertex) + findConvexAngle = + -- TODO: Avoid calling allConsecutiveEdgePairs which needs to scan the whole polygon. + find isLeftTurn . allConsecutiveEdgePairs + + ||| a b c must be a counterclockwise triangle. Returns Nothing if d isn't in + ||| the triangle a b c. Otherwise returns a number representing the distance + ||| from the a-c side to d. This number can be used to compare different + ||| `d`s. + getBaselineDistanceInTriangle : (a, b, c, d : vertex) -> Maybe Double + getBaselineDistanceInTriangle a b c d = + do let [aPos, bPos, cPos, dPos] = map pos $ the (Vect 4 vertex) [a, b, c, d] + let distanceToAB = cross [aPos, bPos] ^.^ dPos + guard (distanceToAB >= 0) + -- TODO: I can probably compute one less cross here. + let distanceToBC = cross [bPos, cPos] ^.^ dPos + guard (distanceToBC >= 0) + let distanceToCA = cross [cPos, aPos] ^.^ dPos + guard (distanceToCA >= 0) + pure distanceToCA + + -- Find a triangle `t` that can be part of a triangulation of `polygon`. + -- Return `t` and the remaining polygon, or return `Nothing` if `polygon` + -- is empty. + covering + takeTriangle : Ord vertex => Show vertex => PolygonAsPaths vertex -> + MaybeT (Either String) (Triangle vertex, PolygonAsPaths vertex) + takeTriangle polygon = + do (a, b, c) <- maybeToMaybeT $ findConvexAngle polygon + let verticesExceptABC = SortedSet.toList $ delete a $ delete b $ delete c $ allVertices polygon + let verticesInABCByDistanceToAC = + do d <- verticesExceptABC + distance <- maybeToList $ getBaselineDistanceInTriangle a b c d + pure (distance, d) + case maximum verticesInABCByDistanceToAC of + Nothing => do polygon' <- shortcut polygon a b c + pure (MkTriangle a b c, polygon') + Just (_, d) => -- Insert a detour from b to d and back. + do Just (wasAfterD, wasBeforeD) <- + pure (leftmost pos b d [(after, before) | (before, d', after) <- toTriples polygon, d' == d]) + | Nothing => throwError "takeTriangle: leftmost failed" + let polygon' = + setNextVertex a b d $ + setNextVertex b d wasAfterD $ + setNextVertex wasBeforeD d b $ + setNextVertex d b c $ + polygon + takeTriangle polygon' + + ||| Builds the given polygon out of triangles. The edges in the polygon + ||| should be oriented the same way as the triangle [1, 0, 0], [0, 1, 0], [0, + ||| 1, 0]. (If all coordinates have the form [x, y, 1], where x is right and + ||| y is up, that means counter-clockwise.) The triangles returned will also + ||| be oriented that way. + export covering + triangulate : Ord vertex => Show vertex => EdgeNetwork vertex -> Either String (List (Triangle vertex)) + triangulate polygonAsNet = + do polygonAsPaths <- netToPaths pos polygonAsNet + unfoldrM (runMaybeT . takeTriangle) polygonAsPaths blob - /dev/null blob + b4a9d17bf51621a21aa71a2f64266b61223d7edc (mode 644) --- /dev/null +++ src/SolidGeometry/PolygonAsPaths.idr @@ -0,0 +1,174 @@ +-- Represent the edges of a polygon as a map that says what vertex comes after each edge. +-- +-- E.g. if the polygon consists of the loops A->B->C->D and E->F->G, that's represented by this mapping: +-- +-- (A, B) => C +-- (B, C) => D +-- (C, D) => A +-- (D, A) => B +-- (E, F) => G +-- (F, G) => E +-- (G, E) => F + +module SolidGeometry.PolygonAsPaths + +import public Control.Monad.Error.Interface +import Control.Order.Fold +import Data.List +import Data.SortedMap +import Data.SortedSet +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.EdgeNetwork + +%default total + +export +data PolygonAsPaths vertex = MkPolygonAsPaths (SortedMap (vertex, vertex) vertex) + +export +Show vertex => Show (PolygonAsPaths vertex) where + show (MkPolygonAsPaths m) = show m + +export +toTriples : PolygonAsPaths vertex -> List (vertex, vertex, vertex) +toTriples (MkPolygonAsPaths m) = [(u, v, w) | ((u, v), w) <- SortedMap.toList m] + +empty : Ord vertex => PolygonAsPaths vertex +empty = MkPolygonAsPaths empty + +addNextVertex : PolygonAsPaths vertex -> vertex -> vertex -> vertex -> Either String (PolygonAsPaths vertex) +addNextVertex (MkPolygonAsPaths m) u v w = + case lookup (u, v) m of + Nothing => pure $ MkPolygonAsPaths $ insert (u, v) w m + Just _ => Left "Next vertex already assigned" + +export +setNextVertex : vertex -> vertex -> vertex -> PolygonAsPaths vertex -> PolygonAsPaths vertex +setNextVertex a b c (MkPolygonAsPaths m) = MkPolygonAsPaths $ insert (a, b) c m + +-- TODO: Make this fast or avoid calling it. +before : Eq vertex => Show vertex => MonadError String m => PolygonAsPaths vertex -> vertex -> vertex -> m vertex +before polygon@(MkPolygonAsPaths m) v w = + case [u' | ((u', v'), w') <- SortedMap.toList m, v == v', w == w'] of + [u] => pure u + [] => throwError ("nothing before " <+> show (w, v) <+> " in " <+> show polygon) + (_::_::_) => throwError "multiple things come before (polygon is inconsistent)" + +||| shortcut polygon a b c replaces a->b->c with a->c. +export +shortcut : Eq vertex => Show vertex => MonadError String m => PolygonAsPaths vertex -> vertex -> vertex -> vertex -> m (PolygonAsPaths vertex) +shortcut polygon@(MkPolygonAsPaths m) a b c = + do z <- before polygon a b + Just d <- pure (lookup (b, c) m) + | Nothing => throwError ("shortcut " <+> show (polygon, a, b, c) <+> ": nothing after") + -- TODO: Remove the lookup once we're fairly confident in correctness. + case lookup (a, b) m of + Nothing => throwError "Attempt to shortcut path piece that doesn't exist." + Just c' => + if c' /= c + then throwError "Attempt to shortcut a path piece that doesn't exist (it goes somewhere else)." + else if z == c && d == a + then -- The whole triangle was in the polygon. Just delete it. + pure $ MkPolygonAsPaths $ delete (c, a) $ delete (b, c) $ delete (a, b) m + else if z == c || d == a + then throwError "shortcut: inconsistent polygon" + else pure $ MkPolygonAsPaths $ insert (a, c) d $ insert (z, a) c $ delete (b, c) $ delete (a, b) m + +export +allConsecutiveEdgePairs : PolygonAsPaths vertex -> List (vertex, vertex, vertex) +allConsecutiveEdgePairs (MkPolygonAsPaths m) = + [(u, v, w) | ((u, v), w) <- SortedMap.toList m] + +||| Returns the set of all vertices in the polygon. Assumes the polygon is well-formed. +export +allVertices : Ord vertex => PolygonAsPaths vertex -> SortedSet vertex +allVertices (MkPolygonAsPaths m) = fromList [u | ((u, _), _) <- SortedMap.toList m] + +parameters (pos : vertex -> Vector 3 Double) + ||| Given an edge (u, v) and v's forward neighbours (out edges), choose the + ||| one that makes the tightest "left" turn after the edge (u, v). Returns + ||| Nothing if the list is empty. + ||| + ||| A "left" turn is one with the same orientation as [1, 0, 0], [0, 1, 0], + ||| [0, 0, 1]. If the coordinates have the form [x, y, 1] where the x axis + ||| points right and the y axis points up, these will be left turns. + ||| + ||| The extra values of type a are ignored, except the one corresponding to + ||| the chosen vertex is returned. + ||| + ||| pos u and pos v must not be colinear. + export + leftmost : {0 a : Type} -> Ord vertex => Ord a => vertex -> vertex -> List (vertex, a) -> Maybe (vertex, a) + leftmost {a} u v ws = + let uPos : Vector 3 Double + uPos = pos u + vPos : Vector 3 Double + vPos = pos v + uv : Vector 3 Double + uv = vPos ^-^ uPos + leftDirection : Vector 3 Double + leftDirection = cross [uPos, vPos] + isOnLeft : Vector 3 Double -> Bool + isOnLeft w = w ^.^ leftDirection >= 0 + wsAndPositions : List ((vertex, a), Vector 3 Double) + wsAndPositions = [(w, pos (fst w)) | w <- ws] + -- Identify the elements of ws which are a left turn rather than a right turn. + onLeft : List ((vertex, a), Vector 3 Double) + onLeft = filter (isOnLeft . snd) wsAndPositions + in + let (side, posIfLeft) = if isNil onLeft + then (wsAndPositions, -1.0) + else (onLeft, 1.0) + in + let directionOfIncreasingAngle : Vector 3 Double + directionOfIncreasingAngle = posIfLeft *^ cross [leftDirection, vPos] + sideRanked : List (Double, (vertex, a)) + sideRanked = + [ (directionOfIncreasingAngle ^.^ normalize wPos, w) + | (w, wPos) <- side + ] + in + map snd $ minimum sideRanked + + ||| Helper for takeOneLoop. (u, v) is the next edge to be added. + covering + finishTakingLoop : Ord vertex => Show vertex => PolygonAsPaths vertex -> EdgeNetwork vertex -> vertex -> vertex -> vertex -> vertex -> Either String (EdgeNetwork vertex, PolygonAsPaths vertex) + finishTakingLoop polygonSoFar net start afterStart u v = + do let net' = net <+> makeSignedEdge v u -- Subtract (u, v) + False <- pure (v == start) + | True => + do polygonSoFar' <- addNextVertex polygonSoFar u v afterStart + pure (net', polygonSoFar') + (w, ()) <- case leftmost u v [(w', ()) | w' <- getForwardNeighbours net' v] of + Nothing => Left ("no next vertex after (" <+> show u <+> ", " <+> show v <+> ") in network " <+> show net') + Just w => pure w + polygonSoFar' <- addNextVertex polygonSoFar u v w + finishTakingLoop polygonSoFar' net' start afterStart v w + + ||| Returns Right Nothing if the edge network is empty. + covering + takeOneLoop : Ord vertex => Show vertex => PolygonAsPaths vertex -> EdgeNetwork vertex -> Either String (Maybe (EdgeNetwork vertex, PolygonAsPaths vertex)) + takeOneLoop polygonSoFar net = + case getAForwardEdge net of + Nothing => pure Nothing + Just (u, v) => case finishTakingLoop polygonSoFar net u v u v of + Left e => + Left ("takeOneLoop from " <+> show net <+> ": " <+> e) + Right v => Right (Just v) + + ||| Like netToPaths, but adds to an existing polygon. + covering + netToPaths' : Ord vertex => Show vertex => PolygonAsPaths vertex -> EdgeNetwork vertex -> Either String (PolygonAsPaths vertex) + netToPaths' polygonSoFar net = + do Just (net', polygonSoFar') <- takeOneLoop polygonSoFar net + | Nothing => pure polygonSoFar + netToPaths' polygonSoFar' net' + + ||| Convert an EdgeNetwork representation of a ploygon to a PolygonAsPaths. Returns Right _ if successful. Fails with an error message if an edge is present more than once in the input, or if the input isn't made of loops. + export covering + netToPaths : Ord vertex => Show vertex => EdgeNetwork vertex -> Either String (PolygonAsPaths vertex) + netToPaths net = + case netToPaths' empty net of + Left e => Left ("netToPaths " <+> show net <+> ": " <+> e) + Right x => Right x blob - /dev/null blob + be5929f32a1c4ef0531c77ac5c66d035a06070b4 (mode 644) --- /dev/null +++ src/SolidGeometry/PolygonTest.idr @@ -0,0 +1,210 @@ +module SolidGeometry.PolygonTest + +import Data.Fin +import Data.List +import Data.Vect +import Math.LinearAlgebra +import SolidGeometry.EdgeNetwork +import SolidGeometry.Polygon +import Test.TestTree + +%default total + +canonicalizeTriangle : Ord vertex => Triangle vertex -> Triangle vertex +canonicalizeTriangle (MkTriangle a b c) = + if a <= b + then if a <= c + then MkTriangle a b c + else MkTriangle c a b + else if b <= c + then MkTriangle b c a + else MkTriangle c a b + +canonicalizeTriangles : Ord vertex => List (Triangle vertex) -> List (Triangle vertex) +canonicalizeTriangles = sort . map canonicalizeTriangle + +record TriangulationTestCase where + constructor MkTriangulationTestCase + name : String + {numPoints : Nat} + -- TODO: Add test cases using the full 3D coordinate space. + points : Vect numPoints (Vector 2 Double) + edges : List (Fin numPoints, Fin numPoints) + correctCanonicalTriangulations : List (List (Triangle (Fin numPoints))) + +covering +triangulationTestCase : TriangulationTestCase -> TestCase +triangulationTestCase testCase = + let vertex : Type + vertex = Fin testCase.numPoints + + pos : vertex -> Vector 3 Double + pos i = vectToVector $ vectorToVect (index i testCase.points) ++ [1] + + edges : EdgeNetwork vertex + edges = concat [makeSignedEdge u v | (u, v) <- testCase.edges] + + result : TestResult + result = + case triangulate pos edges of + Left error => MkFailure $ "triangulate returned error: " <+> error + Right triangles => + if elem (canonicalizeTriangles triangles) testCase.correctCanonicalTriangulations + then MkSuccess + else MkFailure $ "Bad triangulation: " <+> show triangles + in MkTestCase testCase.name (pure result) + +-- 3<-2 +-- | ^ +-- v | +-- 0->1 +square : TriangulationTestCase +square = + let vertex : Type + vertex = Fin 4 + + correctAnswers : List (List (Triangle vertex)) + correctAnswers = + [ [MkTriangle 0 1 2, MkTriangle 0 2 3] + , [MkTriangle 0 1 3, MkTriangle 1 2 3] + ] + in + MkTriangulationTestCase + { name = "triangulate a square" + , points = [[0, 0], [1, 0], [1, 1], [0, 1]] + , edges = [(0, 1), (1, 2), (2, 3), (3, 0)] + , correctCanonicalTriangulations = correctAnswers + } + +-- 3<-----2 +-- | ^ +-- | 5->6 | +-- | ^ | | +-- | | v | +-- | 4<-7 | +-- v | +-- 0----->1 +loop : TriangulationTestCase +loop = + let vertex : Type + vertex = Fin 8 + + correctAnswers : List (List (Triangle vertex)) + correctAnswers = + do south <- + [ [MkTriangle 0 1 4, MkTriangle 1 7 4] + , [MkTriangle 0 1 7, MkTriangle 0 7 4] + ] + east <- + [ [MkTriangle 1 2 7, MkTriangle 2 6 7] + , [MkTriangle 1 2 6, MkTriangle 1 6 7] + ] + north <- + [ [MkTriangle 2 3 6, MkTriangle 3 5 6] + , [MkTriangle 2 3 5, MkTriangle 2 5 6] + ] + west <- + [ [MkTriangle 0 4 3, MkTriangle 3 4 5] + , [MkTriangle 0 4 5, MkTriangle 0 5 3] + ] + pure $ canonicalizeTriangles (south <+> east <+> north <+> west) + in + MkTriangulationTestCase + { name = "triangulate a square loop" + , points = [[-2, -2], [2, -2], [2, 2], [-2, 2], [-1, -1], [-1, 1], [1, 1], [1, -1]] + , edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4)] + , correctCanonicalTriangulations = correctAnswers + } + +-- Similar to loop, but with different coordinates and with the vertices given in a different order (indices 5 and 7 are) swapped. +anotherLoop : TriangulationTestCase +anotherLoop = + let vertex : Type + vertex = Fin 8 + + correctAnswers : List (List (Triangle vertex)) + correctAnswers = + do south <- + [ [MkTriangle 0 1 4, MkTriangle 1 5 4] + , [MkTriangle 0 1 5, MkTriangle 0 5 4] + ] + east <- + [ [MkTriangle 1 2 5, MkTriangle 2 6 5] + , [MkTriangle 1 2 6, MkTriangle 1 6 5] + ] + north <- + [ [MkTriangle 2 3 6, MkTriangle 3 7 6] + , [MkTriangle 2 3 7, MkTriangle 2 7 6] + ] + west <- + [ [MkTriangle 0 4 3, MkTriangle 3 4 7] + , [MkTriangle 0 4 7, MkTriangle 0 7 3] + ] + pure $ canonicalizeTriangles (south <+> east <+> north <+> west) + in + MkTriangulationTestCase + { name = "triangulate a square loop" + , points = [ [0, 0.3], [0.1, 0.3], [0.1, 0.4], [0, 0.4] + , [0.025, 0.325], [0.075, 0.325], [0.075, 0.375], [0.025, 0.375] + ] + , edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 7), (7, 6), (6, 5), (5, 4)] + , correctCanonicalTriangulations = correctAnswers + } + +-- Same as loop, but with a little island at the very middle: +-- +-- 11<-10 +-- | ^ +-- v | +-- 8->9 +loopWithIsland : TriangulationTestCase +loopWithIsland = + let vertex : Type + vertex = Fin 12 + + correctAnswers : List (List (Triangle vertex)) + correctAnswers = + do south <- + [ [MkTriangle 0 1 4, MkTriangle 1 7 4] + , [MkTriangle 0 1 7, MkTriangle 0 7 4] + ] + east <- + [ [MkTriangle 1 2 7, MkTriangle 2 6 7] + , [MkTriangle 1 2 6, MkTriangle 1 6 7] + ] + north <- + [ [MkTriangle 2 3 6, MkTriangle 3 5 6] + , [MkTriangle 2 3 5, MkTriangle 2 5 6] + ] + west <- + [ [MkTriangle 0 4 3, MkTriangle 3 4 5] + , [MkTriangle 0 4 5, MkTriangle 0 5 3] + ] + island <- + [ [MkTriangle 8 9 10, MkTriangle 8 10 11] + , [MkTriangle 8 9 11, MkTriangle 9 10 11] + ] + pure $ canonicalizeTriangles (south <+> east <+> north <+> west <+> island) + in + MkTriangulationTestCase + { name = "triangulate a square loop with an island" + , points = [ [-2, -2], [2, -2], [2, 2], [-2, 2] + , [-1, -1], [-1, 1], [1, 1], [1, -1] + , [-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5] + ] + , edges = [ (0, 1), (1, 2), (2, 3), (3, 0) + , (4, 5), (5, 6), (6, 7), (7, 4) + , (8, 9), (9, 10), (10, 11), (11, 8) + ] + , correctCanonicalTriangulations = correctAnswers + } + +export covering +tests : TestTree +tests = + MkTestGroup "Polygon.trangulate" $ map (MkTestLeaf . triangulationTestCase) + [ square + , loop + , anotherLoop + , loopWithIsland + ] blob - /dev/null blob + 409a6ad58c44ca7b36a383df9614500e1ad5bb84 (mode 644) --- /dev/null +++ src/SolidGeometry/Polyhedron.idr @@ -0,0 +1,426 @@ +module SolidGeometry.Polyhedron + +import Control.Monad.State +import Data.List +import Data.SortedMap +import Data.Vect +import Math.LinearAlgebra +import Math.LinearAlgebra.Normalized +import SolidGeometry.Circle +import SolidGeometry.EdgeNetwork +import SolidGeometry.HalfSpace +import SolidGeometry.Polygon +import SolidGeometry.Sign + +%default total + +||| A polyhedron is defined by a list of half-spaces and a function saying +||| whether a point is inside or outside the polyhedron, given which side of +||| each half-space it's on. +public export +record Polyhedron (numPlanes : Nat) where + constructor MkPolyhedron + ||| Each plane is defined by a linear map from the 2-sphere in R^3 to the + ||| 3-sphere in R^4. (We don't actually require the matrix to preserve + ||| norms.) + ||| + ||| The normal vector of a plane is the cross product of its rows. Another + ||| point is on the "Positive side" of the plane iff its dot product with + ||| that normal vector is positive. + planes : Vect numPlanes (Matrix 3 4 Double) + ||| Given knowledge of which side of each plane a point is on, determines + ||| whether the point is inside this polyhedron. + ||| + ||| For example, if this polyhedron is the intersection of the positive sides + ||| of all the planes, this function should return Inside iff the input is + ||| all Positive. + inside : (Fin numPlanes -> Sign) -> InOutSide + +||| A side of one of the half-spaces used to define a polyhedron. +public export +record FaceId (numPlanes : Nat) where + constructor MkFaceId + ||| Index into the planes field of the original Polyhedron. + planeIndex : Fin numPlanes + ||| Positive iff the "outside" of this face is the positive side of that + ||| plane. + outside : Sign + +Show (FaceId numPlanes) where + show face = show (face.planeIndex, face.outside) + +||| Coordinates of a vertex of the polyhedron. +||| +||| TODO: Rethink adjustment of coordinates. Normalizing them to magnitude 1 +||| probably isn't the right solution because I expect that will mess up +||| interpolation, e.g. making textures look weird (though I haven't confirmed +||| this). Probably a better approach would be to normalize so that +||| faceCoordinates have the form [x,y,+/-1] whenever possible; if we then use +||| x and y for texture coordinates then everything will make sense. +public export +record VertexCoordinates where + constructor MkVertexCoordinates + ||| The vertex's position on the 3-sphere. + absolute : Normalized (Vector 4 Double) + ||| The vertex's position in the face's coordinate space. + faceCoordinates : Normalized (Vector 3 Double) + +export +Show VertexCoordinates where + show vc = fastConcat + [ "MkVertexCoordinates { absolute = " + , show vc.absolute + , ", faceCoordinates = " + , show vc.faceCoordinates + , " }" + ] + +||| A way to triangulate a single face of a Polyhedron. The "face" might not be +||| connected. Really this is all the faces that come from a particular side of +||| a particular one of the half-spaces used to define the polyhedron. +public export +record TriangulatedFace (numPlanes : Nat) where + constructor MkTriangulatedFace + faceId : FaceId numPlanes + {0 numVertices : Nat} + ||| For each vertex: its position on the 3-sphere, and its position in the + ||| parameter space of the face. + vertices : Vect numVertices VertexCoordinates + triangles : List (Triangle (Fin numVertices)) + +export +emptyPolyhedron : Polyhedron 0 +emptyPolyhedron = MkPolyhedron {planes = Nil, inside = const Outside} + +faceIsEmpty : TriangulatedFace n -> Bool +faceIsEmpty f = null f.triangles + +data IsNonZero : (0 n : Nat) -> Type where + IsSuccessorOf : (0 predecessor : Nat) -> IsNonZero (S predecessor) + +finImpliesNonZero : Fin n -> IsNonZero n +finImpliesNonZero FZ = IsSuccessorOf _ +finImpliesNonZero (FS _) = IsSuccessorOf _ + +reverseListToVect' : (n0 : Nat) -> Vect n0 a -> List a -> (n1 : Nat ** Vect n1 a) +reverseListToVect' n0 acc [] = (n0 ** acc) +reverseListToVect' n0 acc (h :: t) = reverseListToVect' (S n0) (h :: acc) t + +reverseListToVect : List a -> (n : Nat ** Vect n a) +reverseListToVect = reverseListToVect' 0 [] + +mapToVect : SortedMap key value -> (n : Nat ** Vect n (key, value)) +mapToVect = reverseListToVect . toList + +vectToIndexMap' : Fin (S n0) -> SortedMap a (Fin (n0 + n1)) -> Vect n1 a -> SortedMap a (Fin (n0 + n1)) +vectToIndexMap' _ acc [] = acc +vectToIndexMap' {n1 = S n1'} nextIndex acc (h :: t) = + let eq : (S n0 + n1' = n0 + S n1') + in + rewrite sym eq in + vectToIndexMap' + (FS nextIndex) + (insert h (weakenN n1' nextIndex) (rewrite eq in acc)) + t + +vectToIndexMap : Ord a => Vect n a -> SortedMap a (Fin n) +vectToIndexMap = vectToIndexMap' {n0=Z} FZ empty + +||| Given a nonzero vector, returns a vector that's not colinear with it. +notColinear : Abs a => Ord a => {n : Nat} -> Vector (S (S n)) a -> Vector (S (S n)) a +notColinear v = + let (a :: b :: t) = vectorToVect v in + vectToVector $ + if abs a >= abs b + then 0 :: 1 :: replicate _ 0 + else 1 :: 0 :: replicate _ 0 + +||| Given a normal vector v, returns a matrix m such that v is a positive +||| scalar times the cross product of the rows of m. +implicitToParameterized : Abs a => Neg a => Ord a => Vector 3 a -> Matrix 2 3 a +implicitToParameterized v = + let start : Vector 3 a + start = notColinear v + x : Vector 3 a + x = cross [start, v] + y : Vector 3 a + y = cross [v, x] + in vectorsToMatrix [x, y] + +||| Given a plane P1 and an implicitly-defined plane P2, returns a +||| parameterized representation of the circle where they intersect, both on +||| the 3-sphere and in P1's parameter space. +||| +||| The resulting circle will be oriented as follows. If [a, b] are the rows of +||| the first matrix returned and c is the normal vector of P2, Then the dot +||| product of cross [a, b, c] and the cross product of the rows of P1's matrix +||| will be positive. +intersectPlanes : Abs a => Neg a => Ord a => Matrix 3 4 a -> Vector 4 a -> (Matrix 2 4 a, Matrix 2 3 a) +intersectPlanes p1 p2 = + let intersectionInP1ParameterSpace : Matrix 2 3 a + intersectionInP1ParameterSpace = implicitToParameterized (p1 !*^ p2) + in (intersectionInP1ParameterSpace !*! p1, intersectionInP1ParameterSpace) + +||| Convert a `Fin n` value to a `Fin (S n)` value by skipping a particular +||| value (`skip`) in the range. +skipFin : Fin (S n) -> Fin n -> Fin (S n) +skipFin skip x = + let x' : Fin (S n) + x' = weaken x + in if x' < skip then x' else FS x + +data OrderingWithProof : Nat -> Nat -> Type where + LTProof : (0 prf : x `LT` y) -> OrderingWithProof x y + EQProof : (0 prf : x = y) -> OrderingWithProof x y + GTProof : (0 prf : x `GT` y) -> OrderingWithProof x y + +0 sucPred : (n : Nat) -> m `LT` n -> S (pred n) = n +sucPred (S n) (LTESucc _) = Refl + +0 ltToLT : (x, y : Nat) -> LT = compare x y -> x `LT` y +ltToLT Z (S _) _ = LTESucc LTEZero +ltToLT (S x') (S y') eq = LTESucc $ ltToLT x' y' eq + +0 eqToEQ : (x, y : Nat) -> EQ = compare x y -> x = y +eqToEQ Z Z _ = Refl +eqToEQ (S x') (S y') eq = cong S $ eqToEQ x' y' eq + +0 gtToGT : (x, y : Nat) -> GT = compare x y -> x `GT` y +gtToGT Z Z _ impossible +gtToGT Z (S _) _ impossible +gtToGT (S _) Z _ = LTESucc LTEZero +gtToGT (S x') (S y') eq = LTESucc $ gtToGT x' y' eq + +compareWithProof' : (x, y : Nat) -> (o : Ordering) -> (o = compare x y) -> OrderingWithProof x y +compareWithProof' x y LT e = LTProof $ ltToLT x y e +compareWithProof' x y EQ e = EQProof $ eqToEQ x y e +compareWithProof' x y GT e = GTProof $ gtToGT x y e + +compareWithProof : (x, y : Nat) -> OrderingWithProof x y +compareWithProof x y = compareWithProof' x y (compare x y) Refl + +0 finToNatLTE : (n : Nat) -> (x : Fin (S n)) -> finToNat x `LTE` n +finToNatLTE _ FZ = LTEZero +finToNatLTE (S n') (FS x') = LTESucc $ finToNatLTE n' x' + +||| The opposite of `skipFin`. Convert a `Fin (S n)` value to a `Fin n` value +||| by skipping a particular element of the domain (`skipped`). If the input is +||| `skipped` itself, return `Nothing`. +unskipFin : (skipped : Fin (S n)) -> Fin (S n) -> Maybe (Fin n) +unskipFin skipped x = + let xNat : Nat + xNat = finToNat x + 0 xLteN : xNat `LTE` n + xLteN = finToNatLTE n x + in case compareWithProof (finToNat skipped) xNat of + LTProof lt => + let 0 sucPredX : (S (pred xNat) = xNat) + sucPredX = sucPred xNat lt + 0 predXLtN : (pred xNat) `LT` n + predXLtN = rewrite sucPredX in xLteN + in Just $ natToFinLT (pred xNat) {prf=predXLtN} + EQProof _ => Nothing + GTProof xLtSkipped => + let 0 skippedLteN : finToNat skipped `LTE` n + skippedLteN = finToNatLTE n skipped + 0 xLtN : xNat `LT` n + xLtN = transitive xLtSkipped skippedLteN + in Just $ natToFinLT xNat {prf=xLtN} + +sequenceLazy' : Monad f => List a -> List (Lazy (f a)) -> f (List a) +sequenceLazy' acc [] = pure $ reverse acc +sequenceLazy' acc (h :: t) = + do x <- h + sequenceLazy' (x :: acc) t + +sequenceLazy : Monad f => List (Lazy (f a)) -> f (List a) +sequenceLazy = sequenceLazy' [] + +perturbPolyhedron : {numPlanes : Nat} -> Polyhedron numPlanes -> Polyhedron numPlanes +perturbPolyhedron polyhedron = + let eps : Double + eps = 0.0001 + perturbation : Double -> Matrix 3 4 Double + perturbation i = + let v : Vect 4 Double + -- Awful hack. Arbitrary numbers. + v = [eps, 1.37 * eps, i / cast numPlanes * eps, 0.92 * eps] + in vectsToMatrix (replicate 3 v) + in { planes := + [ p !+! perturbation (cast $ finToNat i) + | (i, p) <- zip range polyhedron.planes + ] + } + polyhedron + +export partial +triangulate : {numPlanes : Nat} -> (polyhedron : Polyhedron numPlanes) -> Either String (List (TriangulatedFace numPlanes)) +triangulate {numPlanes = Z} _ = pure [] +triangulate {numPlanes = numPlanes@(S numOtherPlanes)} originalPolyhedron = + let -- Perturb originalPolyhedron to avoid degenerate cases. + polyhedron : Polyhedron numPlanes + polyhedron = perturbPolyhedron originalPolyhedron + planeNormals : Vect numPlanes (Vector 4 Double) + planeNormals = map (cross . matrixToVectors) polyhedron.planes + partial + triangulateFace : + FaceId numPlanes -> + Either String (TriangulatedFace numPlanes) + triangulateFace faceId = + let 0 OtherPlaneIndex : Type + OtherPlaneIndex = Fin numOtherPlanes + ||| An intersection of two other planes with the plane of this + ||| face. Swapping the indices gives the opposite vertex. + 0 Vertex : Type + Vertex = (OtherPlaneIndex, OtherPlaneIndex) + ||| A point on the plane corresponding to this face is on the + ||| face itself if the faceId.outside side of this plane is + ||| Outside the polyhedron and the other side of this plane is + ||| Inside. + isOnFace : (Fin numOtherPlanes -> Sign) -> Bool + isOnFace f = + let insertMySide : Sign -> Fin numPlanes -> Sign + insertMySide sign planeIndex = + case unskipFin faceId.planeIndex planeIndex of + Nothing => sign + Just otherPlaneIndex => f otherPlaneIndex + in + polyhedron.inside (insertMySide faceId.outside) == Outside && + polyhedron.inside (insertMySide (flipSign faceId.outside)) == Inside + thisPlane : Matrix 3 4 Double + thisPlane = + index faceId.planeIndex polyhedron.planes + otherPlaneNormals : Vect numOtherPlanes (Vector 4 Double) + otherPlaneNormals = deleteAt faceId.planeIndex planeNormals + ||| Given the index of another plane, returns the edges of this + ||| face where it meets that other plane. Also fills in the + ||| locations of any vertices that got used. + edgesFromOtherPlane : + {numOtherPlanes' : Nat} -> + {auto 0 isNumOtherPlanes : numOtherPlanes' = numOtherPlanes} -> + Fin numOtherPlanes' -> + State (SortedMap Vertex VertexCoordinates) (List (Vertex, Vertex)) + edgesFromOtherPlane {numOtherPlanes'} otherPlaneIndex with (finImpliesNonZero otherPlaneIndex) + edgesFromOtherPlane {numOtherPlanes' = S numThirdPlanes} otherPlaneIndex | IsSuccessorOf numThirdPlanes = + let (line3d, line2d) = + intersectPlanes + thisPlane $ index otherPlaneIndex $ + rewrite isNumOtherPlanes in + otherPlaneNormals + in + let thirdPlaneToOtherPlaneIndex : Fin numThirdPlanes -> Fin numOtherPlanes + thirdPlaneToOtherPlaneIndex = + rewrite sym isNumOtherPlanes in skipFin otherPlaneIndex + thirdPlaneNormals : Vect numThirdPlanes (Vector 4 Double) + thirdPlaneNormals = + deleteAt otherPlaneIndex $ + rewrite isNumOtherPlanes in otherPlaneNormals + -- Intersections of these two planes with other planes, in + -- line3d's (and line2d's) parameter space. + intersections : Vect numThirdPlanes (Vector 2 Double) + intersections = [negateVector $ cross [line3d !*^ thirdPlane] | thirdPlane <- thirdPlaneNormals] + edgeDirection : (Fin numThirdPlanes -> Sign) -> Maybe Sign + edgeDirection f = + let insertMySide : Sign -> Fin numOtherPlanes -> Sign + insertMySide sign planeIndex = + case unskipFin otherPlaneIndex (rewrite isNumOtherPlanes in planeIndex) of + Nothing => sign + Just thirdPlaneIndex => f thirdPlaneIndex + in + case (isOnFace (insertMySide Positive), isOnFace (insertMySide Negative)) of + (False, True) => Just Positive + (True, False) => Just Negative + _ => Nothing + in + case boundariesToEdges edgeDirection intersections of + -- TODO: handle the AllOneWay case. + Edges edges => + let otherPlaneIndex' : Fin numOtherPlanes + otherPlaneIndex' = + rewrite sym isNumOtherPlanes in otherPlaneIndex + storeVertexLocation : + (Fin numThirdPlanes, Sign) -> Vertex -> + State (SortedMap Vertex VertexCoordinates) () + storeVertexLocation (i, sign) vertex = + do old <- gets (lookup vertex) + case old of + Just _ => pure () + Nothing => + let unsignedPointOnCircle : Vector 2 Double + -- TODO: Use an array to avoid costly lookup. + unsignedPointOnCircle = index i intersections + pointOnCircle : Vector 2 Double + pointOnCircle = + case sign of + Positive => unsignedPointOnCircle + Negative => negateVector unsignedPointOnCircle + in + modify $ insert vertex $ + MkVertexCoordinates + { absolute = makeNormalized (pointOnCircle ^*! line3d) + , faceCoordinates = makeNormalized (pointOnCircle ^*! line2d) + } + convertVertex : (Fin numThirdPlanes, Sign) -> Vertex + convertVertex (index, sign) = + let index' : Fin numOtherPlanes + index' = thirdPlaneToOtherPlaneIndex index + in + case sign of + Positive => (otherPlaneIndex', index') + Negative => (index', otherPlaneIndex') + in + sequence + [ let u' : Vertex + u' = convertVertex u + v' : Vertex + v' = convertVertex v + in + do storeVertexLocation u u' + storeVertexLocation v v' + pure (u', v') + | (u, v) <- edges + ] + edgesFromAllOtherPlanes : + State (SortedMap Vertex VertexCoordinates) (List (Vertex, Vertex)) + edgesFromAllOtherPlanes = map concat $ sequence $ map edgesFromOtherPlane range + positionsAndEdges : ( SortedMap Vertex VertexCoordinates + , List (Vertex, Vertex) + ) + positionsAndEdges = runState empty edgesFromAllOtherPlanes + in + let (positions, edges) = positionsAndEdges in + let (numUsedVertices ** usedVertices) = mapToVect positions in + let vertexPositions : Vect numUsedVertices VertexCoordinates + vertexPositions = map snd usedVertices + vertexToFaceCoordinates : Fin numUsedVertices -> Vector 3 Double + vertexToFaceCoordinates i = + fromNormalized (index i vertexPositions).faceCoordinates + vertexToIndexMap : SortedMap Vertex (Fin numUsedVertices) + vertexToIndexMap = vectToIndexMap (map fst usedVertices) + partial + vertexToIndex : Vertex -> Fin numUsedVertices + vertexToIndex v = case lookup v vertexToIndexMap of Just i => i + edgeNetwork : EdgeNetwork (Fin numUsedVertices) + edgeNetwork = concat [makeSignedEdge (vertexToIndex u) (vertexToIndex v) | (u, v) <- edges] + in + do triangles <- + mapFst + ( \err => + fastConcat + [ "face ", show faceId + , " with vertices ", show usedVertices + , ": ", err + ] + ) + (Polygon.triangulate vertexToFaceCoordinates edgeNetwork) + pure $ + MkTriangulatedFace + { faceId = faceId + , vertices = vertexPositions + , triangles = triangles + } + in do triangulatedFaces <- + sequenceLazy [ delay $ triangulateFace (MkFaceId i s) + | i <- toList range, s <- [Positive, Negative] + ] + pure $ filter (not . faceIsEmpty) $ triangulatedFaces blob - /dev/null blob + 718abb83019daa67ba416e77f6a015be6416b771 (mode 644) --- /dev/null +++ src/SolidGeometry/Sign.idr @@ -0,0 +1,33 @@ +module SolidGeometry.Sign + +%default total + +public export +data Sign = Negative | Positive + +export +Eq Sign where + Negative == Negative = True + Positive == Positive = True + _ == _ = False + +export +Ord Sign where + compare Negative Negative = EQ + compare Negative Positive = LT + compare Positive Negative = GT + compare Positive Positive = EQ + +export +Show Sign where + show Negative = "Negative" + show Positive = "Positive" + +export +flipSign : Sign -> Sign +flipSign Positive = Negative +flipSign Negative = Positive + +export +softSign : Double -> Sign +softSign x = if x >= 0 then Positive else Negative blob - /dev/null blob + 1dcf7dfbcf9c0e2a660ce5601ceecfd3258a3cee (mode 644) --- /dev/null +++ src/Test/TestTree.idr @@ -0,0 +1,52 @@ +module Test.TestTree + +%default total + +public export +data TestResult = MkSuccess | MkFailure String + +public export +record TestCase where + constructor MkTestCase + name : String + run : IO TestResult + +runTestCase : TestCase -> IO TestResult +runTestCase testCase = + do putStrLn $ "Running test case: " <+> testCase.name + result <- testCase.run + case result of + MkSuccess => putStrLn "Success" + MkFailure message => putStrLn $ "Failure: " <+> message + pure result + +public export +data TestTree = MkTestLeaf TestCase | MkTestGroup String (List TestTree) + +mutual + -- Helper. We can't just call use sequence . map runTestTree in runTestTree because the compiler wouldn't be able to tell runTestTree is total. + runTestTrees : List TestResult -> List TestTree -> IO (List TestResult) + runTestTrees acc [] = pure $ reverse acc + runTestTrees acc (tree :: trees) = + do result <- runTestTree tree + runTestTrees (result :: acc) trees + + export + runTestTree : TestTree -> IO TestResult + runTestTree (MkTestLeaf testCase) = runTestCase testCase + runTestTree (MkTestGroup name subTrees) = + let resultToNumFailures : TestResult -> Nat + resultToNumFailures MkSuccess = 0 + resultToNumFailures (MkFailure _) = 1 + in + do putStrLn $ "Starting test group: " <+> name + results <- runTestTrees [] subTrees + let numFailures = sum $ map resultToNumFailures results + putStrLn $ + "Test group finished: " <+> name <+> ". " <+> + show numFailures <+> "/" <+> show (length results) <+> + " failed." + pure $ + if numFailures == 0 + then MkSuccess + else MkFailure name blob - /dev/null blob + 1d11170cf07a0c4840e89654cbe975a8bc5e84a0 (mode 644) --- /dev/null +++ src/TestMain.idr @@ -0,0 +1,20 @@ +module TestMain + +import Math.LinearAlgebraTest +import SolidGeometry.CircleTest +import SolidGeometry.PolygonTest +import Test.TestTree + +%default total + +tests : TestTree + +main : IO () +main = + ignore $ + runTestTree $ + MkTestGroup "all tests" + [ LinearAlgebraTest.tests + , CircleTest.tests + , PolygonTest.tests + ]