diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,10 @@ # Revision history for juicy-gcode +## 0.2.1.0 -- 2022-11-26 + +- The approximation error is now calculated along the radial direction +- The approximation error calculation is now exact instead of sampling-based + ## 0.2.0.2 -- 2022-10-31 - Fix a problem triggered by non-quadratic inflexion point equations diff --git a/juicy-gcode.cabal b/juicy-gcode.cabal --- a/juicy-gcode.cabal +++ b/juicy-gcode.cabal @@ -1,5 +1,5 @@ name: juicy-gcode -version: 0.2.0.2 +version: 0.2.1.0 license: BSD3 license-file: LICENSE author: dlacko @@ -19,7 +19,22 @@ hs-source-dirs: src main-is: Main.hs - other-modules: Approx BiArc CircularArc CubicBezier GCode Line Render SvgArcSegment Transformation Types SVGExt Paths_juicy_gcode + other-modules: Approx.BiArc + Graphics.BiArc + Graphics.CircularArc + Graphics.CubicBezier + Graphics.Curve + Graphics.Line + Graphics.LineSegment + Graphics.Path + Graphics.Point + Graphics.Transformation + GCode + Render + Utils + SvgArcSegment + SVGExt + Paths_juicy_gcode build-depends: base >=4.8 && <5, diff --git a/src/Approx.hs b/src/Approx.hs deleted file mode 100644 --- a/src/Approx.hs +++ /dev/null @@ -1,115 +0,0 @@ -module Approx ( bezier2biarc - ) where - -import qualified CubicBezier as B -import qualified BiArc as BA -import qualified Line as L - -import Data.Bool (bool) -import Linear - -import Types - --- Approximate a bezier curve with biarcs (Left) and line segments (Right) -bezier2biarc :: B.CubicBezier - -> Double - -> [Either BA.BiArc (V2 Double)] -bezier2biarc mbezier resolution - -- Edge case: all points on the same line -> it is a line - | (L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c1 mbezier)) && - (L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c2 mbezier)) - = [Right (B._p2 mbezier)] - -- Edge case: p1 == c1, don't split - | (B._p1 mbezier) == (B._c1 mbezier) - = approxOne mbezier - -- Edge case: p2 == c2, don't split - | (B._p2 mbezier) == (B._c2 mbezier) - = approxOne mbezier - -- Split by the inflexion points (if any) - | otherwise - = byInflection (B.inflectionPoints mbezier) - where - order a b | b < a = (b, a) - | otherwise = (a, b) - - byInflection [t] = approxOne b1 ++ approxOne b2 - where - (b1, b2) = B.bezierSplitAt mbezier t - - byInflection [t1, t2] = approxOne b1 ++ approxOne b2 ++ approxOne b3 - where - (it1, it2') = order t1 t2 - - -- Make the first split and save the first new curve. The second one has to be splitted again - -- at the recalculated t2 (it is on a new curve) - it2 = (1 - it1) * it2' - - (b1, toSplit) = B.bezierSplitAt mbezier it1 - (b2, b3) = B.bezierSplitAt toSplit it2 - - byInflection _ = approxOne mbezier - - -- TODO: make it tail recursive - approxOne :: B.CubicBezier -> [Either BA.BiArc (V2 Double)] - approxOne bezier - -- Approximate bezier length. if smaller than resolution, do not approximate - | (distance (B._p1 bezier) (B._c1 bezier)) + - (distance (B._c1 bezier) (B._c2 bezier)) + - (distance (B._c2 bezier) (B._p2 bezier)) < resolution - = [Right (B._p2 bezier)] - -- Edge case: start- and endpoints are the same - | (B._p1 bezier) == (B._p2 bezier) - = splitAndRecur 0.5 - -- Edge case: control lines are parallel - | (L._m t1) == (L._m t2) || (isNaN (L._m t1) && isNaN (L._m t2)) - = splitAndRecur 0.5 - -- Approximation is not close enough yet, refine - | BA.isStable biarc && maxDistance > resolution - = splitAndRecur maxDistanceAt - -- Desired case: approximation is stable and close enough - | BA.isStable biarc - = [Left biarc] - -- Unstable approximation: split the bezier into half, basically switching to - -- linear approximation mode - | otherwise - = splitAndRecur 0.5 - - where - -- Edge case: P1==C1 or P2==C2 - -- there is no derivative at P1 or P2, use the other control point - c1 = bool (B._c1 bezier) (B._c2 bezier) ((B._p1 bezier) == (B._c1 bezier)) - c2 = bool (B._c2 bezier) (B._c1 bezier) ((B._p2 bezier) == (B._c2 bezier)) - - -- V: Intersection point of tangent lines - t1 = L.fromPoints (B._p1 bezier) c1 - t2 = L.fromPoints (B._p2 bezier) c2 - v = L.intersection t1 t2 - - -- G: incenter point of the triangle (P1, V, P2) - dP2V = distance (B._p2 bezier) v - dP1V = distance (B._p1 bezier) v - dP1P2 = distance (B._p1 bezier) (B._p2 bezier) - g = (dP2V *^ B._p1 bezier + dP1V *^ B._p2 bezier + dP1P2 *^ v) ^/ (dP2V + dP1V + dP1P2) - - -- Calculate the BiArc - biarc = BA.create (B._p1 bezier) (B._p1 bezier - c1) (B._p2 bezier) (B._p2 bezier - c2) g - - -- Calculate the error - -- TODO: we only calculate the distance at 8 points (first and last skipped as - -- they should be precise), seems a resonable approximation as for now - parameterStep = 1 / 10 - - (maxDistance, maxDistanceAt) = maxDistance' 0 0 parameterStep - - maxDistance' m mt t - | t < 1 - = if' (d > m) (maxDistance' d t nt) (maxDistance' m mt nt) - | otherwise - = (m, mt) - where - d = distance (BA.pointAt biarc t) (B.pointAt bezier t) - nt = t + parameterStep - - splitAndRecur t = let (b1, b2) = B.bezierSplitAt bezier t - in approxOne b1 ++ approxOne b2 - diff --git a/src/Approx/BiArc.hs b/src/Approx/BiArc.hs new file mode 100644 --- /dev/null +++ b/src/Approx/BiArc.hs @@ -0,0 +1,193 @@ +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Redundant bracket" #-} +module Approx.BiArc ( + bezier2biarcs +) where + +import qualified Graphics.CubicBezier as B +import qualified Graphics.BiArc as BA +import qualified Graphics.CircularArc as CA +import qualified Graphics.Line as L +import Graphics.Curve +import Graphics.Path +import Graphics.Point +import Utils + +import Data.Bool (bool) +import Control.Lens +import Linear + +eps :: Double +eps = 0.0001 +maxiter :: Double +maxiter = 10 + +-- Approximate a bezier curve with biarcs (Left) and line segments (Right) +bezier2biarcs :: B.CubicBezier + -> Double + -> [PathCommand] +bezier2biarcs mbezier resolution + -- Degenerate curve: all points on the same line -> it is a line + | L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c1 mbezier) && + L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c2 mbezier) + = [LineTo (toPoint (B._p2 mbezier))] + -- Degenerate curve: p1 == c1, don't split + | B._p1 mbezier == B._c1 mbezier + = approxOne mbezier + -- Degenerate curve: p2 == c2, don't split + | B._p2 mbezier == B._c2 mbezier + = approxOne mbezier + -- Split by the inflexion points (if any) + | otherwise + = byInflection (B.inflectionPoints mbezier) + where + order a b | b < a = (b, a) + | otherwise = (a, b) + + byInflection [t] = approxOne b1 ++ approxOne b2 + where + (b1, b2) = B.splitAt mbezier t + + byInflection [t1, t2] = approxOne b1 ++ approxOne b2 ++ approxOne b3 + where + (it1, it2') = order t1 t2 + + -- Make the first split and save the first new curve. The second one has to be splitted again + -- at the recalculated t2 (it is on a new curve) + it2 = (1 - it1) * it2' + + (b1, toSplit) = B.splitAt mbezier it1 + (b2, b3) = B.splitAt toSplit it2 + + byInflection _ = approxOne mbezier + + -- Recursive step (TODO: tail recursive) + approxOne :: B.CubicBezier -> [PathCommand] + approxOne bezier + -- Approximate bezier length. if max length is smaller than resolution, do not approximate + | B.maxArcLength bezier < resolution + = [LineTo (toPoint (B._p2 bezier))] + -- Edge case: start- and endpoints are the same + | B._p1 bezier == B._p2 bezier + = splitAndRecur 0.5 + -- Edge case: control lines are parallel + | L._m t1 == L._m t2 || isNaN (L._m t1) && isNaN (L._m t2) + = splitAndRecur 0.5 + -- Biarc triangle has the wrong orientation + -- Curve looks like this: https://pomax.github.io/bezierinfo/images/chapters/decasteljau/df92f529841f39decf9ad62b0967855a.png + | B.isClockwise bezier /= isClockwise3 (B._p1 bezier) (B._p2 bezier) v + = splitAndRecur 0.5 + -- Unstable approximation: split the bezier into half, it will switch to linear approximation if the segments get too small + | not (isStable biarc) + = splitAndRecur 0.5 + -- Approximation is not close enough yet, refine + | maxDistance > resolution + = splitAndRecur maxDistanceAt + -- Desired case: approximation is stable and close enough + | otherwise + = biarc2path biarc + + where + -- Edge case: P1==C1 or P2==C2 + -- there is no derivative at P1 or P2, use the other control point + c1 = bool (B._c1 bezier) (B._c2 bezier) (B._p1 bezier == B._c1 bezier) + c2 = bool (B._c2 bezier) (B._c1 bezier) (B._p2 bezier == B._c2 bezier) + + -- V: Intersection point of tangent lines + t1 = L.fromPoints (B._p1 bezier) c1 + t2 = L.fromPoints (B._p2 bezier) c2 + v = L.intersection t1 t2 + + -- G: incenter point of the triangle (P1, V, P2) + dP2V = distance (B._p2 bezier) v + dP1V = distance (B._p1 bezier) v + dP1P2 = distance (B._p1 bezier) (B._p2 bezier) + g = (dP2V *^ B._p1 bezier + dP1V *^ B._p2 bezier + dP1P2 *^ v) ^/ (dP2V + dP1V + dP1P2) + + -- Calculate the BiArc + biarc = BA.fromPoints (B._p1 bezier) (B._p1 bezier - c1) (B._p2 bezier) (B._p2 bezier - c2) g + + (maxDistanceAt, maxDistance) = calculateMaxDistance bezier biarc + + splitAndRecur t = let (b1, b2) = B.splitAt bezier t + in approxOne b1 ++ approxOne b2 + +biarc2path :: BA.BiArc -> [PathCommand] +biarc2path biarc = map + (\arc -> ArcTo (toPoint (CA._c arc)) (toPoint (CA._p2 arc)) (CA.isClockwise arc)) + [BA._a1 biarc, BA._a2 biarc] + +-- Heuristics for unstable biarc: the radius of at least one of the arcs +-- is too big or too small. Not too scientific... +isStable :: BA.BiArc -> Bool +isStable biarc + = not (CA._r (BA._a1 biarc) > 99999 || CA._r (BA._a1 biarc) < 0.001 || + CA._r (BA._a2 biarc) > 99999 || CA._r (BA._a2 biarc) < 0.001) + +-- Calculate the maximum approximation error along the radial direction +-- D.J. Walton*, D.S. Meek, Approximation of a planar cubic Bezier spiral by circular arcs (1996) +calculateMaxDistance :: B.CubicBezier -> BA.BiArc -> (Double, Double) +calculateMaxDistance bezier biarc + -- This should not happenm but if, split the bezier at the middle + | tj == -1 = (0.5, 0x7FEFFFFFFFFFFFFF) + | otherwise = bigger (bigger (tj, dj) (t0, d0)) (t1, d1) + where + tj = findRadialIntersection bezier biarc (BA.jointAt biarc) + dj = distance (pointAt bezier tj) (pointAt biarc (BA.jointAt biarc)) + + g arc u = dot (pointAt bezier u - CA._c arc) (B.firstDerivativeAt bezier u) + g' arc u = quadrance (B.firstDerivativeAt bezier u) + + dot (pointAt bezier u - CA._c arc) (B.secondDerivativeAt bezier u) + + bigger f@(_, df) s@(ts, ds) + | ts == -1 = f + | df > ds = f + | otherwise = s + + -- Valid in (0,tj] + t0 = findRoot (g (BA._a1 biarc)) (g' (BA._a1 biarc)) eps tj + d0 = abs ((distance (pointAt bezier t0) (CA._c (BA._a1 biarc))) - (CA._r (BA._a1 biarc))) + -- Valid in [tj,1) + t1 = findRoot (g (BA._a2 biarc)) (g' (BA._a2 biarc)) tj (1 - eps) + d1 = abs ((distance (pointAt bezier t1) (CA._c (BA._a2 biarc))) - (CA._r (BA._a2 biarc))) + +-- Takes a paramater `t` fore the `biarc` and calculates the related parameter fo +-- the `bezier` (which is the intersection point in the radial direction) +findRadialIntersection :: B.CubicBezier -> BA.BiArc -> Double -> Double +findRadialIntersection bezier biarc t + | t == 0 || t == 1 = t + | otherwise = findRoot (\u -> dot (pointAt bezier u - p) h) (\u -> dot (B.firstDerivativeAt bezier u) h) 0 1 + where + p = pointAt biarc t + c = CA._c $ if' (t <= BA.jointAt biarc) (BA._a1 biarc) (BA._a2 biarc) + m = p - c + h = normalize $ V2 (negate (m ^. _y)) (m ^. _x) + +-- Tries to find the root of f in interval [lowerBound,upperBound] using a combination of +-- Newton and bisection methods. +-- It is supposed to have at most one solution. If no solution is found, returns -1 +findRoot :: (Double -> Double) -> (Double -> Double) -> Double -> Double -> Double +findRoot f df lowerBound upperBound + | fl * fu > 0 = -1 + | fl == 0 = lowerBound + | fu == 0 = upperBound + | otherwise = iter maxiter fl fu lowerBound upperBound ((lowerBound + upperBound) / 2) + where + fl = f lowerBound + fu = f upperBound + + iter i fmin fmax lb ub root + -- we're good, or if i==0, we may not reached tolarence yet, but hopefully it is close enough + | abs fx < eps || i <= 0 = root + -- overshoot or undershoot -> switch to bisection + | n < lb || n > ub + = if' (fmin * fx < 0) + (iter (i-1) fmin fx lb root ((lb + root) / 2)) + (iter (i-1) fx fmax root ub ((root + ub) / 2)) + -- Newton step + | otherwise + = iter (i-1) fmin fmax lb ub n + where + fx = f root + h = fx / df root + n = root - h diff --git a/src/BiArc.hs b/src/BiArc.hs deleted file mode 100644 --- a/src/BiArc.hs +++ /dev/null @@ -1,90 +0,0 @@ -module BiArc ( BiArc (..) - , create - , pointAt - , arcLength - , isStable - ) where - -import qualified CircularArc as CA -import qualified Line as L - -import Linear hiding (angle) -import Control.Lens - -data BiArc = BiArc { _a1 :: CA.CircularArc - , _a2 :: CA.CircularArc - } deriving Show - -create :: V2 Double -- Start point - -> V2 Double -- Tangent vector at start point - -> V2 Double -- End point - -> V2 Double -- Tangent vector at end point - -> V2 Double -- Transition point (connection point of the arcs) - -> BiArc -create p1 t1 p2 t2 t - = BiArc (CA.CircularArc c1 r1 startAngle1 sweepAngle1 p1 t) (CA.CircularArc c2 r2 startAngle2 sweepAngle2 t p2) - where - -- Calculate the orientation - osum = (t ^. _x - p1 ^. _x) * (t ^. _y + p1 ^. _y) - + (p2 ^. _x - t ^. _x) * (p2 ^. _y + t ^. _y) - + (p1 ^. _x - p2 ^. _x) * (p1 ^. _y + p2 ^. _y) - cw = osum < 0 - - -- Calculate perpendicular lines to the tangent at P1 and P2 - tl1 = L.createPerpendicularAt p1 (p1 + t1) - tl2 = L.createPerpendicularAt p2 (p2 + t2) - - -- Calculate the perpendicular bisector of P1T and P2T - p1t2 = (p1 + t) ^/ 2 - pb_p1t = L.createPerpendicularAt p1t2 t - - p2t2 = (p2 + t) ^/ 2 - pb_p2t = L.createPerpendicularAt p2t2 t - - -- The origo of the circles are at the intersection points - c1 = L.intersection tl1 pb_p1t - c2 = L.intersection tl2 pb_p2t - - -- Calculate the radii - r1 = distance c1 p1 - r2 = distance c2 p2 - - -- Calculate start and sweep angles - startVector1 = p1 - c1; - endVector1 = t - c1; - startAngle1 = atan2 (startVector1 ^. _y) (startVector1 ^. _x) - sweepAngle1' = (atan2 (endVector1 ^. _y) (endVector1 ^. _x)) - startAngle1 - - startVector2 = t - c2 - endVector2 = p2 - c2 - startAngle2 = atan2 (startVector2 ^. _y) (startVector2 ^. _x) - sweepAngle2' = (atan2 (endVector2 ^. _y) (endVector2 ^. _x)) - startAngle2 - - -- Adjust angles according to the orientation of the curve - sweepAngle1 = adjustSweepAngle cw sweepAngle1' - sweepAngle2 = adjustSweepAngle cw sweepAngle2' - -adjustSweepAngle :: Bool -> Double -> Double -adjustSweepAngle True angle | angle < 0 = 2 * pi + angle -adjustSweepAngle False angle | angle > 0 = angle - 2 * pi -adjustSweepAngle _ angle = angle - -pointAt :: BiArc -> Double -> V2 Double -pointAt arc t - | t <= s - = CA.pointAt (_a1 arc) (t / s) - | otherwise - = CA.pointAt (_a2 arc) ((t - s) / (1 - s)) - where - s = CA.arcLength (_a1 arc) / (arcLength arc) - -arcLength :: BiArc -> Double -arcLength arc = CA.arcLength (_a1 arc) + CA.arcLength (_a2 arc) - --- Heuristics for unstable biarc: the radius of at least one of the arcs --- is too big or too small -isStable :: BiArc -> Bool -isStable biarc - = not (CA._r (_a1 biarc) > 99999 || CA._r (_a1 biarc) < 0.001 || - CA._r (_a2 biarc) > 99999 || CA._r (_a2 biarc) < 0.001) - diff --git a/src/CircularArc.hs b/src/CircularArc.hs deleted file mode 100644 --- a/src/CircularArc.hs +++ /dev/null @@ -1,29 +0,0 @@ -module CircularArc ( CircularArc (..) - , isClockwise - , pointAt - , arcLength - ) where - -import Linear -import Control.Lens - -data CircularArc = CircularArc { _c :: V2 Double - , _r :: Double - , _startAngle :: Double - , _sweepAngle :: Double - , _p1 :: V2 Double - , _p2 :: V2 Double - } deriving Show - -isClockwise :: CircularArc -> Bool -isClockwise arc = _sweepAngle arc > 0 - -pointAt :: CircularArc -> Double -> V2 Double -pointAt arc t = V2 x y - where - x = _c arc ^. _x + _r arc * cos (_startAngle arc + t * _sweepAngle arc) - y = _c arc ^. _y + _r arc * sin (_startAngle arc + t * _sweepAngle arc) - -arcLength :: CircularArc -> Double -arcLength arc = _r arc * abs(_sweepAngle arc) - diff --git a/src/CubicBezier.hs b/src/CubicBezier.hs deleted file mode 100644 --- a/src/CubicBezier.hs +++ /dev/null @@ -1,68 +0,0 @@ -module CubicBezier ( CubicBezier (..) - , pointAt - , bezierSplitAt - , isClockwise - , inflectionPoints - ) where - -import Linear -import Control.Lens -import Data.Complex - -data CubicBezier = CubicBezier { _p1 :: V2 Double - , _c1 :: V2 Double - , _c2 :: V2 Double - , _p2 :: V2 Double - } deriving Show - -pointAt :: CubicBezier -> Double -> V2 Double -pointAt bezier t = ((1 - t) ** 3) *^ _p1 bezier + - ((1 - t) ** 2) * 3 * t *^ _c1 bezier + - (t ** 2) * (1 - t) * 3 *^ _c2 bezier + - (t ** 3) *^ _p2 bezier - -bezierSplitAt :: CubicBezier -> Double -> (CubicBezier, CubicBezier) -bezierSplitAt bezier t = (CubicBezier (_p1 bezier) p0 p01 dp, CubicBezier dp p12 p2 (_p2 bezier)) - where - p0 = _p1 bezier + t *^ (_c1 bezier - _p1 bezier) - p1 = _c1 bezier + t *^ (_c2 bezier - _c1 bezier) - p2 = _c2 bezier + t *^ (_p2 bezier - _c2 bezier) - - p01 = p0 + t *^ (p1 - p0) - p12 = p1 + t *^ (p2 - p1) - - dp = p01 + t *^ (p12 - p01) - -isClockwise :: CubicBezier -> Bool -isClockwise bezier = s < 0 - where - s = (_c1 bezier ^. _x - _p1 bezier ^. _x) * (_c1 bezier ^. _y + _p1 bezier ^. _y) - + (_c2 bezier ^. _x - _c1 bezier ^. _x) * (_c2 bezier ^. _y + _c1 bezier ^. _y) - + (_p2 bezier ^. _x - _c2 bezier ^. _x) * (_p2 bezier ^. _y + _c2 bezier ^. _y) - + (_p1 bezier ^. _x - _p2 bezier ^. _x) * (_p1 bezier ^. _y + _p2 bezier ^. _y) - -inflectionPoints :: CubicBezier -> [Double] -inflectionPoints bezier - | a /= 0 = realInflectionPoints [t1, t2] - | otherwise = realInflectionPoints [t] - where - pa = _c1 bezier - _p1 bezier - pb = _c2 bezier - _c1 bezier - pa - pc = _p2 bezier - _c2 bezier - pa - 2 *^ pb - - a = (pb ^. _x * pc ^. _y - pb ^. _y * pc ^. _x) :+ 0 - b = (pa ^. _x * pc ^. _y - pa ^. _y * pc ^. _x) :+ 0 - c = (pa ^. _x * pb ^. _y - pa ^. _y * pb ^. _x) :+ 0 - - -- linear case - t = -c / b - - -- quadratic case - t1 = (-b + sqrt (b * b - 4 * a * c)) / (2 * a) - t2 = (-b - sqrt (b * b - 4 * a * c)) / (2 * a) - -realInflectionPoints :: [Complex Double] -> [Double] -realInflectionPoints = map realPart . filter isInflectionPoint - -isInflectionPoint :: Complex Double -> Bool -isInflectionPoint c = imagPart c == 0 && realPart c > 0 && realPart c < 1 diff --git a/src/GCode.hs b/src/GCode.hs --- a/src/GCode.hs +++ b/src/GCode.hs @@ -6,7 +6,8 @@ import Data.List import Text.Printf -import Types +import Graphics.Path +import Utils data GCodeFlavor = GCodeFlavor { _begin :: String , _end :: String @@ -17,7 +18,7 @@ defaultFlavor :: GCodeFlavor defaultFlavor = GCodeFlavor "G17\nG90\nG0 Z1\nG0 X0 Y0" "G0 Z1" "G01 Z0 F10.00" "G00 Z1" -toString :: GCodeFlavor -> Int -> [GCodeOp] -> String +toString :: GCodeFlavor -> Int -> [PathCommand] -> String toString (GCodeFlavor begin end on off) dpi gops = begin ++ "\n" ++ @@ -32,15 +33,15 @@ mm :: Double -> Double mm px = (px * 2.54 * 10) / dd - toString' (GMoveTo p@(x,y) : gs) _ False + toString' (MoveTo p@(x,y) : gs) _ False = printf "G00 X%.4f Y%.4f" (mm x) (mm y) : toString' gs p False - toString' (GMoveTo p@(x,y) : gs) _ True + toString' (MoveTo p@(x,y) : gs) _ True = off : printf "G00 X%.4f Y%.4f" (mm x) (mm y) : toString' gs p False toString' gs cp False = on : toString' gs cp True - toString' (GLineTo p@(x,y) : gs) _ True + toString' (LineTo p@(x,y) : gs) _ True = printf "G01 X%.4f Y%.4f" (mm x) (mm y) : toString' gs p True - toString' (GArcTo (ox,oy) p@(x,y) cw : gs) (cx,cy) True + toString' (ArcTo (ox,oy) p@(x,y) cw : gs) (cx,cy) True = arcStr : toString' gs p True where i = ox - cx @@ -49,7 +50,7 @@ cmd = if' cw "G03" "G02" arcStr = printf "%s X%.4f Y%.4f I%.4f J%.4f" cmd (mm x) (mm y) (mm i) (mm j) - toString' (GBezierTo (c1x,c1y) (c2x,c2y) p2@(p2x,p2y) : gs) (p1x,p1y) True + toString' (BezierTo (c1x,c1y) (c2x,c2y) p2@(p2x,p2y) : gs) (p1x,p1y) True = bStr : toString' gs p2 True where i = c1x - p1x diff --git a/src/Graphics/BiArc.hs b/src/Graphics/BiArc.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/BiArc.hs @@ -0,0 +1,87 @@ +module Graphics.BiArc ( + BiArc (..) + , fromPoints + , arcLength + , jointAt +) where + +import qualified Graphics.CircularArc as CA +import qualified Graphics.Line as L + +import Linear hiding (angle) +import Control.Lens + +import Graphics.Curve + +data BiArc = BiArc { _a1 :: CA.CircularArc + , _a2 :: CA.CircularArc + } deriving Show + +instance Curve BiArc where + pointAt arc t + | t <= s + = pointAt (_a1 arc) (t / s) + | otherwise + = pointAt (_a2 arc) ((t - s) / (1 - s)) + where + s = jointAt arc + +fromPoints :: V2 Double -- Start point + -> V2 Double -- Tangent vector at start point + -> V2 Double -- End point + -> V2 Double -- Tangent vector at end point + -> V2 Double -- Transition point (connection point of the arcs) + -> BiArc +fromPoints p1 t1 p2 t2 t + = BiArc (CA.CircularArc c1 r1 startAngle1 sweepAngle1 p1 t) (CA.CircularArc c2 r2 startAngle2 sweepAngle2 t p2) + where + -- Calculate the orientation + osum = (t ^. _x - p1 ^. _x) * (t ^. _y + p1 ^. _y) + + (p2 ^. _x - t ^. _x) * (p2 ^. _y + t ^. _y) + + (p1 ^. _x - p2 ^. _x) * (p1 ^. _y + p2 ^. _y) + cw = osum < 0 + + -- Calculate perpendicular lines to the tangent at P1 and P2 + tl1 = L.createPerpendicularAt p1 (p1 + t1) + tl2 = L.createPerpendicularAt p2 (p2 + t2) + + -- Calculate the perpendicular bisector of P1T and P2T + p1t2 = (p1 + t) ^/ 2 + pb_p1t = L.createPerpendicularAt p1t2 t + + p2t2 = (p2 + t) ^/ 2 + pb_p2t = L.createPerpendicularAt p2t2 t + + -- The origo of the circles are at the intersection points + c1 = L.intersection tl1 pb_p1t + c2 = L.intersection tl2 pb_p2t + + -- Calculate the radii + r1 = distance c1 p1 + r2 = distance c2 p2 + + -- Calculate start and sweep angles + startVector1 = p1 - c1; + endVector1 = t - c1; + startAngle1 = atan2 (startVector1 ^. _y) (startVector1 ^. _x) + sweepAngle1' = atan2 (endVector1 ^. _y) (endVector1 ^. _x) - startAngle1 + + startVector2 = t - c2 + endVector2 = p2 - c2 + startAngle2 = atan2 (startVector2 ^. _y) (startVector2 ^. _x) + sweepAngle2' = atan2 (endVector2 ^. _y) (endVector2 ^. _x) - startAngle2 + + -- Adjust angles according to the orientation of the curve + sweepAngle1 = adjustSweepAngle cw sweepAngle1' + sweepAngle2 = adjustSweepAngle cw sweepAngle2' + +adjustSweepAngle :: Bool -> Double -> Double +adjustSweepAngle True angle | angle < 0 = 2 * pi + angle +adjustSweepAngle False angle | angle > 0 = angle - 2 * pi +adjustSweepAngle _ angle = angle + +arcLength :: BiArc -> Double +arcLength arc = CA.arcLength (_a1 arc) + CA.arcLength (_a2 arc) + +jointAt :: BiArc -> Double +jointAt arc = CA.arcLength (_a1 arc) / arcLength arc diff --git a/src/Graphics/CircularArc.hs b/src/Graphics/CircularArc.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/CircularArc.hs @@ -0,0 +1,34 @@ +module Graphics.CircularArc ( + CircularArc (..) + , isClockwise + , arcLength +) where + +import Linear +import Control.Lens + +import Graphics.Curve + +-- Would be enough one of these sets: +-- 1. c, r, startAngle, sweepAngle +-- 2. c, r, p1, p1, direction +data CircularArc = CircularArc { _c :: V2 Double + , _r :: Double + , _startAngle :: Double + , _sweepAngle :: Double + , _p1 :: V2 Double + , _p2 :: V2 Double + } deriving Show + +instance Curve CircularArc where + pointAt arc t = V2 x y + where + x = _c arc ^. _x + _r arc * cos (_startAngle arc + t * _sweepAngle arc) + y = _c arc ^. _y + _r arc * sin (_startAngle arc + t * _sweepAngle arc) + +isClockwise :: CircularArc -> Bool +isClockwise arc = _sweepAngle arc > 0 + +arcLength :: CircularArc -> Double +arcLength arc = _r arc * abs(_sweepAngle arc) + diff --git a/src/Graphics/CubicBezier.hs b/src/Graphics/CubicBezier.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/CubicBezier.hs @@ -0,0 +1,86 @@ +module Graphics.CubicBezier ( + CubicBezier (..) + , firstDerivativeAt + , secondDerivativeAt + , splitAt + , isClockwise + , inflectionPoints + , maxArcLength +) where + +import Prelude hiding (splitAt) + +import Linear +import Control.Lens +import Data.Complex + +import Graphics.Curve + +data CubicBezier = CubicBezier { _p1 :: V2 Double + , _c1 :: V2 Double + , _c2 :: V2 Double + , _p2 :: V2 Double + } deriving Show + +instance Curve CubicBezier where + pointAt bezier t = (1 - t) ** 3 *^ _p1 bezier + + (1 - t) ** 2 * 3 * t *^ _c1 bezier + + t ** 2 * (1 - t) * 3 *^ _c2 bezier + + t ** 3 *^ _p2 bezier + +firstDerivativeAt :: CubicBezier -> Double -> V2 Double +firstDerivativeAt bezier t = (1 - t) ** 2 * 3 *^ (_c1 bezier - _p1 bezier) + + (1 - t) * t * 6 *^ (_c2 bezier - _c1 bezier) + + t * t * 3 *^ (_p2 bezier - _c2 bezier) + +secondDerivativeAt :: CubicBezier -> Double -> V2 Double +secondDerivativeAt bezier t = (1 - t) * 6 *^ (_c2 bezier - 2 *^ _c1 bezier + _p1 bezier) + + t * 6 *^ (_p2 bezier - 2 *^ _c2 bezier + _c1 bezier) + + +splitAt :: CubicBezier -> Double -> (CubicBezier, CubicBezier) +splitAt bezier t = (CubicBezier (_p1 bezier) p0 p01 dp, CubicBezier dp p12 p2 (_p2 bezier)) + where + p0 = _p1 bezier + t *^ (_c1 bezier - _p1 bezier) + p1 = _c1 bezier + t *^ (_c2 bezier - _c1 bezier) + p2 = _c2 bezier + t *^ (_p2 bezier - _c2 bezier) + + p01 = p0 + t *^ (p1 - p0) + p12 = p1 + t *^ (p2 - p1) + + dp = p01 + t *^ (p12 - p01) + +isClockwise :: CubicBezier -> Bool +isClockwise bezier = isClockwise4 (_p1 bezier) (_c1 bezier) (_c2 bezier) (_p2 bezier) + +inflectionPoints :: CubicBezier -> [Double] +inflectionPoints bezier + | a /= 0 = realInflectionPoints [t1, t2] + | otherwise = realInflectionPoints [t] + where + pa = _c1 bezier - _p1 bezier + pb = _c2 bezier - _c1 bezier - pa + pc = _p2 bezier - _c2 bezier - pa - 2 *^ pb + + a = (pb ^. _x * pc ^. _y - pb ^. _y * pc ^. _x) :+ 0 + b = (pa ^. _x * pc ^. _y - pa ^. _y * pc ^. _x) :+ 0 + c = (pa ^. _x * pb ^. _y - pa ^. _y * pb ^. _x) :+ 0 + + -- linear case + t = -c / b + + -- quadratic case + t1 = (-b + sqrt (b * b - 4 * a * c)) / (2 * a) + t2 = (-b - sqrt (b * b - 4 * a * c)) / (2 * a) + +realInflectionPoints :: [Complex Double] -> [Double] +realInflectionPoints = map realPart . filter isInflectionPoint + +isInflectionPoint :: Complex Double -> Bool +isInflectionPoint c = imagPart c == 0 && realPart c > 0 && realPart c < 1 + +maxArcLength :: CubicBezier -> Double +maxArcLength bezier = + distance (_p1 bezier) (_c1 bezier) + + distance (_c1 bezier) (_c2 bezier) + + distance (_c2 bezier) (_p2 bezier) diff --git a/src/Graphics/Curve.hs b/src/Graphics/Curve.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Curve.hs @@ -0,0 +1,26 @@ +module Graphics.Curve ( + Curve(..), + isClockwise4, + isClockwise3 +) where + +import Linear +import Control.Lens + +class Curve c where + pointAt :: c -> Double -> V2 Double + +isClockwise4 :: V2 Double -> V2 Double -> V2 Double -> V2 Double -> Bool +isClockwise4 p1 p2 p3 p4 = s < 0 + where + s = (p2 ^. _x - p1 ^. _x) * (p2 ^. _y + p1 ^. _y) + + (p3 ^. _x - p2 ^. _x) * (p3 ^. _y + p2 ^. _y) + + (p4 ^. _x - p3 ^. _x) * (p4 ^. _y + p3 ^. _y) + + (p1 ^. _x - p4 ^. _x) * (p1 ^. _y + p4 ^. _y) + +isClockwise3 :: V2 Double -> V2 Double -> V2 Double -> Bool +isClockwise3 p1 p2 p3 = s < 0 + where + s = (p3 ^. _x - p1 ^. _x) * (p3 ^. _y + p1 ^. _y) + + (p2 ^. _x - p3 ^. _x) * (p2 ^. _y + p3 ^. _y) + + (p1 ^. _x - p2 ^. _x) * (p1 ^. _y + p2 ^. _y) diff --git a/src/Graphics/Line.hs b/src/Graphics/Line.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Line.hs @@ -0,0 +1,74 @@ +module Graphics.Line ( + Line (..) + , throughPoint + , fromPoints + , createPerpendicularAt + , slope + , intersection + , isOnLine +) where + +import Linear +import Control.Lens + +-- TODO: letting _p to be NaN is actually a really bad idea +data Line = Line { _m :: Double + , _p :: V2 Double + } deriving Show + +throughPoint :: V2 Double -> Double -> Line +throughPoint p m = Line m p + +fromPoints :: V2 Double -> V2 Double -> Line +fromPoints p1 p2 = throughPoint p1 (slope p1 p2) + +-- Creates a a line which is perpendicular to the line defined by P and P1 and goes through P +createPerpendicularAt :: V2 Double -> V2 Double -> Line +createPerpendicularAt p p1 + | m == 0 + = throughPoint p nan + | isNaN m + = throughPoint p 0 + | otherwise + = throughPoint p (-1 / m) + where + m = slope p p1 + +slope :: V2 Double -> V2 Double -> Double +slope p1 p2 + | p2 ^. _x == p1 ^. _x + = nan + | otherwise + = (p2 ^. _y - p1 ^. _y) / (p2 ^. _x - p1 ^. _x) + +nan :: Double +nan = 0/0 + +-- If the solution is not found it actually returns +/-infinity +intersection :: Line -> Line -> V2 Double +intersection line1 line2 + | isNaN (_m line1) + = verticalIntersection line1 line2 + | isNaN (_m line2) + = verticalIntersection line2 line1 + | otherwise + = V2 x y + where + x = (_m line1 * _p line1 ^. _x - _m line2 * _p line2 ^. _x - _p line1 ^. _y + _p line2 ^. _y) / (_m line1 - _m line2) + y = _m line1 * x - _m line1 * _p line1 ^. _x + _p line1 ^. _y + +-- First line is vertical +verticalIntersection :: Line -> Line -> V2 Double +verticalIntersection vline line = V2 x y + where + x = _p vline ^. _x + y = _m line * (x - _p line ^. _x) + _p line ^. _y + +isOnLine :: Line -> V2 Double -> Bool +isOnLine l p2 + | isNaN (_m l) + = p1 ^. _x == p2 ^. _x + | otherwise + = (p2 ^. _x - p1 ^. _x) * _m l == (p2 ^. _y - p1 ^. _y) + where + p1 = _p l diff --git a/src/Graphics/LineSegment.hs b/src/Graphics/LineSegment.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/LineSegment.hs @@ -0,0 +1,17 @@ +module Graphics.LineSegment ( + fromPoints +) where + +import Linear + +import Graphics.Curve + +data LineSegment = LineSegment { _p1 :: V2 Double + , _p2 :: V2 Double + } deriving Show + +fromPoints :: V2 Double -> V2 Double -> LineSegment +fromPoints = LineSegment + +instance Curve LineSegment where + pointAt ls t = _p1 ls + ((_p2 ls - _p1 ls) ^* t) diff --git a/src/Graphics/Path.hs b/src/Graphics/Path.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Path.hs @@ -0,0 +1,13 @@ +module Graphics.Path ( + PathCommand(..) +) where + +import Graphics.Point + +-- all of them are invariant under affine transformation +data PathCommand + = MoveTo Point + | LineTo Point -- End point + | ArcTo Point Point Bool -- Center point, end point, clockwise + | BezierTo Point Point Point -- Control point1, control point2, end point + deriving Show diff --git a/src/Graphics/Point.hs b/src/Graphics/Point.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Point.hs @@ -0,0 +1,15 @@ +module Graphics.Point ( + Point + , toPoint + , fromPoint +) where + +import Linear + +type Point = (Double, Double) -- A point in the plane, absolute coordinates + +toPoint :: V2 Double -> Point +toPoint (V2 x y) = (x, y) + +fromPoint :: Point -> V2 Double +fromPoint (x, y) = V2 x y diff --git a/src/Graphics/Transformation.hs b/src/Graphics/Transformation.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Transformation.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE FlexibleInstances #-} + +module Graphics.Transformation ( + TransformationMatrix + , fromElements + , identityTransform + , mirrorYTransform + , multiply + , translateTransform + , scaleTransform + , Transform(..) + ) where + +import Data.Matrix as M + +import Graphics.Point +import Graphics.Path + +type TransformationMatrix = Matrix Double + +identityTransform :: TransformationMatrix +identityTransform = identity 3 + +mirrorYTransform :: Double -> Double -> TransformationMatrix +mirrorYTransform _ h = fromElements [1, 0, 0, -1, 0, h] + +translateTransform :: Double -> Double -> TransformationMatrix +translateTransform x y = fromElements [1, 0, 0, 1, x, y] + +scaleTransform :: Double -> Double -> TransformationMatrix +scaleTransform sx sy = fromElements [sx, 0, 0, sy, 0, 0] + +multiply :: TransformationMatrix -> TransformationMatrix -> TransformationMatrix +multiply = multStd + +fromElements :: [Double] -> TransformationMatrix +fromElements [a,b,c,d,e,f] = fromList 3 3 [a,c,e,b,d,f,0,0,1] +fromElements _ = error "Malformed transformation matrix" + +class Transform t where + transform :: TransformationMatrix -> t -> t + +instance Transform Point where + transform m (x,y) = (a * x + c * y + e, b * x + d * y + f) + where + (a:c:e:b:d:f:_) = M.toList m + +instance Transform PathCommand where + transform m (MoveTo p) = MoveTo (transform m p) + transform m (LineTo p) = LineTo (transform m p) + transform m (ArcTo p1 p2 d) = ArcTo (transform m p1) (transform m p2) d + transform m (BezierTo c1 c2 p2) = BezierTo (transform m c1) (transform m c2) (transform m p2) + + diff --git a/src/Line.hs b/src/Line.hs deleted file mode 100644 --- a/src/Line.hs +++ /dev/null @@ -1,73 +0,0 @@ -module Line ( Line (..) - , throughPoint - , fromPoints - , createPerpendicularAt - , slope - , intersection - , isOnLine - ) where - -import Linear -import Control.Lens - --- TODO: letting _p to be NaN is actually a really bad idea -data Line = Line { _m :: Double - , _p :: V2 Double - } deriving Show - -throughPoint :: V2 Double -> Double -> Line -throughPoint p m = Line m p - -fromPoints :: V2 Double -> V2 Double -> Line -fromPoints p1 p2 = throughPoint p1 (slope p1 p2) - --- Creates a a line which is perpendicular to the line defined by P and P1 and goes through P -createPerpendicularAt :: V2 Double -> V2 Double -> Line -createPerpendicularAt p p1 - | m == 0 - = throughPoint p nan - | isNaN m - = throughPoint p 0 - | otherwise - = throughPoint p (-1 / m) - where - m = slope p p1 - -slope :: V2 Double -> V2 Double -> Double -slope p1 p2 - | p2 ^. _x == p1 ^. _x - = nan - | otherwise - = (p2 ^. _y - p1 ^. _y) / (p2 ^. _x - p1 ^. _x) - -nan :: Double -nan = 0/0 - --- If the solution is not found it actually returns +/-infinity -intersection :: Line -> Line -> V2 Double -intersection line1 line2 - | isNaN (_m line1) - = verticalIntersection line1 line2 - | isNaN (_m line2) - = verticalIntersection line2 line1 - | otherwise - = V2 x y - where - x = (_m line1 * _p line1 ^. _x - _m line2 * _p line2 ^. _x - _p line1 ^. _y + _p line2 ^. _y) / (_m line1 - _m line2) - y = _m line1 * x - _m line1 * _p line1 ^. _x + _p line1 ^. _y - --- First line is vertical -verticalIntersection :: Line -> Line -> V2 Double -verticalIntersection vline line = V2 x y - where - x = _p vline ^. _x - y = _m line * (x - _p line ^. _x) + _p line ^. _y - -isOnLine :: Line -> V2 Double -> Bool -isOnLine l p2 - | isNaN (_m l) - = p1 ^. _x == p2 ^. _x - | otherwise - = (p2 ^. _x - p1 ^. _x) * (_m l) == (p2 ^. _y - p1 ^. _y) - where - p1 = _p l diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -1,4 +1,6 @@ {-# LANGUAGE TemplateHaskell #-} +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Use lambda-case" #-} import qualified Graphics.Svg as SVG @@ -26,32 +28,32 @@ <$> argument str ( metavar "SVGFILE" <> help "The SVG file to be converted" ) - <*> (optional $ strOption + <*> optional (strOption ( long "flavor" <> short 'f' <> metavar "CONFIGFILE" <> help "Configuration of G-Code flavor" )) - <*> (optional $ strOption + <*> optional (strOption ( long "output" <> short 'o' <> metavar "OUTPUTFILE" <> help "The output G-Code file (default is standard output)" )) - <*> (option auto + <*> option auto ( long "dpi" <> value 96 <> short 'd' <> metavar "DPI" - <> help "Used to determine the size of the SVG when it does not contain any units; dot per inch (default is 96)" )) - <*> (option auto + <> help "Used to determine the size of the SVG when it does not contain any units; dot per inch (default is 96)" ) + <*> option auto ( long "resolution" <> value 0.1 <> short 'r' <> metavar "RESOLUTION" - <> help "Shorter paths are replaced by line segments; mm (default is 0.1)" )) - <*> (switch + <> help "Shorter paths are replaced by line segments; mm (default is 0.1)" ) + <*> switch ( long "generate-bezier" <> short 'b' - <> help "Generate bezier curves (G5) instead of arcs (G2,G3)" )) + <> help "Generate bezier curves (G5) instead of arcs (G2,G3)" ) runWithOptions :: Options -> IO () runWithOptions (Options svgFile mbCfg mbOut dpi resolution generateBezier) = @@ -62,7 +64,7 @@ (Just doc) -> writer (toString flavor dpi $ renderDoc generateBezier dpi resolution doc) Nothing -> putStrLn "juicy-gcode: error during opening the SVG file" where - writer = maybe putStr (\fn -> writeFile fn) mbOut + writer = maybe putStr writeFile mbOut toLines :: Text -> String toLines t = unpack $ replace (pack ";") (pack "\n") t @@ -77,7 +79,7 @@ return $ GCodeFlavor (toLines begin) (toLines end) (toLines toolon) (toLines tooloff) versionOption :: Parser (a -> a) -versionOption = infoOption +versionOption = infoOption (concat ["juicy-gcode ", showVersion version, ", git revision ", $(gitHash)]) (long "version" <> short 'v' <> help "Show version") diff --git a/src/Render.hs b/src/Render.hs --- a/src/Render.hs +++ b/src/Render.hs @@ -1,19 +1,21 @@ -module Render ( renderDoc - ) where +module Render ( + renderDoc +) where +import Data.Maybe ( fromMaybe ) + import qualified Graphics.Svg as SVG import qualified Graphics.Svg.CssTypes as CSS import qualified Linear -import Types -import Transformation +import Graphics.Path +import Graphics.Point +import Graphics.Transformation +import Approx.BiArc import SvgArcSegment -import Approx import SVGExt -import qualified CircularArc as CA -import qualified BiArc as BA -import qualified CubicBezier as B +import qualified Graphics.CubicBezier as B mapTuple :: (a -> b) -> (a, a) -> (b, b) mapTuple f (a1, a2) = (f a1, f a2) @@ -24,12 +26,6 @@ fromRPoint :: SVG.RPoint -> Point fromRPoint (Linear.V2 x y) = (x, y) -toPoint :: Linear.V2 Double -> Point -toPoint (Linear.V2 x y) = (x, y) - -fromPoint :: Point -> Linear.V2 Double -fromPoint (x, y) = (Linear.V2 x y) - -- TODO: em, percentage fromSvgNumber :: Int -> SVG.Number -> Double fromSvgNumber dpi num = fromNumber' (CSS.toUserUnit dpi num) @@ -42,16 +38,38 @@ mirrorControlPoint (cx, cy) (cpx, cpy) = (cx + cx - cpx, cy + cy - cpy) -- convert a quadratic bezier to a cubic one -bezierQ2C :: Point -> Point -> Point -> DrawOp +bezierQ2C :: Point -> Point -> Point -> PathCommand bezierQ2C (qp0x, qp0y) (qp1x, qp1y) (qp2x, qp2y) - = DBezierTo (qp0x + 2.0 / 3.0 * (qp1x - qp0x), qp0y + 2.0 / 3.0 * (qp1y - qp0y)) - (qp2x + 2.0 / 3.0 * (qp1x - qp2x), qp2y + 2.0 / 3.0 * (qp1y - qp2y)) - (qp2x, qp2y) + = BezierTo (qp0x + 2.0 / 3.0 * (qp1x - qp0x), qp0y + 2.0 / 3.0 * (qp1y - qp0y)) + (qp2x + 2.0 / 3.0 * (qp1x - qp2x), qp2y + 2.0 / 3.0 * (qp1y - qp2y)) + (qp2x, qp2y) toAbsolute :: (Double, Double) -> SVG.Origin -> (Double, Double) -> (Double, Double) toAbsolute _ SVG.OriginAbsolute p = p toAbsolute (cx,cy) SVG.OriginRelative (dx,dy) = (cx+dx, cy+dy) +-- Apply SVG transformations to a TransformationMatrix +applyTransformations :: TransformationMatrix -> Maybe [SVG.Transformation] -> TransformationMatrix +applyTransformations m Nothing = m +applyTransformations m (Just ts) = foldl applyTransformation m ts + +radiansPerDegree :: Double +radiansPerDegree = pi / 180.0 + +-- https://developer.mozilla.org/en/docs/Web/SVG/Attribute/transform +applyTransformation :: TransformationMatrix -> SVG.Transformation -> TransformationMatrix +applyTransformation m (SVG.TransformMatrix a b c d e f) = multiply m (fromElements [a,b,c,d,e,f]) +applyTransformation m (SVG.Translate x y) = multiply m (fromElements [1,0,0,1,x,y]) +applyTransformation m (SVG.Scale sx mbSy) = multiply m (fromElements [sx,0,0,Data.Maybe.fromMaybe sx mbSy,0,0]) +applyTransformation m (SVG.Rotate a Nothing) + = multiply m (fromElements [cos r, sin r, -sin r, cos r , 0, 0]) + where + r = a * radiansPerDegree +applyTransformation m (SVG.Rotate a (Just (x, y))) = applyTransformations m (Just [SVG.Translate x y , SVG.Rotate a Nothing , SVG.Translate (-x) (-y)]) +applyTransformation m (SVG.SkewX a) = multiply m (fromElements [1,0,tan(a*radiansPerDegree),1,0,0]) +applyTransformation m (SVG.SkewY a) = multiply m (fromElements [1,tan(a*radiansPerDegree),0,1,0,0]) +applyTransformation m SVG.TransformUnknown = m + docTransform :: Int -> SVG.Document -> TransformationMatrix docTransform dpi doc = multiply mirrorTransform (viewBoxTransform $ SVG._viewBox doc) where @@ -62,37 +80,33 @@ mirrorTransform = mirrorYTransform w h - (w, h) = (documentSize dpi doc) + (w, h) = documentSize dpi doc -renderDoc :: Bool -> Int -> Double -> SVG.Document -> [GCodeOp] +renderDoc :: Bool -> Int -> Double -> SVG.Document -> [PathCommand] renderDoc generateBezier dpi resolution doc = stage2 $ renderTrees (docTransform dpi doc) (SVG._elements doc) where - pxresolution = (fromIntegral dpi) / 2.45 / 10 * resolution + pxresolution = fromIntegral dpi / 2.45 / 10 * resolution -- TODO: make it tail recursive - stage2 :: [DrawOp] -> [GCodeOp] - stage2 dops = convert dops (Linear.V2 0 0) + stage2 :: [PathCommand] -> [PathCommand] + stage2 dops = approximate dops (Linear.V2 0 0) where - convert [] _ = [] - convert (DMoveTo p:ds) _ = GMoveTo p : convert ds (fromPoint p) - convert (DLineTo p:ds) _ = GLineTo p : convert ds (fromPoint p) - convert (DBezierTo c1 c2 p2:ds) cp - | generateBezier - = [GBezierTo c1 c2 p2] ++ convert ds (fromPoint p2) - | otherwise - = concatMap biarc2garc biarcs ++ convert ds (fromPoint p2) - where - biarcs = bezier2biarc (B.CubicBezier cp (fromPoint c1) (fromPoint c2) (fromPoint p2)) pxresolution - biarc2garc (Left biarc) - = [arc2garc (BA._a1 biarc), arc2garc (BA._a2 biarc)] - biarc2garc (Right (Linear.V2 x y)) - = [GLineTo (x,y)] - arc2garc arc = GArcTo (toPoint (CA._c arc)) (toPoint (CA._p2 arc)) (CA.isClockwise arc) + approximate [] _ = [] + approximate (MoveTo p:ds) _ = MoveTo p : approximate ds (fromPoint p) + approximate (LineTo p:ds) _ = LineTo p : approximate ds (fromPoint p) + approximate (ArcTo p1 p2 d:ds) _ = ArcTo p1 p2 d : approximate ds (fromPoint p2) + approximate (BezierTo c1 c2 p2:ds) cp + | generateBezier + = BezierTo c1 c2 p2 : approximate ds (fromPoint p2) + | otherwise + = bezier2biarcs + (B.CubicBezier cp (fromPoint c1) (fromPoint c2) (fromPoint p2)) pxresolution + ++ approximate ds (fromPoint p2) - renderPathCommands :: Point -> Point -> Maybe Point -> [SVG.PathCommand] -> [DrawOp] + renderPathCommands :: Point -> Point -> Maybe Point -> [SVG.PathCommand] -> [PathCommand] renderPathCommands _ currentp _ (SVG.MoveTo origin (p:ps):ds) - = DMoveTo ap : renderPathCommands ap ap Nothing (cont ps) + = MoveTo ap : renderPathCommands ap ap Nothing (cont ps) where ap = toAbsolute currentp origin (fromRPoint p) @@ -100,7 +114,7 @@ cont ps' = SVG.LineTo origin ps' : ds renderPathCommands firstp currentp _ (SVG.LineTo origin (p:ps):ds) - = DLineTo ap : renderPathCommands firstp ap Nothing (cont ps) + = LineTo ap : renderPathCommands firstp ap Nothing (cont ps) where ap = toAbsolute currentp origin (fromRPoint p) @@ -108,7 +122,7 @@ cont ps' = SVG.LineTo origin ps' : ds renderPathCommands firstp (_, cy) _ (SVG.HorizontalTo SVG.OriginAbsolute (px:pxs):ds) - = DLineTo ap : renderPathCommands firstp ap Nothing (cont pxs) + = LineTo ap : renderPathCommands firstp ap Nothing (cont pxs) where ap = (px,cy) @@ -116,7 +130,7 @@ cont pxs' = SVG.HorizontalTo SVG.OriginAbsolute pxs' : ds renderPathCommands firstp (cx, cy) _ (SVG.HorizontalTo SVG.OriginRelative (dx:dxs):ds) - = DLineTo ap : renderPathCommands firstp ap Nothing (cont dxs) + = LineTo ap : renderPathCommands firstp ap Nothing (cont dxs) where ap = (cx+dx,cy) @@ -124,7 +138,7 @@ cont dxs' = SVG.HorizontalTo SVG.OriginRelative dxs' : ds renderPathCommands firstp (cx, _) _ (SVG.VerticalTo SVG.OriginAbsolute (py:pys):ds) - = DLineTo ap : renderPathCommands firstp ap Nothing (cont pys) + = LineTo ap : renderPathCommands firstp ap Nothing (cont pys) where ap = (cx,py) @@ -132,7 +146,7 @@ cont pys' = SVG.VerticalTo SVG.OriginAbsolute pys' : ds renderPathCommands firstp (cx, cy) _ (SVG.VerticalTo SVG.OriginRelative (dy:dys):ds) - = DLineTo ap : renderPathCommands firstp ap Nothing (cont dys) + = LineTo ap : renderPathCommands firstp ap Nothing (cont dys) where ap = (cx,cy+dy) @@ -140,7 +154,7 @@ cont dys' = SVG.VerticalTo SVG.OriginRelative dys' : ds renderPathCommands firstp currentp _ (SVG.CurveTo origin ((c1,c2,p):ps):ds) - = DBezierTo ac1 ac2 ap : renderPathCommands firstp ap (Just ac2) (cont ps) + = BezierTo ac1 ac2 ap : renderPathCommands firstp ap (Just ac2) (cont ps) where ap = toAbsolute currentp origin (fromRPoint p) ac1 = toAbsolute currentp origin (fromRPoint c1) @@ -150,7 +164,7 @@ cont ps' = SVG.CurveTo origin ps' : ds renderPathCommands firstp currentp mbControlp (SVG.SmoothCurveTo origin ((c2,p):ps):ds) - = DBezierTo ac1 ac2 ap : renderPathCommands firstp ap (Just ac2) (cont ps) + = BezierTo ac1 ac2 ap : renderPathCommands firstp ap (Just ac2) (cont ps) where ap = toAbsolute currentp origin (fromRPoint p) ac1 = maybe ac2 (mirrorControlPoint currentp) mbControlp @@ -191,28 +205,28 @@ renderPathCommands firstp@(fx,fy) (cx,cy) mbControlp (SVG.EndPath:ds) | fx /= cx || fy /= cy - = DLineTo firstp : renderPathCommands firstp firstp mbControlp ds + = LineTo firstp : renderPathCommands firstp firstp mbControlp ds | otherwise = renderPathCommands firstp firstp mbControlp ds renderPathCommands _ _ _ _ = [] - renderTree :: TransformationMatrix -> SVG.Tree -> [DrawOp] + renderTree :: TransformationMatrix -> SVG.Tree -> [PathCommand] renderTree m (SVG.GroupTree g) = renderTrees (applyTransformations m (SVG._transform (SVG._groupDrawAttributes g))) (SVG._groupChildren g) - renderTree m (SVG.PathTree p) = map (transformDrawOp tr) $ renderPathCommands (0,0) (0,0) Nothing (SVG._pathDefinition p) + renderTree m (SVG.PathTree p) = map (transform tr) $ renderPathCommands (0,0) (0,0) Nothing (SVG._pathDefinition p) where tr = applyTransformations m (SVG._transform (SVG._pathDrawAttributes p)) renderTree m (SVG.RectangleTree r) | rx == 0.0 && ry == 0.0 - = map (transformDrawOp tr) [DMoveTo (x,y), DLineTo (x+w,y), DLineTo (x+w,y+h), DLineTo (x,y+h), DLineTo (x,y)] + = map (transform tr) [MoveTo (x,y), LineTo (x+w,y), LineTo (x+w,y+h), LineTo (x,y+h), LineTo (x,y)] | otherwise - = map (transformDrawOp tr) - ([DMoveTo (x,y+ry)] ++ convertSvgArc (x,y+ry) rx ry 0 False True (x+rx, y) ++ - [DLineTo (x+w-rx,y)] ++ convertSvgArc (x+w-rx,y) rx ry 0 False True (x+w, y+ry) ++ - [DLineTo (x+w,y+h-ry)] ++ convertSvgArc (x+w,y+h-ry) rx ry 0 False True (x+w-rx, y+h) ++ - [DLineTo (x+rx,y+h)] ++ convertSvgArc (x+rx, y+h) rx ry 0 False True (x, y+h-ry) ++ - [DLineTo (x,y+ry)]) + = map (transform tr) + ([MoveTo (x,y+ry)] ++ convertSvgArc (x,y+ry) rx ry 0 False True (x+rx, y) ++ + [LineTo (x+w-rx,y)] ++ convertSvgArc (x+w-rx,y) rx ry 0 False True (x+w, y+ry) ++ + [LineTo (x+w,y+h-ry)] ++ convertSvgArc (x+w,y+h-ry) rx ry 0 False True (x+w-rx, y+h) ++ + [LineTo (x+rx,y+h)] ++ convertSvgArc (x+rx, y+h) rx ry 0 False True (x, y+h-ry) ++ + [LineTo (x,y+ry)]) where (x,y) = fromSvgPoint dpi (SVG._rectUpperLeftCorner r) w = fromSvgNumber dpi (SVG._rectWidth r) @@ -220,23 +234,23 @@ (rx, ry) = mapTuple (fromSvgNumber dpi) (SVG._rectCornerRadius r) tr = applyTransformations m (SVG._transform (SVG._rectDrawAttributes r)) - renderTree m (SVG.LineTree l) = [DMoveTo p1, DLineTo p2] + renderTree m (SVG.LineTree l) = [MoveTo p1, LineTo p2] where - p1 = transformPoint tr (fromSvgPoint dpi (SVG._linePoint1 l)) - p2 = transformPoint tr (fromSvgPoint dpi (SVG._linePoint2 l)) + p1 = transform tr (fromSvgPoint dpi (SVG._linePoint1 l)) + p2 = transform tr (fromSvgPoint dpi (SVG._linePoint2 l)) tr = applyTransformations m (SVG._transform (SVG._lineDrawAttributes l)) - renderTree m (SVG.PolyLineTree l) = map (transformDrawOp tr) (DMoveTo p0:map DLineTo ps) + renderTree m (SVG.PolyLineTree l) = map (transform tr) (MoveTo p0:map LineTo ps) where (p0:ps) = map (\(Linear.V2 x y) -> (x,y)) (SVG._polyLinePoints l) tr = applyTransformations m (SVG._transform (SVG._polyLineDrawAttributes l)) - renderTree m (SVG.PolygonTree l) = map (transformDrawOp tr) (DMoveTo p0:map DLineTo (ps ++ [p0])) + renderTree m (SVG.PolygonTree l) = map (transform tr) (MoveTo p0:map LineTo (ps ++ [p0])) where (p0:ps) = map (\(Linear.V2 x y) -> (x,y)) (SVG._polygonPoints l) tr = applyTransformations m (SVG._transform (SVG._polygonDrawAttributes l)) - renderTree m (SVG.EllipseTree e) = map (transformDrawOp tr) (DMoveTo (cx-rx,cy) : bs1++bs2++bs3++bs4) + renderTree m (SVG.EllipseTree e) = map (transform tr) (MoveTo (cx-rx,cy) : bs1++bs2++bs3++bs4) where bs1 = convertSvgArc (cx-rx, cy) rx ry 0 False True (cx, cy-ry) bs2 = convertSvgArc (cx, cy-ry) rx ry 0 False True (cx+rx, cy) @@ -248,7 +262,7 @@ ry = fromSvgNumber dpi (SVG._ellipseYRadius e) tr = applyTransformations m (SVG._transform (SVG._ellipseDrawAttributes e)) - renderTree m (SVG.CircleTree c) = map (transformDrawOp tr) (DMoveTo (cx-r,cy) : bs1++bs2++bs3++bs4) + renderTree m (SVG.CircleTree c) = map (transform tr) (MoveTo (cx-r,cy) : bs1++bs2++bs3++bs4) where bs1 = convertSvgArc (cx-r, cy) r r 0 False True (cx, cy-r) bs2 = convertSvgArc (cx, cy-r) r r 0 False True (cx+r, cy) @@ -262,5 +276,5 @@ {- The rest: None, UseTree, SymbolTree, TextTree, ImageTree -} renderTree _ _ = [] - renderTrees :: TransformationMatrix -> [SVG.Tree] -> [DrawOp] - renderTrees m es = concat $ map (renderTree m) es + renderTrees :: TransformationMatrix -> [SVG.Tree] -> [PathCommand] + renderTrees m es = concatMap (renderTree m) es diff --git a/src/SVGExt.hs b/src/SVGExt.hs --- a/src/SVGExt.hs +++ b/src/SVGExt.hs @@ -18,8 +18,8 @@ documentSize _ SVG.Document { SVG._width = Just (SVG.Num w) , SVG._height = Just (SVG.Num h) } = (w, h) -documentSize dpi doc@(SVG.Document { SVG._width = Just w - , SVG._height = Just h }) = +documentSize dpi doc@SVG.Document { SVG._width = Just w + , SVG._height = Just h } = documentSize dpi $ doc { SVG._width = Just $ SVG.toUserUnit dpi w , SVG._height = Just $ SVG.toUserUnit dpi h } diff --git a/src/SvgArcSegment.hs b/src/SvgArcSegment.hs --- a/src/SvgArcSegment.hs +++ b/src/SvgArcSegment.hs @@ -1,10 +1,12 @@ -module SvgArcSegment ( - convertSvgArc - ) where +module SvgArcSegment ( + convertSvgArc +) where -import Types - -radiansPerDegree :: Double +import Graphics.Path +import Graphics.Point +import Utils + +radiansPerDegree :: Double radiansPerDegree = pi / 180.0 calculateVectorAngle :: Double -> Double -> Double -> Double -> Double @@ -16,15 +18,15 @@ where ta = atan2 uy ux tb = atan2 vy vx - + -- ported from: https://github.com/vvvv/SVG/blob/master/Source/Paths/SvgArcSegment.cs -convertSvgArc :: Point -> Double -> Double -> Double -> Bool -> Bool -> Point -> [DrawOp] +convertSvgArc :: Point -> Double -> Double -> Double -> Bool -> Bool -> Point -> [PathCommand] convertSvgArc (x0,y0) radiusX radiusY angle largeArcFlag sweepFlag (x,y) | x0 == x && y0 == y = [] | radiusX == 0.0 && radiusY == 0.0 - = [DLineTo (x,y)] - | otherwise + = [LineTo (x,y)] + | otherwise = calcSegments x0 y0 theta1' segments' where sinPhi = sin (angle * radiansPerDegree) @@ -38,32 +40,32 @@ s = sqrt(1.0 - numerator / (radiusX * radiusX * radiusY * radiusY)) rx = if' (numerator < 0.0) (radiusX * s) radiusX ry = if' (numerator < 0.0) (radiusY * s) radiusY - root = if' (numerator < 0.0) - (0.0) - ((if' ((largeArcFlag && sweepFlag) || (not largeArcFlag && not sweepFlag)) (-1.0) 1.0) * + root = if' (numerator < 0.0) + 0.0 + (if' ((largeArcFlag && sweepFlag) || (not largeArcFlag && not sweepFlag)) (-1.0) 1.0 * sqrt(numerator / (radiusX * radiusX * y1dash * y1dash + radiusY * radiusY * x1dash * x1dash))) - + cxdash = root * rx * y1dash / ry cydash = -root * ry * x1dash / rx cx = cosPhi * cxdash - sinPhi * cydash + (x0 + x) / 2.0 cy = sinPhi * cxdash + cosPhi * cydash + (y0 + y) / 2.0 - + theta1' = calculateVectorAngle 1.0 0.0 ((x1dash - cxdash) / rx) ((y1dash - cydash) / ry) dtheta' = calculateVectorAngle ((x1dash - cxdash) / rx) ((y1dash - cydash) / ry) ((-x1dash - cxdash) / rx) ((-y1dash - cydash) / ry) - dtheta = if' (not sweepFlag && dtheta' > 0) + dtheta = if' (not sweepFlag && dtheta' > 0) (dtheta' - 2 * pi) (if' (sweepFlag && dtheta' < 0) (dtheta' + 2 * pi) dtheta') - + segments' = ceiling (abs (dtheta / (pi / 2.0))) delta = dtheta / fromInteger segments' t = 8.0 / 3.0 * sin(delta / 4.0) * sin(delta / 4.0) / sin(delta / 2.0) - - calcSegments startX startY theta1 segments + + calcSegments startX startY theta1 segments | segments == 0 = [] | otherwise - = (DBezierTo (startX + dx1, startY + dy1) (endpointX + dxe, endpointY + dye) (endpointX, endpointY) : calcSegments endpointX endpointY theta2 (segments - 1)) + = BezierTo (startX + dx1, startY + dy1) (endpointX + dxe, endpointY + dye) (endpointX, endpointY) : calcSegments endpointX endpointY theta2 (segments - 1) where cosTheta1 = cos theta1 sinTheta1 = sin theta1 @@ -80,4 +82,3 @@ dxe = t * (cosPhi * rx * sinTheta2 + sinPhi * ry * cosTheta2) dye = t * (sinPhi * rx * sinTheta2 - cosPhi * ry * cosTheta2) - diff --git a/src/Transformation.hs b/src/Transformation.hs deleted file mode 100644 --- a/src/Transformation.hs +++ /dev/null @@ -1,66 +0,0 @@ -module Transformation ( TransformationMatrix - , identityTransform - , mirrorYTransform - , translateTransform - , scaleTransform - , transformPoint - , transformDrawOp - , applyTransformations - , multiply - ) where - -import qualified Graphics.Svg as SVG -import Data.Matrix as M -import Types - -type TransformationMatrix = Matrix Double - -identityTransform :: TransformationMatrix -identityTransform = identity 3 - -mirrorYTransform :: Double -> Double -> TransformationMatrix -mirrorYTransform _ h = fromElements [1, 0, 0, -1, 0, h] - -translateTransform :: Double -> Double -> TransformationMatrix -translateTransform x y = fromElements [1, 0, 0, 1, x, y] - -scaleTransform :: Double -> Double -> TransformationMatrix -scaleTransform sx sy = fromElements [sx, 0, 0, sy, 0, 0] - -multiply :: TransformationMatrix -> TransformationMatrix -> TransformationMatrix -multiply a b = multStd a b - -fromElements :: [Double] -> TransformationMatrix -fromElements [a,b,c,d,e,f] = fromList 3 3 [a,c,e,b,d,f,0,0,1] -fromElements _ = error "Malformed transformation matrix" - -transformPoint :: TransformationMatrix -> Point -> Point -transformPoint m (x,y) = (a * x + c * y + e, b * x + d * y + f) - where - (a:c:e:b:d:f:_) = M.toList m - -transformDrawOp :: TransformationMatrix -> DrawOp -> DrawOp -transformDrawOp m (DMoveTo p) = DMoveTo (transformPoint m p) -transformDrawOp m (DLineTo p) = DLineTo (transformPoint m p) -transformDrawOp m (DBezierTo c1 c2 p2) = DBezierTo (transformPoint m c1) (transformPoint m c2) (transformPoint m p2) - -applyTransformations :: TransformationMatrix -> Maybe [SVG.Transformation] -> TransformationMatrix -applyTransformations m Nothing = m -applyTransformations m (Just ts) = foldl applyTransformation m ts - -radiansPerDegree :: Double -radiansPerDegree = pi / 180.0 - --- https://developer.mozilla.org/en/docs/Web/SVG/Attribute/transform -applyTransformation :: Matrix Double -> SVG.Transformation -> Matrix Double -applyTransformation m (SVG.TransformMatrix a b c d e f) = multStd m (fromElements [a,b,c,d,e,f]) -applyTransformation m (SVG.Translate x y) = multStd m (fromElements [1,0,0,1,x,y]) -applyTransformation m (SVG.Scale sx mbSy) = multStd m (fromElements [sx,0,0,maybe sx id mbSy,0,0]) -applyTransformation m (SVG.Rotate a Nothing) - = multStd m (fromElements [cos(r),sin(r),-sin(r),cos(r),0,0]) - where - r = a * radiansPerDegree -applyTransformation m (SVG.Rotate a (Just (x, y))) = applyTransformations m (Just [SVG.Translate x y , SVG.Rotate a Nothing , SVG.Translate (-x) (-y)]) -applyTransformation m (SVG.SkewX a) = multStd m (fromElements [1,0,tan(a*radiansPerDegree),1,0,0]) -applyTransformation m (SVG.SkewY a) = multStd m (fromElements [1,tan(a*radiansPerDegree),0,1,0,0]) -applyTransformation m (SVG.TransformUnknown) = m diff --git a/src/Types.hs b/src/Types.hs deleted file mode 100644 --- a/src/Types.hs +++ /dev/null @@ -1,25 +0,0 @@ -module Types ( Point - , DrawOp (..) - , GCodeOp (..) - , if' - ) where - -type Point = (Double, Double) -- A point in the plane, absolute coordinates - --- all of them are invariant under affine transformation -data DrawOp = DMoveTo Point - | DLineTo Point -- End point - | DBezierTo Point Point Point -- Control point1, control point2, end point - deriving Show - --- this is basically what GCode can do -data GCodeOp = GMoveTo Point - | GLineTo Point -- End point - | GArcTo Point Point Bool -- Center point, end point, clockwise - | GBezierTo Point Point Point -- First and second control points, end point - deriving Show - --- just to make it available everywhere -if' :: Bool -> t -> t -> t -if' True t _ = t -if' False _ f = f diff --git a/src/Utils.hs b/src/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Utils.hs @@ -0,0 +1,6 @@ +module Utils ( if' ) where + +-- just to make it available everywhere +if' :: Bool -> t -> t -> t +if' True t _ = t +if' False _ f = f