;;; -*- Mode: Lisp; Package: DESIGN; Syntax: Ansi-common-lisp -*-

;; Utility functions for calculating axial alignment and rectangles of maximum extent (for
;; openness measure).

;; A line is represented as a normal and distance from origin; a normal is a unit vector of dx
;; and dy.  A point p is on a line if dot product of it and normal for the line gives d for
;; the line; above a line if value is greater than d; below a line if value is less than d.

;; For alignment of two segments, project lines from endpoints of each segment and check that
;; both endpoints of the other segment are not both on the same side of both lines for the
;; other segment.  For axial alignment, also check that the segments are parallel.

;; For rectangle of maximum extent, project lines from each endpoint of given segment.  For
;; every segment s in bounded-regions to be used to form rectangle, check for intersection of
;; each endpoint line with line through s.  If an intersection point exists, then check that
;; the point lies on s.  Collect all the intersection points that lie on segments and find the
;; ones that are closest to the original segment.  This will produce two points, one above the
;; segment, one below.  Construct the other two points and corresponding segments to form
;; rectangle.

;; For intersection of rectangles, assume rectangles will be aligned with x or y axis (so can
;; avoid general polygon intersection problem).  Find ranges for x, ranges for y, and
;; construct rectangle if possible.

;; See rectangles.lisp for toplevel maximum rectangle code and rectangle objects.


;; Other utility functions:  computing rectangle center
;;                           computing change in direction between 2 segments
;;                           checking that point is inside polygon 


(defstruct vectr
  x
  y)

(defmethod x ((v1 vectr))
  (vectr-x v1))

(defmethod y ((v1 vectr))
  (vectr-y v1))
#|
(defmethod (setf x) (value (v vectr))
  (setf (vectr-x v) value))

(defmethod (setf y) (value (v vectr))
  (setf (vectr-y v) value))
|#
(defun dot-product (a b)
  (+ (* (x a) (x b)) (* (y a) (y b))))

(defun cross-product (a b)
  (- (* (x a) (y b)) (* (y a) (x b)))) 

(defun scalar-mult (s vector)
  (setf (vectr-x vector) (* s (vectr-x vector)))
  (setf (vectr-y vector) (* s (vectr-y vector)))
  vector)


(defun vector-between-pts (pt1 pt2)
  ;; pt1 to pt2
  (make-vectr :x (- (x pt2) (x pt1)) :y (- (y pt2) (y pt1))))

(defun unit-vector-between-pts (pt1 pt2)
  ;; pt1 to pt2
  (let* ((dx (- (x pt2) (x pt1)))
	 (dy (- (y pt2) (y pt1)))
	 (length (sqrt (+ (expt dx 2)
			  (expt dy 2)))))
      (make-vectr :x (/ dx length) :y (/ dy length))))

(defstruct dline
  n              ; normal					     
  d)             ; perpendicular distance from origin

(defun multiply-matrix-simple (2x2 1x2)
  ;; 2x2 (by rows): 11, 12, 21, 22, returns 1x2
  (list (+ (* (first 2x2) (first 1x2)) (* (second 2x2) (second 1x2)))
	(+ (* (third 2x2) (first 1x2)) (* (fourth 2x2) (second 1x2)))))

(defun multiply-inverse-and-distances (n1 n2 d1 d2)
  (let* ((a (vectr-x n1))
	(b (vectr-y n1))
	(c (vectr-x n2))
	(d (vectr-y n2))
	(denom (- (* a d) (* b c))))
    (if (zerop denom) nil				; lines are parallel
	(multiply-matrix-simple `(,(/ d denom) ,(/ (- b) denom) ,(/ (- c) denom) ,(/ a denom))
				`(,d1 ,d2)))))

(defun line-intersection (l1 l2)
  (let ((pt (multiply-inverse-and-distances (dline-n l1) (dline-n l2) (dline-d l1)
					    (dline-d l2))))
    (when pt (make-point (first pt)(second pt)))))

(defun lines-for-segment-endpoints (seg)
  (multiple-value-bind (d1 d2 n)
      (distances-and-vector seg)
    (values (make-dline :n n :d d1) (make-dline :n n :d d2))))

(defun line-for-segment (seg)
  (let ((n (normal-to-segment seg)))
    (make-dline :n n :d (dot-product (endpoint1 seg) n))))

(defun line-for-segment* (seg)
  (let ((n (normal-to-segment seg)))
    (values n (dot-product (endpoint1 seg) n))))

(defun point-on-line-p (point line)
  (= (dot-product point (dline-n line)) (dline-d line)))

(defun point-above-line-p (point line)
  (>= (dot-product point (dline-n line)) (dline-d line)))

(defun point-below-line-p (point line)
  (<= (dot-product point (dline-n line)) (dline-d line)))

(defun point-on-segment-p (point segment)
  (and (or (<= (point-x (endpoint1 segment)) (point-x point) (point-x (endpoint2 segment)))
	   (<= (point-x (endpoint2 segment)) (point-x point) (point-x (endpoint1 segment))))
       (or (<= (point-y (endpoint1 segment)) (point-y point) (point-y (endpoint2 segment)))
	   (<= (point-y (endpoint2 segment)) (point-y point) (point-y (endpoint1 segment))))))

(defun distance (x1 y1 x2 y2)
  (sqrt (+ (expt (- x1 x2) 2)
	     (expt (- y1 y2) 2))))

(defun midpoint-x-y (segment)
  (let ((endpoint1 (endpoint1 segment))
	(endpoint2 (endpoint2 segment)))
    (values (avg (point-x endpoint1) (point-x endpoint2))
	    (avg (point-y endpoint1) (point-y endpoint2)))))

(defun midpoint-distance (segment1 segment2)
  (multiple-value-call #'distance (midpoint-x-y segment1) (midpoint-x-y segment2)))

(defun midpoint-point-distance (segment point)
  (multiple-value-call #'distance (midpoint-x-y segment) (point-x-y* point)))

(defun segment-length (segment)
  (when segment
    (let ((endpoint1 (endpoint1 segment))
	  (endpoint2 (endpoint2 segment)))
      (distance (point-x endpoint1) (point-y endpoint1)
		(point-x endpoint2) (point-y endpoint2)))))

(defun edge-length (edge)
  (segment-length edge))

#+ignore
(defun segment-unit-vector (segment)
  ;; from endpoint2 to endpoint1
  (let ((endpoint1 (endpoint1 segment))
	(endpoint2 (endpoint2 segment)))
    (let ((dx (- (point-x endpoint1) (point-x endpoint2)))
	  (dy (- (point-y endpoint1) (point-y endpoint2)))
	  (length (segment-length segment)))
      (make-vectr :x (/ dx length) :y (/ dy length)))))

(defun segment-unit-vector (segment)
  ;; from endpoint2 to endpoint1
  (unit-vector-between-pts (endpoint2 segment) (endpoint1 segment)))

(defun normal-to-segment (segment)
  ;; segment-normal is line through (either) endpoint of line-segment
  ;; return vector as cons
  (let ((endpoint1 (endpoint1 segment))
	(endpoint2 (endpoint2 segment)))
    (let ((dx (- (point-x endpoint1) (point-x endpoint2)))
	  (dy (- (point-y endpoint1) (point-y endpoint2)))
	  (length (segment-length segment)))
      (make-vectr :x (/ (- dy) length) :y (/ dx length)))))

(defun distances-and-vector (segment)
  ;; two lines:  one passing perpendicular to line segment, through each endpoint
  (let ((endpoint1 (endpoint1 segment))
	(endpoint2 (endpoint2 segment)))
    (let ((vector (segment-unit-vector segment)))
      (values (dot-product endpoint1 vector) (dot-product endpoint2 vector) vector))))


;;; Alignment

(defun segment-outside-p (segment1 segment2)
  ;; are both endpoints of segment1 on same side of perpendicular lines through
  ;; endpoints of segment2
  (multiple-value-bind (l1 l2)
      (lines-for-segment-endpoints segment2)
    (let ((pt1 (endpoint1 segment1))
	  (pt2 (endpoint2 segment1)))
      (or (and (point-above-line-p pt1 l1) (point-above-line-p pt2 l1)
	       (point-above-line-p pt1 l2) (point-above-line-p pt2 l2))
	  (and (point-below-line-p pt1 l1) (point-below-line-p pt2 l1)
	       (point-below-line-p pt1 l2) (point-below-line-p pt2 l2))))))


(defun aligned-segments-p (segment1 segment2)
;; don't have to be parallel, any point on one segment lies within bounds of normals to
;; endpoints of other segment and vv
  (not (or (segment-outside-p segment1 segment2) (segment-outside-p segment2 segment1))))

(defvar *sine-parallel-tolerance-angle* .087)			; 5 degrees

(defun parallel-segments-p (segment1 segment2)
  ;; normals are (approx) equal
  ;; cross-product of normals gives sine of angle; check this for epsilon
  (<= (abs (cross-product (segment-unit-vector segment1)
			  (segment-unit-vector segment2)))
      *sine-parallel-tolerance-angle*))

(defun parallel-normals-p (n1 n2)
  (<= (abs (cross-product n1 n2)) *sine-parallel-tolerance-angle*))

(defun axially-aligned-segments-p (segment1 segment2)
;; have to be aligned and parallel
  (and (parallel-segments-p segment1 segment2)
       (aligned-segments-p segment1 segment2)))

(defun colinearp (segment1 segment2)
  ;; can't just check absolute value of d1 and d2 because might be same d but on opposite
  ;; sides of origin (d is + in same direction as n)
  (multiple-value-bind (n1 d1)
    (line-for-segment* segment1)
    (multiple-value-bind (n2 d2)
      (line-for-segment* segment2)
      (and (zerop (cross-product n1 n2))
	   (if (plusp (dot-product n1 n2))		       
	       (= d1 d2)
	       (= d1 (- d2)))))))

#+ignore
(defvar  *colinear-tolerance-d* 1.0)				;?? units
#+ignore
(defun almost-colinearp (segment1 segment2)
  (multiple-value-bind (n1 d1)
    (line-for-segment* segment1)
    (multiple-value-bind (n2 d2)
      (line-for-segment* segment2)
      (and (parallel-normals-p n1 n2)
	   (if (plusp (dot-product n1 n2))
	       (or (= d1 d2) (<= (abs (- d1 d2)) *colinear-tolerance-d*))
	       (or (= d1 (- d2)) (<= (abs (+ d1 d2)) *colinear-tolerance-d*)))))))

(defvar *tan-colinear-tolerance-angle* .052)			; 3 degrees

(defun almost-colinearp (segment1 segment2)
  ;; version that is independent of coordinate system units: find largest distance between
  ;; segments (by checking distance between each segment's endpoints and the line through the
  ;; other segment), normalize that distance by e.g. length of longest segment; check that
  ;; this ratio is less than a tolerance (based on tan of angle)
  (multiple-value-bind (n1 d1)
    (line-for-segment* segment1)
    (multiple-value-bind (n2 d2)
      (line-for-segment* segment2)
      (and (parallel-normals-p n1 n2)
	   (multiple-value-bind (delta1 delta2 delta3 delta4)
	       (endpoint-distances-to-lines segment1 segment2 n1 d1 n2 d2)
	     (<= (/  (max delta1 delta2 delta3 delta4)
		     (max (segment-length segment1) (segment-length segment2)))
		 *tan-colinear-tolerance-angle*))))))

(defun endpoint-distances-to-lines (seg1 seg2 n1 d1 n2 d2)
  ;; endpoints of seg1 to line2 (d2 n2), endpoints of seg2 to line1 (d1 n1)
  (values (distance-from-pt-to-line* (endpoint1 seg1) n2 d2)
	  (distance-from-pt-to-line* (endpoint2 seg1) n2 d2)
	  (distance-from-pt-to-line* (endpoint1 seg2) n1 d1)
	  (distance-from-pt-to-line* (endpoint2 seg2) n1 d1)))

(defun intersecting-segments-p (segment1 segment2)
  (if (colinearp segment1 segment2)
      nil
      (let ((pt (line-intersection (line-for-segment segment1) (line-for-segment segment2))))
	(if (and pt (point-on-segment-p pt segment1) (point-on-segment-p pt segment2))
	    t
	    nil))))

(defun segment2-line-intersects-lines-through-segment1-endpoints (seg1 seg2
							       &optional l1a l1b)
  (unless (and l1a l1b) (multiple-value-setq (l1a l1b) (lines-for-segment-endpoints seg1)))
  (let ((l2 (line-for-segment seg2)))
    (values (line-intersection l1a l2) (line-intersection l1b l2))))

(defun check-for-segment-line-intersections (seg1 segments &optional endpoints-to-avoid)
  ;; which segments intersect lines perpendicular to seg1 endpoints?
  (multiple-value-bind (l1a l1b)
      (lines-for-segment-endpoints seg1)
    (loop for seg2 in segments
	  with endpoint1-intersections
	  with endpoint2-intersections
	  do (unless (or (edges-share-equal-endpoint-p seg1 seg2)
			 (member (endpoint1 seg2) endpoints-to-avoid :test #'point-equal)
			 (member (endpoint2 seg2) endpoints-to-avoid :test #'point-equal))
	       (multiple-value-bind (pt1 pt2)
		   (segment2-line-intersects-lines-through-segment1-endpoints seg1 seg2
									      l1a l1b)
		 (when (and pt1 (point-on-segment-p pt1 seg2))
		   (pushnew pt1 endpoint1-intersections :test #'point-equal))
		 (when (and pt2 (point-on-segment-p pt2 seg2))
		   (pushnew pt2 endpoint2-intersections :test #'point-equal))))
	  finally (return (values endpoint1-intersections endpoint2-intersections
				  l1a l1b)))))


(defun distance-from-pt-to-line (pt line)
  (abs (- (dline-d line) (dot-product pt (dline-n line)))))

(defun distance-from-pt-to-line* (pt n d)
  (abs (- d (dot-product pt n))))

;; Functions for use in finding valid intersection pts

(defun segment-endpt-p (pt segments)
	   (loop for seg in segments
		 when (or (point-equal pt (endpoint1 seg))
			  (point-equal pt (endpoint2 seg)))
		   do (return t)))

(defun shared-endpoint (seg1 seg2)
	     (if (or (point-equal (endpoint1 seg1) (endpoint1 seg2))
		     (point-equal (endpoint1 seg1) (endpoint2 seg2)))
		 (endpoint1 seg1) (endpoint2 seg2)))

(defun intersection-point-on-same-side-p (pt seg test-line other-line)
	     ;; find intersection of line through seg with other-line, then see if pt and
	     ;; intersection pt are on same side of test-line
	     (let ((int-pt (line-intersection (line-for-segment seg) other-line)))
	       (or (and (point-above-line-p pt test-line)
			(point-above-line-p int-pt test-line))
		   (and (point-below-line-p pt test-line)
			(point-below-line-p int-pt test-line)))))

(defun	segment-pts-on-opposite-sides (seg1 seg2 line line2)
	     ;; segs share an endpt so just check other endpts
	     (let* ((pt (shared-endpoint seg1 seg2))
		   (other-seg1-pt (other-endpoint seg1 pt))
		   (other-seg2-pt (other-endpoint seg2 pt)))
	       (or (and (point-above-line-p other-seg1-pt line)
			(point-below-line-p other-seg2-pt line))
		   (and (point-below-line-p other-seg1-pt line)
			(point-above-line-p other-seg2-pt line))
		   (and (point-on-line-p other-seg1-pt line)
			(intersection-point-on-same-side-p other-seg2-pt seg2 line line2))
		   (and (point-on-line-p other-seg2-pt line)
			(intersection-point-on-same-side-p other-seg1-pt seg1 line line2)))))

(defun segments-for-pt-cross-line (pt segments line other-line)
	     (let ((segs (loop for seg in segments
			       when (or (point-equal pt (endpoint1 seg))
					(point-equal pt (endpoint2 seg)))
				 collect seg)))
	     ;; assume at least 2 segments for now
	     (loop for seg1 in segs
		     do  (loop for seg2 in (remove seg1 segs)
			    when (segment-pts-on-opposite-sides seg1 seg2 line other-line)
			      do (return-from segments-for-pt-cross-line t)))))

(defun valid-intersection-pt-p (pt segments line other-line)
  ;; if pt is not segment endpoint, ok  
  ;; if pt is segment endpoint: if the segments for which the pt is an endpoint together
  ;; cross the line on which pt lies (which passes through one of seg1 endpoints) then the 
  ;; pt can be counted as a valid intersection point; if one of the (endpoint) segments lies
  ;; on the line, have to check that the other (endpoint) segment is between the two lines:
  ;; check that its nonshared endpoint and the intersection of its line and other-line are
  ;; on the same side of line
  (or (not (segment-endpt-p pt segments))
      (segments-for-pt-cross-line pt segments line other-line)))

(defun valid-intersection-pts (pts segments line other-line)
  (loop for pt in pts
	when (valid-intersection-pt-p pt segments line other-line)
	  collect pt))

(defun find-minimum-intersection-points (seg1 segments &optional endpoints-to-avoid)
  ;; line1 and line2 are lines perdendicular and passing through seg1 endpoints 
  (multiple-value-bind (endpoint1-pts endpoint2-pts line1 line2)
      (check-for-segment-line-intersections seg1 segments endpoints-to-avoid)
    (setq endpoint1-pts (valid-intersection-pts endpoint1-pts segments line1 line2))
    (setq endpoint2-pts (valid-intersection-pts endpoint2-pts segments line2 line1))
    (flet ((rtn-value (pt pt-d)
	     (list pt pt-d (if (member pt endpoint1-pts) line1
			  line2))))
    (loop for pt in (append endpoint1-pts endpoint2-pts)
	  with seg1-line = (line-for-segment seg1)
	  as pt-d = (distance-from-pt-to-line pt seg1-line)
	  with above-pts
	  with below-pts
	  with min-above-pt  ; dotted pair of pt and pt-d
	  with min-below-pt
	  when (point-above-line-p pt seg1-line)
	    do (progn (push pt above-pts)
		      (if min-above-pt
			  (when (< pt-d (cdr min-above-pt))
			    (setq min-above-pt (cons pt pt-d)))
			  (setq min-above-pt (cons pt pt-d))))
	  when (point-below-line-p pt seg1-line)
	    do (progn (push pt below-pts)
		      (if min-below-pt
			  (when (< pt-d (cdr min-below-pt))
			    (setq min-below-pt (cons pt pt-d)))
			  (setq min-below-pt (cons pt pt-d))))
	  finally (return (values (rtn-value (car min-above-pt) (cdr min-above-pt))
				  (rtn-value (car min-below-pt) (cdr min-below-pt))
				  line1 line2 seg1-line))))))

(defun find-segments-for-rectangle-of-max-extent (seg1 segments &optional endpoints-to-avoid)
  ;; seg1 is parallel to end of rectangle, and its endpoints intersect sides of rectangle
  ;; min-... is list of pt, distance from pt to segment, seg1 endpoint pt is "even" with
  (let (new-pt-above new-pt-below new-pt-above-line-from-seg new-pt-above-line-from-pt
	new-pt-below-line-from-seg new-pt-below-line-from-pt)
  (multiple-value-bind (min-above-pt-data min-below-pt-data line1 line2 seg1-line)
      (find-minimum-intersection-points seg1 segments endpoints-to-avoid)
    ;; construct points for the other end of the segment on each side
    (destructuring-bind (pt-above pt-above-d pt-above-line-from-seg)
	min-above-pt-data
      (setq new-pt-above-line-from-seg (if (eq pt-above-line-from-seg line1) line2 line1))
      ;; new point should be intersection of new-point-above-line-from-seg and line with
      ;; n = seg1-line normal and d = pt-d (i.e. line through both above line points)
      (setq new-pt-above-line-from-pt
	    (make-dline :n (dline-n seg1-line) :d (+ (dline-d seg1-line) pt-above-d)))
      (setq new-pt-above (line-intersection new-pt-above-line-from-seg
					    new-pt-above-line-from-pt))
      (destructuring-bind (pt-below pt-below-d pt-below-line-from-seg)
	  min-below-pt-data
	(setq new-pt-below-line-from-seg (if (eq pt-below-line-from-seg line1) line2 line1))
	(setq new-pt-below-line-from-pt
	    (make-dline :n (dline-n seg1-line) :d (- (dline-d seg1-line) pt-below-d)))
	(setq new-pt-below (line-intersection new-pt-below-line-from-seg
					    new-pt-below-line-from-pt))
	(append (list (make-basic-edge pt-above new-pt-above)
		      (make-basic-edge pt-below new-pt-below))
		(if (eq pt-above-line-from-seg pt-below-line-from-seg)
		    (list (make-basic-edge pt-above pt-below)
			  (make-basic-edge new-pt-above new-pt-below))
		    (list (make-basic-edge pt-above new-pt-below)
			  (make-basic-edge new-pt-above pt-below)))))))))

;; Rectangle centers

(defun rectangle-center (minx miny maxx maxy)
  (values (avg minx maxx) (avg miny maxy)))

(defun rectangle-center+ (minx miny maxx maxy)
  (list (avg minx maxx) (avg miny maxy)))


;;;  Change in direction: angle between extension of one segment and other segment
;;;  (i.e. supplement of angle between)

;;;  Note:  this function would still work if segments don't share an endpoint but can
;;;  be ordered seg1 to seg2; have to worry about sign though; and have to use all 4 endpoints

;;;  assumes units are radians
  
(defun direction-change (x1 y1 x2 y2 x3 y3)
  ;; shared point is x2y2, vectors are from x2y2 to x1y1 and x3y3 to x2y2 
  ;; order of pts along path: x1y1 x2y2 x3y3
  (let ((ax (- x2 x1))
	(ay (- y2 y1))
	(bx (- x3 x2))
	(by (- y3 y2)))
    (atan (- (* ax by) (* ay bx)) (+ (* ax bx) (* ay by)))))

(defun segment-direction-change (seg1 seg2)
  ;; segments are ordered seg1 to seg2, shared endpoint is designated endpoint2
  ;; points in order are designated endpoint1 endpoint2 endpoint3.
    (let* ((endpoint1 (if (endpoint-of-edge-p seg2 (endpoint1 seg1))
			  (endpoint2 seg1) (endpoint1 seg1)))
	   (endpoint2 (other-endpoint seg1 endpoint1))
	   (endpoint3 (other-endpoint seg2 endpoint2)))
      (direction-change (point-x endpoint1) (point-y endpoint1) (point-x endpoint2)
			(point-y endpoint2) (point-x endpoint3) (point-y endpoint3))))

;;;+++ Broken?  get different answer if reverse args

(defun angle-between (seg1 seg2)
  ;; negate one of segments and compute as in change in direction
  ;; also have to reverse sign of angle because it's calculated from seg1 to seg2 and
  ;; supplement should be calculated -seg2 to seg1
    (let* ((endpoint1 (if (endpoint-of-edge-p seg2 (endpoint1 seg1))
			  (endpoint2 seg1) (endpoint1 seg1)))
	   (endpoint2 (other-endpoint seg1 endpoint1))
	   (endpoint3 (other-endpoint seg2 endpoint2))
	   (x1 (point-x endpoint1))
	   (y1 (point-y endpoint1))
	   (x2 (point-x endpoint2))
	   (y2 (point-y endpoint2))
	   (x3 (point-x endpoint3))
	   (y3 (point-y endpoint3))
	   (bx (- x1 x2))
	   (by (- y1 y2))
	   (ax (- x3 x2))
	   (ay (- y3 y2)))
    (atan (- (* ax by) (* ay bx)) (+ (* ax bx) (* ay by)))))

;; Point inside polygon
#+ignore
(defun insidep (x y segments)
  ;; segments are edges of the polygon
  (loop with cnt = 0
	for segment in segments
	as pt1 = (endpoint1 segment)
	as pt2 = (endpoint2 segment)
	as x1 = (- (point-x pt1) x)
	as y1 = (- (point-y pt1) y)
	as x2 = (- (point-x pt2) x)
	as y2 = (- (point-y pt2) y)
	unless (or (and (> x1 0) (> x2 0))
		   (and (< x1 0) (< x2 0))
		   (and (> y1 0) (> y2 0))
		   (< 0 (/ (- (* x2 y1) (* x1 y2)) (- x2 x1))))
	  do (setf cnt (1+ cnt))
	finally (return (oddp cnt))))


;;; Need fancier version to take care of hitting a vertex.  Code courtesy of JGA.

;;; INSIDE-POLYGON-P tests whether a point is enclosed in a face, which is defined by
;;; edges.  It uses the fact that a line from a test point to infinity will cross an odd
;;; number of edges if the point is enclosed.  The test line is in the negative y
;;; direction.  Special logic covers the case of hitting a vertex.

(defun inside-polygon-p (x-tst y-tst edges)
  (let ((encl-flag nil) (half-flag nil) (half-flag-direction nil))
    (macrolet ((intersects (x0 y0 x1 y1)
                 `(and (eq (< ,x0 x-tst) (> ,x1 x-tst))
                       (> y-tst
                          (/ (- (* ,y1 (- x-tst ,x0)) (* ,y0 (- x-tst ,x1)))
                             (- ,x1 ,x0)))))
               (half-edge-bookkeeping (x)
                 `(let ((new-half-edge (> ,x x-tst)))   ; t = goes right
                    (cond ((null half-flag)
                           (setq half-flag T half-flag-direction new-half-edge))
                          ((eq new-half-edge half-flag-direction)
                           (setq half-flag nil))
                          (T
                           (setq half-flag nil)
                           (setq encl-flag (not encl-flag)))))))
      (loop for edge in edges
            for p0 = (endpoint1 edge) for x0 = (point-x p0) for y0 = (point-y p0)
            for p1 = (endpoint2 edge) for x1 = (point-x p1) for y1 = (point-y p1)
            do (cond ((= x0 x1))
                     ((and (= x0 x-tst) (< y0 y-tst)) (half-edge-bookkeeping x1))
                     ((and (= x1 x-tst) (< y1 y-tst)) (half-edge-bookkeeping x0))
                     ((intersects x0 y0 x1 y1) (setq encl-flag (not encl-flag))))
            finally (return encl-flag)))))

#+ignore
(defun inside-or-on-polygon-p (x-tst y-tst edges)
  ;; change intersection test to include =
  (let ((encl-flag nil) (half-flag nil) (half-flag-direction nil))
    (macrolet ((intersects (x0 y0 x1 y1)
                 `(and (eq (<= ,x0 x-tst) (>= ,x1 x-tst))
                       (>= y-tst
                          (/ (- (* ,y1 (- x-tst ,x0)) (* ,y0 (- x-tst ,x1)))
                             (- ,x1 ,x0)))))
               (half-edge-bookkeeping (x)
                 `(let ((new-half-edge (> ,x x-tst)))   ; t = goes right
                    (cond ((null half-flag)
                           (setq half-flag T half-flag-direction new-half-edge))
                          ((eq new-half-edge half-flag-direction)
                           (setq half-flag nil))
                          (T
                           (setq half-flag nil)
                           (setq encl-flag (not encl-flag)))))))
      (loop for edge in edges
            for p0 = (endpoint1 edge) for x0 = (point-x p0) for y0 = (point-y p0)
            for p1 = (endpoint2 edge) for x1 = (point-x p1) for y1 = (point-y p1)
            do (cond ((= x0 x1))
                     ((and (= x0 x-tst) (< y0 y-tst)) (half-edge-bookkeeping x1))
                     ((and (= x1 x-tst) (< y1 y-tst)) (half-edge-bookkeeping x0))
                     ((intersects x0 y0 x1 y1) (setq encl-flag (not encl-flag))))
            finally (return encl-flag)))))

(defun inside-or-on-polygon-p (x-tst y-tst &rest endpoints)
  ;; change intersection test to include =
  (let ((encl-flag nil) (half-flag nil) (half-flag-direction nil))
    (macrolet ((intersects (x0 y0 x1 y1)
                 `(and (eq (<= ,x0 x-tst) (>= ,x1 x-tst))
                       (>= y-tst
                          (/ (- (* ,y1 (- x-tst ,x0)) (* ,y0 (- x-tst ,x1)))
                             (- ,x1 ,x0)))))
               (half-edge-bookkeeping (x)
                 `(let ((new-half-edge (> ,x x-tst)))   ; t = goes right
                    (cond ((null half-flag)
                           (setq half-flag T half-flag-direction new-half-edge))
                          ((eq new-half-edge half-flag-direction)
                           (setq half-flag nil))
                          (T
                           (setq half-flag nil)
                           (setq encl-flag (not encl-flag)))))))
      (loop for (pt . rest) on (append endpoints (list (car endpoints)))
	    until (not rest)
            for p0 = pt for x0 = (point-x p0) for y0 = (point-y p0)
            for p1 = (car rest) for x1 = (point-x p1) for y1 = (point-y p1)
            do (cond ((= x0 x1))
                     ((and (= x0 x-tst) (< y0 y-tst)) (half-edge-bookkeeping x1))
                     ((and (= x1 x-tst) (< y1 y-tst)) (half-edge-bookkeeping x0))
                     ((intersects x0 y0 x1 y1) (setq encl-flag (not encl-flag))))
            finally (return encl-flag)))))

;; seg1 and seg2 are aligned, but their endpoints overlap (in one dimension); find min segment
;; "between" the endpoints by projecting each endpoint onto a line and taking the middle
;; endpoints

(defun min-overlapping-segment (seg1 seg2)
  ;; project endpoints onto direction (so get one dimension), check for overlap, then take
  ;; middle two points
  (let* ((e11 (endpoint1 seg1)) (e21 (endpoint2 seg1))
	 (e12 (endpoint1 seg2)) (e22 (endpoint2 seg2))
	 (v (vector-between-pts e11 e21))
	 (n11 (dot-product e11 v)) (n21 (dot-product e21 v))
	 (n12 (dot-product e12 v)) (n22 (dot-product e22 v))
	 (small1 (min n11 n21))
	 (small1-endpoint (if (= small1 n11) e11 e21))
	 (large1 (max n11 n21))
	 (large1-endpoint (if (= large1 n11) e11 e21))
	 (small2 (min n12 n22))
	 (small2-endpoint (if (= small2 n12) e12 e22))
	 (large2 (max n12 n22))
	 (large2-endpoint (if (= large2 n12) e12 e22))
	 (large (min large1 large2))
	 (large-endpoint (if (= large large1) large1-endpoint large2-endpoint))
	 (small (max small1 small2))
	 (small-endpoint (if (= small small1) small1-endpoint small2-endpoint)))
    (when (and (> large small))
      ;; project large-endpoint onto small-endpoint's segment
      (let* ((seg-small (if (member seg1 (edges small-endpoint)) seg1 seg2))
	     (seg-large (if (eq seg-small seg1) seg2 seg1))
	     (line (line-for-segment seg-small))
	     (vector (segment-unit-vector seg-large))
	     (line-perpendicular-to-large-endpoint
	       (make-dline :n vector :d (dot-product large-endpoint vector)))
	     (pt (line-intersection line line-perpendicular-to-large-endpoint)))
      (make-basic-edge small-endpoint pt)))))