;;; -*- Mode: Lisp; Syntax: Common-lisp; Package: MU; Base: 10.; -*-

;;; Copyright (c) 1987, Massachusetts Institute of Technology
;;; Author: Mike Drumheller

;;;*****************************************************************************
;;; CHANGE HISTORY
;;;
;;; 10/12/87  Changed the package to mu.  Recompiled for Release 7.  (W. Gillett)
;;;*****************************************************************************

(in-package :mu)

;;; This is my memory-efficient representation of points.
(defsubst MAKE-POINT (x y)
  (let ((point 0.))
    (setq point (dpb x (byte 16 0) point))
    (setq point (dpb y (byte 16 16) point))
    point))

;;;  Accessors for points.
(defsubst PTX (point)
  (ldb (byte 16 0) point))

(defsubst PTY (point)
  (ldb (byte 16 16) point))

;;;  INDICES-OK? just tells you whether X and Y are alright to use as arguments to
;;;  RASTER-AREF for a raster having width XDIM and height YDIM.
(defun INDICES-OK? (x y xdim ydim)
  (not (or (minusp x)
	   (minusp y)
	   (>= x xdim)
	   (>= y ydim))))

(defvar *POINTS* nil)
;;;  What the hell does this do?
(defun COLLECT-REGION-SUBROUTINE
       (array x y &optional (show-every-other-time--this-arg-can-be-nil 20)
	      (xdim (array-dimension array 0)) (ydim (array-dimension array 1)))
  (if show-every-other-time--this-arg-can-be-nil
      (if (= 1 (random show-every-other-time--this-arg-can-be-nil))
	  (show-1b-raster array 0 0 *mike-display-pane* nil nil)))
  (setf (aref array x y) 0)
  (cond ((and     (indices-ok? (1+ x) y xdim ydim)
	      (= 1 (aref array (1+ x) y)))
	 (setf (aref array (1+ x) y) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1+ x) y  xdim ydim)))
  (cond ((and     (indices-ok? (1+ x) (1+ y) xdim ydim)
	      (= 1 (aref array (1+ x) (1+ y))))
	 (setf (aref array (1+ x) (1+ y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1+ x) (1+ y)  xdim ydim)))
  (cond ((and     (indices-ok? x (1+ y) xdim ydim)
	      (= 1 (aref array x (1+ y))))
	 (setf (aref array x (1+ y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array x (1+ y) xdim ydim)))
  (cond ((and     (indices-ok? (1- x) (1+ y) xdim ydim)
	      (= 1 (aref array (1- x) (1+ y))))
	 (setf (aref array (1- x) (1+ y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1- x) (1+ y)  xdim ydim)))
  (cond ((and     (indices-ok? (1- x) y xdim ydim)
	      (= 1 (aref array (1- x) y)))
	 (setf (aref array (1- x) y) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1- x) y xdim ydim)))
  (cond ((and     (indices-ok? (1- x) (1- y) xdim ydim)
	      (= 1 (aref array (1- x) (1- y))))
	 (setf (aref array (1- x) (1- y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1- x) (1- y)  xdim ydim)))
  (cond ((and     (indices-ok? x (1- y) xdim ydim)
	      (= 1 (aref array x (1- y))))
	 (setf (aref array x (1- y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array x (1- y)  xdim ydim)))
  (cond ((and     (indices-ok? (1+ x) (1- y) xdim ydim)
	      (= 1 (aref array (1+ x) (1- y))))
	 (setf (aref array (1+ x) (1- y)) 0)
	 (push (make-point x y) *points*)
	 (collect-region-subroutine array (1+ x) (1- y) xdim ydim))))

(defun COLLECT-REGION
       (array x y &optional (show-every-other-time--this-arg-can-be-nil 20))
  (setq *points* (list (make-point x y)))
  (collect-region-subroutine array x y show-every-other-time--this-arg-can-be-nil)
  *points*)

(defun CENTROID-OF-LIST-OF-POINTS (points)
  (let ((total (length points))
	(xsum (loop for p in points summing
		    (ptx p)))
	(ysum (loop for p in points summing
				      (pty p))))
    (values (floor xsum total) (floor ysum total))))

(defun COLLECT-REGIONS
       (array &optional (show-every-other-time--this-arg-can-be-nil 20)
	      (scratch-array (make-array (array-dimensions array) :element-type '(mod 2))))
  (move-array array scratch-array)
  (let ((blobs nil))
    (loop for x from 0 below (array-dimension scratch-array 0) do
      (loop for y from 0 below (array-dimension scratch-array 1) do
	(if (= 1 (aref scratch-array x y))
	    (push (collect-region
		    scratch-array x y show-every-other-time--this-arg-can-be-nil) blobs))))
    blobs))

;;; ARRAY is left untouched, SCRATCH-ARRAY is returned.
(defun REPLACE-REGIONS-WITH-THEIR-CENTROIDS
       (array
	 &optional (blob-size-threshold 10) (show-every-other-time--this-arg-can-be-nil 20)
	 (scratch-array (make-array (array-dimensions array) :element-type '(mod 2))))
  (move-array array scratch-array)
  (let ((blobs (collect-regions
		 scratch-array show-every-other-time--this-arg-can-be-nil scratch-array)))
    (loop for blob in blobs do
      (if (>= (length blob) blob-size-threshold)
	  (multiple-value-bind (x y)
	      (centroid-of-list-of-points blob)
	    (setf (aref scratch-array x y) 1)))))
  scratch-array)

(defun DRAW-POINTS-INTO-ARRAY (list-of-points array &optional (show? t))
  (loop for point in list-of-points do
    (if show?
	(show-1b-raster array))
    (setf (aref array (ptx point) (pty point)) 1)))

(defun DRAW-REGIONS-INTO-ARRAY
       (list-of-lists-of-points array &optional (show? t) (zero-raster-first? t))
  (if zero-raster-first?
      #+symbolics
      (vu:zero-raster array)
      #-symbolics
      (error "Unimplemented, complain to PAO")
      )
  (loop for region in list-of-lists-of-points do
    (draw-points-into-array region array show?)))
  
;;; THESE ARE OLD VERSIONS OF ENDPOINT-FIT AND BREAK-CONTOUR WHICH BROKE THE CONTOUR AT THE 
;;; POINT FARTHEST AWAY (IN THE PERPENDICULAR DIRECTION) FROM THE LINE THROUGH THE
;;; TWO ENDPOINTS OF THE CONTOUR.
;;;(DEFUN ENDPOINT-FIT 
;;;        (MIN-PTS MAX-DEV LIST-OF-CONTOURS &AUX FIRST-HALF SECOND-HALF
;;;              DEV-AT-BREAKPOINT)
;;;  (IF (NULL LIST-OF-CONTOURS)
;;;      NIL
;;;      (IF (< (LENGTH (CAR LIST-OF-CONTOURS)) MIN-PTS)
;;;	  (ENDPOINT-FIT MIN-PTS MAX-DEV (CDR LIST-OF-CONTOURS))
;;;	  (MULTIPLE-VALUE (FIRST-HALF SECOND-HALF DEV-AT-BREAKPOINT)
;;;	    (BREAK-CONTOUR (CAR LIST-OF-CONTOURS)))
;;;	  (IF (> DEV-AT-BREAKPOINT MAX-DEV)
;;;	      (ENDPOINT-FIT MIN-PTS MAX-DEV
;;;              (APPEND (LIST FIRST-HALF) (LIST SECOND-HALF) (CDR LIST-OF-CONTOURS)))
;;;	      (COND (*SHOWP*
;;;		     (DRAW-CONTOURS-INTO-ARRAY (LIST (CAR LIST-OF-CONTOURS)))
;;;		     (SHOW-1B-RASTER *DRAW-ARRAY*)
;;;		     (IF *PAUSEP* (PAUSE))))
;;;	      (APPEND (LIST (CAR LIST-OF-CONTOURS))
;;;		      (ENDPOINT-FIT MIN-PTS MAX-DEV (CDR LIST-OF-CONTOURS)))))))
;;;
;;;(DEFUN BREAK-CONTOUR (CONTOUR &AUX (FIRST (CAR CONTOUR)) (LAST (CAR (LAST CONTOUR)))
;;;		      (CURR-DEV 0) (MAX-DEV-ON-THIS-CONTOUR 0) (NTH 1))
;;;  (DO ((POINT (CDR CONTOUR) (CDR POINT))
;;;       (N 2 (1+ N)))
;;;      ((NULL (CDR POINT))
;;;       (VALUES (FIRSTN NTH CONTOUR) (NTHCDR (1- NTH) CONTOUR) MAX-DEV-ON-THIS-CONTOUR))
;;;    (SETQ CURR-DEV (ABS (DIST-FROM-POINT-TO-LINE (CAR POINT) FIRST LAST)))
;;;    (COND ((> CURR-DEV MAX-DEV-ON-THIS-CONTOUR)
;;;	   (SETQ MAX-DEV-ON-THIS-CONTOUR CURR-DEV
;;;		 NTH N)))))

#+symbolics
(defun DRAW-REGION
       (region
	 &optional
	 double-thick?
	 (start-point 0)
	 (end-point (1- (length region)))
	 (window-or-array *mike-display-pane*)
	 (alu tv:alu-seta))
  (move-cursor-if-nec 260)
  (let ((c (nthcdr start-point region)))
    (loop for i from 0 to (- end-point start-point)
	  for point in c do
      (let ((x (ptx point))
	    (y (pty point)))
	(cond ((arrayp window-or-array)
	       (if double-thick?
		   (si:%draw-rectangle 2 2 x y alu window-or-array)
		   (si:%draw-rectangle 1 1 x y alu window-or-array)))
	      (t
	       (zl-user:send window-or-array :draw-point x y)
	       (cond (double-thick?
		      (zl-user::send window-or-array :draw-rectangle x y 2 2)))))))))

(defvar *SELECTED-CONTOURS* nil)
#+symbolics
(defun DRAW-REGIONS (list-of-regions  &optional pause? double-thick?
				      (window-or-array *mike-display-pane*))
  (setq *selected-contours* nil)
  (loop for region in list-of-regions do
    (draw-region region double-thick? 0 (1- (length region)) window-or-array)
    (cond (pause?
	   (pause)
	   (if (tv:key-state (char-code #\s))
	       (push region *selected-contours*))))))

(defvar *SHOWP* nil)
(defun SHOWP () (setq *showp* (not *showp*)))

(defvar *PAUSEP* nil)

;;; THE FUNCTION ENDPOINT-FIT SHOULD BE USED WITH EXTREME CAUTION.  IT HAS BEEN KNOWN TO CRASH
;;; LISP MACHINES FOR NO APPARENT REASON.
;;;
;;; This newer, better version of ENDPOINT-FIT (besides breaking Lispms) breaks the contours
;;; at the point A where the sum of the lengths of the two segments from endpoint1 to A and
;;; from endpoint2 to A is greatest.
;;;
;;; This is important because a contour with the following shape might not be broken by the
;;; previous method, because even though it is a very sharply bent contour, the maximum
;;; deviation of any point on the contour from the line between e1 and e2 is still very
;;; small.
;;;
;;;         **
;;;        *  *
;;;        *   *
;;;        *   *
;;;        *   *
;;;        *   *
;;;      e1 *   *
;;;            *
;;;          *
;;;         * e2
;;;
;;;
;;; Low RATIO implies that it will try to find very straight lines.
;;;
;;; RATIO cannot be less than 1 because it is impossible to draw a straight line between
;;; the two endpoints of a contour and have it be less than the sum of two lines from
;;; <endpoint-a> to <point-x-on-the-contour> and from <point-x-on-the-contour> to
;;; <endpoint-b>.
;;;
;;; 1.01 is usually a good value for RATIO.
(defun ENDPOINT-FIT (list-of-contours &optional (min-pts 10) (ratio 1.01))
  (let (first-half second-half dist-sum-at-breakpoint dist-between-endpts)
    (if (not (> ratio 1.0))
	(error
	  "RATIO (~a) must be greater than 1.
           See the documentation of this function." ratio))
    (if (null list-of-contours)
	nil
	(if (< (length (car list-of-contours)) min-pts)
	    (endpoint-fit (cdr list-of-contours) min-pts ratio)
	    (setq dist-between-endpts
		  (let* ((p1 (caar list-of-contours))
			 (p2 (car (last (car list-of-contours))))
			 (x1 (ptx p1))
			 (y1 (pty p1))
			 (x2 (ptx p2))
			 (y2 (pty p2)))
		    (dist x1 y1 x2 y2)))
	    (multiple-value-setq
	      (first-half second-half dist-sum-at-breakpoint)
	      (break-contour (car list-of-contours)))
	    (if (> (/ dist-sum-at-breakpoint dist-between-endpts) ratio) ;if contour is
                                                                      ; not straight enough
		(endpoint-fit
		  (append (list first-half)
			  (list second-half)
			  (cdr list-of-contours))
		  min-pts ratio)
		(cond (*showp*
		       (draw-region (car list-of-contours))
		       (if *pausep* (pause))))
		(append (list (car list-of-contours))
			(endpoint-fit (cdr list-of-contours) min-pts ratio)))))))

(defun BREAK-CONTOUR (contour)
  (let ((first (car contour))
	(last (car (last contour)))
	(curr-dist-sum 0)
	(max-dist-sum-on-this-contour 0)
	(nth 1))
    (do ((point (cdr contour) (cdr point))
	 (n 2 (1+ n)))
	((null (cdr point))
	 (values (zl:firstn nth contour)
		 (nthcdr (1- nth) contour)
		 max-dist-sum-on-this-contour))
      (let* ((p (car point))
	     (xp (ptx p))
	     (yp (pty p))
	     (xf (ptx first))
	     (yf (pty first))
	     (xl (ptx last))
	     (yl (pty last)))
	(setq curr-dist-sum (+ (dist xp yp xf yf) (dist xp yp xl yl))))
      (cond ((> curr-dist-sum max-dist-sum-on-this-contour)
	     (setq max-dist-sum-on-this-contour curr-dist-sum nth n))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; This following code and global variables are for tracing contours.

;;; DESTRUCTIVE-CURVE-TRACE-ONE-DIRECTION: With regard to the contour length, this function
;;; assumes that the STARTING POINT (x .  y) is a contour point, i.e., it gets counted in the
;;; length.  Given a 1b-array that is all zeros except for some contours (which are 1's), this
;;; function conses up a list of all the points on the contour passing through (x .  y).  Not
;;; really.  Only those points found on the contour by moving out from x-y in ONE DIRECTION.
;;; You have to go back to (x . y) later and trace it in the other direction.
;;;
;;; This particular curve-tracing function prefers vertical contours (given a forked
;;; contour), since it was designed for use on the vertical edges stereo project.  See the
;;; brief comments to the right of the code.
(defun DESTRUCTIVE-CURVE-TRACE-ONE-DIRECTION
       (array x y xdim ydim &optional points (length 1))
  (if *showp* (show-1b-raster array))
  (if *pausep* (pause))
  (cond ((and (within-raster-bounds x (1+ y) xdim ydim)
	      (= 1 (aref array x (1+ y))))  ; first search straight down...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array x (1+ y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1+ x) (1+ y) xdim ydim)
	      (= 1 (aref array (1+ x) (1+ y))))  ; then down and to the right...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1+ x) (1+ y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1- x) (1+ y) xdim ydim)
	      (= 1 (aref array (1- x) (1+ y))))  ; then down and to the left...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1- x) (1+ y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1+ x) y xdim ydim)
	      (= 1 (aref array (1+ x) y)))  ; then to the right...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1+ x) y xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1- x) y xdim ydim)
	      (= 1 (aref array (1- x) y)))  ; then to the left...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1- x) y xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds x (1- y) xdim ydim)
	      (= 1 (aref array x (1- y))))  ; then search straight up...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array x (1- y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1+ x) (1- y) xdim ydim)
	      (= 1 (aref array (1+ x) (1- y))))   ; then up and to the right...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1+ x) (1- y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((and (within-raster-bounds (1- x) (1- y) xdim ydim)
	      (= 1 (aref array (1- x) (1- y))))   ; then up and to the left...
	 (setf (aref array x y) 0)
	 (destructive-curve-trace-one-direction
	   array (1- x) (1- y) xdim ydim (cons (make-point x y) points) (1+ length)))
	((= 1 (aref array x y))
	 (setf (aref array x y) 0)
	 (values (cons (make-point x y) points) length))
	(t (values points (1- length)))))

;(defvar TEST-ARRAY (make-array '(256 256) :element-type '(mod 2)))
;(defun TEST (&optional (diameter 5))
;  (vu:zero-raster test-array)
;  (loop for x from 100 to (+ 100 diameter -1) do
;    (loop for y from 100 to (+ 100 diameter -1) do
;      (setf (aref test-array x y) 1)))
;  (collect-all-contours test-array))
  
(defun DESTRUCTIVE-CURVE-TRACE (array x y xdim ydim)
  (let (curve1 length1 curve2 length2)
    (multiple-value-setq (curve1 length1)
      (destructive-curve-trace-one-direction array x y xdim ydim))
    (multiple-value-setq (curve2 length2)
      (destructive-curve-trace-one-direction array x y xdim ydim))
    (values (append curve1 (cdr (reverse curve2))) (+ length1 (1- length2)))))

(defvar *TRACING-ARRAY* (make-array '(256 256) :element-type '(mod 2)))

(defvar *CONTOURS* nil)

;;; COLLECT-ALL-CONTOURS is designed to prefer vertical contours (since it was written as part
;;; of the vertical edges stereo project).  What makes it prefer vertical contours is the fact
;;; that the areffing is column-major, i.e., the main loop scans the array like a normal TV
;;; scan (left to right, not top to bottom).  This areffing pattern prevents
;;; DESTRUCTIVE-CURVE-TRACE from breaking vertical contours in the middle.
(defun COLLECT-ALL-CONTOURS (array &optional (min-length 2))
  (let ((xdim (array-dimension array 0))
	(ydim (array-dimension array 1))
	contour
	length
	contours)
    (if (or (neq (array-dimension *tracing-array* 0) xdim)
	    (neq (array-dimension *tracing-array* 1) ydim))
	(setq *tracing-array* (make-array (list xdim ydim) :element-type '(mod 2))))
    (move-array array *tracing-array*)
    (loop for y from 0 below ydim 
	  do  
      ;; column-major areffing (left-to-right scanning) is essential 
      ;; for detecting vert. cont.
      (loop for x from 0 below xdim do
	(cond ((= 1 (aref *tracing-array* x y))
        (multiple-value-setq (contour length)
	  (destructive-curve-trace *tracing-array* x y xdim ydim))
	(if (>= length min-length)
	    (push contour contours))))))
    (setq *contours* contours))
  nil)

;;; SPECIAL-NEIGHBORS takes a 1b-array that has a bunch of one-pixel-wide contours in it, and
;;; it returns an array that has much the same contours.  The difference between the old
;;; contours and the new ones is that in the new array, there are no right-anlge turns in the
;;; contours.  It also attempts to do a little bit of vertical contour extension if you ask it
;;; to.  Summary:
;;;
;;;      X       X            X          X
;;;   X X  =>   X   ,    X X X  =>  X X   ,  etc,  i.e., it makes sure that all the 
;;;   X       X                              new contours are entirely "diagonal" in
;;;                                          character.
;;;
;;;  Why?
;;;
;;;  Because the contour follower below barfs on arrays that don't have this property.  It is
;;;  recursive and it just spins forever.  Therefore is is important to run this function on
;;;  the contour-image before running COLLECT-ALL-CONTOURS.
(defun SPECIAL-NEIGHBORS
       (contour-image
	 &optional (connect-vert? t) 
	 (out-arr (make-array (array-dimensions contour-image) :element-type '(mod 2)))
	 (xdim (array-dimension contour-image 0))
	 (ydim (array-dimension contour-image 1)))
  (let (pixel left right above below)
    (copy-array-contents contour-image out-arr)
    (loop for x from 0 below xdim do  
      (setf (aref out-arr x 0) 0)
      (setf (aref out-arr x 1) 0)
      (setf (aref out-arr x (1- ydim)) 0)
      (setf (aref out-arr x (- ydim 2)) 0))
    (loop for y from 0 below ydim do 
      (setf (aref out-arr 0 y) 0)
      (setf (aref out-arr 1 y) 0)
      (setf (aref out-arr (1- xdim) y) 0)
      (setf (aref out-arr (- xdim 2) y) 0))
    (loop for x from 1 below (1- xdim) do
      (loop for y from 1 below (1- ydim) do
	(setq pixel (= 1 (aref out-arr x y))
	      left (= 1 (aref out-arr (1- x) y))
	      right (= 1 (aref out-arr (1+ x) y))
	      above (= 1 (aref out-arr x (1- y)))
	      below (= 1 (aref out-arr x (1+ y))))
	(if (not pixel)
	    nil
	    (if (or (and left (or above below))
		    (and right (or above below))
		    (and above (or left right))
		    (and below (or left right)))
		(setf (aref out-arr x y) 0)))
	(if connect-vert?
	    (if (and (not pixel)
		     (not left) (not right) above below
		     (= 1 (aref out-arr x (- y 2)))
		     (= 1 (aref out-arr x (+ y 2))))
		(setf (aref out-arr x y) 1))))))
  out-arr)