ray tracing and finding the normal vector to the surface at the intersection point

Question

When doing a ray trace with rayTraceP, I can find the point where a ray intersects with a diagram.

> rayTraceP (p2 (0, 0)) (r2 (1, 0)) ((p2 (1,-1) ~~ p2 (1,1))
Just (p2 (1.0, 0.0))

I want to use this to find not only the "collision point", but also the collision time and the normal vector to the surface at that point.

-- A Collision has a time, a contact point, and a normal vector.
-- The normal vector is perpendicular to the surface at the contact
-- point.
data Collision v n = Collision n (Point v n) (v n)
  deriving (Show)

Given a start point for the ray and a velocity vector along the ray, I can find the contact point end using rayTraceP:

end <- rayTraceP start vel dia

And I can find the collision time using the distance between start and end:

time = distance start end / norm vel

But I'm stuck on finding the normal vector. I'm working within this function:

rayTraceC :: (Metric v, OrderedField n)
             => Point v n -> v n -> QDiagram B v n Any -> Maybe (Collision v n)
-- Takes a starting position for the ray, a velocity vector for the
-- ray, and a diagram to trace the ray to. If the ray intersects with
-- the diagram, it returns a Collision containing:
--  * The amount of time it takes for a point along the ray going at
--    the given velocity to intersect with the diagram.
--  * The point at which it intersects with the diagram.
--  * The normal vector to the surface at that point (which will be
--    perpendicular to the surface there).
-- If the ray does not intersect with the diagram, it returns Nothing.
rayTraceC start vel dia =
  do
    end <- rayTraceP start vel dia
    let time = distance start end / norm vel
    -- This is where I'm getting stuck. 
    -- How do I find the normal vector?
    let normalV = ???
    return (Collision time end normalV)

Some examples of what I want it to do:

> -- colliding straight on:
> rayTraceC (p2 (0, 0)) (r2 (1, 0)) (p2 (1,-1) ~~ p2 (1,1))
Just (Collision 1 (p2 (1, 0)) (r2 (-1, 0)))
> -- colliding from a diagonal:
> rayTraceC (p2 (0, 0)) (r2 (1, 1)) (p2 (1,0) ~~ p2 (1,2))
Just (Collision 1 (p2 (1, 1)) (r2 (-1, 0))
> -- colliding onto a diagonal:
> rayTraceC (p2 (0, 0)) (r2 (1, 0)) (p2 (0,-1) ~~ p2 (2,1))
Just (Collision 1 (p2 (1, 0)) (r2 (-√2/2, √2/2)))
> -- no collision
> rayTraceC (p2 (0, 0)) (r2 (1, 0)) (p2 (1,1) ~~ p2 (1,2))
Nothing

It is correct on everything in these examples except for the normal vector.

I have looked in the documentation for both Diagrams.Trace and Diagrams.Core.Trace, but maybe I'm looking in the wrong places.

Answer 1

There is no way to do this in general; it depends on what exactly you hit. There is a module Diagrams.Tangent for computing tangents of trails, but to compute the tangent at a given point you have to know its parameter with respect to the trail; and one thing we are missing at the moment is a way to convert from a given point to the parameter of the closest point on a given segment/trail/path (it's been on the to-do list for a while).

Dreaming even bigger, perhaps traces themselves ought to return something more informative---not just parameters telling you how far along the ray the hit are, but also information about what you hit (from which one could more easily do things like compute a normal vector).

What kinds of things are you computing traces of? There might be a way to take advantage of the particular details of your use case to get the normals you want in a not-too-terrible way.

Answer 2

Brent Yorgey's answer points out the Diagrams.Tangent module, and in particular normalAtParam, which works on Parameteric functions, including trails, but not all Diagrams.

Fortunately, many 2D diagram functions, like circle, square, rect, ~~, etc. can actually return any TrailLike type, including Trail V2 n. So a function with the type

rayTraceTrailC :: forall n . (RealFloat n, Epsilon n)
                  => 
                  Point V2 n 
                  -> V2 n
                  -> Located (Trail V2 n)
                  -> Maybe (Collision V2 n)

Would actually work on the values returned by circle, square, rect, ~~, etc. if it could be defined:

> rayTraceTrailC 
    (p2 (0, 0))
    (r2 (1, 0))
    (circle 1 # moveTo (p2 (2,0)))
Just (Collision 1 (p2 (1, 0)) (r2 (-1, 0)))

And this function can be defined by breaking the trail up into a list of fixed segments which are either linear or bezier curves, using the fixTrail function. That reduces the problem to the simpler rayTraceFixedSegmentC.

rayTraceTrailC start vel trail =
  combine (mapMaybe (rayTraceFixedSegmentC start vel) (fixTrail trail))
  where
    combine [] = Nothing
    combine cs = Just (minimumBy (\(Collision a _ _) (Collision b _ _) -> compare a b) cs)

The rayTraceFixedSegmentC can use rayTraceP to calculate the contact point, but we can't find the normal vector right away because we don't know what the parameter is at that contact point. So punt further and add fixedSegmentNormalV helper function to the wish list:

rayTraceFixedSegmentC :: forall n . (RealFloat n, Epsilon n)
                         =>
                         Point V2 n
                         -> V2 n
                         -> FixedSegment V2 n
                         -> Maybe (Collision V2 n)
rayTraceFixedSegmentC start vel seg =
  do
    end <- rayTraceP start vel (unfixTrail [seg])
    let time = distance start end / norm vel
    let normalV = normalize (project (fixedSegmentNormalV seg end) (negated vel))
    return (Collision time end normalV)

This fixedSegmentNormalV function just has to return a normal vector for a single segment going through a single point, without worrying about the vel direction. It can destruct the FixedSegment type, and if it's linear, that's easy:

fixedSegmentNormalV :: forall n . (OrderedField n)
                       => 
                       FixedSegment V2 n -> Point V2 n -> V2 n
fixedSegmentNormalV seg pt =
  case seg of
    FLinear a b -> perp (b .-. a)
    FCubic a b c d ->
      ???

In the FCubic case, to calculate the parameter where the curve goes through pt, I'm not sure what to do, but if you don't mind approximations here we can just take a bunch of points along it and find the one closest to pt. After that we can call normalAtParam as Brent Yorgey suggested.

fixedSegmentNormalV seg pt =
  case seg of
    FLinear a b -> perp (b .-. a)
    FCubic a b c d ->
      -- APPROXIMATION: find the closest parameter value t
      let ts = map ((/100) . fromIntegral) [0..100]
          dist t = distance (seg `atParam` t) pt
          t = minimumBy (\a b -> compare (dist a) (dist b)) ts
      -- once we have that parameter value we can call a built-in function
      in normalAtParam seg t

With this, the rayTraceTrailC function is working with this approximation. However, it doesn't work for Diagrams, only Located Trails.

It can work on the values returned by functions like circle and rect, but not on combined diagrams. So you have to keep those building blocks of diagrams separate, as trails, for as long as you need this collision ray tracing.

Using the normal vectors to reflect the rays (the outgoing ray has an equal angle from the normal vector) looks like this: