arc.lisp

   1 ;;; Copyright (c) 2008 Zachary Beane, All Rights Reserved
   2 ;;;
   3 ;;; Redistribution and use in source and binary forms, with or without
   4 ;;; modification, are permitted provided that the following conditions
   5 ;;; are met:
   6 ;;;
   7 ;;;   * Redistributions of source code must retain the above copyright
   8 ;;;     notice, this list of conditions and the following disclaimer.
   9 ;;;
  10 ;;;   * Redistributions in binary form must reproduce the above
  11 ;;;     copyright notice, this list of conditions and the following
  12 ;;;     disclaimer in the documentation and/or other materials
  13 ;;;     provided with the distribution.
  14 ;;;
  15 ;;; THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESSED
  16 ;;; OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  17 ;;; WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  18 ;;; ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
  19 ;;; DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  20 ;;; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  21 ;;; GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  22 ;;; INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
  23 ;;; WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  24 ;;; NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  25 ;;; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26 ;;;
  27
  28 ;;; arc.lisp
  29
  30 (in-package #:vecto)
  31
  32 ;;; Adapted from Ben Deane's com.elbeno.curve 0.1 library, with
  33 ;;; permission. See http://www.elbeno.com/lisp/ for the original.
  34
  35 (defparameter +cubic-error-coeffs-0+
  36   (make-array '(2 4 4) :initial-contents
  37               '((( 3.85268  -21.229      -0.330434   0.0127842)
  38                  (-1.61486    0.706564    0.225945   0.263682)
  39                  (-0.910164   0.388383    0.00551445 0.00671814)
  40                  (-0.630184   0.192402    0.0098871  0.0102527))
  41                 ((-0.162211   9.94329     0.13723    0.0124084)
  42                  (-0.253135   0.00187735  0.0230286  0.01264)
  43                  (-0.0695069 -0.0437594   0.0120636  0.0163087)
  44                  (-0.0328856 -0.00926032 -0.00173573 0.00527385)))))
  45
  46 (defparameter +cubic-error-coeffs-1+
  47   (make-array '(2 4 4) :initial-contents
  48               '((( 0.0899116 -19.2349    -4.11711     0.183362)
  49                  ( 0.138148   -1.45804    1.32044     1.38474)
  50                  ( 0.230903   -0.450262   0.219963    0.414038)
  51                  ( 0.0590565  -0.101062   0.0430592   0.0204699))
  52                 (( 0.0164649   9.89394    0.0919496   0.00760802)
  53                  ( 0.0191603  -0.0322058  0.0134667  -0.0825018)
  54                  ( 0.0156192  -0.017535   0.00326508 -0.228157)
  55                  (-0.0236752   0.0405821 -0.0173086   0.176187)))))
  56
  57 ;;; [elbeno:]
  58 ;;; compute the error of a cubic bezier
  59 ;;; that approximates an elliptical arc
  60 ;;; with radii a, b
  61 ;;; between angles eta1 and eta2
  62
  63 (defun calc-c-term (i b/a etasum arr)
  64   (loop
  65      for j from 0 to 3
  66      sum (* (/ (+ (* (aref arr i j 0) b/a b/a)
  67                   (* (aref arr i j 1) b/a)
  68                   (aref arr i j 2))
  69                (+ (aref arr i j 3) b/a))
  70             (cos (* j etasum)))))
  71
  72 (defun bezier-error (a b eta1 eta2)
  73   (let* ((b/a (/ b a))
  74          (etadiff (- eta2 eta1))
  75          (etasum (+ eta2 eta1))
  76          (arr (if (< b/a 0.25)
  77                   +cubic-error-coeffs-0+
  78                   +cubic-error-coeffs-1+)))
  79     (* (/ (+ (* 0.001 b/a b/a) (* 4.98 b/a) 0.207)
  80           (+ b/a 0.0067))
  81        a
  82        (exp (+ (calc-c-term 0 b/a etasum arr)
  83                (* (calc-c-term 1 b/a etasum arr) etadiff))))))
  84
  85 (defun ellipse-val (cx cy a b theta eta)
  86   (values
  87    (+ cx
  88       (* a (cos theta) (cos eta))
  89       (* (- b) (sin theta) (sin eta)))
  90    (+ cy
  91       (* a (sin theta) (cos eta))
  92       (* b (cos theta) (sin eta)))))
  93
  94 (defun ellipse-deriv-val (a b theta eta)
  95   (values
  96    (+ (* (- a) (cos theta) (sin eta))
  97       (* (- b) (sin theta) (cos eta)))
  98    (+ (* (- a) (sin theta) (sin eta))
  99       (* b (cos theta) (cos eta)))))
 100
 101 ;;; FIXME: The original elbeno code used real abstraction to manage
 102 ;;; points and curves and splines. I ripped it out and replaced it
 103 ;;; with the following ugly mess. Fix, fix, fix. For example, it might
 104 ;;; be possible to use cl-vectors to accumulate the path instead of
 105 ;;; using lists of conses.
 106
 107 (defun approximate-arc-single (cx cy a b theta eta1 eta2)
 108   (let* ((etadiff (- eta2 eta1))
 109          (k (tan (/ etadiff 2)))
 110          (alpha (* (sin etadiff)
 111                    (/ (1- (sqrt (+ 4 (* 3 k k)))) 3)))
 112          px1 py1
 113          px2 py2
 114          qx1 qy1
 115          qx2 qy2
 116          sx1 sy1
 117          sx2 sy2)
 118     (setf (values px1 py1) (ellipse-val cx cy a b theta eta1)
 119           (values px2 py2) (ellipse-val cx cy a b theta eta2)
 120           (values sx1 sy1) (ellipse-deriv-val a b theta eta1)
 121           (values sx2 sy2) (ellipse-deriv-val a b theta eta2)
 122           qx1 (+ px1 (* alpha sx1))
 123           qy1 (+ py1 (* alpha sy1))
 124           qx2 (- px2 (* alpha sx2))
 125           qy2 (- py2 (* alpha sy2)))
 126     (list (cons px1 py1)
 127           (cons qx1 qy1)
 128           (cons qx2 qy2)
 129           (cons px2 py2))))
 130
 131 (defun approximate-arc (cx cy a b theta eta1 eta2 err)
 132   (cond ((< eta2 eta1)
 133          (error "approximate-arc: eta2 must be bigger than eta1"))
 134         ((> eta2 (+ eta1 (/ pi 2) (* eta2 long-float-epsilon)))
 135          (let ((etamid (+ eta1 (/ pi 2) (* eta2 long-float-epsilon))))
 136            (nconc
 137             (approximate-arc cx cy a b theta eta1 etamid err)
 138             (approximate-arc cx cy a b theta etamid eta2 err))))
 139         (t (if (> err (bezier-error a b eta1 eta2))
 140                (list (approximate-arc-single cx cy a b theta eta1 eta2))
 141                (let ((etamid (/ (+ eta1 eta2) 2)))
 142                  (nconc
 143                   (approximate-arc cx cy a b theta eta1 etamid err)
 144                   (approximate-arc cx cy a b theta etamid eta2 err)))))))
 145
 146 (defun approximate-elliptical-arc (cx cy a b theta eta1 eta2
 147                                    &optional (err 0.5))
 148   "Approximate an elliptical arc with a cubic bezier spline."
 149   (if (> b a)
 150       (approximate-arc cx cy b a
 151                        (+ theta (/ pi 2))
 152                        (- eta1 (/ pi 2))
 153                        (- eta2 (/ pi 2)) err)
 154       (approximate-arc cx cy a b theta eta1 eta2 err)))