Merge pull request #10 from phoe-trash/master
[vecto.git] / arc.lisp
blob6291c5db1169fe9c560d4cf1d8e0ac333fc37d40
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 ;;;
28 ;;; arc.lisp
30 (in-package #:vecto)
32 ;;; Adapted from Ben Deane's com.elbeno.curve 0.1 library, with
33 ;;; permission. See http://www.elbeno.com/lisp/ for the original.
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)))))
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)))))
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
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)))))
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))
82 (exp (+ (calc-c-term 0 b/a etasum arr)
83 (* (calc-c-term 1 b/a etasum arr) etadiff))))))
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)))))
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)))))
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.
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))))
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)))))))
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)))