Added various 2D geometry routines and a small McCLIM app to test a proportion of them.

(defsystem :geometry

  :depends-on ()

  :serial t

  :components ((:file "geometry")))

;; References: http://local.wasp.uwa.edu.au/~pbourke/geometry/

(defpackage :geometry

  (:use :cl)

  (:export

   ;; Dimensions

   :lerp

   :lerp2

   :clamp

   ;:square

   ;:determinant

   ;; Points

   :distance

   :distance-squared

   ;; Lines/Vectors

   :scalar-product

   :point-is-on-line-p

   :point-is-on-line-segment-p

   :lines-are-parallel-p

   :lines-are-identical-p

   :line-intersection

   :segment-intersection

   ;; Circles

   :line-circle-intersection

   :circle-circle-intersection

   ))

(in-package :geometry)

;; Dimensions

(defun lerp (a b u)

  "Linearly Interpolate from A to B by U.

The return value when U = 0 is A, U = 0.5 is halfway from A to B, U = 1 is B."

  (+ a (* u (- b a))))

(defun lerp2 (a1 b1 a2 b2 u)

  "Linearly Interpolate from (A1,A1) to (B1,B2) by U."

  (values (lerp a1 a2 u)

	  (lerp b1 b2 u)))

(defun clamp (min x max)

  "Clamp x between min and max."

  (cond ((< x min) min)

	((> x max) max)

	(t x)))

(defun square (x) (* x x))

;;; Compute the determinant

;;; |a b|

;;; |   | = ad - bc

;;; |c d|

(defun determinant (a b c d)

  (- (* a d) (* b c)))

;; Points

(defun distance (x1 y1 x2 y2)

  "Return the straight line distance between two points, uses Pythagoras."

  (sqrt (+ (square (- x1 x2))

	   (square (- y1 y2)))))

(defun distance-squared (x1 y1 x2 y2)

  "Return the straight line distance between two points squared, uses Pythagoras."

  (+ (square (- x1 x2))

     (square (- y1 y2))))

;; Lines/Vectors

(defun scalar-product (x1 y1 x2 y2)

  "return the scalar product beween two vectors"

  (+ (* x1 x2) (* y1 y2)))

(defun point-is-on-line-p (x1 y1 x2 y2 xp yp)

  "return true iff (xp,yp) is on the line through (x1,y1) and (x2,y2)"

  (= (* (- xp x1) (- y2 y1))

     (* (- yp y1) (- x2 x1))))

(defun point-is-on-line-segment-p (x1 y1 x2 y2 xp yp)

  "return true iff (xp,yp) is on the line segment between (x1,y1) and (x2,y2)"

  (and (point-is-on-line-p x1 y1 x2 y2 xp yp)

       (<= x1 xp x2)

       (<= y1 yp y2)))

(defun lines-are-parallel-p (x11 y11 x12 y12 x21 y21 x22 y22)

  "return true iff the lines are parallel"

  (= (* (- x12 x11) (- y22 y21))

     (* (- x22 x21) (- y12 y11))))

(defun lines-are-identical-p (x11 y11 x12 y12 x21 y21 x22 y22)

  (and (lines-are-parallel-p x11 y11 x12 y12 x21 y21 x22 y22)

       (point-is-on-line-p x11 y11 x12 y12 x21 y21)))

;;; Return the point of intersection of two lines, one through

;;; (x11,y11) and (x12,y12) and the other one through (x21,y21)

;;; and (x22,y22), as two values, x and y.

;;;

;;; If the lines are the same, return T instead.

;;;

;;; If the lines are parallel but not the same, return NIL instead.

;;;

;;; Mathworld says that to find the intersection of two lines,

;;; one going through (x11,y11) and (x12,y12) and the other one

;;; going through (x21,y21) and (x22,y22) you compute:

;;;

;;;

;;;     || x11 y11 |             |

;;;     ||         |  x11 - x12  |

;;;     || x12 y12 |             |

;;;     |                        |

;;;     || x21 y21 |             |

;;;     ||         |  x21 - x22  |

;;;     || x22 y22 |             |

;;; x = --------------------------

;;;      | x11 - x12  y11 - y12 |

;;;      |                      |

;;;      | x21 - x22  y21 - y22 |

;;;  

;;;

;;;     || x11 y11 |             |

;;;     ||         |  y11 - y12  |

;;;     || x12 y12 |             |

;;;     |                        |

;;;     || x21 y21 |             |

;;;     ||         |  y21 - y22  |

;;;     || x22 y22 |             |

;;; y = --------------------------

;;;      | x11 - x12  y11 - y12 |

;;;      |                      |

;;;      | x21 - x22  y21 - y22 |

;;;  

;;;

;;; This will of course fail when

;;;

;;;      | x11 - x12  y11 - y12 |

;;;      |                      | = 0

;;;      | x21 - x22  y21 - y22 |

;;;

;;; If we assume that we are not given degenerate lines with

;;; two identical points, this happens when the lines are parallel.

(defun  line-intersection (x11 y11 x12 y12 x21 y21 x22 y22)

  (cond ((lines-are-identical-p x11 y11 x12 y12 x21 y21 x22 y22) t)

        ((lines-are-parallel-p x11 y11 x12 y12 x21 y21 x22 y22) nil)

        (t (let* ((dx1 (- x11 x12))

                  (dx2 (- x21 x22))

                  (dy1 (- y11 y12))

                  (dy2 (- y21 y22))

                  (a (determinant x11 y11 x12 y12))

                  (b (determinant x21 y21 x22 y22))

                  (c (determinant dx1 dy1 dx2 dy2)))

             (values (/ (determinant a dx1 b dx2) c)

                     (/ (determinant a dy1 b dy2) c))))))

;;; Return the point of intersection of two line segments,

;;; (x11,y11) -- (x12,y12) and (x21,y21) -- (x22,y22),

;;; as two values, x and y.

;;;

;;; If there are several points of intersection, return T instead.

;;;

;;; If there are no points of intersection, return NIL instead.

;;;

(defun segment-intersection (x11 y11 x12 y12 x21 y21 x22 y22)

  (multiple-value-bind (x y)

      (line-intersection x11 y11 x12 y12 x21 y21 x22 y22)

    (case x

      ((nil) nil)

      ((t) (cond ((and (= x11 x21)

                       (= y11 y21)

                       (not (point-is-on-line-segment-p x11 y11 x12 y12 x22 y22))

                       (not (point-is-on-line-segment-p x21 y21 x22 y22 x12 y12)))

                  (values x11 y11))

                 ((and (= x11 x22)

                       (= y11 y22)

                       (not (point-is-on-line-segment-p x11 y11 x12 y12 x21 y21))

                       (not (point-is-on-line-segment-p x21 y21 x22 y22 x12 y12)))

                  (values x11 y11))

                 ((and (= x12 x21)

                       (= y12 y21)

                       (not (point-is-on-line-segment-p x11 y11 x12 y12 x22 y22))

                       (not (point-is-on-line-segment-p x21 y21 x22 y22 x11 y11)))

                  (values x12 y12))

                 ((and (= x12 x22)

                       (= y12 y22)

                       (not (point-is-on-line-segment-p x11 y11 x12 y12 x21 y21))

                       (not (point-is-on-line-segment-p x21 y21 x22 y22 x11 y11)))

                  (values x12 y12))

                 (t t)))

      (t (cond ((and (or (<= x11 x x12)

                         (<= x12 x x11))

                     (or (<= x21 x x22)

                         (<= x22 x x21))

                     (or (<= y11 y y12)

                         (<= y12 y y11))

                     (or (<= y21 y y22)

                         (<= y22 y y21)))

                (values x y))

               (t nil))))))

;; Circles

(defun line-circle-intersection (x1 y1 x2 y2 cx cy r)

  "Return either zero one or two points of intersection between the line (x1,y1) (x2,y2) and the circle centered at (cx,cy) with radius r"

  (let* ((a (+ (expt (- x2 x1) 2) (expt (- y2 y1) 2)))

	 (b (* 2 (+ (* (- x2 x1) (- x1 cx))

		    (* (- y2 y1) (- y1 cy)))))

	 (c (- (+ (expt cx 2) (expt cy 2)

		  (expt x1 2) (expt y1 2))

	       (* 2 (+ (* cx x1)

		       (* cy y1)))

	       (* r r)))

	 (discriminant (- (* b b) (* 4 a c))))

    (cond ((= discriminant 0)

	   (multiple-value-call #'lerp2 x1 y1 x2 y2

				(/ (- b)

				   2 a)))

	  ((> discriminant 0)

	   (let* ((d-sqrt (sqrt discriminant))

		  (u1 (/ (+ (- b) d-sqrt) 2 a))

		  (u2 (/ (- (- b) d-sqrt) 2 a)))

	     (multiple-value-bind (ix1 iy1)

		 (lerp2 x1 y1 x2 y2 u1)

	       (multiple-value-bind (ix2 iy2)

		   (lerp2 x1 y1 x2 y2 u2)

		 (values ix1 iy1 ix2 iy2))))))))

(defun circle-circle-intersection (x1 y1 r1 x2 y2 r2)

  "Return the two points of intersection of the two circles circumferences, or nil"

  (let ((d (distance x1 y1 x2 y2)))

    (if (<= (abs (- r1 r2)) d (+ r1 r2))

	(let* ((a (/ (+ (- (* r1 r1) (* r2 r2)) (* d d)) 2 d))

	       (h (sqrt (- (* r1 r1) (* a a))))

	       (x3 (+ x1 (/ (* a (- x2 x1)) d)))

	       (y3 (+ y1 (/ (* a (- y2 y1)) d))))

	  (values

	   (+ x3 (/ (* h (- y2 y1)) d))

	   (- y3 (/ (* h (- x2 x1)) d))

	   (- x3 (/ (* h (- y2 y1)) d))

	   (+ y3 (/ (* h (- x2 x1)) d))))

	(values))))

(defsystem :mcclim-geom-tests

  :depends-on (:geometry :mcclim)

  :serial t

  :components ((:file "mcclim-geom-tests")))

(defpackage :mcclim-geom-tests

  (:use :geometry :clim :clim-lisp))

(in-package :mcclim-geom-tests)

(defun geom-test-1 ()

  ;; line intersection

  (let ((frame (make-application-frame 'geom-test))

	p1 p2 p3 p4

	l1 l2

	i1)

    (push (setf p1 (make-instance 'geom-point-x-y :x 50         :y 50)) (construction-of frame))

    (push (setf p2 (make-instance 'geom-point-x-y :x (- 256 50) :y 50)) (construction-of frame))

    (push (setf p3 (make-instance 'geom-point-x-y :x (- 256 50) :y (- 256 50))) (construction-of frame))

    (push (setf p4 (make-instance 'geom-point-x-y :x 50         :y (- 256 50))) (construction-of frame))

    (push (setf l1 (make-instance 'geom-line-a-b :a p1 :b p3)) (construction-of frame))

    (push (setf l2 (make-instance 'geom-line-a-b :a p2 :b p4)) (construction-of frame))

    (push (setf i1 (make-instance 'geom-point-line-line-intersection :a l1 :b l2))

	  (construction-of frame))

    (run-frame-top-level frame)))

(defun geom-test-2 ()

  ;; circle and line intersection

  (let ((frame (make-application-frame 'geom-test))

	p1 p2 p3 p4

	c1

	l1

	i1 i2)

    (push (setf p1 (make-instance 'geom-point-x-y :x 128 :y 128)) (construction-of frame))

    (push (setf p2 (make-instance 'geom-point-x-y :x 128 :y (+ 128 60))) (construction-of frame))

    (push (setf p3 (make-instance 'geom-point-x-y :x 80 :y 80)) (construction-of frame))

    (push (setf p4 (make-instance 'geom-point-x-y :x 200 :y 150)) (construction-of frame))

    ;;

    (push (setf c1 (make-instance 'geom-circle-center-radius :center p1 :point-on-circumference p2))

	  (construction-of frame))

    (push (setf l1 (make-instance 'geom-line-a-b :a p3 :b p4))

	  (construction-of frame))

    ;;

    (push (setf i1 (make-instance 'geom-point-circle-line-intersection :circle c1 :line l1))

	  (construction-of frame))

    (push (setf i2 (make-instance 'geom-point-circle-line-intersection :circle c1 :line l1 :other t))

	  (construction-of frame))

    (run-frame-top-level frame)))

(defun geom-test-3 ()

  ;; circle and circle intersection

  (let ((frame (make-application-frame 'geom-test))

	p1 p2 p3 p4

	c1 c2

	i1 i2)

    (push (setf p1 (make-instance 'geom-point-x-y :x (- 128 50) :y 128)) (construction-of frame))

    (push (setf p2 (make-instance 'geom-point-x-y :x 128 :y 128)) (construction-of frame))

    (push (setf p3 (make-instance 'geom-point-x-y :x (+ 128 50) :y 128)) (construction-of frame))

    (push (setf p4 (make-instance 'geom-point-x-y :x 128 :y (+ 128 50))) (construction-of frame))

    ;;

    (push (setf c1 (make-instance 'geom-circle-center-radius :center p1 :point-on-circumference p2))

	  (construction-of frame))

    (push (setf c2 (make-instance 'geom-circle-center-radius :center p3 :point-on-circumference p4))

	  (construction-of frame))

    ;;

    (push (setf i1 (make-instance 'geom-point-circle-circle-intersection :left c1 :right c2))

	  (construction-of frame))

    (push (setf i2 (make-instance 'geom-point-circle-circle-intersection :left c1 :right c2 :other t))

	  (construction-of frame))

    (run-frame-top-level frame)))

(defun geom-test-4 ()

  ;; circle from 3 points

  (let ((frame (make-application-frame 'geom-test))

	p1 p2 p3

	l1 l2

	b1 b2

	n1 n2

	i1

	c1)

    (push (setf p1 (make-instance 'geom-point-x-y :x 20 :y 120)) (construction-of frame))

    (push (setf p2 (make-instance 'geom-point-x-y :x 100 :y 40)) (construction-of frame))

    (push (setf p3 (make-instance 'geom-point-x-y :x 180 :y 90)) (construction-of frame))

    ;;

    (push (setf l1 (make-instance 'geom-line-a-b :a p1 :b p2)) (construction-of frame))

    (push (setf l2 (make-instance 'geom-line-a-b :a p2 :b p3)) (construction-of frame))

    ;;

    (push (setf b1 (make-instance 'geom-point-bisector :left p1 :right p2)) (construction-of frame))

    (push (setf b2 (make-instance 'geom-point-bisector :left p2 :right p3)) (construction-of frame))

    ;;

    (push (setf n1 (make-instance 'geom-line-normal :at b1 :to p1)) (construction-of frame))

    (push (setf n2 (make-instance 'geom-line-normal :at b2 :to p2)) (construction-of frame))

    ;;

    (push (setf i1 (make-instance 'geom-point-line-line-intersection :a n1 :b n2))

	  (construction-of frame))

    ;;

    (push (setf c1 (make-instance 'geom-circle-center-radius :center i1 :point-on-circumference p1))

	  (construction-of frame))

    (run-frame-top-level frame)))

(defclass geom-point () ())

(defclass geom-point-x-y (geom-point)

     ((x :initarg :x :accessor x-of)

      (y :initarg :y :accessor y-of)))

(defclass geom-point-line-line-intersection (geom-point)

     ((a :initarg :a :accessor a-of)

      (b :initarg :b :accessor b-of)))

(defclass geom-point-circle-line-intersection (geom-point)

     ((circle :initarg :circle :accessor circle-of)

      (line :initarg :line :accessor line-of)

      (other :initform nil :initarg :other :accessor other)))

(defclass geom-point-circle-circle-intersection (geom-point)

     ((left :initarg :left :accessor left-of)

      (right :initarg :right :accessor right-of)

      (other :initform nil :initarg :other :accessor other)))

(defclass geom-point-bisector (geom-point)

     ((left :initarg :left :accessor left-of)

      (right :initarg :right :accessor right-of)))

(defgeneric point-x-y (geom-point))

(defmethod point-x-y ((point geom-point-x-y))

  (values (x-of point) (y-of point)))

(defmethod point-x-y ((point geom-point-line-line-intersection))

  (multiple-value-bind (x1 y1) (point-x-y (a-of (a-of point)))

    (multiple-value-bind (x2 y2) (point-x-y (b-of (a-of point)))

      (multiple-value-bind (x3 y3) (point-x-y (a-of (b-of point)))

	(multiple-value-bind (x4 y4) (point-x-y (b-of (b-of point)))

	  (multiple-value-call #'line-intersection

	    x1 y1 x2 y2 x3 y3 x4 y4))))))

(defmethod point-x-y ((point geom-point-circle-line-intersection))

  (multiple-value-bind (x1 y1) (point-x-y (a-of (line-of point)))

    (multiple-value-bind (x2 y2) (point-x-y (b-of (line-of point)))

      (multiple-value-bind (cx cy cr) (center-and-radius (circle-of point))

	(let ((points (multiple-value-list (line-circle-intersection x1 y1 x2 y2 cx cy cr))))

	  (if (other point)

	      (when (= (length points) 4)

		(values (elt points 2) (elt points 3)))

	      (when (>= (length points) 2)

		(values (elt points 0) (elt points 1)))))))))

(defmethod point-x-y ((point geom-point-circle-circle-intersection))

  (multiple-value-bind (c1x c1y c1r) (center-and-radius (left-of point))

    (multiple-value-bind (c2x c2y c2r) (center-and-radius (right-of point))

      (let ((points (multiple-value-list (circle-circle-intersection c1x c1y c1r c2x c2y c2r))))

	(when points

	  (if (other point)

	      (values (elt points 2) (elt points 3))

	      (values (elt points 0) (elt points 1))))))))

(defmethod point-x-y ((point geom-point-bisector))

  (multiple-value-bind (x1 y1) (point-x-y (left-of point))

    (multiple-value-bind (x2 y2) (point-x-y (right-of point))

      (multiple-value-call #'lerp2 x1 y1 x2 y2 0.5))))

(defclass geom-line () ())

(defclass geom-line-a-b (geom-line)

     ((a :initarg :a :accessor a-of)

      (b :initarg :b :accessor b-of)))

(defclass geom-line-normal (geom-line)

     ((at :initarg :at :accessor normal-start)

      (to :initarg :to :accessor normal-reference)))

(defmethod a-of ((line geom-line-normal))

  (normal-start line))

(defmethod b-of ((line geom-line-normal))

 (multiple-value-bind (x1 y1) (point-x-y (normal-start line))

   (multiple-value-bind (x2 y2) (point-x-y (normal-reference line))

     (let ((vx (- x2 x1))

	   (vy (- y2 y1)))

       (make-instance 'geom-point-x-y :x (+ x1 vy) :y (+ y1 (- vx)))))))

(defclass geom-circle () ())

(defclass geom-circle-center-radius (geom-circle)

     ((center-point :initarg :center :accessor center-of)

      (radius-point :initarg :point-on-circumference :accessor circumference-point)))

(defgeneric center-and-radius (circle))

(defmethod center-and-radius ((circle geom-circle-center-radius))

  (multiple-value-bind (x y) (point-x-y (center-of circle))

    (multiple-value-bind (cx cy) (point-x-y (circumference-point circle))

      (values x y (distance x y cx cy)))))

(defgeneric display (object pane))

(defmethod display ((point geom-point) pane)

  (multiple-value-bind (x y) (point-x-y point)

    (ignore-errors

      (draw-circle* pane x y 2)

      (when (typep point 'geom-point-x-y)

	(draw-circle* pane x y 4 :filled nil)))))

(defmethod display ((line geom-line) pane)

  (multiple-value-bind (x1 y1) (point-x-y (a-of line))

    (multiple-value-bind (x2 y2) (point-x-y (b-of line))

      (ignore-errors

	(draw-line* pane x1 y1 x2 y2)))))

(defmethod display ((circle geom-circle) pane)

  (multiple-value-bind (x y r) (center-and-radius circle)

    (ignore-errors

      (draw-circle* pane x y r :filled nil))))

(defclass geom-pane (basic-gadget) ())

(define-application-frame geom-test ()

  ((construction :initform '() :accessor construction-of)

   (selected-point :initform '() :accessor selected-point))

  (:panes (canvas (make-pane 'geom-pane)))

  (:layouts (default (vertically (:width 256 :height 256)

		       canvas))))

(defmethod handle-repaint ((pane geom-pane) region)

  (declare (ignore region))

  (with-bounding-rectangle* (x0 y0 x1 y1) (sheet-region pane)

    (draw-rectangle* pane x0 y0 x1 y1 :filled t :ink (make-rgb-color 1 1 1))

    (mapc #'(lambda (bit) (display bit pane))

	  (construction-of *application-frame*))))

(defmethod handle-event ((pane geom-pane) (event pointer-button-press-event))

  (let ((mouse-x (pointer-event-x event))

	(mouse-y (pointer-event-y event)))

    (dolist (bit (construction-of *application-frame*))

      (when (typep bit 'geom-point-x-y)

	(multiple-value-bind (point-x point-y) (point-x-y bit)

	  (when (< (distance mouse-x mouse-y point-x point-y) 10)

	    (setf (selected-point *application-frame*) bit)

	    (return)))))))

(defmethod handle-event ((pane geom-pane) (event pointer-button-release-event))

  (setf (selected-point *application-frame*) nil))

(defmethod handle-event ((pane geom-pane) (event pointer-motion-event))

  (when (selected-point *application-frame*)

    (let ((mouse-x (pointer-event-x event))

	  (mouse-y (pointer-event-y event)))

      (setf (x-of (selected-point *application-frame*)) mouse-x)

      (setf (y-of (selected-point *application-frame*)) mouse-y))

    (repaint-sheet pane (sheet-region pane))))