File: bezier.lisp

package info (click to toggle)
cl-mcclim 0.9.6.dfsg.cvs20100315-2
  • links: PTS
  • area: main
  • in suites: jessie, jessie-kfreebsd, stretch
  • size: 6,624 kB
  • sloc: lisp: 106,015; makefile: 47; sh: 15
file content (812 lines) | stat: -rw-r--r-- 30,281 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
(in-package :clim-internals)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Utilities

(defun point-to-complex (point)
  "convert a point to a complex number"
  (complex (point-x point) (point-y point)))

(defun complex-to-point (complex)
  "convert a complex number to a point"
  (make-point (realpart complex) (imagpart complex)))

(defun distance (p0 p1)
  "return the euclidian distance between two points"
  (multiple-value-bind (x0 y0) (point-position p0)
    (multiple-value-bind (x1 y1) (point-position p1)
      (let* ((dx (- x1 x0))
	     (dx2 (* dx dx))
	     (dy (- y1 y0))
	     (dy2 (* dy dy)))
	(sqrt (+ dx2 dy2))))))

(defun part-way (p0 p1 alpha)
  "return a point that is part way between two other points"
  (multiple-value-bind (x0 y0) (point-position p0)
    (multiple-value-bind (x1 y1) (point-position p1)
      (make-point (+ (* (- 1 alpha) x0) (* alpha x1))
		  (+ (* (- 1 alpha) y0) (* alpha y1))))))

(defun dot-dist (p p0 p1)
  "dot distance between a point and a line"
  (let ((dx (- (point-x p1) (point-x p0)))
	(dy (- (point-y p1) (point-y p0))))
    (- (* (point-x p) dy)
       (* (point-y p) dx))))

(defun solve-quadratic (a2 a1 a0 &key complex-roots multiple-roots)
  (when (zerop a2)
    (return-from solve-quadratic (- (/ a0 a1))))
  (unless (= a2 1)
    (setf a1 (/ a1 a2)
	  a0 (/ a0 a2)))
  (let* ((-a1/2 (- (/ a1 2.0)))
	 (r (- (* -a1/2 -a1/2) a0)))
    (cond ((zerop r)
	   (if multiple-roots (values -a1/2 -a1/2) -a1/2))
	  ((minusp r)
	   (if complex-roots (values (+ -a1/2 (sqrt r)) (- -a1/2 (sqrt r))) (values)))
	  (t
	   (values (+ -a1/2 (sqrt r)) (- -a1/2 (sqrt r)))))))
  
(defun dist (v z)
  "compute the distance between a point and a vector represented as a complex number"
  (- (* (realpart z) (point-y v))
     (* (imagpart z) (point-x v))))

(defclass bezier-design (design) 
  ((%or :accessor original-region :initform nil)))

(defgeneric medium-draw-bezier-design* (stream design))

(defclass bezier-design-output-record (standard-graphics-displayed-output-record)
  ((stream :initarg :stream)
   (design :initarg :design)))

(defmethod initialize-instance :after ((record bezier-design-output-record) &key)
  (with-slots (design) record
    (setf (rectangle-edges* record)
	  (bounding-rectangle* design))))

(defmethod medium-draw-bezier-design* :around ((stream output-recording-stream) design)
  (with-sheet-medium (medium stream)
    (let ((transformed-design (transform-region (medium-transformation medium) design)))
      (when (stream-recording-p stream)
	(let ((record (make-instance 'bezier-design-output-record
				     :stream stream
				     :design transformed-design)))
	  (stream-add-output-record stream record)))
      (when (stream-drawing-p stream)
	(medium-draw-bezier-design* medium design)))))

(defmethod medium-draw-bezier-design* :around 
    ((medium transform-coordinates-mixin) design)
  (let* ((tr (medium-transformation medium))
         (design (transform-region tr design)))
    (call-next-method medium design)))

(defmethod replay-output-record ((record bezier-design-output-record) stream &optional
				 (region +everywhere+) (x-offset 0) (y-offset 0))
  (declare (ignore x-offset y-offset region))
  (with-slots (design) record
    (medium-draw-bezier-design* (sheet-medium stream) design)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Bezier curves and areas

(defclass bezier-segment ()
  ((p0 :initarg :p0)
   (p1 :initarg :p1)
   (p2 :initarg :p2)
   (p3 :initarg :p3)))

(defun make-bezier-segment (p0 p1 p2 p3)
  (make-instance 'bezier-segment
		 :p0 p0 :p1 p1 :p2 p2 :p3 p3))

(defclass bounding-rectangle-mixin ()
  ((min-x) (min-y) (max-x) (max-y)))

(defmethod bounding-rectangle* ((region bounding-rectangle-mixin))
  (with-slots (min-x min-y max-x max-y) region
    (values min-x min-y max-x max-y)))

(defclass segments-mixin (bounding-rectangle-mixin)
  ((%segments :initarg :segments :initform '() :reader %segments)))

(defmethod compute-bounding-rectangle* ((segments-mixin segments-mixin))
  (multiple-value-bind (final-min-x final-min-y final-max-x final-max-y)
      (segment-bounding-rectangle (car (%segments segments-mixin)))
    (loop for segment in (cdr (%segments segments-mixin))
	  do (multiple-value-bind (min-x min-y max-x max-y)
		 (segment-bounding-rectangle segment)
	       (setf final-min-x (min final-min-x min-x)
		     final-min-y (min final-min-y min-y)
		     final-max-x (max final-max-x max-x)
		     final-max-y (max final-max-y max-y))))
    (values final-min-x final-min-y final-max-x final-max-y)))

(defmethod initialize-instance :after ((region segments-mixin) &rest args)
  (declare (ignore args))
  (multiple-value-bind (computed-min-x computed-min-y computed-max-x computed-max-y)
      (compute-bounding-rectangle* region)
    (with-slots (min-x min-y max-x max-y) region
      (setf min-x computed-min-x
	    min-y computed-min-y
	    max-x computed-max-x
	    max-y computed-max-y))))

;;; a path defined as a sequence of Bezier curve segments
(defclass bezier-curve (path segments-mixin bounding-rectangle-mixin) ())

(defun make-bezier-thing (class point-seq)
  (assert (= (mod (length point-seq) 3) 1))
  (make-instance class
		 :segments (loop for (p0 p1 p2 p3) on point-seq by #'cdddr
				 until (null p1)
				 collect (make-bezier-segment p0 p1 p2 p3))))

(defun make-bezier-thing* (class coord-seq)
  (assert (= (mod (length coord-seq) 6) 2))
  (make-instance class
		 :segments (loop for (x0 y0 x1 y1 x2 y2 x3 y3 x4 y4)
				 on coord-seq by #'(lambda (x) (nthcdr 6 x))
				 until (null x1)
				 collect (make-bezier-segment
					  (make-point x0 y0)
					  (make-point x1 y1)
					  (make-point x2 y2)
					  (make-point x3 y3)))))

(defun make-bezier-curve (point-seq)
  (make-bezier-thing 'bezier-curve point-seq))

(defun make-bezier-curve* (coord-seq)
  (make-bezier-thing* 'bezier-curve coord-seq))

(defun transform-segment (transformation segment)
  (with-slots (p0 p1 p2 p3) segment
    (make-bezier-segment (transform-region transformation p0)
			 (transform-region transformation p1)
			 (transform-region transformation p2)
			 (transform-region transformation p3))))

(defmethod transform-region (transformation (path bezier-curve))
  (make-instance 'bezier-curve
		 :segments (mapcar (lambda (segment)
				     (transform-segment transformation segment))
				   (%segments path))))

(defmethod region-equal ((p1 point) (p2 point))
  (let ((coordinate-epsilon (* #.(expt 2 10) double-float-epsilon)))
    (and (<= (abs (- (point-x p1) (point-x p2))) coordinate-epsilon)
	 (<= (abs (- (point-y p1) (point-y p2))) coordinate-epsilon))))	 

(defmethod region-union ((r1 bezier-curve) (r2 bezier-curve))
  (let ((p (slot-value (car (last (%segments r1))) 'p3))
	(seg (car (%segments r2))))
    (if (region-equal p (slot-value seg 'p0))
	(with-slots (p1 p2 p3) seg
	  (make-instance 'bezier-curve
			 :segments (append (%segments r1)
					   (cons (make-bezier-segment p p1 p2 p3)
						 (cdr (%segments r2))))))
	(call-next-method))))

;;; an area defined as a closed path of Bezier curve segments
(defclass bezier-area (area bezier-design segments-mixin bounding-rectangle-mixin) 
  ((%trans :initarg :transformation :reader transformation :initform +identity-transformation+)))

(defgeneric close-path (path))

(defmethod close-path ((path bezier-curve))
  (let ((segments (%segments path)))
    (assert (region-equal (slot-value (car segments) 'p0)
			  (slot-value (car (last segments)) 'p3)))
    (make-instance 'bezier-area :segments segments)))

(defun path-start (path)
  (slot-value (car (%segments path)) 'p0))

(defun path-end (path)
  (slot-value (car (last (%segments path))) 'p3))

(defun make-bezier-area (point-seq)
  (assert (region-equal (car point-seq) (car (last point-seq))))
  (make-bezier-thing 'bezier-area point-seq))

(defun make-bezier-area* (coord-seq)
  (assert (and (coordinate= (car coord-seq) (car (last coord-seq 2)))
	       (coordinate= (cadr coord-seq) (car (last coord-seq)))))
  (make-bezier-thing* 'bezier-area coord-seq))

(defmethod segments ((area bezier-area))
  (let ((tr (transformation area)))
    (mapcar (lambda (s) (transform-segment tr s)) (%segments area))))

(defmethod transform-region (transformation (area bezier-area))
  (let* ((tr (transformation area))
         (result (if (translation-transformation-p transformation)
                     (make-instance 'bezier-area :segments (%segments area)
                                    :transformation 
                                    (compose-transformations transformation tr))
                     (make-instance 'bezier-area 
                                    :segments (mapcar (lambda (s) (transform-segment transformation s)) (segments area))))))
    (when (translation-transformation-p transformation)
      (setf (original-region result) (or (original-region area) area)))
    result))

(defmethod compute-bounding-rectangle* ((area bezier-area))
  (multiple-value-bind (lx ly ux uy) (call-next-method)
    (let ((tr (transformation area)))
      (transform-rectangle* tr lx ly ux uy))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Special cases of combined Bezier areas

;;; A union of bezier areas.  This is not itself a bezier area.
(defclass bezier-union (area bezier-design)
  ((%trans :initarg :transformation :reader transformation :initform +identity-transformation+)
   (%areas :initarg :areas :initform '() :reader areas)))

(defmethod transform-region (transformation (union bezier-union))
  (let* ((tr (transformation union))
         (new-tr (compose-transformations transformation tr))
         (result (if (translation-transformation-p transformation)
                     (make-instance 'bezier-union :areas (areas union)
                                    :transformation new-tr)
                     (make-instance 'bezier-union
                                    :areas (loop for area in (areas union) collect (transform-region new-tr area))))))
    (when (translation-transformation-p transformation)
      (setf (original-region result) (or (original-region union) union)))
    result))

(defun bounding-rectangle-of-areas (areas)
  (multiple-value-bind (final-min-x final-min-y final-max-x final-max-y)
      (bounding-rectangle* (car areas))
    (loop for area in (cdr areas)
	  do (multiple-value-bind (min-x min-y max-x max-y)
		 (bounding-rectangle* area)
	       (setf final-min-x (min final-min-x min-x)
		     final-min-y (min final-min-y min-y)
		     final-max-x (max final-max-x max-x)
		     final-max-y (max final-max-y max-y))))
    (values final-min-x final-min-y final-max-x final-max-y)))

(defmethod bounding-rectangle* ((design bezier-union))
  (multiple-value-bind (lx ly ux uy)
      (bounding-rectangle-of-areas (areas design))
    (transform-rectangle* (transformation design) lx ly ux uy)))

(defmethod region-union ((r1 bezier-area) (r2 bezier-area))
  (make-instance 'bezier-union :areas (list r1 r2)))

(defmethod region-union ((r1 bezier-union) (r2 bezier-area))
  (let ((tr (transformation r1)))
    (make-instance 'bezier-union 
                   :areas (cons (untransform-region tr r2) (areas r1))
                   :transformation tr)))

(defmethod region-union ((r1 bezier-area) (r2 bezier-union))
  (let ((tr (transformation r2)))
    (make-instance 'bezier-union 
                   :areas (cons (untransform-region tr r1) (areas r2))
                   :transformation tr)))

(defmethod region-union ((r1 bezier-union) (r2 bezier-union))
  (let ((tr1 (transformation r1))
        (tr2 (transformation r2)))
    (if (transformation-equal tr1 tr2)
        (make-instance 'bezier-union 
                       :areas (append (areas r1) (areas r2))
                       :transformation tr1)
        (let ((len1 (length (areas r1)))
              (len2 (length (areas r2))))
          (if (> len2 len1)
              (make-instance 'bezier-union
                             :areas (append (mapcar (lambda (r) (untransform-region tr2 (transform-region tr1 r))) (areas r1)) (areas r2))
                             :transformation tr2)
              (make-instance 'bezier-union
                             :areas (append (mapcar (lambda (r) (untransform-region tr1 (transform-region tr2 r))) (areas r2)) (areas r1))
                             :transformation tr1))))))

(defclass bezier-difference (area bezier-design)
  ((%positive-areas :initarg :positive-areas :initform '() :reader positive-areas)
   (%negative-areas :initarg :negative-areas :initform '() :reader negative-areas)))

(defmethod transform-region (transformation (area bezier-difference))
  (let* ((pareas (loop for area in (positive-areas area)
                       collect (transform-region transformation area)))
         (nareas (loop for area in (negative-areas area)
                       collect (transform-region transformation area)))
         (result (make-instance 'bezier-difference
                                :positive-areas pareas
                                :negative-areas nareas)))
    (when (translation-transformation-p transformation)
      (setf (original-region result) (or (original-region area) area)))
    result))

(defmethod bounding-rectangle* ((design bezier-difference))
  (bounding-rectangle-of-areas (positive-areas design)))

(defmethod region-difference ((r1 bezier-area) (r2 bezier-area))
  (make-instance 'bezier-difference
		 :positive-areas (list r1)
		 :negative-areas (list r2)))

(defmethod region-difference ((r1 bezier-area) (r2 bezier-union))
  (let ((tr (transformation r2)))
    (make-instance 'bezier-difference
                   :positive-areas (list r1)
                   :negative-areas (mapcar (lambda (r) (transform-region tr r)) (areas r2)))))

(defmethod region-difference ((r1 bezier-union) (r2 bezier-area))
  (let ((tr (transformation r1)))
    (make-instance 'bezier-difference
                   :positive-areas (mapcar (lambda (r) (transform-region tr r)) (areas r1))
                   :negative-areas (list r2))))

(defmethod region-difference ((r1 bezier-union) (r2 bezier-union))
  (let ((tr1 (transformation r1))
        (tr2 (transformation r2)))
    (make-instance 'bezier-difference
                   :positive-areas (mapcar (lambda (r) (transform-region tr1 r)) (areas r1))
                   :negative-areas (mapcar (lambda (r) (transform-region tr2 r)) (areas r2)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Converting a path to a polyline or an area to a polygon

;;; convert a cubic bezier segment to a list of 
;;; line segments
(defun %polygonalize (p0 p1 p2 p3 &key (precision 0.01))
  (if (< (- (+ (distance p0 p1)
	       (distance p1 p2)
	       (distance p2 p3))
	    (distance p0 p3))
	 precision)
      (list p3)
      (let* ((p01 (part-way p0 p1 0.5))
	     (p12 (part-way p1 p2 0.5))
	     (p23 (part-way p2 p3 0.5))
	     (p012 (part-way p01 p12 0.5))
	     (p123 (part-way p12 p23 0.5))
	     (p0123 (part-way p012 p123 0.5)))
	(nconc (%polygonalize p0 p01 p012 p0123 :precision precision)
	       (%polygonalize p0123 p123 p23 p3 :precision precision)))))

(defgeneric polygonalize (thing))

(defmethod polygonalize ((segment bezier-segment))
  (with-slots (p0 p1 p2 p3) segment
    (%polygonalize p0 p1 p2 p3)))

(defmethod polygonalize ((path bezier-curve))
  (let ((segments (%segments path)))
    (make-polyline
     (cons (slot-value (car segments) 'p0)
	   (mapcan #'polygonalize segments)))))

(defmethod polygonalize ((area bezier-area))
  (let ((segments (segments area)))
    (make-polygon (mapcan #'polygonalize segments))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Reversing a path

(defgeneric reverse-path (path))

(defun reverse-segment (bezier-segment)
  (with-slots (p0 p1 p2 p3) bezier-segment
    (make-bezier-segment p3 p2 p1 p0)))

(defmethod reverse-path ((path bezier-curve))
  (make-instance 'bezier-curve
		 :segments (reverse (mapcar #'reverse-segment (%segments path)))))

(defmethod reverse-path ((path bezier-area))
  (make-instance 'bezier-area
		 :segments (reverse (mapcar #'reverse-segment (%segments path)))
                 :transformation (transformation path)))

;;; slanting transformation are used by Metafont
(defun make-slanting-transformation (slant)
  (make-transformation 1.0 slant 0.0 1.0 0.0 0.0))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Bounding rectangle

(defun evaluate-bezier (w0 w1 w2 w3 a)
  (let ((1-a (- 1.0 a)))
    (+ (* 1-a 1-a 1-a w0)
       (* 3.0 1-a 1-a a w1)
       (* 3.0 1-a a a w2)
       (* a a a w3))))

(defun dimension-min-max (w0 w1 w2 w3)
  (when (> w0 w3)
    (rotatef w0 w3)
    (rotatef w1 w2))
  (when (and (<= w0 w1 w3)
	     (<= w0 w2 w3))
    (return-from dimension-min-max
      (values w0 w3)))
  (let ((a (+ (- w0) (* 3 w1) (* -3 w2) w3))
	(b (+ (* 2 w0) (* -4 w1) (* 2 w2)))
	(c (- w1 w0)))
    (if (zerop a)
	(if (zerop b)
	    (values w0 w3)
	    (let ((candidate (/ (- c) b)))
	      (if (or (<= candidate 0.0)
		      (>= candidate 1.0))
		  (values w0 w3)
		  (let ((w (evaluate-bezier w0 w1 w2 w3 candidate)))
		    (values (min w w0) (max w w3))))))
	(multiple-value-bind (candidate0 candidate1)
	    (solve-quadratic a b c :multiple-roots t)
	  (if (null candidate0)
	      (values w0 w3)
	      (let ((wa (evaluate-bezier w0 w1 w2 w3 candidate0))
		    (wb (evaluate-bezier w0 w1 w2 w3 candidate1)))
		(if (or (<= candidate0 0.0) (>= candidate0 1.0))
		    (if (or (<= candidate1 0.0) (>= candidate1 1.0))
			(values w0 w3)
			(values (min wb w0) (max wb w3)))
		    (if (or (<= candidate1 0.0) (>= candidate1 1.0))
			(values (min wa w0) (max wa w3))
			(values (min wa wb w0) (max wa wb w3))))))))))

(defun segment-bounding-rectangle (segment)
  (with-slots (p0 p1 p2 p3) segment
    (let ((x0 (point-x p0))
	  (x1 (point-x p1))
	  (x2 (point-x p2))
	  (x3 (point-x p3))
	  (y0 (point-y p0))
	  (y1 (point-y p1))
	  (y2 (point-y p2))
	  (y3 (point-y p3)))
    (multiple-value-bind (min-x max-x)
	(dimension-min-max x0 x1 x2 x3)
      (multiple-value-bind (min-y max-y)
	  (dimension-min-max y0 y1 y2 y3)
	(values min-x min-y max-x max-y))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Convolution

(defun find-split-points-for-side (aa bb cc)
  (let ((roots '()))
    (multiple-value-bind (r1 r2)
	(solve-quadratic aa bb cc)
      (unless (or (null r1) (<= r1 0.0) (>= r1 1.0)) (push r1 roots))
      (unless (or (null r2) (<= r2 0.0) (>= r2 1.0)) (push r2 roots))
      roots)))

(defun find-split-points (sides segment)
  (let ((split-points '()))
    (with-slots (p0 p1 p2 p3) segment
      (let ((x0 (point-x p0)) (y0 (point-y p0))
	    (x1 (point-x p1)) (y1 (point-y p1))
	    (x2 (point-x p2)) (y2 (point-y p2))
	    (x3 (point-x p3)) (y3 (point-y p3)))
	(let ((xa (+ (- x0) (* 3 x1) (* -3 x2) x3))
	      (ya (+ (- y0) (* 3 y1) (* -3 y2) y3))
	      (xb (* 2 (+ x0 (* -2 x1) x2)))
	      (yb (* 2 (+ y0 (* -2 y1) y2)))
	      (xc (- x1 x0))
	      (yc (- y1 y0)))
	  (loop for side in sides
		do (let* ((sr (realpart side))
			  (si (imagpart side))
			  (aa (- (* xa si)
				 (* ya sr)))
			  (bb (- (* xb si)
				 (* yb sr)))
			  (cc (- (* xc si)
				 (* yc sr))))
		     (setf split-points
			   (append (find-split-points-for-side aa bb cc) split-points))))))
      (sort (remove-duplicates split-points) #'<))))

(defun split-segment (segment split-points)
  (if (null split-points)
      (list segment)
      (with-slots (p0 p1 p2 p3) segment
	(let* ((n (floor (length split-points) 2))
	       (pivot (nth n split-points))
	       (left (mapcar (lambda (x) (/ x pivot))
			     (subseq split-points 0 n)))
	       (right (mapcar (lambda (x) (/ (- x pivot) (- 1.0 pivot)))
			      (subseq split-points (1+ n))))
	       (p01 (part-way p0 p1 pivot))
	       (p12 (part-way p1 p2 pivot))
	       (p23 (part-way p2 p3 pivot))
	       (p012 (part-way p01 p12 pivot))
	       (p123 (part-way p12 p23 pivot))
	       (p0123 (part-way p012 p123 pivot)))
	  (append (split-segment (make-bezier-segment p0 p01 p012 p0123) left)
		  (split-segment (make-bezier-segment p0123 p123 p23 p3) right))))))

(defun mid-derivative (p0 p1 p2 p3)
  (setf p0 (point-to-complex p0)
	p1 (point-to-complex p1)
	p2 (point-to-complex p2)
	p3 (point-to-complex p3))
  (let ((a 0.5))
    (+ (* a a (+ (- p0) (* 3 p1) (* -3 p2) p3))
       (* 2 a (+ p0 (* -2 p1) p2))
       (- p1 p0))))

(defun make-line-segment (p0 p1)
  (make-bezier-segment p0 (part-way p0 p1 1/3) (part-way p0 p1 2/3) p1))

(defun add-points (p0 p1)
  (make-point (+ (point-x p0) (point-x p1)) (+ (point-y p0) (point-y p1))))

(defun convert-primitive-segment-to-bezier-area (polygon segment)
  (with-slots (p0 p1 p2 p3) segment
    (let* ((m (mid-derivative p0 p1 p2 p3))
	   (right (reduce (lambda (a b) (if (> (dist a m) (dist b m)) a b))
			  polygon))
	   (left (reduce (lambda (a b) (if (< (dist a m) (dist b m)) a b))
			 polygon)))
      (make-instance 'bezier-area
        :segments 
	(list (make-bezier-segment (add-points p0 right) (add-points p1 right)
				   (add-points p2 right) (add-points p3 right))
	      (make-line-segment (add-points p3 right) (add-points p3 left))
	      (make-bezier-segment (add-points p3 left) (add-points p2 left)
				   (add-points p1 left) (add-points p0 left))
	      (make-line-segment (add-points p0 left) (add-points p0 right)))))))

(defun area-at-point (area point)
  (let ((transformation 
	 (make-translation-transformation (point-x point) (point-y point))))
    (transform-region transformation area)))

(defun convolve-polygon-and-segment (area polygon segment first)
  (declare (optimize debug))
  (let* ((points (polygon-points polygon))
	 (sides (loop for (p0 p1) on (append (last points) points)
		      until (null p1)
		      collect (- (point-to-complex p1) (point-to-complex p0))))
	 (split-points (find-split-points sides segment))
	 (segments (split-segment segment split-points)))
    (loop for segment in segments 
	  if first collect (area-at-point area (slot-value segment 'p0))
	  collect (convert-primitive-segment-to-bezier-area 
		   (polygon-points polygon) segment)
	  collect (area-at-point area (slot-value segment 'p3)))))

(defgeneric convolve-regions (area path))

(defmethod convolve-regions ((area bezier-area) (path bezier-curve))
  (let ((polygon (polygonalize area)))
    (make-instance 
     'bezier-union :areas 
     (loop for segment in (%segments path)
	   for first = t then nil
	   append (convolve-polygon-and-segment area polygon segment first)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Rendering

(defclass scanlines ()
  ((%first-line :initform 0 :accessor first-line)
   (%chain :initform (make-instance 'flexichain:standard-flexichain) :reader chain)))

(defun nb-lines (lines)
  (flexichain:nb-elements (chain lines)))

(defun crossings (lines i)
  (flexichain:element* (chain lines) (- i (first-line lines))))

(defun line-number-to-index (lines line-number)
  (let* ((chain (chain lines))
	 (size (flexichain:nb-elements chain)))
    ;; make sure there is an element corresponding to the line number
    (cond ((zerop size)
	   (flexichain:insert* chain 0 '())
	   (setf (first-line lines) line-number))
	  ((< line-number (first-line lines))
	   (loop for i from line-number below (first-line lines)
		 do (flexichain:insert* chain 0 '()))
	   (setf (first-line lines) line-number))
	  ((>= line-number (+ (first-line lines) size))
	   (loop for i from (+ (first-line lines) size) to line-number
		 do (flexichain:insert* chain size '()))))
    (- line-number (first-line lines))))

;;; insert a single crossing into LINES
(defun insert-crossing (lines line-number x inverse-p)
  (let ((chain (chain lines))
	(index (line-number-to-index lines line-number)))
    (setf (flexichain:element* chain index)
	  (merge 'list
		 (flexichain:element* chain index)
		 (list (cons x inverse-p)) #'< :key #'car))))

;;; compute the crossings of a line segment and insert
;;; them into LINES
(defun compute-crossings (lines p0 p1)
  (let ((inverse-p nil))
    (when (< (point-y p1) (point-y p0))
      (rotatef p0 p1)
      (setf inverse-p t))
    (let ((x0 (point-x p0)) (y0 (point-y p0))
	  (x1 (point-x p1)) (y1 (point-y p1)))
      (loop for y from (round y0) below (round y1)
	    for x = (+ x0 (* (- x1 x0) (/ (- (+ y 0.5) y0) (- y1 y0))))
	    do (insert-crossing lines y x inverse-p)))))

(defun scan-lines (polygon)
  (let ((lines (make-instance 'scanlines))
	(points (polygon-points polygon)))
    (loop for (p0 p1) on (append (last points) points)
	  until (null p1)
	  do (compute-crossings lines p0 p1))
    lines))

(defun render-scan-lines (array pixel-value line crossings min-x min-y)
  (let ((level 0)
	(start nil)
	(height (array-dimension array 0))
	(width (array-dimension array 1)))
    (loop for (x . inverse-p) in crossings
	  do (when (zerop level)
	       (setf start x))
	  do (setf level (if inverse-p (1+ level) (1- level)))
	  do (when (zerop level)
	       (loop for c from (round start) below (round x)
		     do (when (and (<= 0 (round (- line min-y)) (1- height))
				   (<= 0 (- c min-x) (1- width)))
			  (setf (aref array (round (- line min-y)) (- c min-x))
				pixel-value)))))))

(defun render-polygon (array polygon pixel-value min-x min-y)
  (let ((lines (scan-lines polygon)))
    (loop for i from (first-line lines)
	  repeat (nb-lines lines)
	  do (render-scan-lines array pixel-value i (crossings lines i) min-x min-y))))

(defgeneric positive-negative-areas (design))

(defmethod positive-negative-areas ((design bezier-area))
  (values (list design) '()))

(defmethod positive-negative-areas ((design bezier-union))
  (values (areas design) '()))

(defmethod positive-negative-areas ((design bezier-difference))
  (values (positive-areas design) (negative-areas design)))

(defun render-to-array (design)
  (multiple-value-bind (positive-areas negative-areas)
      (positive-negative-areas design)
    (multiple-value-bind (min-x min-y max-x max-y)
	(bounding-rectangle-of-areas positive-areas)
      (setf min-x (* 4 (floor min-x))
	    min-y (* 4 (floor min-y))
	    max-x (* 4 (ceiling max-x))
	    max-y (* 4 (ceiling max-y)))
      (let ((result (make-array (list (- max-y min-y) (- max-x min-x))
				:element-type 'bit :initial-element 1))
	    (transformation (make-scaling-transformation* 4 4)))
	(loop for area in positive-areas
	      do (let* ((transformed-area (transform-region transformation area))
			(polygon (polygonalize transformed-area)))
		   (render-polygon result polygon 0 min-x min-y)))
	(loop for area in negative-areas
	      do (let* ((transformed-area (transform-region transformation area))
			(polygon (polygonalize transformed-area)))
		   (render-polygon result polygon 1 min-x min-y)))
	result))))

(defparameter *pixmaps* (make-hash-table :test #'equal))

(defun resolve-ink (medium)
  (if (eq (medium-ink medium) +foreground-ink+)
      (medium-foreground medium)
      (medium-ink medium)))

(defun make-ink (medium transparency)
  (let* ((a (/ transparency 16.0))
	 (1-a (- 1.0 a)))
    (multiple-value-bind (r g b) (color-rgb (resolve-ink medium))
      (make-rgb-color (+ (* a 1.0) (* 1-a r))
		      (+ (* a 1.0) (* 1-a g))
		      (+ (* a 1.0) (* 1-a b))))))

(defgeneric ensure-pixmap (medium design))

(defmethod ensure-pixmap (medium rdesign)
  (let* ((design (or (original-region rdesign) rdesign))
         (pixmap (gethash (list (medium-sheet medium) (resolve-ink medium) design)
			 *pixmaps*)))
    (when (null pixmap)
      (let* ((picture (render-to-array design))
	     (height (array-dimension picture 0))
	     (width (array-dimension picture 1))
	     (reduced-picture (make-array (list (/ height 4) (/ width 4)) :initial-element 16)))
	(loop for l from 0 below height
	      do (loop for c from 0 below width
		       do (when (zerop (aref picture l c))
			    (decf (aref reduced-picture (floor l 4) (floor c 4))))))
	(setf pixmap
	      (with-output-to-pixmap (pixmap-medium
				      (medium-sheet medium)
				      :width (/ width 4) :height (/ height 4))
		(loop for l from 0 below (/ height 4)
		      do (loop for c from 0 below (/ width 4)
			       do (draw-point*
				   pixmap-medium c l
				   :ink (make-ink
					 medium
					 (aref reduced-picture l c)))))))
	(setf (gethash (list (medium-sheet medium) (resolve-ink medium) design)
		       *pixmaps*)
	      pixmap)))
    pixmap))

(defun render-through-pixmap (design medium)
  (multiple-value-bind (min-x min-y)
      (bounding-rectangle* design)
    ;; the design we've got has already been transformed by the
    ;; medium/user transformation, and COPY-FROM-PIXMAP is in user
    ;; coordinates.  So we need to transform back (or set the medium's
    ;; transformation to be +IDENTITY-TRANSFORMATION+ temporarily, but
    ;; that's even uglier)
    (multiple-value-bind (utmin-x utmin-y)
        (untransform-position (medium-transformation medium) min-x min-y)
      (setf min-x (floor utmin-x)
            min-y (floor utmin-y))
      (let ((pixmap (ensure-pixmap medium design)))
        (copy-from-pixmap pixmap 0 0 (pixmap-width pixmap) (pixmap-height pixmap)
                          (medium-sheet medium) min-x min-y)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Generic drawing

(defun draw-bezier-design* (sheet design &rest options)
  (climi::with-medium-options (sheet options)
    (medium-draw-bezier-design* sheet design)))

(defmethod draw-design (medium (design bezier-design) &rest options &key &allow-other-keys)
  (apply #'draw-bezier-design* medium design options))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Drawing bezier designs to screen


;;; Fallback method (suitable for CLX)

(defmethod medium-draw-bezier-design* (medium design)
  (render-through-pixmap design medium))

#|
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Test cases

(defparameter *r1* (make-bezier-area* '(10 10 20 20 30 20 40 10 30 5 20 5 10 10)))

(defparameter *r2* (make-bezier-area* '(15 10 20 12 30 15 35 10 30 8 20 8 15 10)))

(defparameter *r3* (region-difference *r1* *r2*))

(defparameter *r4* (make-bezier-curve* '(100 100 120 150 160 160 170 160)))

(defparameter *r5* (convolute-regions *r2* *r4*))
|#