-- |
-- Module: M.Collision.Internal.March2
-- Description: Grid-based ray marching implementation using digital differential analyzer
-- Copyright: (c) axionbuster, 2025
-- License: BSD-3-Clause
--
-- Alternative implementation of ray marching algorithm focused on performance.
module M.Collision.Internal.March2
  ( VHit (..),
    isnotahit,
    isahit,
    march,
  )
where

import Control.Lens
import Data.Int
import Debug.Trace
import Linear
import Text.Printf

-- | voxel hit data structure. may encode hit or no hit.
data VHit i a = VHit
  { -- | finite nonnegative float on hit; otherwise
    -- will be positive infinity and there will be
    -- no hit (and the other fields will be
    -- defined, but meaningless and unspecified)
    forall i a. VHit i a -> a
vhittim :: !a,
    -- | integer coordinate location of where it hit
    forall i a. VHit i a -> V3 i
vhitloc :: !(V3 i),
    -- | the normal vectors (signum; opposite to
    -- the displacement)
    forall i a. VHit i a -> V3 i
vhitnor :: !(V3 i)
  }
  deriving (VHit i a -> VHit i a -> Bool
(VHit i a -> VHit i a -> Bool)
-> (VHit i a -> VHit i a -> Bool) -> Eq (VHit i a)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall i a. (Eq a, Eq i) => VHit i a -> VHit i a -> Bool
$c== :: forall i a. (Eq a, Eq i) => VHit i a -> VHit i a -> Bool
== :: VHit i a -> VHit i a -> Bool
$c/= :: forall i a. (Eq a, Eq i) => VHit i a -> VHit i a -> Bool
/= :: VHit i a -> VHit i a -> Bool
Eq, Int -> VHit i a -> String -> String
[VHit i a] -> String -> String
VHit i a -> String
(Int -> VHit i a -> String -> String)
-> (VHit i a -> String)
-> ([VHit i a] -> String -> String)
-> Show (VHit i a)
forall a.
(Int -> a -> String -> String)
-> (a -> String) -> ([a] -> String -> String) -> Show a
forall i a. (Show a, Show i) => Int -> VHit i a -> String -> String
forall i a. (Show a, Show i) => [VHit i a] -> String -> String
forall i a. (Show a, Show i) => VHit i a -> String
$cshowsPrec :: forall i a. (Show a, Show i) => Int -> VHit i a -> String -> String
showsPrec :: Int -> VHit i a -> String -> String
$cshow :: forall i a. (Show a, Show i) => VHit i a -> String
show :: VHit i a -> String
$cshowList :: forall i a. (Show a, Show i) => [VHit i a] -> String -> String
showList :: [VHit i a] -> String -> String
Show)

-- | decide if something is not a hit
isnotahit :: (RealFloat a) => VHit i a -> Bool
isnotahit :: forall a i. RealFloat a => VHit i a -> Bool
isnotahit = a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite (a -> Bool) -> (VHit i a -> a) -> VHit i a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VHit i a -> a
forall i a. VHit i a -> a
vhittim

-- | decide if something is a hit
isahit :: (RealFloat a) => VHit i a -> Bool
isahit :: forall a i. RealFloat a => VHit i a -> Bool
isahit = Bool -> Bool
not (Bool -> Bool) -> (VHit i a -> Bool) -> VHit i a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VHit i a -> Bool
forall a i. RealFloat a => VHit i a -> Bool
isnotahit

-- | march along an integer grid using a digital differential
-- analyzer (DDA)-based algorithm
march ::
  (RealFloat a, Integral i, Monad m) =>
  -- | test for stoppage ('True' to stop)
  (V3 i -> m Bool) ->
  -- | direction, any physical dimension
  V3 a ->
  -- | initial position, dimensionless
  V3 a ->
  -- | number of iterations (maximum)
  Int ->
  -- | hit information (success or failure)
  m (VHit i a)
march :: forall a i (m :: * -> *).
(RealFloat a, Integral i, Monad m) =>
(V3 i -> m Bool) -> V3 a -> V3 a -> Int -> m (VHit i a)
march V3 i -> m Bool
test V3 a
ray V3 a
pos0 = V3 a -> V3 i -> V3 i -> a -> Int -> m (VHit i a)
go V3 a
dis0 (a -> i
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (a -> i) -> V3 a -> V3 i
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> V3 a
pos0) V3 i
0 a
0
  where
    -- prevent NaNs. for the comparison in finding 'closest'
    -- NaNs need to have no influence and therefore actually
    -- must be Infinity.
    mul1 :: m b -> m b -> m b
mul1 m b
p m b
q =
      [ if b -> Bool
forall a. RealFloat a => a -> Bool
isInfinite b
a Bool -> Bool -> Bool
|| b -> Bool
forall a. RealFloat a => a -> Bool
isInfinite b
b
          then b
1 b -> b -> b
forall a. Fractional a => a -> a -> a
/ b
0
          else b
a b -> b -> b
forall a. Num a => a -> a -> a
* b
b
      | b
a <- m b
p
      | b
b <- m b
q
      ]
    -- ////// dimensional analysis //////
    --
    -- . dimensions of "pos" and "pos0"
    --
    -- by comparing these three expressions:
    --  dis0 = rcp `mul1` ... pos0 ...
    --       = [X] * [L]
    --  rcp  = [X]
    --  dis' = dis + [1] `mul1` rcp
    --       = [X] * [L] + [X]
    -- we can deduce that
    --  [X] = [X] * [L]
    -- which has two solutions:
    --  [X] = [0] (meaningless)
    --  [L] = [1]
    -- which makes sense given that our algorithm
    -- works on 1x1x1 voxels.
    --
    -- . dimensions of "ray"
    --
    -- if [L] = [1] then we are free in our choice
    -- of [X], the dimensions of rcp. and the definition
    -- of rcp is
    --  rcp = recip ray
    -- so ray has [1]/[X] dimensions. thus ray is
    -- compatible with any physical dimensions.
    dis0 :: V3 a
dis0 =
      V3 a
rcp
        V3 a -> V3 a -> V3 a
forall {b} {m :: * -> *}.
(RealFloat b, MonadZip m) =>
m b -> m b -> m b
`mul1` V3 a -> V3 a
forall a. Num a => a -> a
abs
          ( (a -> a) -> V3 a -> V3 a
forall a b. (a -> b) -> V3 a -> V3 b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. (Integral a, Num b) => a -> b
fromIntegral @Int64 (Int64 -> a) -> (a -> Int64) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Int64
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor) V3 a
pos0
              V3 a -> V3 a -> V3 a
forall a. Num a => a -> a -> a
- V3 a
pos0
              V3 a -> V3 a -> V3 a
forall a. Num a => a -> a -> a
+ V3 a -> V3 a
forall a. Num a => a -> a
signum V3 a
ray
          )
    rcp :: V3 a
rcp = V3 a -> V3 a
forall a. Fractional a => a -> a
recip V3 a
ray
    sgn :: V3 i
sgn = a -> i
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (a -> i) -> V3 a -> V3 i
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> V3 a -> V3 a
forall a. Num a => a -> a
signum V3 a
ray
    go :: V3 a -> V3 i -> V3 i -> a -> Int -> m (VHit i a)
go V3 a
dis V3 i
pos V3 i
closest0 a
time0 Int
iter
      | Int
iter Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 = do
          t <- V3 i -> m Bool
test V3 i
pos
          if t
            then pure $ VHit time0 pos (-sgn * closest0)
            else
              -- for each of x, y, and z coordinates,
              -- decide if it's the "closest" one among
              -- the other coordinates to the respective
              -- next voxel boundary. there could be many
              -- simultaneous advancements, though, which
              -- is handled smoothly with no issue.
              let closest =
                    [ Int -> i
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> i) -> Int -> i
forall a b. (a -> b) -> a -> b
$ Bool -> Int
forall a. Enum a => a -> Int
fromEnum (a
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
b Bool -> Bool -> Bool
&& a
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
c)
                    | a
a <- V3 a
dis
                    | a
b <- V3 a
dis V3 a -> Getting (V3 a) (V3 a) (V3 a) -> V3 a
forall s a. s -> Getting a s a -> a
^. Getting (V3 a) (V3 a) (V3 a)
Lens' (V3 a) (V3 a)
forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_yzx
                    | a
c <- V3 a
dis V3 a -> Getting (V3 a) (V3 a) (V3 a) -> V3 a
forall s a. s -> Getting a s a -> a
^. Getting (V3 a) (V3 a) (V3 a)
Lens' (V3 a) (V3 a)
forall (t :: * -> *) a. R3 t => Lens' (t a) (V3 a)
_zxy
                    ]
                  fclosest = i -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (i -> a) -> V3 i -> V3 a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> V3 i
closest
               in go
                    (dis + fclosest `mul1` rcp)
                    (pos + closest * sgn)
                    closest
                    (time0 + minimum dis)
                    (iter - 1)
      | Bool
otherwise = VHit i a -> m (VHit i a)
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (VHit i a -> m (VHit i a)) -> VHit i a -> m (VHit i a)
forall a b. (a -> b) -> a -> b
$ a -> V3 i -> V3 i -> VHit i a
forall i a. a -> V3 i -> V3 i -> VHit i a
VHit (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) V3 i
0 V3 i
0

-- debug

_march2tester :: V3 Int -> IO Bool
_march2tester :: V3 Int -> IO Bool
_march2tester V3 Int
v = Bool
False Bool -> IO () -> IO Bool
forall a b. a -> IO b -> IO a
forall (f :: * -> *) a b. Functor f => a -> f b -> f a
<$ String -> IO ()
traceIO (String -> String -> String
forall r. PrintfType r => String -> r
printf String
"v = %s" (V3 Int -> String
forall a. Show a => a -> String
show V3 Int
v))

_march2tester_distance :: V3 Int -> Double -> V3 Int -> IO Bool
_march2tester_distance :: V3 Int -> Double -> V3 Int -> IO Bool
_march2tester_distance V3 Int
home Double
dst V3 Int
v =
  let h :: V3 Double
h = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> V3 Int -> V3 Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> V3 Int
home
      w :: V3 Double
w = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> V3 Int -> V3 Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> V3 Int
v
   in Bool -> IO Bool
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (V3 Double -> V3 Double -> Double
forall a. Floating a => V3 a -> V3 a -> a
forall (f :: * -> *) a. (Metric f, Floating a) => f a -> f a -> a
distance V3 Double
h V3 Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
dst)