Last active
December 21, 2024 03:39
-
-
Save axionbuster/8b7017bc89a3e03aaa751f0ba88c6ee1 to your computer and use it in GitHub Desktop.
Grid ray marching with floating-point error compensation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| -- | march along a line segment, finding all intersections with grid points | |
| module Lib (march) where | |
| import Control.Monad.Fix | |
| import Control.Monad.ST.Lazy | |
| import Data.Functor | |
| import Data.Functor.Rep | |
| import Data.STRef.Lazy | |
| import Prelude hiding (read) | |
| -- | march along a line segment, finding all intersections | |
| -- with grid points | |
| -- | |
| -- the starting point is NOT included in the output | |
| -- | |
| -- a compensated sum is used to reduce floating point error | |
| -- | |
| -- the returned list being infinite, it is recommended to | |
| -- use 'take' to limit the number of points to be computed | |
| march :: | |
| forall f a. | |
| ( Applicative f, | |
| Foldable f, | |
| Representable f, | |
| RealFloat a | |
| ) => | |
| -- | starting point | |
| f a -> | |
| -- | direction (no need to be normalized) | |
| f a -> | |
| -- | list of (delta time, point) pairs | |
| [(a, f a)] | |
| march start direction = runST do | |
| let fi = fromIntegral :: Int -> a | |
| new = newSTRef | |
| read = readSTRef | |
| write = writeSTRef | |
| minnonan a b | |
| | isNaN a = b | |
| | isNaN b = a | |
| | otherwise = min a b -- if both are NaN, then pick either | |
| minimum_ = foldr1 minnonan | |
| cur <- new start | |
| com <- new $ pure 0 -- Kahan sum compensator | |
| let sig = floor . signum <$> direction | |
| fix \this -> do | |
| -- mechanism: | |
| -- using the parametric equation of the line segment | |
| -- find the closest intersection with the grid -> get 'time' value | |
| -- then use the 'time' to get the coordinates of the intersection | |
| let (!) = index | |
| t cur' i = | |
| -- solve for time to next intersection in dimension i | |
| let s = fi (truncate (cur' ! i) + sig ! i) - cur' ! i | |
| in s / direction ! i | |
| add c x y = | |
| -- Kahan's compensated sum (x += y) | |
| let y' = y - c | |
| s = x + y' | |
| c' = (s - x) - y' | |
| in (s, c') | |
| com' <- read com | |
| cur' <- read cur | |
| let tim = minimum_ $ tabulate @f $ t cur' | |
| vadd v w = tabulate @f \i -> | |
| -- elementwise error-compensated vector addition | |
| add (com' ! i) (v ! i) (w ! i) | |
| s = vadd cur' $ direction <&> (* tim) | |
| newcur = fst <$> s | |
| newcom = snd <$> s | |
| write cur newcur | |
| write com newcom | |
| ((tim, newcur) :) <$> this | |
| {-# INLINEABLE march #-} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
the only dependency besides 'base' is 'adjunctions'