/[gsharp]/gsharp/mf.lisp
ViewVC logotype

Contents of /gsharp/mf.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (show annotations)
Sat Nov 15 17:19:38 2008 UTC (5 years, 5 months ago) by rstrandh
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +4 -2 lines
Make it possible to use cycle after --.
1 (in-package :mf)
2
3 (defclass basic-path-join () ())
4
5 (defclass concatenate-path-join (basic-path-join) ())
6
7 (defparameter *concatenate-path-join* (make-instance 'concatenate-path-join))
8
9 (defclass controls (basic-path-join)
10 ((a :initarg :a)
11 (b :initarg :b)))
12
13 (defclass tensions (basic-path-join)
14 ((a :initarg :a)
15 (b :initarg :b)))
16
17 (defclass direction-specifier () ())
18
19 (defclass curl (direction-specifier)
20 ((curl :initarg :curl)))
21
22 (defclass direction (direction-specifier)
23 ((direction :initarg :direction)))
24
25 (defclass cycle () ())
26
27 (defparameter *cycle* (make-instance 'cycle))
28
29 (defun flatten (list)
30 (loop for elem in list
31 append (if (consp elem)
32 (flatten elem)
33 (list elem))))
34
35 (defclass context ()
36 ((%neighbor :initform nil :initarg :neighbor :accessor neighbor)
37 (%curl :initform nil :initarg :curl :accessor curl)
38 (%direction :initform nil :initarg :direction :accessor direction)
39 (%tension :initform nil :initarg :tension :accessor tension)
40 (%control :initform nil :initarg :control :accessor control)
41 (%angle :initform nil :initarg :angle :accessor angle)))
42
43 (defclass point ()
44 ((%point :initarg :point :reader point)
45 (%rank :initform nil :accessor rank)))
46
47 (defclass left-context-point (point)
48 ((%left-context :initform (make-instance 'context)
49 :initarg :left-context :reader left-context)))
50
51 (defclass right-context-point (point)
52 ((%right-context :initform (make-instance 'context)
53 :initarg :right-context :reader right-context)))
54
55 (defclass left-endpoint (right-context-point) ())
56
57 (defclass right-endpoint (left-context-point) ())
58
59 (defclass corner-point (left-context-point right-context-point) ())
60
61 (defclass interior-point (left-context-point right-context-point) ())
62
63 (defun remove-concatenates (path)
64 (loop until (null path)
65 if (and (not (null (cddr path))) ; path has at least 3 elements
66 (eq (second path) *concatenate-path-join*))
67 collect (prog1 (make-instance 'corner-point :point (car path))
68 (setf path (cdddr path)))
69 else collect (pop path)))
70
71 (defun check-start-end (path)
72 (assert (numberp (car path))
73 ()
74 "the path must start with a point, but ~s was found" (car path))
75 (assert (or (numberp (car (last path))) (eq (car (last path)) *cycle*))
76 ()
77 "the path must end with a point, but ~s was found" (car (last path))))
78
79 (defun check-cycle (path)
80 ;; check that there is no `cycle' other than at
81 ;; the end of the path.
82 (assert (not (member-if (lambda (x) (eq x *cycle*)) (butlast path)))
83 ()
84 "only the last element of a path can be CYCLE"))
85
86 (defun check-syntax (path)
87 (check-start-end path)
88 (check-cycle path)
89 (loop for (x y z) on path
90 ;; check that each direction specifier is surrounded
91 ;; by a point on one side and a tensions object on the other
92 ;; or possibly a cycle object on the right
93 do (when (typep y 'direction-specifier)
94 (assert (or (and (typep x 'tensions)
95 (or (numberp z) (eq z *cycle*)))
96 (and (numberp x)
97 (typep z 'tensions)))
98 ()
99 "a direction specifier must have a point on one side and a tensions object on the other, but ~a and ~a were found" x z))
100 ;; check that each controls object and each
101 ;; concatenate path join is surrounded by points
102 ;; or possibly a cycle object on the right
103 do (when (or (eq y *concatenate-path-join*)
104 (typep y 'controls))
105 (assert (numberp x)
106 ()
107 "a concatenate path join or a controls object must have a point to the left, but ~a was found" x)
108 (assert (or (numberp z) (eq z *cycle*))
109 ()
110 "a concatenate path join or a controls object must have a point or `cycle' to the right, but ~a was found" z))
111 ;; check that each tensions object is
112 ;; surrounded by a direction specifier or a point
113 ;; or possibly a cycle on the right
114 do (when (typep y 'tensions)
115 (assert (typep x '(or number direction-specifier))
116 ()
117 "a tension object must have a direction specifier or a point to the left, but ~a was found" x)
118 (assert (or (typep z '(or number direction-specifier)) (eq z *cycle*))
119 ()
120 "a tension object must have a direction specifier, a point, or `cycle' to the right, but ~a was found" z))
121 ;; check that each point or cycle object is surrounded by a
122 ;; direction specifier or a basic path join
123 do (when (or (numberp y) (eq y *cycle*))
124 (assert (typep x '(or direction-specifier basic-path-join))
125 ()
126 "a point or a cycle object must have a direction specifier or a basic path joint to the left, but ~a was found" x)
127 (assert (or (null z) (typep z '(or direction-specifier basic-path-join)))
128 ()
129 "a point must have a direction specifier or a basic path joint to the right, but ~a was found" z))))
130
131 (defun propagate-direction-specifiers (path)
132 (loop for (x y z) on path
133 do (when (typep y 'direction-specifier)
134 (if (typep x 'point)
135 (if (typep y 'curl)
136 (setf (curl (right-context x)) (slot-value y 'curl))
137 (setf (direction (right-context x)) (slot-value y 'direction)))
138 (if (typep y 'curl)
139 (setf (curl (left-context (if (eq z *cycle*) (car path) z)))
140 (slot-value y 'curl))
141 (setf (direction (left-context (if (eq z *cycle*) (car path) z)))
142 (slot-value y 'direction)))))))
143
144 (defun propagate-tensions-controls (path)
145 (loop for (x y z) on path
146 do (typecase y
147 (tensions (setf (tension (right-context x)) (slot-value y 'a))
148 (setf (tension (left-context (if (eq z *cycle*)
149 (car path)
150 z)))
151 (slot-value y 'b)))
152 (controls (setf (control (right-context x)) (slot-value y 'a))
153 (setf (control (left-context (if (eq z *cycle*)
154 (car path)
155 z)))
156 (slot-value y 'b))))))
157
158 (defun link-and-rank-points (path)
159 (when (eq (car (last path)) *cycle*)
160 (nbutlast path)
161 (setf (neighbor (right-context (car (last path)))) (car path))
162 (setf (neighbor (left-context (car path))) (car (last path))))
163 (loop for (x y) on path
164 for i from 0
165 do (setf (rank x) i)
166 until (null y)
167 do (setf (neighbor (right-context x)) y
168 (neighbor (left-context y)) x)))
169
170 (defun propagate-directions (path)
171 (flet ((possibly-fill-in-curl (context)
172 (when (and (null (control context))
173 (null (direction context))
174 (null (curl context)))
175 (setf (curl context) 1.0))))
176 (loop for point in path
177 do (typecase point
178 (left-endpoint
179 (possibly-fill-in-curl (right-context point)))
180 (right-endpoint
181 (possibly-fill-in-curl (left-context point)))
182 (corner-point
183 (possibly-fill-in-curl (right-context point))
184 (possibly-fill-in-curl (left-context point)))
185 (interior-point
186 (let ((lc (left-context point))
187 (rc (right-context point)))
188 (when (and (null (curl lc))
189 (null (direction lc))
190 (null (control lc)))
191 (cond ((not (null (curl rc)))
192 (setf (curl lc)
193 (curl rc)))
194 ((not (null (direction rc)))
195 (setf (direction lc)
196 (direction rc)))
197 ((not (null (control rc)))
198 (setf (direction lc)
199 (- (control rc) (point point))))))
200 (when (and (null (curl rc))
201 (null (direction rc))
202 (null (control rc)))
203 (cond ((not (null (curl lc)))
204 (setf (curl rc)
205 (curl lc)))
206 ((not (null (direction lc)))
207 (setf (direction rc)
208 (direction lc)))
209 ((not (null (control lc)))
210 (setf (direction rc)
211 (- (point point) (control lc))))))))))))
212
213 (defun solve (system)
214 (let* ((width (array-dimension system 1))
215 (height (array-dimension system 0))
216 (rows (loop for i from 0 below height collect i)))
217 (flet ((eliminate (rows column)
218 (flet ((eliminate-row (row1 row2)
219 (let ((factor (/ (aref system row2 column)
220 (aref system row1 column))))
221 (loop for i from column below width
222 do (decf (aref system row2 i)
223 (* factor (aref system row1 i)))))))
224 (let ((pivot-row (member-if (lambda (row)
225 (not (zerop (aref system row column))))
226 rows)))
227 (rotatef (car pivot-row) (car rows))
228 (loop for row in (cdr rows)
229 do (eliminate-row (car rows) row))))))
230 (loop for column from 0 below (- width 2)
231 for remaining-rows on rows
232 do (eliminate remaining-rows column))
233 (let ((last-column (1- width)))
234 (loop for rev-rows on (reverse rows)
235 for column downfrom (- width 2)
236 do (loop with row1 = (car rev-rows)
237 for row2 in (cdr rev-rows)
238 do (decf (aref system row2 last-column)
239 (* (/ (aref system row2 column)
240 (aref system row1 column))
241 (aref system row1 last-column)))
242 do (setf (aref system row2 column) 0.0)))
243 (loop for column from 0
244 for row in rows
245 do (setf (aref system row last-column)
246 (/ (aref system row last-column)
247 (aref system row column))))))
248 (let ((solution (make-array height)))
249 (loop for i from 0 below height
250 do (setf (aref solution i) (aref system (elt rows i) (1- width))))
251 solution)))
252
253 (defun solve-angles (path)
254 (let* ((open-p (typep (car path) 'left-endpoint))
255 (length (length path))
256 (nb-variables (- (* 2 (length path))
257 (if open-p 2 0)))
258 (width (1+ nb-variables))
259 (matrix (make-array (list nb-variables width) :initial-element 0.0))
260 (equation-number -1))
261 (labels ((out (i) (* 2 i))
262 (in (i) (1- (* 2 (if (zerop i) length i))))
263 (handle-right-context (point context)
264 (cond ((not (null (control context)))
265 (let ((out-angle (phase (/ (- (control context)
266 (point point))
267 (- (point (neighbor context))
268 (point point))))))
269 (setf (aref matrix (incf equation-number) (out (rank point)))
270 1.0)
271 (setf (aref matrix equation-number (1- width))
272 out-angle)))
273 ((not (null (direction context)))
274 (let ((out-angle (phase (/ (direction context)
275 (- (point (neighbor context))
276 (point point))))))
277 (setf (aref matrix (incf equation-number) (out (rank point)))
278 1.0)
279 (setf (aref matrix equation-number (1- width))
280 out-angle)))
281 ((not (null (curl context)))
282 (let* ((a0 (tension context))
283 (b1 (tension (left-context (neighbor context))))
284 (g0 (curl context))
285 (c1 (- (* a0 a0 a0 (- 1 (* 3.0 b1)))
286 (* g0 b1 b1 b1)))
287 (c2 (+ (* a0 a0 a0)
288 (- (* g0 b1 b1 b1))
289 (* 3.0 a0))))
290 (setf (aref matrix (incf equation-number) (out (rank point)))
291 c1)
292 (setf (aref matrix equation-number (in (rank (neighbor context))))
293 c2)))))
294 (handle-left-context (point context)
295 (cond ((not (null (control context)))
296 (let ((in-angle (phase (/ (- (point point)
297 (point (neighbor context)))
298 (- (point point)
299 (control context))))))
300 (setf (aref matrix (incf equation-number) (in (rank point)))
301 1.0)
302 (setf (aref matrix equation-number (1- width))
303 in-angle)))
304 ((not (null (direction context)))
305 (let ((in-angle (phase (/ (- (point point)
306 (point (neighbor context)))
307 (direction context)))))
308 (setf (aref matrix (incf equation-number) (in (rank point)))
309 1.0)
310 (setf (aref matrix equation-number (1- width))
311 in-angle)))
312 ((not (null (curl context)))
313 (let* ((bn (tension context))
314 (an-1 (tension (right-context (neighbor context))))
315 (gn (curl context))
316 (c1 (+ (* bn bn bn)
317 (- (* gn an-1 an-1 an-1))
318 (* 3.0 bn)))
319 (c2 (- (* bn bn bn (- 1 (* 3.0 an-1)))
320 (* gn an-1 an-1 an-1))))
321 (setf (aref matrix (incf equation-number) (out (rank (neighbor context))))
322 c1)
323 (setf (aref matrix equation-number (in (rank point)))
324 c2))))))
325 (loop for point in path
326 do (typecase point
327 (left-endpoint
328 (handle-right-context point (right-context point)))
329 (right-endpoint
330 (handle-left-context point (left-context point)))
331 (corner-point
332 (handle-right-context point (right-context point))
333 (handle-left-context point (left-context point)))
334 (interior-point
335 (let ((lc (left-context point))
336 (rc (right-context point)))
337 (if (and (null (curl lc)) (null (direction lc)) (null (control lc))
338 (null (curl rc)) (null (direction rc)) (null (control rc)))
339 (let* ((ak-1 (tension (right-context (neighbor lc))))
340 (bk (tension lc))
341 (ak (tension rc))
342 (bk+1 (tension (left-context (neighbor rc))))
343 (lk (abs (- (point point) (point (neighbor lc)))))
344 (lk+1 (abs (- (point (neighbor rc)) (point point))))
345 (c1 (* bk bk bk+1 lk))
346 (c2 (* bk bk bk+1 lk+1 (- 1.0 (* 3.0 ak-1))))
347 (c3 (- (* ak ak ak-1 lk (- 1.0 (* 3.0 bk+1)))))
348 (c4 (- (* ak ak ak-1 lk))))
349 (setf (aref matrix
350 (incf equation-number)
351 (out (rank (neighbor lc))))
352 c1)
353 (setf (aref matrix equation-number (in (rank point)))
354 c2)
355 (setf (aref matrix equation-number (out (rank point)))
356 c3)
357 (setf (aref matrix
358 equation-number
359 (in (rank (neighbor rc))))
360 c4)
361 (setf (aref matrix (incf equation-number) (out (rank point)))
362 1.0)
363 (setf (aref matrix equation-number (in (rank point)))
364 1.0)
365 (setf (aref matrix equation-number (1- width))
366 (- (phase (/ (- (point (neighbor rc)) (point point))
367 (- (point point) (point (neighbor lc))))))))
368 (progn (handle-left-context point (left-context point))
369 (handle-right-context point (right-context point))))))))
370 (let ((solution (solve matrix)))
371 (loop for point in path
372 do (when (typep point 'left-context-point)
373 (setf (angle (left-context point)) (aref solution (in (rank point)))))
374 do (when (typep point 'right-context-point)
375 (setf (angle (right-context point)) (aref solution (out (rank point))))))))))
376
377 (defun hobby (theta phi)
378 (/ (+ 2.0
379 (* #.(sqrt 2.0)
380 (- (sin theta) (* 1/16 (sin phi)))
381 (- (sin phi) (* 1/16 (sin theta)))
382 (- (cos theta) (cos phi))))
383 (* 3.0
384 (+ 1.0
385 (* #.(* 0.5 (- (sqrt 5.0) 1.0))
386 (cos theta))
387 (* #.(* 0.5 (- 3.0 (sqrt 5.0)))
388 (cos phi))))))
389
390 (defun handle-point-pair (p0 p1 tr tl theta phi)
391 (values
392 (+ p0
393 (/ (* (exp (* #c(0.0 1.0) theta))
394 (- p1 p0)
395 (hobby theta phi))
396 tr))
397 (- p1
398 (/ (* (exp (* #c(0.0 -1.0) phi))
399 (- p1 p0)
400 (hobby phi theta))
401 tl))))
402
403 (defun assign-control-points (path)
404 (loop for (p0 p1) on path
405 until (null p1)
406 do (let* ((rc (right-context p0))
407 (lc (left-context p1))
408 (theta (angle rc))
409 (phi (angle lc)))
410 (when (null (control rc))
411 (multiple-value-bind (c0 c1)
412 (handle-point-pair (point p0) (point p1)
413 (tension rc) (tension lc) theta phi)
414 (setf (control rc) c0
415 (control lc) c1)))))
416 (unless (typep (car path) 'left-endpoint)
417 (let* ((p0 (car (last path)))
418 (p1 (car path))
419 (rc (right-context p0))
420 (lc (left-context p1))
421 (theta (angle rc))
422 (phi (angle lc)))
423 (when (null (control rc))
424 (multiple-value-bind (c0 c1)
425 (handle-point-pair (point p0) (point p1)
426 (tension rc) (tension lc) theta phi)
427 (setf (control rc) c0
428 (control lc) c1))))))
429
430 (defun point-to-complex (point)
431 "convert a point to a complex number"
432 (complex (clim:point-x point) (clim:point-y point)))
433
434 (defun complex-to-point (complex)
435 "convert a complex number to a point"
436 (clim:make-point (realpart complex) (imagpart complex)))
437
438 (defun make-mf-path (&rest body)
439 (let ((path (mapcar (lambda (x)
440 (if (clim:pointp x)
441 (point-to-complex x)
442 x))
443 (flatten body))))
444 (check-syntax path)
445 ;; replace each sequence of type `p & p' by a corner point
446 (setf path (remove-concatenates path))
447 ;; replace the end points if path is not a cycle
448 (unless (eq (car (last path)) *cycle*)
449 (setf (car path) (make-instance 'left-endpoint :point (car path)))
450 (setf (car (last path)) (make-instance 'right-endpoint :point (car (last path)))))
451 ;; replace all other points by interior points
452 (setf path (loop for element in path
453 collect (if (numberp element)
454 (make-instance 'interior-point :point element)
455 element)))
456 ;; propagate direction specifiers to their respective points
457 (propagate-direction-specifiers path)
458 ;; remove all direction specifiers
459 (setf path (remove-if (lambda (x) (typep x 'direction-specifier)) path))
460 ;; propagate tensions and controls to their respective points
461 (propagate-tensions-controls path)
462 ;; remove all tensions and controls objects
463 (setf path (remove-if (lambda (x) (typep x '(or tensions controls))) path))
464 ;; link and rank the points of the path, remove the cycle object
465 (link-and-rank-points path)
466 ;; now the path contains only point objects
467 (propagate-directions path)
468 (solve-angles path)
469 (assign-control-points path)
470 (if (typep (car path) 'left-endpoint)
471 (let ((segments (loop for point in (butlast path)
472 collect (let ((rc (right-context point)))
473 (climi::make-bezier-segment
474 (complex-to-point (point point))
475 (complex-to-point (control rc))
476 (complex-to-point (control (left-context (neighbor rc))))
477 (complex-to-point (point (neighbor rc))))))))
478 (make-instance 'climi::bezier-curve :segments segments))
479 (let ((segments (loop for point in path
480 collect (let ((rc (right-context point)))
481 (climi::make-bezier-segment
482 (complex-to-point (point point))
483 (complex-to-point (control rc))
484 (complex-to-point (control (left-context (neighbor rc))))
485 (complex-to-point (point (neighbor rc))))))))
486 (make-instance 'climi::bezier-area :segments segments)))))
487
488 (defparameter *infinity* 4095.99998) ;see the MF book
489
490 (defmacro mf (&body body)
491 `(flet ((control (a)
492 (make-instance 'controls :a a :b a))
493 (controls (a b)
494 (make-instance 'controls :a a :b b))
495 (tension (a)
496 (assert (>= a 0.75)
497 ()
498 "tension values must be greater than 0.75: ~a"
499 a)
500 (make-instance 'tensions :a a :b a))
501 (tensions (a b)
502 (assert (>= (min a b) 0.75)
503 ()
504 "tension values must be greater than 0.75: ~a"
505 (min a b))
506 (make-instance 'tensions :a a :b b))
507 (direction (d)
508 (make-instance 'direction :direction d))
509 (curl (c)
510 (make-instance 'curl :curl c)))
511 (declare (ignorable (function control)
512 (function controls)
513 (function tension)
514 (function tensions)
515 (function direction)
516 (function curl)))
517 (let* ((++ (tension 1.0))
518 (+++ (tension 1.0)) ; this is not right
519 (& *concatenate-path-join*)
520 (--- (tension *infinity*))
521 (-- (list (make-instance 'curl :curl 1)
522 (make-instance 'tensions :a 2 :b 2) ; should be 1 rather than 2
523 (make-instance 'curl :curl 1)))
524 (cycle *cycle*)
525 (up (direction #c(0 1)))
526 (down (direction #c(0 -1)))
527 (left (direction #c(-1 0)))
528 (right (direction #c(1 0))))
529 (declare (ignorable +++ & --- -- cycle up down left right))
530 (make-mf-path
531 ,@body))))
532
533 (defun part-way (p0 p1 alpha)
534 (+ (* (- 1 alpha) p0) (* alpha p1)))
535
536 ;;; some standard paths
537 (defparameter +quarter-circle+
538 (let* ((a (* 0.5 (- (sqrt 2) 1)))
539 (q0 (clim:make-point 0.5 0.0))
540 (q1 (clim:make-point 0.0 0.5))
541 (p0 (clim:make-point (/ 0.5 (sqrt 2)) (/ 0.5 (sqrt 2))))
542 (p1 (clim:make-point 0.5 a))
543 (p2 (clim:make-point a 0.5))
544 (alpha 0.7))
545 (climi::make-bezier-curve
546 (list q0
547 (climi::part-way q0 p1 alpha)
548 (climi::part-way p0 p1 alpha)
549 p0
550 (climi::part-way p0 p2 alpha)
551 (climi::part-way q1 p2 alpha)
552 q1))))
553
554 (defparameter +half-circle+
555 (let* ((tr (clim:make-rotation-transformation (/ pi 2)))
556 (rotated-quarter-circle (clim:transform-region tr +quarter-circle+)))
557 (clim:region-union +quarter-circle+ rotated-quarter-circle)))
558
559 (defparameter +full-circle+
560 (let* ((tr (clim:make-rotation-transformation pi))
561 (rotated-half-circle (clim:transform-region tr +half-circle+)))
562 (climi::close-path (clim:region-union +half-circle+ rotated-half-circle))))
563
564 (defparameter +unit-square+
565 (climi::close-path
566 (mf #c(0.5 0.5) -- #c(-0.5 0.5) -- #c(-0.5 -0.5) -- #c(0.5 -0.5) -- #c(0.5 0.5))))
567
568 (defun superellipse (r top l bot superness)
569 (let ((xtr (part-way (realpart top) (realpart r) superness))
570 (yrt (part-way (imagpart r) (imagpart top) superness))
571 (xtl (part-way (realpart top) (realpart l) superness))
572 (ylt (part-way (imagpart l) (imagpart top) superness))
573 (xbl (part-way (realpart bot) (realpart l) superness))
574 (ylb (part-way (imagpart l) (imagpart bot) superness))
575 (xbr (part-way (realpart bot) (realpart r) superness))
576 (yrb (part-way (imagpart r) (imagpart bot) superness)))
577 (mf r up +++ (complex xtr yrt) (direction (- top r)) +++
578 top left +++ (complex xtl ylt) (direction (- l top)) +++
579 l down +++ (complex xbl ylb) (direction (- bot l)) +++
580 bot right +++ (complex xbr yrb) (direction (- r bot)) +++ cycle)))
581
582 (defparameter +razor+
583 (climi::close-path (mf #c(-0.5 0) -- #c(0.5 0) -- #c(-0.5 0))))
584
585 ;;; pen drawing
586
587 (defvar *pen* nil)
588
589 (defmacro with-pen (pen &body body)
590 `(let ((*pen* ,pen))
591 ,@body))
592
593 (defun draw-path (path)
594 (climi::convolve-regions *pen* path))

  ViewVC Help
Powered by ViewVC 1.1.5