;;; This software is being provided to you, the LICENSEE, by the
;;; Massachusetts Institute of Technology (M.I.T.), and the NEC Research
;;; Institute (NECI) under the following license.  By obtaining, using
;;; and/or copying this software, you agree that you have read,
;;; understood, and will comply with these terms and conditions:
;;; 
;;; Permission to use, copy, modify and distribute, including the right to
;;; grant others rights to distribute at any tier, this software and its
;;; documentation for research purposes only and without fee or royalty is
;;; hereby granted, provided that you agree to comply with the following
;;; copyright notice and statements, including the disclaimer, and that
;;; the same appear on ALL copies of the software and documentation,
;;; including modifications that you make for internal use or for
;;; distribution:
;;; 
;;; Copyright 1995 by the Massachusetts Institute of Technology.  All
;;; rights reserved.  Copyright 1995 by NEC Research Institute Inc. All
;;; rights reserved
;;; 
;;; THIS SOFTWARE IS PROVIDED "AS IS", AND M.I.T. and NECI MAKE NO
;;; REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By way of example,
;;; but not limitation, M.I.T. and NECI MAKE NO REPRESENTATIONS OR
;;; WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
;;; THAT THE USE OF THE LICENSED SOFTWARE OR DOCUMENTATION WILL NOT
;;; INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER
;;; RIGHTS.
;;; 
;;; The name of the Massachusetts Institute of Technology or M.I.T., or
;;; NEC Research Institute, or NECI may NOT be used in advertising or
;;; publicity pertaining to distribution of the software.  Title to
;;; copyright in this software and any associated documentation shall at
;;; all times remain with M.I.T. and NECI, and USER agrees to preserve
;;; same.
;;; 

;;; -*- Mode: LISP; Base: 10; Package: 2DG; Syntax: Common-Lisp -*-
(in-package '2DG)

(in-package 'g)

(export '(*sqrt2*))

(in-package '2dg)

(export '(theta-s-length theta-s-thetas theta-s-ss make-theta-s theta-s-p
          set-theta-s-theta set-theta-s-s theta-s-theta-list theta-s-s theta-s-theta print-theta-s
	  theta-s-from-string string-from-theta-s))

(defstruct (THETA-S
	    :named
	    (:constructor make-theta-s (length &aux (thetas (make-array (list length) :element-type 'double-float))
					            (ss (make-array (list length) :element-type 'double-float))))
	    :predicate
	    (:print-function (lambda (theta-s stream &rest ignore)
			      (format stream "#<THETA-S-GRAPH>"))))
  "A theta-s space graph"
  length
  thetas
  ss)

;;; Accessors for theta-s-graph

(defun SET-THETA-S-THETA (ts i val)
  (setf (aref (theta-s-thetas ts) i) val))

(defun SET-THETA-S-S (ts i val)
  (setf (aref (theta-s-ss ts) i) val))

(defun THETA-S-THETA-LIST (ts)
  (let ((thetas (theta-s-thetas ts)))
    (loop for i from 0 to (1- (theta-s-length ts))
	  collecting
	  (aref thetas i))))

(defun THETA-S-S (ts i)
  (aref (theta-s-ss ts) i))

(defun THETA-S-THETA (ts i)
  (aref (theta-s-thetas ts) i))

(defun PRINT-THETA-S (ts &optional (stream t) max-length (start 0) (new-lines t))
  (when new-lines (format stream "~%"))
  (format stream ":length = ~A" (theta-s-length ts))
  (loop for i from start to (1- (if max-length (min (+ start max-length) (theta-s-length ts))
				    (theta-s-length ts)))
	doing
    (when (and new-lines (= 0 (mod i 5))) (format stream "~%"))
    (if (integerp (theta-s-s ts i)) (format stream "(~3D " (theta-s-s ts i))
	(format stream "(~3F " (theta-s-s ts i)))
    (format stream "~5F) " (theta-s-theta ts i))))
  
;;; This routine describes the direction between two neighboring points with a number, as follows:
;;;
;;;   3  2  1
;;;   4  P  0
;;;   5  6  7
;;; where P shows the point.
(defun THETA-S-FROM-INTEGER-STRING (string)
  "Construct a theta-s space graph from a list of connected points."
  (flet ((neighbor-direction (point neighbor)
	   (cond ((= (point-x point) (point-x neighbor))
		  (if (< (point-y point) (point-y neighbor)) 2 6))
		 ((< (point-x point) (point-x neighbor))
		  (cond ((= (point-y point) (point-y neighbor)) 0)
			((< (point-y point) (point-y neighbor)) 1)
			(t 7)))
		 (t (cond ((= (point-y point) (point-y neighbor)) 4)
			  ((< (point-y point) (point-y neighbor)) 3)
			  (t 5))))))
    (let* ((string-circular? (points-neighbors? (car string) (car (last string))))
	   (excess-distance 0)			; The length of diagonal connections in excess of 1 is accumulated so
						; that some points are counted twice.
	   (thetas nil))
      (loop for point in string
	    for next-point in (if string-circular? (append (cdr string) (list (car string))) (cdr string))
	    for i from 0
	    doing
	(let ((dir (neighbor-direction point next-point)))
	  (when (oddp dir)
	    (setf excess-distance (+ excess-distance (- *sqrt2* 1))))
	  (push (* dir *PI/4*) thetas)
	  (when (<= 1 excess-distance)
	    (push (* dir *PI/4*) thetas)
	    (decf excess-distance))))
      (let* ((theta-s-length (length thetas))
	     (theta-s-graph (make-theta-s theta-s-length)))
	(loop for theta in thetas
	      for i from (1- theta-s-length) downto 0	; The thetas were put in in reverse order
	      doing
	  (set-theta-s-theta theta-s-graph i theta)
	  (set-theta-s-s theta-s-graph i i))
	theta-s-graph))))

(defun THETA-S-FROM-STRING (string &optional string-circular?)
  "Construct a theta-s space graph from a list of connected points."
  (flet ((neighbor-direction (point neighbor)
	   (let ((dx (- (point-x neighbor) (point-x point)))
		 (dy (- (point-y neighbor) (point-y point))))
	     (values (atan-0-2pi dy dx) (sqrt (+ (sq dx) (sq dy))))))
	 (direction-difference (d1 d2)  ; A direction is from 0 to 2pi.  This gives d1 - d2, ranging from -pi to pi.
	   (if (< (abs (- d1 d2)) pi)
	       (- d1 d2)
	       (if (< d1 d2)
		   (- d1 (- d2 *2pi*))
		   (- (- d1 *2pi*) d2)))))
    (let* ((excess-distance 0)			; The length of diagonal connections in excess of 1 is accumulated so
						; that some points are counted twice.
	   (thetas nil)
	   (last-theta 0)
	   (last-direction 0))  ; The first theta will be between -PI and PI
      (loop for point in string
	    for next-point in (if string-circular? (append (cdr string) (list (car string))) (cdr string))
	    for i from 0
	    doing
	(multiple-value-bind (direction distance)
	    (neighbor-direction point next-point)
	  (let ((theta (+ last-theta (direction-difference direction last-direction))))
	    (setf last-theta theta)
	    (setf last-direction direction)
	    (setf excess-distance (+ excess-distance distance))
	    (loop until (< excess-distance 1)
		  doing
	      (push theta thetas)
	      (decf excess-distance)))))
      ;; We essentially sample the string at intervals of 1 pixel, finding the theta at those points.
      ;; Excess-distance keeps track of how far "point" is from the last integer distance.
      (let* ((theta-s-length (length thetas))
	     (theta-s-graph (make-theta-s theta-s-length)))
	(loop for theta in thetas
	      for i from (1- theta-s-length) downto 0	; The thetas were put in in reverse order
	      doing
	  (set-theta-s-theta theta-s-graph i theta)
	  (set-theta-s-s theta-s-graph i i))
	theta-s-graph))))

(defun STRING-FROM-THETA-S (ts &optional (start-point (make-point 0 0)))
  "Returns a string of connected points that begin at (0 0).  If two
   adjacent entries in the theta-s space graph map to non-neighbor points
   in the plane, we interpolate between them with a line."
  (let* ((last-point start-point)
	 (round-last-point last-point)
	 (last-theta (theta-s-theta ts 0))
	 (last-s (theta-s-s ts 0)))
    (cons start-point
	  (loop for i from 1 to (1- (theta-s-length ts))
		appending
		  (let* ((ds (- (theta-s-s ts i) last-s))
			 (new-theta (theta-s-theta ts i))
			 (new-point (add last-point (make-point (* ds (cos last-theta)) (* ds (sin last-theta)))))
			 (round-new-point (round-point new-point))
			 (string-points
			   (if (point-equal round-last-point round-new-point) nil
			       (if (points-neighbors? round-last-point round-new-point) (list round-new-point)
				   (cdr (points-between round-last-point round-new-point))))))
		    ;; Cdr to prevent repitition of last point.
		    (setf last-s (+ last-s ds))
		    (setf last-point new-point)
		    (setf round-last-point round-new-point)
		    (setf last-theta new-theta)
		    string-points)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; These go into other files in the geometry package.

(export '(round-point points-neighbors? mutate-point-x mutate-point-y geometric-object? grow-strings))

(defun ROUND-POINT (pt)
  (make-point (round (point-x pt)) (round (point-y pt))))
						
(defun POINTS-NEIGHBORS? (point1 point2)
  (and (<= (abs (- (point-x point1) (point-x point2))) 1)
       (<= (abs (- (point-y point1) (point-y point2))) 1)))

(defun MUTATE-POINT-X (pt new-x)
  (setf (point-x pt) new-x))

(defun MUTATE-POINT-Y (pt new-y)
  (setf (point-y pt) new-y))

(defun GEOMETRIC-OBJECT? (object)
  (case (geometric-type object)
    ((point line line-segment rectangle polygon) t)
    (t nil)))

(defun ATAN-0-2PI (y x)
  (let ((ans (atan y x)))
    (if (< ans 0) (+ *2pi* ans) ans)))

(defun GROW-STRINGS (strings raster &optional (max-dist 2.9))
  (let* ((sq-max-dist (square max-dist))
	 (last-points (loop for s in strings collecting (car (last s))))
	 (potential-matches
	   (loop for s1 in strings
		 for p1 = (car s1)
		 for p1l in last-points
		 for rest-strings on (cdr strings)
		 for rest-last on (cdr last-points)
		 appending
		 (when (not (points-neighbors? p1 p1l))
		   (loop for s2 in rest-strings
			 for p2 = (car s2)
			 for p2l in rest-last
			 appending
			 (when (not (points-neighbors? p2 p2l))
			   (append (when (and (< (distance-squared p1 p2) sq-max-dist)
					      (point-endpoint? p1 raster) (point-endpoint? p2 raster))
				     (list (list s1 s2 nil nil)))
				   (when (and (< (distance-squared p1l p2) sq-max-dist)
					      (point-endpoint? p1l raster) (point-endpoint? p2 raster))
					      (list (list s1 s2 t nil)))
				   (when (and (< (distance-squared p1 p2l) sq-max-dist)
					      (point-endpoint? p1 raster) (point-endpoint? p2l raster))
				     (list (list s1 s2 nil t)))
				   (when (and (< (distance-squared p1l p2l) sq-max-dist)
					      (point-endpoint? p1l raster) (point-endpoint? p2l raster))
				     (list (list s1 s2 t t)))))))))
	 (matches (loop for match in potential-matches
			appending
			  (when (and (= 1 (loop for match2 in potential-matches
						summing
						  (if (and (eq (car match) (car match2)) (eq (nth 2 match) (nth 2 match2)))
						      1 0)))
				     (= 1 (loop for match2 in potential-matches
						summing
						  (if (and (eq (cadr match) (cadr match2)) (eq (nth 3 match) (nth 3 match2)))
						      1 0))))
			    (list match))))
	 (matches-left matches))
    (let* ((strings-to-remove nil)
	   (new-strings
	     (loop until (null matches-left)
		   appending
	       (let* ((match (pop matches-left))
		      (s1 (if (not (nth 2 match)) (reverse (car match)) (car match))) 
		      ;; The first string, with its last point the one that's near the second string.
		      (rev-s1 (if (nth 2 match) (reverse (car match)) (car match)))	; reverse of s1
		      (s2 (if (nth 3 match) (reverse (cadr match)) (cadr match)))	; first point is the one near s1.
		      (connecting-points (string-connecting-points (car rev-s1) (car s2) raster)))
		 (when connecting-points
		   (let ((new-string (append s1 connecting-points s2)))
		     (push (car match) strings-to-remove)
		     (push (cadr match) strings-to-remove)
		     (loop for match2 in matches-left
			   doing
		       (when (eq (car match) (car match2))
			 (progn (rplaca match2 new-string)
				(rplaca (nth 2 match2) nil)))	; The nearby point in this match must be the end of the new string
		       ;; We don't have to compare (car match) and (cadr match2), because of the order of the matches 
		       (if (eq (cadr match) (car match2))
			   (progn (rplaca match2 new-string)
				  (rplaca (cdddr match2) t))
			   (if (eq (cadr match) (cadr match2))
			       (progn (rplaca (cdr match2) new-string)
				      (rplaca (cdddr match2) t)))))
		     (list new-string)))))))
      (append new-strings (loop for string in strings appending (when (not (member string strings-to-remove)) (list string)))))))

(defun STRING-CONNECTING-POINTS (p1 p2 raster)
  "The first points of the two strings are nearby.  We return a list of points that bridges the gap.  This must
   done with points that aren't adjacent to other 
   For now, they must be separated by only a single pixel."
  (loop for candidate in (intersection (point-neighbors p1 raster) (point-neighbors p2 raster) :test #'point-equal)
	doing
    (when (loop for candidate-neighbor in (remove p1 (remove p2 (point-neighbors candidate raster) :test #'point-equal)
						  :test #'point-equal)
		always
		(= 0 (raster-aref raster (point-x candidate-neighbor) (point-y candidate-neighbor))))
      (return (list candidate)))))

(defun POINT-ENDPOINT? (point raster)
  (= 1 (loop for pt in (point-neighbors point raster)
	     summing
	     (raster-aref raster (point-x pt) (point-y pt)))))

(defun POINT-NEIGHBORS (point raster-image)
  (multiple-value-bind (max-x max-y)
      (decode-raster-array raster-image)
    (loop for dx in (list -1 0 1)
	  with x = (point-x point)
	  with y = (point-y point)
	  appending
	    (loop for dy in (if (= dx 0) (list -1 1) (list -1 0 1))
		  for nx = (+ x dx)
		  for ny = (+ y dy)
		  appending
		    (when (and (<= 0 nx) (<= 0 ny) (< nx max-x) (< ny max-y))
		      (list (make-point nx ny)))))))

(defun STRING-CIRCULAR? (s)
  (and (null (cdr s)) (points-neighbors? (car s) (car (last s)))))

