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

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

;;;*****************************************************************************
;;; CHANGE HISTORY
;;;
;;; 9/25/87  WEG  Changed the package to mu.  Recompiled for Release 7.
;;; 92MAR09  PAO Ported to Lucid.
;;;*****************************************************************************

(in-package :mu)

;;; A simple calibration scheme using a test pattern composed of dots.

;;; now included in the system mike-utils7
;;;(load "b:>mike>graphics>draw-line")
;;;(load "b:>mike>image>1b-region-label")

(proclaim '(special *left-256* *right-256*))

(defvar *MIN-X-WARP* 0)
(defvar *MAX-X-WARP* 0)

(defvar *MIN-Y-WARP* 0)
(defvar *MAX-Y-WARP* 0)

(defvar *X-WARP-ARRAY* (make-array '(256 256) :element-type '(mod 256.)))
(defvar *Y-WARP-ARRAY* (make-array '(256 256) :element-type '(mod 256.)))

;;; IMAGE!! is a FIELD-PVAR.
;;;
;;; This function is analogous to the "expand" operation found in morphological image
;;; processing.  It takes a pvar (IMAGE!!) and fills in all pixels that are zero.  It does so
;;; by repeatedly setting each pixel that is zero to be the average of its NONZERO neighbors.
;;; If EIGHT-CONNECTED? is T, then zero-pixels looks at all eight of their neighbors.  If
;;; EIGHT-CONNECTED? is nil, then zero-pixels look only at their 4 orthogonal neighbors.
;;;
;;; This function is used for filling in after warping.  (Sometimes "rips and tears" develop
;;; after using a naive warping function.) Therefore, the IMAGE!! that is being warped should
;;; first have a constant added to it to make sure none of its pixels are zero, THEN it
;;; should be warped.  This will insure that the result of warping has zeros only where tears
;;; occured.
(*defunc FILL-IN-HOLES!! (image!! &optional (times 1))
  (format t "~%Filling in holes...")
  (dotimes (i times)
    (*let ((result!! image!!)
	   (sum-of-neighbors!! (!! 0))
	   (num-of-nonzero-neighbors!! (!! 0)))
      (declare (type (field-pvar (pvar-length image!!)) result!!)
	       (type (field-pvar (+ (pvar-length image!!) 3)) sum-of-neighbors!!)
	       ;; (type boolean-pvar a-neighbor-was-nonzero)
	       )
      (loop for n in '(0 1 2 3 4 5 6 7) do
	(*let ((neighbor!! (nearest-neighbor!! result!! n)))
	  (*incf sum-of-neighbors!! neighbor!!)
	  (*if (and!! (zerop!! result!!) (nzerop!! neighbor!!))
	       (*incf num-of-nonzero-neighbors!!))))
      (*if (nzerop!! num-of-nonzero-neighbors!!)
	   (*set result!! (floor!! sum-of-neighbors!! num-of-nonzero-neighbors!!)))
      (format t "done.")
      (return result!!))))

;;; Used only by warp-with-router!! (see below).  Given an unsigned array, read it into an unsigned pvar,
;;; then superimpose a signed pvar.  The unsigned array holds signed numbers.
(*defunc UNSIGNED-TO-SIGNED!! (array)
  (cmv::with-temp-pvars
    (let* ((unsigned-pvar (cmv::write-raster-to-cm-grid!! array))
	   (signed-pvar (*lx:allocate-temp-pvar
			  :length (pvar-length unsigned-pvar)
			  :location (pvar-location unsigned-pvar)
			  :type :signed)))
      signed-pvar)))

;;; IMAGE!! is an unsigned pvar.
;;; X-WARP-FUNC is a pvar of integers telling each pixel in the right image how much it should be displaced
;;; in the X-direction.
;;; Y-WARP-FUNC is a pvar of integers telling each pixel in the right image how much it should be displaced
;;; in the Y-direction.
;;; X-WARP-FUNC and Y-WARP-FUNC can be signed pvars or unsigned arrays.  The unsigned arrays hold signed
;;; quantities, so we have to be sure to interpret the contents correctly.  (WEG 4/5/88)
;;;
;;; The number one is added to the image!! before warping so that result is guaranteed to
;;; contain zeros only in places where the image has been "ripped" by the warping process.
;;;
;;; Note: after warping and filling in holes, the right image is zero at locations which were left empty by the
;;; warping.  The stereo code takes advantage of this feature, so don't change this routine without examining
;;; the effect on the stereo code.  The relevant parameter in the stereo code is
;;; *suppress-matches-outside-of-right-image-boundary?*.  (WEG 2/18/88)

(*defunc WARP-WITH-ROUTER!!
	 (image!! &optional x-warp-func y-warp-func)
  (format t "~%Warping with router...")
  (cmv::with-temp-pvars
    (let ((x-warp-func!! (if (pvarp x-warp-func) x-warp-func (unsigned-to-signed!! x-warp-func)))
	  (y-warp-func!! (if (pvarp y-warp-func) y-warp-func (unsigned-to-signed!! y-warp-func))))
      (*let* ((1+image!! (1+!! image!!))   ;;; <== Make sure there are no zeros in the image to be warped.
	      (new-x!! (clamp!! (+!! (x!!) x-warp-func!!) (!! 0) (!! (1- (dimension-size 0)))))
	      (new-y!! (clamp!! (+!! (y!!) y-warp-func!!) (!! 0) (!! (1- (dimension-size 1)))))
	      (new-cube!! (cube-from-grid-address!! new-x!! new-y!!))
	      (shifted-1+image!! (ash!! 1+image!! (!! 5)))
	      (temp!! (!! 1))
	      (message-with-0001-in-low-bit!!
		(deposit-byte!! shifted-1+image!! (!! 0) (!! 5) temp!!))
	      (inbox!! (!! 0)))
	;;; MESSAGE-WITH-0001-IN-LOW-BIT!! looks like this:
	;;; bits:  MSB . . . . . . . . . . . . . LSB
	;;;        <  1+IMAGE!!  ASHED  BY  5> 0 0 0 1
	;;; 0001 because I assume a maximum of (2^5 - 1) collisions
	(declare (type (field-pvar (1+ (pvar-length image!!))) 1+image!!)
		 (type (field-pvar 1) temp!!)
		 (type (field-pvar (+ 5 (pvar-length 1+image!!))) shifted-1+image!! message-with-0001-in-low-bit!!)
		 (type (field-pvar (+ 5 (pvar-length message-with-0001-in-low-bit!!))) inbox!!)
		 (type (field-pvar (dimension-address-length 0)) new-x!!)
		 (type (field-pvar (dimension-address-length 1)) new-y!!)
		 (type (field-pvar *current-send-address-length*) new-cube!!))
	(*PSET :ADD MESSAGE-WITH-0001-IN-LOW-BIT!! INBOX!! NEW-CUBE!!)
	(*let ((no-of-collisions!! (load-byte!! inbox!! (!! 0) (!! 5)))
	       (warped-image!! (load-byte!! inbox!! (!! 5) (!! (- (pvar-length inbox!!) 5)))))
	  (declare (type (field-pvar 5) no-of-collisions!!)
		   (type (field-pvar (- (pvar-length inbox!!) 5)) warped-image!!))
	  (*if (and!! (nzerop!! warped-image!!) (nzerop!! no-of-collisions!!))
	       (*set warped-image!! (round!! warped-image!! no-of-collisions!!)))
	  ;; Use this instead of (*decf warped-image!!) to make sure you get no negative number problems.
	  (*set warped-image!! (max!! (!! 0) (1-!! warped-image!!)))
	  (format t "done...")
	  (load-byte!! (fill-in-holes!! warped-image!! 1) (!! 0) (!! (pvar-length image!!))))))))

;;; commented out until I fix problem: array-pvar!! is obsolete
;(*defunc *GRAB-AND-REGISTER (&optional left!! right!! show? (show-x 0) (show-y 0))
;  (let ((size (dimension-size 0)))
;    (grab-stereo-pair-twiddle)
;    (move-frames size nil nil)
;    (*set left!! (transpose!! (array-pvar!! (left-array-of-size size))))
;    (*set right!! (warp-with-router!! (transpose!! (array-pvar!! (right-array-of-size size)))
;				      (transpose!! (array-pvar!! *x-warp-array* t))
;				      (transpose!! (array-pvar!! *y-warp-array* t))))
;    (cond (show?
;	   (*show-pvar left!! :x show-x :y show-y :transpose t)
;	   (*show-pvar right!! :x (+ show-x size) :y show-y :transpose t)))))

;;;  This file depends on a:>mike>graphics>draw-line.lisp

(defvar *LEFT-WARP-1B* (make-array '(256 256) :element-type '(mod 2)))
(defvar *RIGHT-WARP-1B* (make-array '(256 256) :element-type '(mod 2)))

(defvar *X-WARP-ART-Q* (make-array '(256 256)))
(defvar *Y-WARP-ART-Q* (make-array '(256 256)))

;;;  This function does a square-spiral search around each dot in the warping image until it
;;;  finds a dot in the stationary image.  Whenever a warping-dot finds its nearest-neighboring stationary-dot,
;;;  it writes the coords of the stationary-dot into X-ARRAY and Y-ARRAY at the
;;;  position of the warping-dot (NOT at the position of the stationary dot that was just found
;;;  by the warping dot).
(defun MATCH-UP-WITH-NEAREST-DOTS
       (stationary-dots warping-dots x-array y-array max-search-radius)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array stationary-dots)
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(if (= 1 (vu:raster-aref warping-dots x y))
	    (multiple-value-bind (lx ly)
		(square-spiral-search stationary-dots x y max-search-radius)
	      (setf (vu:raster-aref x-array x y) lx (vu:raster-aref y-array x y) ly)))))))

;;; Given the point (START-X, START-Y) in the array DOT-IMAGE, this function returns the RELATIVE coordinates of the nearest
;;; element containing a 1, via a spiral search.  Returns NIL if no such neighbor is found.
(defun SQUARE-SPIRAL-SEARCH (dot-image start-x start-y max-radius)
  (if (= 1 (vu:raster-aref dot-image start-x start-y))
      (values 0 0)
      (let (answer)
	(dotimes (i max-radius)
	  (if (setq answer (sq-sp-search-at-particular-radius dot-image start-x start-y (1+ i)))
	      (return answer)))
	(values (car answer) (cadr answer)))))

(defun SQ-SP-SEARCH-AT-PARTICULAR-RADIUS (dot-image cx cy radius)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array dot-image)
    (or (let ((y (+ cy radius)))
	  (dotimes (i (1+ (* 2 radius)))
	    (let ((realx (+ cx radius (- i))))
	      (if (and (not (or (minusp realx) (minusp y) (>= realx xdim) (>= y ydim)))
		       (= 1 (vu:raster-aref dot-image realx y)))
		  (return (list (- realx cx) (- y cy)))))))
	(let ((x (- cx radius)))
	  (dotimes (i (* 2 radius))
	    (let ((realy (+ cy radius -1 (- i))))
	      (if (and (not (or (minusp x) (minusp realy) (>= x xdim) (>= realy ydim)))
		       (= 1 (vu:raster-aref dot-image x realy)))
		  (return (list (- x cx) (- realy cy)))))))
	(let ((y (- cy radius)))
	  (dotimes (i (* 2 radius))
	    (let ((realx (+ cx (- radius) 1 i)))
	      (if (and (not (or (minusp realx) (minusp y) (>= realx xdim) (>= y ydim)))
		       (= 1 (vu:raster-aref dot-image realx y)))
		  (return (list (- realx cx) (- y cy)))))))
	(let ((x (+ cx radius)))
	  (dotimes (i (1- (* 2 radius)))
	    (let ((realy (+ cy (- radius) 1 i)))
	      (if (and (not (or (minusp x) (minusp realy) (>= x xdim) (>= realy ydim)))
		       (= 1 (vu:raster-aref dot-image x realy)))
		  (return (list (- x cx) (- realy cy))))))))))

;;;  This function starts at the point (CX,CY) in the array DOT-IMAGE and searches for the nearest non-zero element in the
;;;  positive y-direction in the array.  It searchs in an y-ward-opening triangular region having half-angle SEARCH-HALF-ANGLE.
;;;  It returns either a list (x,y) of the coordinates of the non-zero point that it or found, or NIL if it went off the edge
;;;  of the array during the search.
;;;  
;;;  Example:  If you called (FIND-NEAREST-NEIGHBOR-POSY <array> 6 3 10 45), it would return (7,6), because there is a 1 at the
;;;  point (7,6) in <array>.  (The *'s mark the area search by the function, i.e., a 45-degree-half-angle cone.)
;;;
;;;   y^             
;;;    | 0 * * * * * * * * * * * 
;;;    | 0 0 * * * * * * * * * 0 
;;;  6 | 0 0 0 * * * * 1 * * 0 0 
;;;    | 0 0 0 0 * * * * * 0 0 0 
;;;    | 0 0 0 0 0 * * * 0 0 0 0   <---  <array>
;;;  3 | 0 0 0 0 0 0 1 0 0 0 0 0 
;;;    | 0 0 0 0 0 0 0 0 0 0 0 0 
;;;    | 0 0 0 0 0 0 0 0 0 0 0 0 
;;;    | 0 0 0 0 0 0 0 0 0 0 0 0 
;;;  (o,o)----------------------->x
;;;                  6 7
(defun FIND-NEAREST-NEIGHBOR-POSY (dot-image cx cy search-radial-extent search-half-angle)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array dot-image)
    (let ((ratio (tan (/ search-half-angle (/ 360 2 3.1415927)))))
      (dotimes (dy search-radial-extent)
	(let* ((finish-dx (round (* ratio (1+ dy))))
	       (start-dx (- finish-dx))
	       found-x)
	  (if (dotimes (dx (1+ (- finish-dx start-dx)))
		(let ((realx (+ cx start-dx dx))
		      (realy (+ cy 1 dy)))
		  (cond ((or (minusp realx) (minusp realy)
			     (>= realx xdim) (>= realy ydim))  ;;; Went off the edge, return nil.
			 (setq found-x nil)) 
			((not (zerop (vu:raster-aref dot-image realx realy)))
			 (return (setq found-x realx))))))
	      (return (if (numberp found-x)
			  (list found-x (+ cy 1 dy))
			  nil))))))))

;;;  This function is completely analogous to FIND-NEAREST-NEIGHBOR-POSY, except it searches
;;;  in the POSITIVE X-DIRECTION instead of in the postive y-direction.
(defun FIND-NEAREST-NEIGHBOR-POSX (dot-image cx cy search-radial-extent search-half-angle)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array dot-image)
    (let ((ratio (tan (/ search-half-angle (/ 360 2 3.1415927)))))
      (dotimes (dx search-radial-extent)
	(let* ((finish-dy (round (* ratio (1+ dx))))
	       (start-dy (- finish-dy))
	       found-y)
	  (if (dotimes (dy (1+ (- finish-dy start-dy)))
		(let ((realx (+ cx 1 dx))
		      (realy (+ cy start-dy dy)))
		  (cond ((or (minusp realx) (minusp realy)
			     (>= realx xdim) (>= realy ydim))  ;;; Went off the edge, return nil.
			 (setq found-y nil))
			((not (zerop (vu:raster-aref dot-image realx realy)))
			 (return (setq found-y realy))))))
	      (return (if (numberp found-y)
			  (list (+ cx 1 dx) found-y)
			  nil))))))))

;;;  This function is completely analogous to FIND-NEAREST-NEIGHBOR-POSY, except it searches
;;;  in the NEGATIVE y-direction instead of in the positive y-direction.
(defun FIND-NEAREST-NEIGHBOR-NEGY (dot-image cx cy search-radial-extent search-half-angle)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array dot-image)
    (let ((ratio (tan (/ search-half-angle (/ 360 2 3.1415927)))))
      (dotimes (dy search-radial-extent)
	(let* ((finish-dx (round (* ratio (1+ dy))))
	       (start-dx (- finish-dx))
	       found-x)
	  (if (dotimes (dx (1+ (- finish-dx start-dx)))
		(let ((realx (+ cx start-dx dx))
		      (realy (- cy 1 dy)))
		  (cond ((or (minusp realx) (minusp realy)
			     (>= realx xdim) (>= realy ydim))  ;;; Went off the edge, return nil.
			 (setq found-x nil)) 
			((not (zerop (vu:raster-aref dot-image realx realy)))
			 (return (setq found-x realx))))))
	      (return (if (numberp found-x)
			  (list found-x (- cy 1 dy))
			  nil))))))))

;;;  This function is completely analogous to FIND-NEAREST-NEIGHBOR-POSX, except it searches
;;;  in the NEGATIVE x-direction instead of in the postive x-direction.
(defun FIND-NEAREST-NEIGHBOR-NEGX (dot-image cx cy search-radial-extent search-half-angle)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array dot-image)
    (let ((ratio (tan (/ search-half-angle (/ 360 2 3.1415927)))))
      (dotimes (dx search-radial-extent)
	(let* ((finish-dy (round (* ratio (1+ dx))))
	       (start-dy (- finish-dy))
	       found-y)
	  (if (dotimes (dy (1+ (- finish-dy start-dy)))
		(let ((realx (- cx 1 dx))
		      (realy (+ cy start-dy dy)))
		  (cond ((or (minusp realx) (minusp realy)
			     (>= realx xdim) (>= realy ydim))  ;;; Went off the edge, return nil.
			 (setq found-y nil))
			((not (zerop (vu:raster-aref dot-image realx realy)))
			 (return (setq found-y realy))))))
	      (return (if (numberp found-y)
			  (list (- cx 1 dx) found-y)
			  nil))))))))

(defvar *X-ARRAY-FOR-LINEAR-INTERPOLATION* (make-array 200.))
(defvar *Y-ARRAY-FOR-LINEAR-INTERPOLATION* (make-array 200.))
(defun LINEARLY-INTERPOLATE-SEGMENT-BETWEEN-TWO-POINTS (array x1 y1 x2 y2 &optional (round? t))
  (let ((start-val (vu:raster-aref array x1 y1))
	(top-index (draw-line-coordinates-into-1d-arrays
		     *x-array-for-linear-interpolation*
		     *y-array-for-linear-interpolation*
		     x1 y1 x2 y2))
	(slope (/ (- (vu:raster-aref array x2 y2) (vu:raster-aref array x1 y1))
		  (sqrt (+ (square (- x2 x1)) (square (- y2 y1)))))))
    (loop for i from 1 to (1- top-index) do
      (let ((ix (aref *x-array-for-linear-interpolation* i))
	    (iy (aref *y-array-for-linear-interpolation* i)))
	(setf (vu:raster-aref array ix iy)
	      (let ((val (+ start-val (* slope (sqrt (+ (square (- ix x1)) (square (- iy y1))))))))
		(if round?
		    (round val)
		    val)))))))
  
(defvar *TEMP-1B-ARRAY* (make-array '(256 256) :element-type '(mod 2)))
(defun CONNECT-THE-DOTS-WITH-LINEARLY-INTERPOLATED-SEGMENTS
       (1b-guide-dots art-q-array search-radial-extent search-half-angle output-1b)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array art-q-array)
    #+symbolics
    (zl-user:bitblt tv:alu-ior xdim ydim 1b-guide-dots 0 0 output-1b 0 0)  ;;; must be ior!
    #-symbolics
    (error "Conversion from Symbolics-specific code not complete.")
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(cond ((= 1 (vu:raster-aref 1b-guide-dots x y))
	       (let ((posy-neighb (find-nearest-neighbor-posy 1b-guide-dots x y search-radial-extent search-half-angle))
		     (posx-neighb (find-nearest-neighbor-posx 1b-guide-dots x y search-radial-extent search-half-angle)))
		 (cond (posy-neighb
			(linearly-interpolate-segment-between-two-points
			  output-1b x y (car posy-neighb) (cadr posy-neighb) t)
			(linearly-interpolate-segment-between-two-points
			  art-q-array x y (car posy-neighb) (cadr posy-neighb) t)))
		 (cond (posx-neighb
			(linearly-interpolate-segment-between-two-points
			  output-1b x y (car posx-neighb) (cadr posx-neighb) t)
			(linearly-interpolate-segment-between-two-points
			  art-q-array x y (car posx-neighb) (cadr posx-neighb) t)))))))))
  output-1b)

(defun LIN-INTERP-TWO-PTS-OUT-TO-EDGE-OF-ARRAY (array xa ya xb yb)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array array)
    (let* ((xdim-1 (1- xdim))
	   (ydim-1 (1- ydim))
	   (xdim2 (floor xdim 2))
	   (ydim2 (floor ydim 2))
	   (pa-to-ctr (sqrt (+ (square (- xa xdim2)) (square (- ya ydim2)))))
	   (pb-to-ctr (sqrt (+ (square (- xb xdim2)) (square (- yb ydim2)))))
	   (pb-is-closer-to-edge (> pb-to-ctr pa-to-ctr))
	   (x1 (if pb-is-closer-to-edge xa xb))
	   (y1 (if pb-is-closer-to-edge ya yb))
	   (x2 (if pb-is-closer-to-edge xb xa))
	   (y2 (if pb-is-closer-to-edge yb ya))
	   x3 y3)
      (cond ((and (= y1 y2) (> x2 x1))
	     (setq x3 xdim-1 y3 y2))
	    ((and (= y1 y2) (< x2 x1))
	     (setq x3 0 y3 y2))
	    ((and (= x1 x2) (> y2 y1))
	     (setq x3 x2 y3 ydim-1))
	    ((and (= x1 x2) (< y2 y1))
	     (setq x3 x2 y3 0))
	    ((and (> x2 x1) (> y2 y1))
	     (let ((s (/ (- y2 y1) (- x2 x1))))
	       (if (> s 1)
		   (setq x3 (+ x2 (round (/ (- ydim-1 y2) s))) y3 ydim-1)
		   (setq x3 xdim-1 y3 (+ y2 (round (* (- xdim-1 x2) s)))))))
	    ((and (> x2 x1) (< y2 y1))
	     (let ((s (/ (- y1 y2) (- x2 x1))))
	       (if (> s 1)
		   (setq x3 (+ x2 (round (/ y2 s))) y3 0)
		   (setq x3 xdim-1 y3 (- y2 (round (* (- xdim-1 x2) s)))))))
	    ((and (< x2 x1) (> y2 y1))
	     (let ((s (/ (- y2 y1) (- x1 x2))))
	       (if (> s 1)
		   (setq x3 (- x2 (round (/ (- ydim-1 y2) s))) y3 ydim-1)
		   (setq x3 0 y3 (+ y2 (round (* x2 s)))))))
	    ((and (< x2 x1) (< y2 y1))
	     (let ((s (/ (- y1 y2) (- x1 x2))))
	       (if (> s 1)
		   (setq x3 (- x2 (round (/ y2 s))) y3 0)
		   (setq x3 0 y3 (- y2 (round (* x2 s))))))))
      (let* ((d12 (sqrt (+ (square (- x1 x2)) (square (- y1 y2)))))
	     (d23 (sqrt (+ (square (- x2 x3)) (square (- y2 y3)))))
	     (v1 (vu:raster-aref array x1 y1))
	     (v2 (vu:raster-aref array x2 y2))
	     (v3 (round (+ v2 (* d23 (/ (float (- v2 v1)) d12))))))
	(setf (vu:raster-aref array x3 y3) v3)
	(linearly-interpolate-segment-between-two-points array x2 y2 x3 y3)))))

(defun EXTEND-LINES-TOWARD-NEGX-EDGES
       (1b-guide-dots art-q-array search-radial-extent search-half-angle output-1b)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array 1b-guide-dots)
    #+symbolics
    (zl-user:bitblt tv:alu-ior xdim ydim 1b-guide-dots 0 0 output-1b 0 0)  ;;; must be ior!
    #-symbolics
    (error "Conversion from Symbolics-specific code not complete...")
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(if (= 1 (vu:raster-aref 1b-guide-dots x y))
	    (let ((negx-neighb (find-nearest-neighbor-negx 1b-guide-dots x y search-radial-extent search-half-angle))
		  (posx-neighb (find-nearest-neighbor-posx 1b-guide-dots x y search-radial-extent search-half-angle)))
	      (cond ((and (not negx-neighb)
			  posx-neighb)
		     (lin-interp-two-pts-out-to-edge-of-array
		       output-1b x y (car posx-neighb) (cadr posx-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       art-q-array x y (car posx-neighb) (cadr posx-neighb))))))))))
		     
(defun EXTEND-LINES-TOWARD-POSX-EDGES
       (1b-guide-dots art-q-array search-radial-extent search-half-angle output-1b)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array 1b-guide-dots)
    #+symbolics
    (zl-user:bitblt tv:alu-ior xdim ydim 1b-guide-dots 0 0 output-1b 0 0)  ;;; must be ior!
    #-symbolics
    (error "Conversion from Symbolics-specific code not complete...")
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(if (= 1 (vu:raster-aref 1b-guide-dots x y))
	    (let ((negx-neighb (find-nearest-neighbor-negx 1b-guide-dots x y search-radial-extent search-half-angle))
		  (posx-neighb (find-nearest-neighbor-posx 1b-guide-dots x y search-radial-extent search-half-angle)))
	      (cond ((and negx-neighb
			  (not posx-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       output-1b x y (car negx-neighb) (cadr negx-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       art-q-array x y (car negx-neighb) (cadr negx-neighb))))))))))

(defun EXTEND-LINES-TOWARD-NEGY-EDGES
       (1b-guide-dots art-q-array search-radial-extent search-half-angle output-1b)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array 1b-guide-dots)
    #+symbolics
    (zl-user:bitblt tv:alu-ior xdim ydim 1b-guide-dots 0 0 output-1b 0 0)  ;;; must be ior!
    #-symbolics
    (error "Conversion from Symbolics-specific code not complete...")
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(if (= 1 (vu:raster-aref 1b-guide-dots x y))
	    (let ((negy-neighb (find-nearest-neighbor-negy 1b-guide-dots x y search-radial-extent search-half-angle))
		  (posy-neighb (find-nearest-neighbor-posy 1b-guide-dots x y search-radial-extent search-half-angle)))
	      (cond ((and (not negy-neighb)
			  posy-neighb)
		     (lin-interp-two-pts-out-to-edge-of-array
		       output-1b x y (car posy-neighb) (cadr posy-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       art-q-array x y (car posy-neighb) (cadr posy-neighb))))))))))

(defun EXTEND-LINES-TOWARD-POSY-EDGES
       (1b-guide-dots art-q-array search-radial-extent search-half-angle output-1b)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array 1b-guide-dots)
    #+symbolics
    (zl-user:bitblt tv:alu-ior xdim ydim 1b-guide-dots 0 0 output-1b 0 0)  ;;; must be ior!
    #-symbolics
    (error "Conversion from Symbolics-specific code not complete...")
    (loop for y from 0 below ydim do
      (loop for x from 0 below xdim do
	(if (= 1 (vu:raster-aref 1b-guide-dots x y))
	    (let ((negy-neighb (find-nearest-neighbor-negy 1b-guide-dots x y search-radial-extent search-half-angle))
		  (posy-neighb (find-nearest-neighbor-posy 1b-guide-dots x y search-radial-extent search-half-angle)))
	      (cond ((and negy-neighb
			  (not posy-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       output-1b x y (car negy-neighb) (cadr negy-neighb))
		     (lin-interp-two-pts-out-to-edge-of-array
		       art-q-array x y (car negy-neighb) (cadr negy-neighb))))))))))

(defun EXTEND-ALL-EXTERIOR-LINES
       (1b-warping-dots &optional (x-array *x-warp-art-q*) (y-array *y-warp-art-q*)
	(search-radial-extent 80) (search-half-angle 20) (output-1b *temp-1b-array*))
  ;;
  (extend-lines-toward-negy-edges 1b-warping-dots x-array search-radial-extent search-half-angle output-1b)
  (extend-lines-toward-negy-edges 1b-warping-dots y-array search-radial-extent search-half-angle output-1b)
  ;;
  (extend-lines-toward-posy-edges 1b-warping-dots x-array search-radial-extent search-half-angle output-1b)
  (extend-lines-toward-posy-edges 1b-warping-dots y-array search-radial-extent search-half-angle output-1b)
  ;;
  (extend-lines-toward-negx-edges 1b-warping-dots x-array search-radial-extent search-half-angle output-1b)
  (extend-lines-toward-negx-edges 1b-warping-dots y-array search-radial-extent search-half-angle output-1b)
  ;;
  (extend-lines-toward-posx-edges 1b-warping-dots x-array search-radial-extent search-half-angle output-1b)
  (extend-lines-toward-posx-edges 1b-warping-dots y-array search-radial-extent search-half-angle output-1b)
  ;;
  output-1b)

(defun SET-CORNER-0-0 (array)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array array)
    (let ((a nil) (b nil))
      (let* ((first-top (dotimes (x xdim)
			  (cond ((numberp (vu:raster-aref array x 0))
				 (setq a x)
				 (return (vu:raster-aref array x 0))))))
	     (second-top (dotimes (x (- xdim a))
			   (cond ((numberp (vu:raster-aref array (+ a x 1) 0))
				  (setq b (+ a x 1))
				  (return (vu:raster-aref array (+ a x 1) 0))))))
	     (value-due-to-top (+ first-top (* (- (/ (- second-top first-top)
						     (- b a)))
					       a)))
	     (first-left (dotimes (y ydim)
			   (cond ((numberp (vu:raster-aref array 0 y))
				  (setq a y)
				  (return (vu:raster-aref array 0 y))))))
	     (second-left (dotimes (y (- ydim a))
			    (cond ((numberp (vu:raster-aref array 0 (+ a y 1)))
				   (setq b (+ a y 1))
				   (return (vu:raster-aref array 0 (+ a y 1)))))))
	     (value-due-to-left (+ first-left (* (- (/ (- second-left first-left)
						       (- b a)))
						 a))))
	(setf (vu:raster-aref array 0 0) (round (+ value-due-to-top value-due-to-left) 2))))))
  
(defun SET-CORNER-0-1 (array)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array array)
    (let* ((ydim-1 (1- ydim))
	   a b)
      (let* ((first-bottom (dotimes (x xdim)
			     (cond ((numberp (vu:raster-aref array x ydim-1))
				    (setq a x)
				    (return (vu:raster-aref array x ydim-1))))))
	     (second-bottom (dotimes (x (- xdim a))
			      (cond ((numberp (vu:raster-aref array (+ a x 1) ydim-1))
				     (setq b (+ a x 1))
				     (return (vu:raster-aref array (+ a x 1) ydim-1))))))
	     (value-due-to-bottom (+ first-bottom (* (- (/ (- second-bottom first-bottom)
							   (- b a)))
						     a)))
	     (last-left (dotimes (i ydim)
			  (cond ((numberp (vu:raster-aref array 0 (- ydim-1 i)))
				 (setq a i)
				 (return (vu:raster-aref array 0 (- ydim-1 i)))))))
	     (next-to-last-left (dotimes (i (- ydim a))
				  (cond ((numberp (vu:raster-aref array 0 (- ydim-1 a i 1)))
					 (setq b (+ a i 1))
					 (return (vu:raster-aref array 0 (- ydim-1 a i 1)))))))
	     (value-due-to-left (+ last-left (* (- (/ (- next-to-last-left last-left)
						      (- b a)))
						a))))
	(setf (vu:raster-aref array 0 ydim-1) (round (+ value-due-to-bottom value-due-to-left) 2))))))

(defun SET-CORNER-1-0 (array)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array array)
    (let* ((xdim-1 (1- xdim))
	   a b)
      (let* ((last-top (dotimes (i xdim)
			 (cond ((numberp (vu:raster-aref array (- xdim-1 i) 0))
				(setq a i)
				(return (vu:raster-aref array (- xdim-1 i) 0))))))
	     (next-to-last-top (dotimes (i (- xdim a))
				 (cond ((numberp (vu:raster-aref array (- xdim-1 a i 1) 0))
					(setq b (+ a i 1))
					(return (vu:raster-aref array (- xdim-1 a i 1) 0))))))
	     (value-due-to-top (+ last-top (* (- (/ (- next-to-last-top last-top)
						    (- b a)))
					      a)))
	     (first-right (dotimes (y ydim)
			    (cond ((numberp (vu:raster-aref array xdim-1 y))
				   (setq a y)
				   (return (vu:raster-aref array xdim-1 y))))))
	     (second-right (dotimes (y (- ydim a))
			     (cond ((numberp (vu:raster-aref array xdim-1 (+ a y 1)))
				    (setq b (+ a y 1))
				    (return (vu:raster-aref array xdim-1 (+ a y 1)))))))
	     (value-due-to-right (+ first-right (* (- (/ (- second-right first-right)
							 (- b a)))
						   a))))
	(setf (vu:raster-aref array xdim-1 0) (round (+ value-due-to-top value-due-to-right) 2))))))



(defun SET-CORNER-1-1 (array)
  (multiple-value-bind (xdim ydim) (vu:decode-raster-array array)
    (let* ((xdim-1 (1- xdim))
	   (ydim-1 (1- ydim))
	   a b)
      (let* ((last-bottom (dotimes (i xdim)
			    (cond ((numberp (vu:raster-aref array (- xdim-1 i) ydim-1))
				   (setq a i)
				   (return (vu:raster-aref array (- xdim-1 i) ydim-1))))))
	     (next-to-last-bottom (dotimes (i (- xdim a))
				    (cond ((numberp (vu:raster-aref array (- xdim-1 a i 1) ydim-1))
					   (setq b (+ a i 1))
					   (return (vu:raster-aref array (- xdim-1 a i 1) ydim-1))))))
	     (value-due-to-bottom (+ last-bottom (* (- (/ (- next-to-last-bottom last-bottom)
							  (- b a)))
						    a)))
	     (last-right (dotimes (i ydim)
			   (cond ((numberp (vu:raster-aref array xdim-1 (- ydim-1 i)))
				  (setq a i)
				  (return (vu:raster-aref array xdim-1 (- ydim-1 i)))))))
	     (next-to-last-right (dotimes (i (- ydim a))
				   (cond ((numberp (vu:raster-aref array xdim-1 (- ydim-1 a i 1)))
					  (setq b (+ a i 1))
					  (return (vu:raster-aref array xdim-1 (- ydim-1 a i 1)))))))
	     (value-due-to-right (+ last-right (* (- (/ (- next-to-last-right last-right)
							(- b a)))
						  a))))
	(setf (vu:raster-aref array xdim-1 ydim-1) (round (+ value-due-to-bottom value-due-to-right) 2))))))

(defun SET-ALL-CORNERS (array)
  (set-corner-0-0 array)
  (set-corner-1-0 array)
  (set-corner-0-1 array)
  (set-corner-1-1 array))

(defun SET-TOP-EDGES (array)
  (horizontal-linear-interpolate-single-row array 0)
  (horizontal-linear-interpolate-single-row array (1- (vu:raster-height array))))

(defun SET-SIDE-EDGES (array)
  (vertical-linear-interpolate-single-column array 0)
  (vertical-linear-interpolate-single-column array (1- (vu:raster-width array))))

;;;  Linearly interpolates a column of an array containing numbers and nil's.  The
;;;  top and bottom elements of the column MUST be numbers, they cannot be nil.
(defun VERTICAL-LINEAR-INTERPOLATE-SINGLE-COLUMN (array x)
  (let* ((ydim (vu:raster-height array))
	 (ydim-1 (1- ydim))
	 (last-y 0)
	 (current-y (1+ last-y)))
    (loop do
      (when (= current-y ydim-1)
	(return (values)))
      (setq current-y (do ((y (1+ last-y) (1+ y)))
			  ((or (= y ydim-1) (vu:raster-aref array x y))
			   y)))
      (let* ((slope (/ (- (float (vu:raster-aref array x current-y))
			  (vu:raster-aref array x last-y))
		       (- current-y last-y)))
	     (last-value (vu:raster-aref array x last-y)))
	(loop for y from last-y to current-y do
	  (setf (vu:raster-aref array x y)
		(+ last-value (round (* slope (- y last-y)))))))
      (setq last-y current-y))))

;;;  Linearly interpolates a row of an array containing numbers and nil's.  The
;;;  leftmost and rightmost elements of the row MUST
;;;  be numbers, they cannot be nil.
(defun HORIZONTAL-LINEAR-INTERPOLATE-SINGLE-ROW (array y)
  (let* ((xdim (vu:raster-width array))
	 (xdim-1 (1- xdim))
	 (last-x 0)
	 (current-x (1+ last-x)))
    (loop do
      (when (= current-x xdim-1)
	(return (values)))
      (setq current-x (do ((x (1+ last-x) (1+ x)))
			  ((or (= x xdim-1) (vu:raster-aref array x y))
			   x)))
      (let* ((slope (/ (- (float (vu:raster-aref array current-x y)) (vu:raster-aref array last-x y))
		       (- current-x last-x)))
	     (last-value (vu:raster-aref array last-x y)))
	(loop for x from last-x to current-x do
	  (setf (vu:raster-aref array x y) (+ last-value (round (* slope (- x last-x)))))))
      (setq last-x current-x))))

(defun VERTICAL-INTERPOLATE (array)
  (loop for x from 0 below (vu:raster-width array) do
    (vertical-linear-interpolate-single-column array x))
  array)

(defun HORIZONTAL-INTERPOLATE (array)
  (loop for y from 0 below (vu:raster-height array) do
    (horizontal-linear-interpolate-single-row array y))
  array)

(defun MAKE-SURE-EVERY-POINT-IN-RIGHT-IMAGE-GOT-A-MATCH-IN-LEFT (1b-array art-q-array)
  #-symbolics (declare (ignore 1b-array art-q-array))
  #-symbolics (error "Conversion from Symbolics-specific code not complete...")
  #+symbolics
  (let ((xdim (vu:raster-width 1b-array))
	(ydim (vu:raster-height 1b-array)))
    (loop for x from 0 below xdim do
      (loop for y from 0 below ydim do
	(let ((1b-el (vu:raster-aref 1b-array x y))
	      (art-q-el (vu:raster-aref art-q-array x y)))
	  (cond ((not (or (and (zerop 1b-el)
			       (not art-q-el))
			  (and (= 1 1b-el)
			       art-q-el)))
		 (vu::zero-raster *temp-1b-array*)
		 (let ((x1 (max (- x 2) 0))
		       (y1 (max (- y 2) 0))
		       (x2 (min (+ x 2) xdim))
		       (y2 (min (+ y 2) ydim)))
		   (si:%draw-line x1 y1 x2 y2 tv:alu-seta t *temp-1b-array*)
		   (si:%draw-line x1 y2 x2 y1 tv:alu-seta t *temp-1b-array*))
		 (bitblt tv:alu-ior xdim ydim 1b-array 0 0 *temp-1b-array* 0 0)
		 (show-1b-raster *temp-1b-array* 768 0 *mike-display-pane* t t)
		 (error "Some point in the right image did not get a match in the left image!! (~a,~a), 1b=~a art-q=~a"
			x y 1b-el art-q-el))))))))

#-symbolics
(defun LINEARLY-INTERPOLATE-BETWEEN-DOTS-IN-CALIBRATION-TEST-PATTERN
       (&optional (left-1b-array-of-dots *left-warp-1b*)
		  (right-1b-array-of-dots *right-warp-1b*)
		  (spiral-search-radius 20) (cone-search-radius 80)
		  (cone-half-angle 20))
  (declare (ignore left-1b-array-of-dots right-1b-array-of-dots
		   spiral-search-radius cone-half-angle cone-search-radius))
  (error "Conversion from Symbolics-specific code not complete..."))

#+symbolics
(defun LINEARLY-INTERPOLATE-BETWEEN-DOTS-IN-CALIBRATION-TEST-PATTERN
       (&optional (left-1b-array-of-dots *left-warp-1b*)
		  (right-1b-array-of-dots *right-warp-1b*)
		  (spiral-search-radius 20) (cone-search-radius 80)
		  (cone-half-angle 20))
  (format t "~%Clearing arrays...")
  (vu::zero-raster *temp-1b-array*)
  ;; Important: clear-raster sets art-q arrays to nil.  This info is used later on to tell
  ;; whether entries in *x-warp-art-q* and *y-warp-art-q* have been set.
  (clear-raster *x-warp-art-q*)
  (clear-raster *y-warp-art-q*)
  (format t "done.")
  ;;
  (show-1b-raster left-1b-array-of-dots 0 0 *mike-display-pane* t t)
  (show-1b-raster right-1b-array-of-dots (vu::raster-width left-1b-array-of-dots) 0 *mike-display-pane* t t)
  ;;
  (format t "~%Matching...")
  (match-up-with-nearest-dots left-1b-array-of-dots right-1b-array-of-dots *x-warp-art-q* *y-warp-art-q* spiral-search-radius)
  (format t "done.")
  ;;
  (format t "~%Checking...")
  (make-sure-every-point-in-right-image-got-a-match-in-left *right-warp-1b* *x-warp-art-q*)
  (make-sure-every-point-in-right-image-got-a-match-in-left *right-warp-1b* *y-warp-art-q*)
  (format t "done.")
  ;;
  (format t "~%Extending exterior lines...")
  (extend-all-exterior-lines *right-warp-1b* *x-warp-art-q* *y-warp-art-q* cone-search-radius cone-half-angle *temp-1b-array*)
  (show-1b-raster *temp-1b-array* 512)
  (format t "done.")
  ;;
  (format t "~%Connecting interior lines...")
  (connect-the-dots-with-linearly-interpolated-segments
    *right-warp-1b* *x-warp-art-q* cone-search-radius cone-half-angle *temp-1b-array*)
  (connect-the-dots-with-linearly-interpolated-segments
    *right-warp-1b* *y-warp-art-q* cone-search-radius cone-half-angle *temp-1b-array*)
  (show-1b-raster *temp-1b-array* 768)
  (format t "done.")
  ;;
  (format t "~%Setting corners...")
  (set-all-corners *x-warp-art-q*)
  (set-all-corners *y-warp-art-q*)
  (format t "done.")
  ;;
  (format t "~%Horizontally interpolating *x-warp-art-q*...")
  (set-side-edges *x-warp-art-q*)
  (horizontal-interpolate *x-warp-art-q*)
  (format t "done.")
  ;;
  (format t "~%Vertically interpolating *y-warp-art-q*...")
  (set-top-edges *y-warp-art-q*)
  (vertical-interpolate *y-warp-art-q*)
  (format t "done.")
  ;;
  (format t "~%Copying into the 8-bit arrays *x-warp-array* and *y-warp-array* ...")
  (loop for x from 0 below (vu::raster-width *x-warp-array*) do
    (loop for y from 0 below (vu::raster-height *y-warp-array*) do
      (setf (vu:raster-aref *x-warp-array* x y) (vu:raster-aref *x-warp-art-q* x y)
	    (vu:raster-aref *y-warp-array* x y) (vu:raster-aref *y-warp-art-q* x y))))
  (format t "done."))

(defun QUIT? (input-object)
  (and (symbolp input-object)
       (string-equal "Q" (symbol-name input-object))))

;;; ************************************************** CONVERT

#+symbolics
(*defunc THRESHOLD-PVAR-MANUALLY!! (pvar &optional (invert? t) (initial-threshold 60) (window *mike-display-pane*))
  (let ((threshold initial-threshold)
	cursor-y (input t))
    (zl-user:send window :clear-window)
    (setq cursor-y (+ 150 (dimension-size 1)))
    (loop until (quit? input) do
      (*show-bit (if invert? (not!! (>!! pvar (!! threshold))) (>!! pvar (!! threshold)))
		 0 0 0 window)
      (zl-user:send window :set-cursorpos 0 cursor-y)
      (zl-user:send window :clear-rest-of-line)
      (setq input (zl-user:prompt-and-read :expression "Type a threshold increment (threshold currently ~a) or type Q to quit: "
					   threshold))
      (if (numberp input) (setq threshold (+ threshold input))))
    (if invert? (not!! (>!! pvar (!! threshold))) (>!! pvar (!! threshold)))))

;;; ************************************************** CONVERT

#+symbolics
(*defunc CUT-OFF-EDGES-MANUALLY!! (boolean!! &optional (window *mike-display-pane*))
  (let ((input t) (bcut 0) (rcut 0) (lcut 0) (tcut 0) cursor-y)
    (*let ((x-addr!! (x!!))
	   (y-addr!! (y!!))
	   (temp!! boolean!!))
      (zl-user:send window :clear-window)
      (setq cursor-y (+ 150 (dimension-size 1)))
      (loop until (quit? input) do
	(*show-bit (and!! temp!! (<!! y-addr!! (!! (- (dimension-size 1) bcut))))
		   0 0 0 window)
	(zl-user:send window :set-cursorpos 0 cursor-y)
	(zl-user:send window :clear-rest-of-line)
	(setq input (zl-user:prompt-and-read :expression "Cut the right edge by (currently ~a) or type Q to go on: "
					     bcut))
	(if (numberp input) (setq bcut (+ bcut input))))
      (*set temp!! (and!! temp!! (<!! y-addr!! (!! (- (dimension-size 1) bcut)))))
      (setq input t)
      (loop until (quit? input) do
	(*show-bit (and!! temp!! (>!! x-addr!! (!! lcut)))
		   0 0 0 window)
	(zl-user:send window :set-cursorpos 0 cursor-y)
	(zl-user:send window :clear-rest-of-line)
	(setq input (zl-user:prompt-and-read :expression "Cut the top edge by (currently ~a) or type Q to go on: "
					     lcut))
	(if (numberp input) (setq lcut (+ lcut input))))
      (*set temp!! (and!! temp!! (>!! x-addr!! (!! lcut))))
      (setq input t)
      (loop until (quit? input) do
	(*show-bit (and!! temp!! (<!! x-addr!! (!! (- (dimension-size 0) rcut))))
		   0 0 0 window)
	(zl-user:send window :set-cursorpos 0 cursor-y)
	(zl-user:send window :clear-rest-of-line)
	(setq input (zl-user:prompt-and-read :expression "Cut the bottom edge by (currently ~a) or type Q to go on: "
					     rcut))
	(if (numberp input) (setq rcut (+ rcut input))))
      (*set temp!! (and!! temp!! (<!! x-addr!! (!! (- (dimension-size 0) rcut)))))
      (setq input t)
      (loop until (quit? input) do
	(*show-bit (and!! temp!! (>!! y-addr!! (!! tcut)))
		   0 0 0 window)
	(zl-user:send window :set-cursorpos 0 cursor-y)
	(zl-user:send window :clear-rest-of-line)
	(setq input (zl-user:prompt-and-read :expression "Cut the left edge by (currently ~a) or type Q to quit: "
					     tcut))
	(if (numberp input) (setq tcut (+ tcut input))))
      (and!! temp!! (>!! y-addr!! (!! tcut))))))

(*defunc MAKE-TEST-PATTERN-ARRAY
	(8b-array dest-1b-array &optional (window *mike-display-pane*))
  #-symbolics (declare (ignore 8b-array dest-1b-array window))
  #-symbolics (error "Conversion from Symbolics-specific code not complete...")
  #+symbolics
  (*let* ((pvar!! (cmv:write-raster-to-cm-grid!! 8b-array))
	  (thresholded!! (threshold-pvar-manually!! pvar!! t 60 window))
	  (with-trimmed-edges!! (cut-off-edges-manually!! thresholded!! window)))
    (declare (type boolean-pvar thresholded!! with-trimmed-edges!!))
    (cmv:*read-raster-from-cm-grid *temp-1b-array* with-trimmed-edges!!)
    (replace-regions-with-their-centroids *temp-1b-array* 5 5 dest-1b-array)
    (show-1b-raster dest-1b-array)
    dest-1b-array))

(*defunc FILL-CALIBRATION-ARRAYS (&key (window *mike-display-pane*) already-grabbed)
  (declare (special *bilinear-subsampling?*))
  (if (not already-grabbed)
      (grab-stereo-pair-twiddle))
  ;; Experimental feature: downsampling using bilinear interpolation to map the 454x454 central section
  ;; of a 576x454 frame into a 256x256 array.
  (if *bilinear-subsampling?*
      (progn (cmv:*read-raster-from-cm-grid *left-256* (cmv:subsample-frame!! *left-array*))
	     (cmv:*read-raster-from-cm-grid *right-256* (cmv:subsample-frame!! *right-array*)))
      (move-frames 256))
  (make-test-pattern-array *left-256* *left-warp-1b* window)
  (make-test-pattern-array *right-256* *right-warp-1b* window)
  (linearly-interpolate-between-dots-in-calibration-test-pattern *left-warp-1b* *right-warp-1b* 30 80 20))
						; function.  (peo, 19 feb 88)
;;; commented out until problem in *grab-and-register is fixed
;(defun TEST-CALIBRATION-MANUALLY ()
;  (*let ((left!! (!! 0))
;	 (right!! (!! 0)))
;    (declare (type (field-pvar 8) left!! right!!))
;    (*grab-and-register left!! right!! t)
;    (cmv::*read-raster-from-cm-grid *left-256* left!!)
;    (cmv::*read-raster-from-cm-grid *right-256* right!!)
;    (format t "~%Hit space to quit...")
;;    (grey:clear)
;    (send spe:*color-window* :clear-window)
;    (loop until (tv:key-state #\space) do  ;;; hit space bar
;      (spe:show-color spe:*color-window* *left-256*)
;      (sleep .5)
;      (spe:show-color spe:*color-window* *right-256*)
;      (sleep .5))))

(defvar WARPING-ARRAYS nil)
(defun SAVE-WARPING-ARRAYS (&optional (filename "warping-arrays"))
  (setq warping-arrays (list *x-warp-array* *y-warp-array*))
  (save-data warping-arrays filename *calibration-data-directory*))

(defun LOAD-WARPING-ARRAYS (&optional (filename "warping-arrays"))
  (declare (special foo))
  (load-data filename 'foo *calibration-data-directory*)
  (setq *x-warp-array* (first foo)
	*y-warp-array* (second foo)))


;;; ***stale code***
;;;   Not tested for a long time.  Do not use if you have any self respect.
;;;
;;; This piece of code is lying around just in case the router is not working.  It is functionally exactly equivalent to
;;; WARP-WITH-ROUTER!!.
;;;
;;; Here the dx and dy are divided by 2 because all the calibration software calculated the disparities     
;;; to be twice as big as they really are, to be used by half-pixel-precision functions
;;;
;;; PLEASE NOTE that this function is slow, it does an awful lot of redundant newsing-around.
;(*defunc WARP-WITH-GRID!!
;	(image!! &optional x-warp-func-pvar y-warp-func-pvar)
;  (format t "~%Warping with grid...")
;  (*let ((1+image!! (1+!! image!!)) ;;; <-- Make sure there are no zeros in the image to be warped.
;	 (new-x!! (max!! (!! 0) (+!! (x!!) x-warp-func-pvar)))
;	 (new-y!! (max!! (!! 0) (+!! (y!!) y-warp-func-pvar)))
;	 (collision-count!! (!! 0))
;	 (selfx!! (x!!))
;	 (selfy!! (y!!))
;	 (warped-image!! (!! 0))
;	 (receive-new-x!! (!! 0))
;	 (receive-new-y!! (!! 0))
;	 (receive-image!! (!! 0)))
;    (declare (type field-pvar new-x!! new-y!! collision-count!! selfx!! selfy!! receive-new-x!! receive-new-y!!)
;	     (type signed-pvar warped-image!! receive-image!!))
;    ;;  This let is tricky.  Notice that because we're GETTING FROM the image-to-be-warped,
;    ;;  the sign of the amount-to-warp changes.
;    (let ((minx (1- (- (*max x-warp-func-pvar))))
;	  (maxx (1+ (- (*min x-warp-func-pvar))))
;	  (miny (1- (- (*max y-warp-func-pvar))))
;	  (maxy (1+ (- (*min y-warp-func-pvar)))))
;      (loop for y from miny to maxy do
;	(loop for x from minx to maxx do
;	  (*set receive-new-x!! (get-grid!! new-x!! x y))
;	  (*set receive-new-y!! (get-grid!! new-y!! x y))
;	  (*set receive-image!! (get-grid!! 1+image!! x y))
;	  (pref selfx!! 0)
;	  (*cond ((and!! (=!! receive-new-x!! selfx!!) (=!! receive-new-y!! selfy!!))
;		  (*incf warped-image!! receive-image!!)
;		  (*incf collision-count!!))))))
;    (*if (nzerop!! warped-image!!)
;	 (*set warped-image!! (floor!! warped-image!! collision-count!!)))
;    (format t "done.")
;    (max!! (!! 0) (1-!! (fill-in-holes!! warped-image!! 1)))))	; The MAX!! here is done because some of the borders receive negative garbage.