Merge pull request #67 from quchen/various

Various

Author: quchen

src/Draw/Plotting.hs

@@ -821,11 +821,11 @@ instance Plotting Polygon where
         lineTo r
 
 -- | FluidNC doesn’t support G05, so we approximate Bezier curves with line pieces.
--- We use the naive Bezier interpolation 'bezierSubdivideT', because it just so
+-- We use the naive Bezier interpolation 'bezierSubdivideEquiparametric', because it just so
 -- happens to put more points in places with more curvature.
 instance Plotting Bezier where
     plot bezier = commented "Bezier (cubic)" $ do
-        let points = bezierSubdivideT 32 bezier
+        let points = bezierSubdivideEquiparametric 32 bezier
         let p:ointsToPlot = points
         repositionTo p
         traverse_ lineTo ointsToPlot

src/Geometry/Bezier.hs

@@ -9,8 +9,9 @@ module Geometry.Bezier
     , bezierS
 
     -- * Subdividing
-    , bezierSubdivideT
-    , bezierSubdivideS
+    , bezierSubdivideEquiparametric
+    , bezierSubdivideEquidistant
+    , bezierSubdivideCasteljau
 
     -- * Interpolation
     , bezierSmoothen
@@ -41,6 +42,26 @@ import Numerics.LinearEquationSystem
 
 
 
+-- $references
+--
+-- == Arc length parameterization
+--
+-- * Moving Along a Curve with Specified Speed (2019)
+--   by David Eberly
+--   https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
+--
+-- == Smoothening
+--
+-- * Cubic Bézier Splines
+--   by Michael Joost
+--   https://www.michael-joost.de/bezierfit.pdf
+--
+-- * Building Smooth Paths Using Bézier Curves
+--   by Stuart Kent
+--   https://www.stkent.com/2015/07/03/building-smooth-paths-using-bezier-curves.html
+
+
+
 -- | Cubic Bezier curve, defined by start, first/second control points, and end.
 data Bezier = Bezier !Vec2 !Vec2 !Vec2 !Vec2 deriving (Eq, Ord, Show)
 
@@ -144,20 +165,20 @@ bezierLength bezier = retryExponentiallyUntilPrecision (integrateSimpson13 f 0 1
 -- | Trace a 'Bezier' curve with a number of points, using the polynomial curve
 -- parameterization. This is very fast, but leads to unevenly spaced points.
 --
--- For subdivision by arc length, use 'bezierSubdivideS'.
-bezierSubdivideT
+-- For subdivision by arc length, use 'bezierSubdivideEquidistant'.
+bezierSubdivideEquiparametric
     :: Int
     -> Bezier
     -> [Vec2]
-bezierSubdivideT n bz = map (bezierT bz) points
+bezierSubdivideEquiparametric n bz = map (bezierT bz) points
   where
     points = map (\x -> fromIntegral x / fromIntegral (n-1)) [0..n-1]
 
 -- | Trace a 'Bezier' curve with a number of evenly spaced points by arc length.
--- This is much more expensive than 'bezierSubdivideT', but may be desirable for
+-- This is much more expensive than 'bezierSubdivideEquiparametric', but may be desirable for
 -- aesthetic purposes.
 --
--- Here it is alongside 'bezierSubdivideT':
+-- Here it is alongside 'bezierSubdivideEquiparametric':
 --
 -- <<docs/haddock/Geometry/Bezier/subdivide_s_t_comparison.svg>>
 --
@@ -167,8 +188,8 @@ bezierSubdivideT n bz = map (bezierT bz) points
 --     let curve = let curveRaw = transform (rotate (deg (-30))) (Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3))
 --                     fitToBox = transform (transformBoundingBox curveRaw (Vec2 10 10, Vec2 290 90) (TransformBBSettings FitWidthHeight IgnoreAspect FitAlignCenter))
 --                 in fitToBox curveRaw
---         evenlySpaced = bezierSubdivideS 16 curve
---         unevenlySpaced = bezierSubdivideT 16 curve
+--         evenlySpaced = bezierSubdivideEquidistant 16 curve
+--         unevenlySpaced = bezierSubdivideEquiparametric 16 curve
 --         offsetBelow :: Transform geo => geo -> geo
 --         offsetBelow = transform (translate (Vec2 0 50))
 --     cairoScope $ do
@@ -189,8 +210,8 @@ bezierSubdivideT n bz = map (bezierT bz) points
 --         cairoScope (setColor (black `withOpacity` 0.2) >> connect e u)
 -- :}
 -- Generated file: size 17KB, crc32: 0x7c147951
-bezierSubdivideS :: Int -> Bezier -> [Vec2]
-bezierSubdivideS n bz = map bezier distances
+bezierSubdivideEquidistant :: Int -> Bezier -> [Vec2]
+bezierSubdivideEquidistant n bz = map bezier distances
   where
 
     -- The step width should correlate with the length of the curve to get a decent
@@ -246,23 +267,65 @@ s_to_t_lut_ode bz ds = LookupTable1 (sol_to_vec sol)
     t0 = 0
     s0 = 0
 
--- $references
+-- | Approximage a Bezier curve with line segments up to a certain precision, using
+-- relatively few points.
 --
--- == Arc length parameterization
---
--- * Moving Along a Curve with Specified Speed (2019)
---   by David Eberly
---   https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
+-- The idea behind Casteljau subdivision is that each Bézier curve can be exactly
+-- subdivided into two Bézier curves (of same degree). This is done recursively, in
+-- this implementation (and commonly) in the middle of the curve. Once a curve
+-- segment is flat enough (given by the tolerance parameter), it is simply rendered
+-- as a line.
 --
--- == Smoothening
+-- <<docs/haddock/Geometry/Bezier/bezierSubdivideCasteljau.svg>>
 --
--- * Cubic Bézier Splines
---   by Michael Joost
---   https://www.michael-joost.de/bezierfit.pdf
---
--- * Building Smooth Paths Using Bézier Curves
---   by Stuart Kent
---   https://www.stkent.com/2015/07/03/building-smooth-paths-using-bezier-curves.html
+-- === __(image code)__
+-- >>> :{
+-- haddockRender "Geometry/Bezier/bezierSubdivideCasteljau.svg" 500 330 $ do
+--     let curve = let curveRaw = transform (rotate (deg (-30))) (Bezier (Vec2 0 0) (Vec2 1 5) (Vec2 2.5 (-1)) (Vec2 3 3))
+--                     fitToBox = transform (transformBoundingBox curveRaw (shrinkBoundingBox 10 [zero, Vec2 500 200]) (TransformBBSettings FitWidthHeight IgnoreAspect FitAlignCenter))
+--                 in fitToBox curveRaw
+--         paintOffset = Vec2 0 30
+--     for_ (zip [0..] [50,25,10,2]) $ \(i, tolerance) -> cairoScope $ do
+--         let points = bezierSubdivideCasteljau tolerance (transform (translate (fromIntegral i *. paintOffset)) curve)
+--         setColor (mathematica97 i)
+--         C.setLineWidth 2
+--         sketch (Polyline points) >> C.stroke
+--         for_ points $ \p -> sketch (Circle p 3) >> C.fill
+--     cairoScope $ do
+--         C.setLineWidth 3
+--         setColor black
+--         sketch (transform (translate (4*.paintOffset)) curve)
+--         C.stroke
+-- :}
+-- Generated file: size 20KB, crc32: 0x679b311c
+bezierSubdivideCasteljau :: Double -> Bezier -> [Vec2]
+bezierSubdivideCasteljau tolerance curve@(Bezier pFirst _ _ _) = pFirst : go curve
+  where
+    go (Bezier p1 p2@(Vec2 x2 y2) p3@(Vec2 x3 y3) p4@(Vec2 x4 y4)) =
+        let
+            p12   = (p1   +. p2  ) /. 2
+            p23   = (p2   +. p3  ) /. 2
+            p34   = (p3   +. p4  ) /. 2
+            p123  = (p12  +. p23 ) /. 2
+            p234  = (p23  +. p34 ) /. 2
+            p1234 = (p123 +. p234) /. 2
+
+            dp@(Vec2 dx dy) = p4 -. p1
+
+            -- d2, d3 are the distance from p2, p3 from the line
+            -- connecting p1 and p4. A curve is flat when those
+            -- two are short together.
+            d2 = abs ((x2-x4)*dy - (y2-y4)*dx)
+            d3 = abs ((x3-x4)*dy - (y3-y4)*dx)
+            curveIsFlat = (d2 + d3)*(d2 + d3) < tolerance^2 * normSquare dp
+        in if curveIsFlat
+            then -- We return only the last point so we don’t get duplicate
+                 -- points for each start+end of adjacent curves.
+                 -- The very first point is forgotten by the
+                [p4]
+            else go (Bezier p1 p12 p123 p1234)
+                 ++
+                 go (Bezier p1234 p234 p34 p4)
 
 -- | Smoothen a number of points by putting a Bezier curve between each pair.
 -- Useful to e.g. make a sketch nicer, or interpolate between points of a crude

src/Geometry/Core.hs

@@ -1204,10 +1204,11 @@ resizeLineSymmetric f line@(Line start end) = (centerLine . resizeLine f . trans
 --
 -- Useful for painting lines going through a point symmetrically.
 centerLine :: Line -> Line
-centerLine line@(Line start end) = transform (translate delta) line
-  where
-    middle = 0.5 *. (start +. end)
-    delta = start -. middle
+centerLine (Line start end) =
+    let middle = 0.5 *. (start +. end)
+        end' = middle
+        start' = start -. end +. middle
+    in Line start' end'
 
 -- | Move the end point of the line so that it has length 1.
 normalizeLine :: Line -> Line
@@ -1291,6 +1292,8 @@ intersectionPoint _                              = Nothing
 --
 -- Returns the point of the intersection, and whether it is inside both, one, or
 -- none of the provided finite line segments.
+--
+-- See 'intersectInfiniteLines' for a more performant, but less nuanced result.
 intersectionLL :: Line -> Line -> LLIntersection
 intersectionLL lineL lineR
     = intersectionType
@@ -2004,51 +2007,76 @@ isConvex (Polygon ps) = allSameSign angleDotProducts
 --
 -- === __(image code)__
 -- >>> :{
--- haddockRender "Geometry/Core/perpendicular_bisector.svg" 150 150 $ do
---     let line = Line (Vec2 10 60) (Vec2 140 100)
+-- haddockRender "Geometry/Core/perpendicular_bisector.svg" 200 160 $ do
+--     let line = Line (Vec2 20 20) (Vec2 190 90)
 --         bisector = perpendicularBisector line
 --     C.setLineWidth 2
 --     sketch line >> C.stroke
 --     sketch bisector >> setColor (mathematica97 1) >> C.stroke
 -- :}
--- Generated file: size 2KB, crc32: 0x1f7d2821
+-- Generated file: size 2KB, crc32: 0x9940c32d
 perpendicularBisector :: Line -> Line
-perpendicularBisector line@(Line start end) = perpendicularLineThrough middle line
+perpendicularBisector (Line start end) =
+    let middle = 0.5 *. (end +. start)
+    in rotateLine90 (Line middle end)
+
+rotateLine90 :: Line -> Line
+rotateLine90 (Line start end) = Line start end'
   where
-    middle = 0.5 *. (start +. end)
+    end' = rotate90 (end -. start) +. start
+
+rotate90 :: Vec2 -> Vec2
+rotate90 (Vec2 x y) = Vec2 (-y) x
 
--- | Line perpendicular to a given line through a point.
+-- | Line perpendicular to a given line through a point, starting at the
+-- intersection point.
+--
+-- If the point is on the line directly, fall back to a perpendicular line through
+-- the point of half the length of the input.
 --
--- The result has the same length as the input, point in its center, and points
--- to the left (90° turned CCW) relative to the input.
+-- This is also known as the vector projection. For vectors \(\mathbf a\) (pointing
+-- from the start of the line to \(\mathbf p\)) and \(\mathbf b\) (pointing from
+-- the start of the line to its end),
+--
+-- \[
+-- \mathrm{proj}_{\mathbf b}(\mathbf a)
+-- = \frac{\mathbf a\cdot\mathbf b}{\mathbf b\cdot\mathbf b}\mathbf b
+-- \]
 --
 -- <<docs/haddock/Geometry/Core/perpendicular_line_through.svg>>
 --
 -- === __(image code)__
 -- >>> :{
--- haddockRender "Geometry/Core/perpendicular_line_through.svg" 150 140 $ do
---     let line = Line (Vec2 20 20) (Vec2 140 70)
---         point = Vec2 70 70
---         perpendicular = perpendicularLineThrough point line
+-- haddockRender "Geometry/Core/perpendicular_line_through.svg" 170 170 $ do
+--     let line = transform (translate (Vec2 20 20))
+--                          (Line zero (Vec2 (3*40) (4*20)))
+--         points =
+--             [ Vec2 20 110 -- above
+--             , Vec2 70 90 -- above
+--             , Vec2 130 20 -- below
+--             , Vec2 110 80 -- directly on
+--             , Vec2 130 150 -- beyond
+--             ]
 --     C.setLineWidth 2
 --     sketch line >> C.stroke
---     setColor (mathematica97 1)
---     sketch (Circle point 5) >> C.fill
---     sketch perpendicular >> C.stroke
+--     for_ (zip [1..] points) $ \(i, p) -> do
+--         setColor (mathematica97 i)
+--         cairoScope $ sketch (perpendicularLineThrough p line) >> C.setDash [3,5] 0 >> C.stroke
+--         cairoScope $ sketch (Circle p 5) >> C.fill
 -- :}
--- Generated file: size 2KB, crc32: 0xb8e5d1d9
+-- Generated file: size 4KB, crc32: 0xac73c163
 perpendicularLineThrough :: Vec2 -> Line -> Line
-perpendicularLineThrough p line@(Line start _) = centerLine line'
-  where
-    -- Move line so it starts at the origin
-    Line start0 end0 = transform (translate (negateV start)) line
-    -- Rotate end point 90° CCW
-    end0' = let Vec2 x y  = end0
-            in Vec2 (-y) x
-    -- Construct rotated line
-    lineAt0' = Line start0 end0'
-    -- Move line back so it goes through the point
-    line' = transform (translate p) lineAt0'
+perpendicularLineThrough p (Line start end) =
+    let a = p -. start
+        b = end -. start
+        proj = (dotProduct a b/dotProduct b b) *. b +. start
+    in if normSquare (p -. proj) <= 0.01^2
+        then -- p is too close to proj, so a 'Line' does not make
+             -- any sense. Fall back to a line of half the input
+             -- length perpendicular to the original one.
+            let p' = proj +. rotate90 ((end -. start) /. 2)
+            in Line proj p'
+        else Line proj p
 
 -- | Optical reflection of a ray on a mirror. Note that the outgoing line has
 -- reversed direction like light rays would. The second result element is the

test/testsuite/Test/Uncategorized/Bezier.hs

@@ -167,7 +167,7 @@ subdivideBezierCurveTest = testVisual "Subdivide" 300 300 "docs/interpolation/4_
         moveTo 200 70
         showText (show (length beziers) ++ " curves")
 
-    let subpoints = beziers >>= (V.fromList . bezierSubdivideT 10)
+    let subpoints = beziers >>= (V.fromList . bezierSubdivideEquiparametric 10)
     let simplified = simplifyTrajectoryRdp 0.05 subpoints
     cairoScope $ do
         let fit :: Transform geo => geo -> geo