; (c) 2008 Dmitri Hrapof <hrapof@common-lisp.ru>
; This file is distributed under the terms of the LLGPL.

(defpackage oci
   (:use common-lisp coeden)
   (:export
    #:MKLIST
    #:INTERLEAVE
    #:GROUP
    #:GRAPHICALIZE
    #:GRAPH->LINES
    #:IMAGE->LINES
    #:IMAGE->BLOCKS
    #:BLOCKS->GRAPH
    #:SHRINK-LST
    #:INSTANTIATE-RECTANGLE
    #:MERGE-BLOCKS
    #:HAF-TRANSFORM
    #:MOUNT-BLOCKS
    #:BRESENHAM-LINE
    #:DITER
    #:MAKE-GRAPH
    #:DIJKSTRA-SHORTEST-PATH
    #:WIDEN-LINES
    #:RENDER-SVG-PATH
    #:LINES->SVG-PATH
    #:TREE->LINES
    #:REDRAW
    #:WORSEN-VISIBILITY
    #:ASSIGNMENT
    #:MAKE-COST-MATRIX
    #:GRAPH-COST
    #:GRAPH43
    #:ADD23
    #:GRAFS->TREE
    #:GRAFS->FILE
    #:GATHER-FRUIT
    #:ALIST->HASH
    #:OBSERVE-WORD
    #:BOXIFY
    #:CARVE-INTO-BOXES
    #:COPY-RECT
    #:write-netpbm
    #:read-netpbm
    #:fourier-search
    #:negate-grayscale-image
    #:logarithm-filter
    #:shrink-image
    #:fourier-mellin-search
    #:mark-image
    #:join-centres
    #:filter-centres
    #:restore-rect
    #:image-mean
    #:hough+sobel-circle-transform
    #:hough+sobel-line-transform
    #:grayscale-image
    #:average
    #:euclidean-distance
    #:treshold-fourier!
    #:pad-image
    #:img<-path
    #:remember-img
    #:prune
    #:reap-meat
    #:dilate-points
    #:cut-into-words
    #:grayscale->bitmap!
    #:new-point
    #:point-i
    #:point-j
    #:power-of-2-ceiling
    #:chamfer-transform!
    #:image-not
    #:list-boxes
    #:hr-distance
    *this-file*
    *lng*
    *omg*
    *frg*
    *oak*
    *acorn*
    *scripts*
))
(in-package oci)

(declaim (inline cref))
(defun cref (image i j &optional (v 0)) 
  (if (array-in-bounds-p image i j) (aref image i j) v))

(defun graph-node (g i) (car (svref g i)))
(defun (setf graph-node) (val g i) (setf (car (svref g i)) val))

(defun graph-edge (g i j) (aref (cdr (svref g i)) j))
(defun (setf graph-edge) (val g i j) (setf (aref (cdr (svref g i)) j) val))

;(defun make-graph (ld lst)
;  (let* ((len (length lst)) (g (make-array len)))
;    (loop for l in lst for i from 0 do
;	 (setf (svref g i) (cons l (make-array len))))
;    (dotimes (i len g)
;      (dotimes (j len)
;	(setf (graph-edge g i j) 
;	      (funcall ld (graph-node g i) (graph-node g j)))))))

(defun make-graph (ld lst)
  (let ((len (length lst)))
    (let ((va (apply #'vector lst))
	  (ga (make-array (list len len))))
      (dotimes (i len)
	(dotimes (j len)
	  (setf (aref ga i j) (funcall ld (svref va i) (svref va j)))))
      (values ga va))))

;For sparse graphs, that is, graphs with much less than n^2 edges,
;Dijkstra's algorithm can be implemented more efficiently by storing the
;graph in the form of adjacency lists and using a binary heap or Fibonacci
;heap as a priority queue to implement the extract-min function. With a
;binary heap, the algorithm requires O((m+n)log n) time, and the Fibonacci
;heap improves this to O(m + n log n).
  
(defun dijkstra-shortest-path 
    (graph &key (elt 0) (end (1- (array-dimension graph 0))))
  (labels
   ((extract-min 
     (s d)
     (let ((z (sort s #'(lambda (v1 v2) (< (aref d v1) (aref d v2))))))
       (values (car z) (cdr z))))
    (dijkstra 
     (tg e s q d p)
     (if (not q) (reverse s)
       (multiple-value-bind 
	(u k) (extract-min q d)
	(if (and e (= u e)) (reverse (push (cons u (aref d u)) s))
	  (progn
	    (dotimes (v (array-dimension tg 0))
	      (let ((w (aref #|graph-edge|# tg u v)))
		(if (and w (> (aref d v) (+ w (aref d u))))
		    (setf (aref d v) (+ w (aref d u))
			  (aref p v) u))))
	    (dijkstra tg e (push (cons u (aref d u)) s) k d p)))))))
   (let ((d (make-array (array-dimension graph 0) 
			:initial-element most-positive-fixnum))
	 (p (make-array (array-dimension graph 0)
			:initial-element nil)))
     (setf (aref d elt) 0)
     (values
      p
      (dijkstra graph end '() 
		(loop for i below (array-dimension graph 0) collect i) d p)))))

(declaim (inline new-point))
(defun new-point (i j) (cons i j))
(declaim (inline point-i))
(defun point-i (p) (car p))
(declaim (inline point-j))
(defun point-j (p) (cdr p))

(defun square (num) (* num num))
(defun average (&rest args) 
  (let ((len (length args)))
    (if (zerop len) 0 (/ (apply #'+ args) len))))
(defun euclidean-distance (p1 p2) 
  ;(sqrt (apply #'+ (mapcar #'square (mapcar #'- p1 p2)))))
  (let ((dy (- (point-i p1) (point-i p2)))
	(dx (- (point-j p1) (point-j p2))))
    (sqrt (+ (* dy dy) (* dx dx)))))

(defun diter (v &optional s (u (1- (array-dimension v 0))))
  (if (not u) s (diter v (push u s) (aref v u))))

(defun read-netpbm (fname)
  (labels ((bread (f) (case (read-char f) (#\0 0) (#\1 1) (t (bread f)))))
    (with-open-file (f fname)
      (let ((rt (copy-readtable)))
	(multiple-value-bind (fun mid) (get-macro-character #\;)
	  (set-macro-character #\# fun mid rt)
	  (let* ((*readtable* rt)
		 (magic (read f)) (width (read f)) (height (read f))
		 (max (if (string= "P1" (symbol-name magic)) nil (read f)))
		 (rid (if (string= "P1" (symbol-name magic)) #'bread #'read))
		 (image (make-array `(,height ,width) :element-type '(unsigned-byte 8))))
	    (dotimes (i height (values image max))
	      (dotimes (j width) (setf (aref image i j) (funcall rid f))))))))))

(defun write-netpbm (f image)
  (with-open-file                                                               
   (fd f :direction :output :if-exists :supersede :if-does-not-exist :create)   
   (let ((magic (if (search "pbm" f) "P1" "P2"))
	 (height (array-dimension image 0))                                     
         (width (array-dimension image 1))                                      
         (count 1))                                                             
    (format fd "~a~%~a ~a~%" magic width height)
    (if (string= magic "P2") (format fd "255"))
    (dotimes (i height)
      (dotimes (j width)
	(if (string= magic "P2") (format fd "~%~a" (aref image i j))
	    (progn
	      (if (> count 70) (progn (setf count 1) (format fd "~%")))
	      (format fd "~a" (aref image i j))
	      (incf count)))))
    (format fd "~%"))))

(defun point< (a b)
  (or (< (point-i a) (point-i b)) 
      (and (= (point-i a) (point-i b)) (< (point-j a) (point-j b)))))

(defun hor-point< (a b)
  (or (< (point-j a) (point-j b)) 
      (and (= (point-j a) (point-j b)) (< (point-i a) (point-i b)))))

(defun /-point< (a b)
  (or (< (point-i a) (point-i b)) 
      (and (= (point-i a) (point-i b)) (> (point-j a) (point-j b)))))

(defparameter saxons
  (let ((s #(:hor
	     :dia011
	     :dia022
	     :dia033
	     :dia045
	     :dia056
	     :dia067
	     :dia078
	     :ver
	     :dia101
	     :dia112
	     :dia123
	     :dia135
	     :dia146
	     :dia157
	     :dia168)))
    (setf (get :curve-q1 :angle) 45)
    (setf (get :curve-q2 :angle) 135)
    (setf (get :curve-q3 :angle) 45)
    (setf (get :curve-q4 :angle) 135)
    (setf (get :dot :angle) 0)
    (dotimes (a (length s) s)
      (setf (get (svref s a) :angle) (* a 11.25)))))

(defun angle->tmp (a) (svref saxons (floor a 11.25)))

(defun angle~ (diff sym1 sym2)
  (let ((a1 (if (symbolp sym1) (get sym1 :angle) sym1))
	(a2 (if (symbolp sym2) (get sym2 :angle) sym2)))
    (if (> diff (abs (- a1 a2))) t
	(cond
	  ((= 0 a1) (> diff (abs (- 180 a2))))
	  ((= 0 a2) (> diff (abs (- a1 180))))))))

(defun %rotate-point (c s p)
  (new-point (round (+ (* (point-i p) c) (* (point-j p) s)))  
	     (round (- (* (point-j p) c) (* (point-i p) s)))))

(defun rotate-point (a p)
  (let* ((r (/ (* a pi) 180)) (c (cos r)) (s (sin r))) (%rotate-point c s p)))

(defun image->blacks (image)
  (let (blacks 
	(bl 0) 
	(height (array-dimension image 0))
	(width (array-dimension image 1)))
    (dotimes (i height)
      (dotimes (j width)
	(if (not (zerop (aref image i j)))
	    (progn (incf bl) (push (new-point i j) blacks)))))
    (values (nreverse blacks) bl height width)))

(defun move-and-clip-rectangle (lim y x height width rect)
  (let (new (mil 0))
    (dolist (p rect)
      (let ((nx (+ x (point-j p))) (ny (+ y (point-i p))))
	(if (and (>= ny 0) (< ny height) (>= nx 0) (< nx width))
	    (push (new-point ny nx) new)
	    (if (= (incf mil) lim) (return-from move-and-clip-rectangle nil)))))
    (nreverse new)))

(defun instantiate-rectangle (rect)
  (let ((i-max (point-i (car rect))) (i-min (point-i (car rect))) 
	(j-max (point-j (car rect))) (j-min (point-j (car rect))))
    (dolist (p (cdr rect))
      (setf i-max (max i-max (point-i p)))
      (setf j-max (max j-max (point-j p)))
      (setf i-min (min i-min (point-i p)))
      (setf j-min (min j-min (point-j p))))
    (let ((a (make-array (list (1+ (- i-max i-min)) (1+ (- j-max j-min)))
			 :initial-element 0)))
      (dolist (p rect)
	(setf (aref a (- (point-i p) i-min) (- (point-j p) j-min)) 1))
      (values a i-min j-min))))

(eval-when (:execute :load-toplevel :compile-toplevel) 
(defun condlet-binds (vars cl)                                                  
  "Paul Graham. On Lisp, p.146"                                                 
  (mapcar #'(lambda (bindform)                                                  
              (if (consp bindform)                                              
                  (cons (cdr (assoc (car bindform) vars))                       
                        (cdr bindform))))                                       
          (cdr cl))))
                                                            
(eval-when (:execute :load-toplevel :compile-toplevel)
(defun condlet-clause (vars cl bodfn)                                           
  "Paul Graham. On Lisp, p.146"                                                 
  `(,(car cl) (let ,(mapcar #'cdr vars)                                         
                (let ,(condlet-binds vars cl)                                   
                  (,bodfn ,@(mapcar #'cdr vars)))))))

(defmacro condlet (clauses &body body)                                          
  "Paul Graham. On Lisp, p.146"                                                 
  (let ((bodfn (gensym))                                                        
        (vars (mapcar #'(lambda (v) (cons v (gensym)))                          
                      (remove-duplicates                                        
                       (mapcar #'car                                            
                               (apply #'append (mapcar #'cdr clauses)))))))     
    `(labels ((,bodfn ,(mapcar #'car vars)                                      
                      ,@body))                                                  
             (cond ,@(mapcar #'(lambda (cl)                                     
                                 (condlet-clause vars cl bodfn))                
                             clauses)))))  

;REM Bresenham algorithm for a line in an arbitrary octant in pseudo Basic
;dx = xend-xstart
;dy = yend-ystart
;
;REM Initialisations
;adx = ABS(dx): ady = ABS(dy) ' Absolute values of distances
;sdx = SGN(dx): sdy = SGN(dy) ' Signum of distances
;
;IF adx > ady THEN
;  ' x is fast direction
;  pdx = sdx: pdy = 0   ' pd. is parallel step
;  ddx = sdx: ddy = sdy ' dd. is diagonal step
;  ef  = ady: es  = adx ' error steps fast, slow
;             ELSE
;  ' y is fast direction
;  pdx = 0  : pdy = sdy ' pd. is parallel step
;  ddx = sdx: ddy = sdy ' dd. is diagonal step
;  ef  = adx: es  = ady ' error steps fast, slow
;  END IF
;
;x = xstart
;y = ystart
;SETPIXEL x,y
;error = es/2
;
;REM Pixel loop: always a step in fast direction, every now and then also one in slow direction
;FOR i=1 TO es          ' es also is the count of pixels zo be drawn
;   REM update error term
;   error = error - ef
;   IF error < 0 THEN
;      error = error + es ' make error term positive (>=0) again
;      REM step in both slow and fast direction
;      x = x + ddx: y = y + ddy ' Diagonal step
;                ELSE
;      REM step in fast direction
;      x = x + pdx: y = y + pdy ' Parallel step
;      END IF
;   SETPIXEL x,y
;   NEXT

(defun bresenham-line (xstart xend ystart yend)
  (let ((dx (- xend xstart)) (dy (- yend ystart)) line)
    (let ((adx (abs dx)) (ady (abs dy)) (sdx (signum dx)) (sdy (signum dy)))
      (condlet
       (((> adx ady) (pdx sdx) (pdy 0) (ddx sdx) (ddy sdy) (ef ady) (es adx))
	(t (pdx 0) (pdy sdy) (ddx sdx) (ddy sdy) (ef adx) (es ady)))
       (let ((x xstart) (y ystart) (error (/ es 2)))
	 (push (new-point y x) line)
	 (dotimes (i es)
	   (decf error ef)
	   (if (minusp error)
	       (progn (incf error es) (incf x ddx) (incf y ddy))
	       (progn (incf x pdx) (incf y pdy)))
	   (push (new-point y x) line)))))
    (nreverse line)))

(defun build-rectangle (x y w l a)
  (labels
      ((fill-rect (lb rect vertices contour)
	 (if (null contour)
	     (values
	      (nreverse
	       (cond
		 ((null (car vertices)) rect)
		 ((null (cdr vertices)) (incf lb) (push (car vertices) rect))
		 (t (loop for w from (point-j (car (last vertices))) 
		       to (point-j (car vertices)) do
			 (incf lb)
			 (push (new-point (point-i (car vertices)) w) rect))
		    rect)))
	      lb)
	     (if (or (null vertices) 
		     (= (point-i (car contour)) (point-i (car vertices))))
		 (fill-rect lb rect (push (car contour) vertices) (cdr contour))
		 (progn
		   (loop for w from (point-j (car (last vertices))) 
		      to (point-j (car vertices)) do
			(incf lb)
			(push (new-point (point-i (car vertices)) w) rect))
		   (fill-rect lb rect (list (car contour)) (cdr contour)))))))
    (let* ((r (/ (* a pi) 180)) (c (cos r)) (s (sin r)))
      (let ((e0 (new-point 0 0))
	    (e1 (%rotate-point c s (new-point 0 (1- l))))
	    (e2 (%rotate-point c s (new-point (1- w) 0)))
	    (e3 (%rotate-point c s (new-point (1- w) (1- l)))))
	(let ((x0 (+ (point-j e0) x)) (y0 (+ (point-i e0) y))
	      (x1 (+ (point-j e1) x)) (y1 (+ (point-i e1) y))
	      (x2 (+ (point-j e2) x)) (y2 (+ (point-i e2) y))
	      (x3 (+ (point-j e3) x)) (y3 (+ (point-i e3) y)))
	  (fill-rect 
	   0 nil nil
	   (sort
	    (remove-duplicates
	     (append
	      (bresenham-line x0 x1 y0 y1)
	      (bresenham-line x0 x2 y0 y2)
	      (bresenham-line x1 x3 y1 y3)
	      (bresenham-line x2 x3 y2 y3)) :test #'equal) #'point<)))))))

(defparameter *lim1* 0.80)
(defparameter *lim2* 0.99)
(defparameter *lim3* 1/24)
(defparameter *ugol* 11.25) ;22.5) ;11.25) ;5.625)

(defun limit (w l a)
  (declare (ignore a))
  (1+ (floor (* w l *lim3*))))
;  (if (member a '(0 90) :test #'=)
;      (1+ (floor (* w l 1/24)))
;      (1+ (floor (* w l 1/12)))))      

; интересно, а если находить контуры и сравнивать их?

(defun sort-rect (type rect)
  (let ((r (sort 
	    rect 
	    (case type 
	      (:hor #'hor-point<)
	      ((:dia101 :dia112 :dia123 :dia135 :dia146 :dia157 :dia168) 
	       #'/-point<)
	      (t #'point<)))))
    (cons 
     (car r) 
     (sort (cdr r) #'< 
	   :key #'(lambda (p) (euclidean-distance p (car r)))))))

(defun rectal-testi 
    (blacks blocks covered rectangle lb lim bl w l a height width y x)
  (declare (optimize (speed 3) (safety 0)))
  (declare (type (unsigned-byte 12) w l height width))
  (declare (type (unsigned-byte 24) lb bl))
  (declare (type (simple-array bit (* *)) blacks covered))
  (labels
      ((set-n-return (rect blacks covered blocks w l a lb bl la)
	 (dolist (p rect)
	   (setf (aref blacks (point-i p) (point-j p)) 0
		 (aref covered (point-i p) (point-j p)) 1))
	 (values (cons (list (if (= w l) :dot (angle->tmp a))
			     w lb (sort rect #'point<)) blocks) 
		 (- bl la))))
    (let ((la 0)
	  (rect (move-and-clip-rectangle lim y x height width rectangle)))
      (declare (type (unsigned-byte 12) la))
      (if (null rect) (values blocks bl)
	  (progn
	    (dolist (p rect) (incf la (aref blacks (point-i p) (point-j p))))
	    (if (> la 0)
		(if (> lim (- lb la))
		    (set-n-return rect blacks covered blocks w l a lb bl la)
		    (let ((lc 0))
		      (declare (type (unsigned-byte 12) lc))
		      (dolist (p rect) 
			(incf lc (aref covered (point-i p) (point-j p))))
		      (if (> lim (- lb la lc))
			  (set-n-return 
			   rect blacks covered blocks w l a lb bl la)
			  (values blocks bl))))
		(values blocks bl)))))))

(defun copy-image (image)
  (declare (optimize (speed 3) (safety 0)))
  (let ((copy (make-array (array-dimensions image) :element-type 'bit)))
    (declare (type (simple-array bit (* *)) copy))
    (dotimes (i (array-dimension image 0) copy)
      (dotimes (j (array-dimension image 1))
	(setf (aref copy i j) (aref image i j))))))

(defun duplicate-image (image)
  (make-array 
   (array-dimensions image)
   :displaced-to 
   (copy-seq (make-array (array-total-size image) :displaced-to image))))

(defparameter *points* t)

(defun %lcd (ww ll bl blacks covered schwarz height width)
  (declare (optimize (speed 3) (safety 0)))
  (declare (type (unsigned-byte 12) ww ll))
  (declare (type (unsigned-byte 24) bl))
  (declare (type (simple-array bit (* *)) blacks covered))
  (let (blocks (blin bl))
    (declare (type (unsigned-byte 24) blin))
    (labels 
	((looper (points)
	    (loop for w from ww above 0 do
		 (loop for l fixnum from ll above (- w 1) do
		      (let ((lb (* w l)))
			(if (and (> (/ bl lb) *lim1*)
				 (if points
				     (= w l)
				     (and (> lb 3) (>= (/ l 1.0 w) 3.0))))
			    (loop for a single-float below 
				 (if points *ugol* 180.0) by *ugol* do
				 (multiple-value-bind
				  (rr lb) (build-rectangle 0 0 w l a)
				  (let ((lim (limit w l a)))
				    (dolist (p schwarz)
				      (multiple-value-setq 
					  (blocks blin)
					(rectal-testi 
					 blacks blocks covered rr lb lim blin 
					 w l a height width 
					 (point-i p) (point-j p)))
				      (if (zerop blin) 
					  (return-from %lcd blocks))))))))))))
      (looper nil)
      (if *points* (looper t))
      blocks)))

(defun haf-transform (image)
  (multiple-value-bind (schwarz bl height width)
      (image->blacks image)
    (sort (%lcd (min height width) (max height width) bl (copy-image image)
		(make-array 
		 (list height width) :element-type 'bit :initial-element 0)
		schwarz height width)
	  #'> :key #'third)))

(defun sb-length (sb) (euclidean-distance (car sb) (car (last sb))))

(defun sb-straight-line-middle (sb)
  ;(mapcar #'average (car sb) (car (last sb)))
  (new-point (/ (+ (point-i (car sb)) (point-i (car (last sb)))) 2)
	     (/ (+ (point-j (car sb)) (point-j (car (last sb)))) 2)))

(defun adjacent-p (x y) ;(> 2 (euclidean-distance x y)))
  (and (>= 1 (abs (- (point-i x) (point-i y))))
       (>= 1 (abs (- (point-j x) (point-j y))))))

(defun sb-average-centre (sb)
  (let ((x 0) (y 0) (l 0))
    (dolist (p sb)
      (incf l)
      (incf x (point-i p))
      (incf y (point-j p)))
    (new-point (/ x l 1.0) (/ y l 1.0))))

(defun sb-position (sb i)
  (let ((sp (member (car sb) i :test #'adjacent-p))
	(ep (member (car (last sb)) i :test #'adjacent-p)))
  (cond
    ((and sp ep) (values :e :s)) 
    (sp :s)
    (ep :e)
    (t
     (let ((pos (/ (euclidean-distance (sb-straight-line-middle i) (car sb)) 
		   (sb-length sb))))
       (cond    
	 ((< pos 0.33) :s)
	 ((<= 0.33 pos 0.66) :m)
	 ((> pos 0.66) :e))))
)))

(defun same-line-p (v a b r)
  (let ((s (car (sort-rect v (list (car a) (car b)))))
	(e (cadr (sort-rect v (list (car (last a)) (car (last b)))))))
    (let ((q1 (- (point-j s) (point-j e))) 
	  (q2 (- (point-i e) (point-i s)))
	  (q3 (- (* (point-i s) (point-j e)) (* (point-i e) (point-j s))))
	  (q4 (euclidean-distance s e))
	  (points (append a b)))
      (let ((vals 
	     (mapcar
	      #'(lambda (p) 
		  (abs (/ (+ q3 (* q1 (point-i p)) (* q2 (point-j p))) q4)))
	      points)))
	(let ((rat (/ (count-if #'(lambda (v) (>= r v)) vals) (length points))))
	  (> rat 0.9))))))

(defun intersecsion (a b) (intersection a b :test #'adjacent-p))
;  (let ((i (intersection a b :test #'equal)))
;    (or i (intersection a b :test #'adjacent-p))))

(defun cross-relation (a b)
  (let ((fa (fourth a)) (fb (fourth b)))
    (let ((i (intersecsion fa fb)))
      (if (and i (not (eq (car a) (car b))) ; забыл, почему нот эк :(
	       (>= (length i) (min (second a) (second b))))
	  (let ((si (sort i #'point<)))
	    (list (sb-position fa si) (sb-position fb si)))))))

(defun block-extension-p (a b)
  (when (intersecsion (fourth a) (fourth b))
    (if (eq :dot (car b)) t
	(unless (eq :dot (car a))
	  (same-line-p 
	   (car a) (fourth a) (fourth b) (min (cadr a) (cadr b)))))))

(defun block-joint (a b)
  (let ((i (sort (intersecsion (fourth a) (fourth b)) #'point<)))
    (when i
      (values 
       (multiple-value-list (sb-position (fourth a) i))
       (multiple-value-list (sb-position (fourth b) i)) i))))

(defparameter *curves*
  '((165 :s 20 :s)
    (20 :e 45 :s)
    (45 :e 90 :s)
    (90 :e 120 :s)
    (120 :e 135 :s)
    (135 :e 0 :e)
    (165 :e 120 :s)
    (120 :e 90 :s)
    (90 :e 45 :s)
    (45 :e 0 :s)
    (70 :e 45 :s)
    (45 :e 0 :s)
    (0 :s 135 :s)
    (135 :e 90 :s)
    (90 :e 135 :s)
    (135 :e 20 :e)
    (0 :e 20 :s)
    (20 :e 70 :s)))

(defparameter *not-curves*
  '((0 :e 90 :s)
    (0 :s 90 :s)
    (112 :s 67 :s)
    (112 :s 90 :s)
    (90 :s 67 :s)
    (90 :e 112 :e)
    (67 :e 112 :e)
    (0 :e 135 :s)
    (135 :e 0 :s)))

(defun curve? (a1 p1 a2 p2)
  (when (> a1 a2) (rotatef a1 a2) (rotatef p1 p2))
  (cond
    ((and (eq p1 :e) (eq p2 :s))
;    (format t "1 ~a ~a ~a ~a ~a~%" a1 p1 a2 p2 (< 90 (- 180 (- a2 a1))))
     (< 90 (- 180 (- a2 a1)))) ; тупой угол
    ((and (eq p1 :s) (eq p2 :e))
;    (format t "2 ~a ~a ~a ~a ~a~%" a1 p1 a2 p2 (< 90 (- 180 (- a2 a1))))
     (< 90 (- 180 (- a2 a1)))) ; тупой угол
    ((and (eq p1 :s) (eq p2 :s))
;     (format t "3 ~a ~a ~a ~a ~a~%" a1 p1 a2 p2 (< 90 (- a2 a1)))
     (< 90 (- a2 a1))) ; тупой угол
    ((and (eq p1 :e) (eq p2 :e))
;     (format t "4 ~a ~a ~a ~a ~a~%" a1 p1 a2 p2 (< 90 (- a2 a1)))
     (< 90 (- a2 a1))) ; тупой угол
;    (t (format t "0 ~a ~a ~a ~a ~a~%" a1 p1 a2 p2 'nol))
))

(defun block-same-p (a b)
  (multiple-value-bind (p1 p2 i) (block-joint a b)
    (let ((a1 (get (car a) :angle)) (a2 (get (car b) :angle)))
      (cond
	((null i) nil)
	((eq :dot (car b)) t)
	((eq :dot (car a)) nil)
	(t (some
	    #'identity
	    (mapcar 
	     #'(lambda (p) (curve? a1 (car p) a2 (cadr p)))
	     (interleave #'list nil p1 p2))))))))

(defun group (pred lst &key (key #'identity) every sorter rem not-every)
  "groups list elements according to pred[icate]: (1 1 2 1) -> ((1 1 1) (2)); if EVERY is true, ensures that PRED is true for any pair within a group; if EVERY is funcallable, calls it with group so far and a potential member, then ANDs the result with that of predicate, which is called for every pair within a group, unless NOT-EVERY is true; REM filters potential group members."
  (labels
      ((grpr (pred test key lst group rest)
	 (if (null lst) (nreverse group) ;(values (nreverse group) (nreverse rest))
	     (if (null group) 
		 (grpr pred test key (cdr lst) (list (car lst)) rest)
		 (if (cond
		       ((or (functionp test) (fboundp test))
			(and
			 (some 
			  #'(lambda (e)
			      (funcall 
			       pred (funcall key e) (funcall key (car lst))))
			  group)
			 (funcall test group (car lst))))
		       (test
			(every
			 #'(lambda (e)
			     (funcall 
			      pred (funcall key e) (funcall key (car lst))))
			 group))
		       (t
			(some
			 #'(lambda (e)
			     (funcall 
			      pred (funcall key e) (funcall key (car lst))))
			 group)))
		     (if (and test (not not-every))
			 (grpr pred test key (cdr lst)
			       (cons (car lst) group) rest)
			 (grpr pred test key (append (cdr lst) rest) 
			       (cons (car lst) group) nil))
		     (grpr 
		      pred test key (cdr lst) group (cons (car lst) rest))))))
       (groupr (pred lst key every rem groups)
	 (let ((ls (if rem (funcall rem (car lst) (cdr lst))
		       (cdr lst))))
	   (let ((group (grpr pred every key ls (list (car lst)) nil)))
	     (let ((rest (set-difference lst group)))
	       (if (null rest) 
		   (nreverse (if (null group) groups (cons group groups)))
		   (groupr 
		    pred (if sorter (sort rest sorter) (nreverse rest)) ; nreverse??
		    key every rem
		    (if (null group) groups (cons group groups)))))))))
    (groupr pred lst key every rem nil)))

(defun segment (blocks)
  (let ((min 360) (max 0) sec)
    (dolist (b blocks)
      (let ((a (get (car b) :angle)))
	(when a 
	  (setf max (max max a))
	  (when (> min a)
	      (setf sec (if (null sec) a min))
	      (setf min a)))))
    (if (and (zerop min) (> max 90))
      (setf min sec max (max 180 max)))
    (list min max)))

(defun curve-p (seg) (< 30 (- (cadr seg) (car seg))))

(defun curve-ends (blocks)
  (multiple-value-bind (g v) (blocks->graph blocks)
    (let (e (h (worsen-visibility g)) (d (array-dimension g 0)))
;       (format t "~a~%" h)
      (dotimes (i d)
	(let ((term
	       (remove-duplicates
		(remove-if 
		 #'null
		 (loop for j below d collect 
		      (if (and (listp (aref h i j))) (car (aref h i j))))))))
	  (if (= 1 (length term)) (push (list (car term) (svref v i)) e))))
      (let* ((se (sort e #'point< :key 
		       #'(lambda (b) (sb-average-centre (fourth (cadr b))))))
	     (n (car se)) (k (car (last se))))
;	(format t "~a~%" se)
	(list ; тут бы, конечно, как-то учесть, что могут оба конца пересекаться
	 (if (eq :s (car n)) (reverse (fourth (cadr n))) (fourth (cadr n)))
	 (if (eq :e (car k)) (reverse (fourth (cadr k))) (fourth (cadr k))))))))

(defun flatten-blocks (blocks)
  (reduce #'(lambda (a b) (union a b :test #'equal)) blocks :key #'fourth))

(defun flatten-curve (blocks)
  (remove-duplicates 
   (let ((term (curve-ends blocks)) (points (flatten-blocks blocks)))
     (append
      (car term)
      (set-difference 
       points (union (car term) (cadr term) :test #'equal) :test #'equal)
      (cadr term)))
   :test #'equal :from-end t))

(defun make-block (type width points) (list type width (length points) points))

(defun compare (e a b)
  (cond
    ((= a b) '=)
    ((> e (abs (- a b))) '=)
    ((> a b) '>)
    ((< a b) '<)))

(defun point-compare (e a b)
  (list
   (compare e (point-i a) (point-i b))
   (compare e (point-j a) (point-j b))))

(defun curve-type (points)
  (let* ((s (car points)) (e (car (last points)))
	 (d (new-point (/ (+ (point-i s) (point-i e)) 2.0)
		       (/ (+ (point-j s) (point-j e)) 2.0)))
	 (q1 (length 
	      (remove-if-not 
	       #'(lambda (p) (and (> (point-i d) (point-i p))
				  (> (point-j p) (point-j d))))
	       points)))
	 (q2 (length
	      (remove-if-not 
	       #'(lambda (p) (and (> (point-i p) (point-i d))
				  (> (point-j p) (point-j d))))
	       points)))
	 (q3 (length
	      (remove-if-not 
	       #'(lambda (p) (and (> (point-i p) (point-i d))
				  (> (point-j d) (point-j p))))
	       points)))
	 (q4 (length 
	      (remove-if-not 
	       #'(lambda (p) (and (> (point-i d) (point-i p))
				  (> (point-j d) (point-j p))))
	       points))))
;    (format t "~a ~a ~a q ~a ~a ~a ~a~%" s d e q1 q2 q3 q4)
    (if (> (point-j s) (point-j e))
	(if (> q4 q2) :curve-q4 :curve-q2)
	(if (> q3 q1) :curve-q3 :curve-q1))))

(defun quarter-curve (blocks) ; может, вообще сортировать блоки и следовать им?
  (group 
   #'block-same-p
   blocks
   :every 
   #'(lambda (g e)
;       (format t "~%~a~%" (mapcar #'car g))
       (let* ((a1 (get (car e) :angle)) (s (signum (- 90 a1))))
;	      (x (segment (cons e g))))
;	 (and
;	  (> 91 (- (cadr x) (car x)))
	  (every
	   #'(lambda (j)
	       (let* ((a2 (get (car j) :angle)) (z (signum (- 90 a2))))
		 (and
		  (not (= a1 a2))
		  (or (zerop s) (zerop z) (= z s)))))
	   g)))))
;)
   

(defun curve-blocks (blocks)
  (sort
   (apply
    #'append
    (mapcar
     #'(lambda (g)
;	 (format t "~%GROUP:~%~a~%" g)
	 (if (or (curve-p (segment g)) 
		 (curve-p (segment (subseq g 0 (floor (length g) 2)))) ; d.h.
		 (curve-p (segment (subseq g (floor (length g) 2)))))
	     (mapcar
	      #'(lambda (j)
;		  (format t "~%j ~a~%" j)
		  (let ((points (flatten-curve j)))
		    (when (point< (car (last points)) (car points))
		      (setf points (nreverse points))) ; dirty hack
;		    (format t "t ~a ~a~%" (curve-type points) (cadar j))
		    (make-block (curve-type points) (cadar j) points)))
	      (quarter-curve g))
	     (list 
	      (make-block 
	       (caar g) (cadar g) (sort-rect (caar g) (flatten-blocks g))))))
     (group #'block-same-p blocks)))
   #'> :key #'third))

(defun merge-blocks (blocks)
   (mapcar
    #'(lambda (g) 
	(make-block (caar g) (cadar g) (sort-rect (caar g) (flatten-blocks g))))
    (group 
     #'block-extension-p blocks 
     :sorter 
     #'(lambda (a b) (> (third a) (third b)))
     :every  
     #'(lambda (g e) 
	 (let ((s (segment (cons e g)))) (> 90 (- (cadr s) (car s))))))))

(defun sort-blocks (blocks)
  (sort blocks 
	#'(lambda (a b)
;	    (let ((aa (get (car a) :angle))
;		  (ba (get (car b) :angle)))
	      (cond
		((eq :dot (car a)) nil)
		((eq :dot (car b)) t)
		(t (point< (first (fourth a)) (first (fourth b))))))));)
		;((= aa ba) (point< (first (fourth a)) (first (fourth b))))
		;(t (< aa ba)))))))
      
(defun mount-blocks (blocks &optional skcolb occ)
  (if (null blocks) (sort-blocks skcolb)
      (if (if (and (eq :dot (car (car blocks)))
		   (= 1 (second (car blocks))))
	      (subsetp (fourth (car blocks)) occ :test #'equal)
	      (> 3 ; поставить параметр
		 (length
		  (set-difference (fourth (car blocks)) occ :test #'equal))))
	  (mount-blocks (cdr blocks) skcolb occ)
	  (mount-blocks 
	   (cdr blocks) 
	   (push (car blocks) skcolb)
	   (union (fourth (car blocks)) occ :test #'equal)))))

(defun blk-rel (a b)
  (let ((ep (/ (min (sb-length (fourth a)) (sb-length (fourth b))) 5)))
    (list 
     (if (eq a b) (car a) (cross-relation a b))
     (point-compare ep (first (fourth a)) (first (fourth b)))
     (point-compare ep (first (fourth a)) (car (last (fourth b))))
     (point-compare ep (car (last (fourth a))) (first (fourth b)))
     (point-compare ep (car (last (fourth a))) (car (last (fourth b)))))))

(defun image->blocks (img) 
  (curve-blocks ;mount-blocks 
   (mount-blocks ;curve-blocks 
    (merge-blocks 
     (haf-transform img)))))

(defun blocks->graph (blocks) (make-graph #'blk-rel blocks))

(defun graphicalize (image) (blocks->graph (image->blocks image)))

;(defparameter *subst-table* (make-hash-table))
(defparameter *overlay* nil)

(defun mklist (obj) (if (listp obj) obj (list obj)))

(defun interleave (fun args prolist &rest prolists)
  "такая идиома уже есть, называется alexandria:map-product"
  (if (null prolist)
      (list (apply fun (reverse args)))
      (apply
       #'append
       (mapcar
	#'(lambda (elt)
	    (apply 
	     #'interleave fun (cons elt args) (car prolists) (cdr prolists)))
	prolist))))

(defun worsen-visibility (graph)
  (let ((tmp (make-array (array-dimensions graph))))
    (dotimes (i (array-dimension graph 0))
      (dotimes (j (array-dimension graph 1))
	(setf (aref tmp i j) (car (aref graph i j)))))
    tmp))
      
(defun shrink-lst (f lst)
  (remove-duplicates 
   (mapcar 
    #'(lambda (p) (new-point (round (point-i p) f) (round (point-j p) f))) lst)
   :test #'equal))

(defun shrink (f i) (instantiate-rectangle (shrink-lst f (image->blacks i))))

; http://homepages.inf.ed.ac.uk/rbf/HIPR2/
; http://homepages.inf.ed.ac.uk/cgi/rbf/CVONLINE/

(defparameter octogonal-skeleton
  #(#2A((0   0   0  )
	(nil 1   nil)
	(1   1   1  ))
    #2A((nil 0   0  )
	(1   1   0  )
	(nil 1   nil))
    #2A((1   nil 0  )
	(1   1   0  )
	(1   nil 0  ))
    #2A((nil 1   nil)
	(1   1   0  )
	(nil 0   0  ))
    #2A((1   1   1  )
	(nil 1   nil)
	(0   0   0  ))
    #2A((nil 1   nil)
	(0   1   1  )
	(0   0   nil))
    #2A((0   nil 1  )
	(0   1   1  )
	(0   nil 1  ))
    #2A((0   0   nil)
	(0   1   1  )
	(nil 1   nil))))

(defparameter pruning-set
  #(#2A((0   0   0  )
	(0   1   0  )
	(0   nil nil))
    #2A((0   0   0  )
	(0   1   0  )
	(nil nil 0  ))
    #2A((0   0   0  )
	(nil 1   0  )
	(nil 0   0  ))
    #2A((nil 0   0  )
	(nil 1   0  )
	(0   0   0  ))
    #2A((nil nil 0  )
	(0   1   0  )
	(0   0   0  ))
    #2A((0   nil nil)
	(0   1   0  )
	(0   0   0  ))
    #2A((0   0   nil)
	(0   1   nil)
	(0   0   0  ))
    #2A((0   0   0  )
	(0   1   nil)
	(0   0   nil))))

(defparameter horizontal-lines
  #2A((-1 -1 -1)
      ( 2  2  2)
      (-1 -1 -1)))

(defparameter vertical-lines
  #2A((-1  2 -1)
      (-1  2 -1)
      (-1  2 -1)))

(defparameter sobel-gx
  #2A((-1  0  1)
      (-2  0  2)
      (-1  0  1)))

(defparameter sobel-gy
  #2A(( 1  2  1)
      ( 0  0  0)
      (-1 -2 -1)))

(defun check-pixel (image i j selt y g fixel)
  (if (not (aref selt y g)) t
      (if (array-in-bounds-p image i j) (= (aref image i j) (aref selt y g))
	  (if fixel (= fixel (aref selt y g))))))

(defun check-structuring-element (image i j selt fixel y g)
  (if (= y (array-dimension selt 0)) t
      (if (check-pixel image (+ i y) (+ j g) selt y g fixel)
	  (if (= (+ 1 g) (array-dimension selt 1))
	      (check-structuring-element image i j selt fixel (+ 1 y) 0)
	      (check-structuring-element image i j selt fixel y (+ 1 g))))))

(defun apply-structuring-element (selt image &key origin (fixel 0) (nvalue 0))
  (let ((oi (if origin (car origin) 
		(1- (ceiling (/ (array-dimension selt 0) 2.0)))))
	(oj (if origin (cdr origin) 
		(1- (ceiling (/ (array-dimension selt 1) 2.0)))))
	(tmp (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) tmp)
      (dotimes (j (array-dimension image 1))
	(if (check-structuring-element image (- i oi) (- j oj) selt fixel 0 0)
	    (setf (aref tmp i j) nvalue)
	    (setf (aref tmp i j) (aref image i j)))))))

(defun isometric-images-equal (i1 i2 i j)
  (if (= i (array-dimension i1 0)) t
      (if (= (aref i1 i j) (aref i2 i j))
	  (if (= (+ 1 j) (array-dimension i1 1))
	      (isometric-images-equal i1 i2 (+ 1 i) 0)
	      (isometric-images-equal i1 i2 i (+ 1 j))))))

(defun image-equal (i1 i2)
  (and (equal (array-dimensions i1) (array-dimensions i2))
       (isometric-images-equal i1 i2 0 0)))

(defun apply-se-set (sset image &key origin (fixel 0) (nvalue 0) (si 0))
  (if (>= si (array-dimension sset 0)) image
      (apply-se-set
       sset
       (apply-structuring-element 
	(aref sset si) image
	:origin origin :fixel fixel :nvalue nvalue)
       :origin origin :fixel fixel :nvalue nvalue :si (+ 1 si))))

(defun apply-sse-to-convergence (sset image &key origin (fixel 0) (nvalue 0))
  (let ((tmp (apply-se-set 
	      sset image :origin origin :fixel fixel :nvalue nvalue)))
    (if (image-equal image tmp) tmp
	(apply-sse-to-convergence 
	 sset tmp :origin origin :fixel fixel :nvalue nvalue))))
  
(defun apply-se-to-convergence (selt image &key origin (fixel 0) (nvalue 0))
  (let ((tmp (apply-structuring-element
	      selt image :origin origin :fixel fixel :nvalue nvalue)))
    (if (image-equal image tmp) tmp
	(apply-se-to-convergence 
	 selt tmp :origin origin :fixel fixel :nvalue nvalue))))

(defun prune (image) (apply-se-set pruning-set image))

(defun reap-meat (image) (apply-sse-to-convergence octogonal-skeleton image))

(defun convolve (kernel image)
  (let ((tmp (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) tmp)
      (dotimes (j (array-dimension image 1))
	(dotimes (k (array-dimension kernel 0))
	  (dotimes (l (array-dimension kernel 1))
	    (incf (aref tmp i j) 
		  (* (cref image (+ i k) (+ j l)) (aref kernel k l)))))))))

(defun render-svg-path (svg &optional path paths)
  (if (null svg) (apply #'append paths)
      (ecase (caar svg)
	(:M (render-svg-path 
	    (cdr svg) (list (cdr (car svg)))
	    (if (null (cdr path)) 
		paths 
		(push (bresenham-line 
		       (car (second path)) (car (first path))
		       (cadr (second path)) (cadr (first path)))
		      paths))))
	(:L (render-svg-path 
	    (cdr svg) (list (cdr (car svg)))
	    (if (null path)
		paths 
		(push (bresenham-line 
		       (car (car path)) (car (cdr (car svg)))
		       (cadr (car path)) (cadr (cdr (car svg))))
		      paths)))))))

(defun surround-point (radius point)
  (let (environs)
    (loop for i from (- radius) to radius do
	 (loop for j from (- radius) to radius do
	      (push (new-point (+ i (point-i point)) 
				    (+ j (point-j point))) 
		    environs)))
    environs))

(defun widen-lines (radius points)
  (apply #'append (mapcar #'(lambda (p) (surround-point radius p)) points)))

(defun unify-type (type)
  (case type
    ((:dia011
      :dia022
      :dia033
      :dia045
      :dia056
      :dia067
      :dia078) :dia045)
    ((:dia101
      :dia112
      :dia123
      :dia135
      :dia146
      :dia157
      :dia168) :dia135)
    (t type)))

(defun sort-vertices (points graph)
  (labels ((get-op (a b graph)
	     (nth (first a) (nth (+ (* 2 (third a)) (third b)) 
				 (cdr (aref graph (second a) (second b)))))))
    (sort (group #'(lambda (a b) (eq '= (get-op a b graph))) points :every t)
	  #'(lambda (a b) (eq '< (get-op a b graph))) :key #'car)))

(defun realize-vertices (groups)
  (let (res)
    (loop for g in groups for n from 0 do
	 (dolist (p g)
	   (push (cons p n) res)))
    (nreverse res)))

(defun graph->lines (graph)
  (let (lines xpoints ypoints)
    (dotimes (i (array-dimension graph 0))
      (let ((y0 (list 0 i 0))
	    (y1 (list 0 i 1))
	    (x0 (list 1 i 0))
	    (x1 (list 1 i 1)))
	(push y0 ypoints)
	(push y1 ypoints)
	(push x0 xpoints)
	(push x1 xpoints)
	(push `((,y0 ,x0) (,y1 ,x1)) lines)))
    (let ((points (append (realize-vertices (sort-vertices xpoints graph))
			  (realize-vertices (sort-vertices ypoints graph)))))
      (sort
       (mapcar
	#'(lambda (line)
	    (mapcar
	     #'(lambda (end)
		 (mapcar 
		  #'(lambda (a) (cdr (assoc a points :test #'equal)))
		  end))
	     line))
	lines)
       #'(lambda (a b) ; заменить на point< ?
	   (or (< (car a) (car b)) 
	       (and (= (car a) (car b)) (< (cadr a) (cadr b)))))
       :key #'car))))
  
(defun lines->svg-path (zoom lines)
  (let (segs) ; поставить сюда T для кривых?
    (dolist (line lines)
      (push `(:M ,(* zoom (cadar line)) ,(* zoom (caar line))) segs)
      (push `(:L ,(* zoom (cadadr line)) ,(* zoom (caadr line))) segs))
    (nreverse segs)))

(defun image->lines (image) (graph->lines (graphicalize image)))

(defun lines->image (lines)
  (instantiate-rectangle (render-svg-path (lines->svg-path 5 lines))))

(defun redraw (image)
  (instantiate-rectangle 
   (render-svg-path (lines->svg-path 5 (image->lines image)))))

; спер из CLOCC
(defun assignment (cost-mx &key (out *standard-output*))
  "Solve the assignment problem using Munkres' algorithm
\(AKA Hungarian Algorithm, AKA Bipartite Minimum-Weight Matching).
Returns the total cost and two assignment vectors: X->Y and Y->X."
  (loop :with x-count = (array-dimension cost-mx 0)
    :and y-count = (array-dimension cost-mx 1)
    ;; tentative least weights ("D label"):
    :with x-tlv = (make-array x-count :initial-element nil)
    :and y-tlv = (make-array y-count :initial-element nil)
    :and x-matching = (make-array x-count :initial-element nil)
    :and y-matching = (make-array y-count :initial-element nil)
    :and s-x = (make-array x-count :element-type 'bit)
    :and s-y = (make-array y-count :element-type 'bit)
    :and x-next = (make-array x-count :initial-element nil)
    :and y-next = (make-array y-count :initial-element nil)
    :and cost = 0
    ;;:and eps = (loop :with min :and max :for i :from 0 :to (1- x-count) :do
    ;;             (loop :for j :from 0 :to (1- y-count)
    ;;               :for cost = (abs (aref cost-mx i j)) :do
    ;;               (when (or (null min) (< cost min)) (setq min cost))
    ;;               (when (or (null max) (> cost max)) (setq max cost)))
    ;;             :finally (return (* (/ (- max min) (+ max min))
    ;;                                 )))
    :for iteration-count :upfrom 1
    :do (when out
          (format out "~&~S(~:D): ~S~%   ~S~%   ~S~%" 'assignment
                  iteration-count cost x-matching y-matching))
    ;; s-x == all free X vertices
    (dotimes (i x-count)
      (if (svref x-matching i)
          (setf (sbit s-x i) 0
                (svref x-tlv i) nil)
          (setf (sbit s-x i) 1
                (svref x-tlv i) 0)))
    (fill y-tlv nil)
    (fill x-next nil) (fill y-next nil)
    (loop :while (find 1 s-x) :do
      (fill s-y 0)
      (dotimes (i x-count)           ; for all X
        (unless (zerop (sbit s-x i)) ; in Sx
          (let ((dx (svref x-tlv i)))
            (dotimes (j y-count)     ; for all Y
              (unless (eql j (svref x-matching i)) ; (x,y) not in M
                (let ((dy (svref y-tlv j))
                      (ndy (+ dx (aref cost-mx i j))))
                  (when (or (null dy) (< ndy dy))
                    (setf (svref y-tlv j) ndy
                          (sbit s-y j) 1
                          (svref y-next j) i))))))))
      (fill s-x 0)
      (dotimes (j y-count)           ; for all Y
        (unless (zerop (sbit s-y j)) ; in Sy
          (let ((i (svref y-matching j)))
            (when i             ; (x,y) in M
              (let ((dx (svref x-tlv i))
                    (ndx (- (svref y-tlv j) (aref cost-mx i j))))
                (when (or (null dx) (< ndx dx))
                  (setf (svref x-tlv i) ndx
                        (sbit s-x i) 1
                        (svref x-next i) j))))))))
    ;; free vertex with min TLV
    (let ((min nil) (pos nil) (weight 0) (len 0))
      (dotimes (j y-count)
        (unless (svref y-matching j)
          (let ((dy (svref y-tlv j)))
            (when (or (null min) (< dy min))
              (setq min dy pos j)))))
      (unless min
        (return-from assignment (values cost x-matching y-matching)))
      ;; augment M with the path associated with POS
      (loop :for i = (svref y-next pos) :for j = (svref x-next i)
        :do (setf (svref y-matching pos) i
                  (svref x-matching i) pos)
        (incf weight (aref cost-mx i pos))
        (incf len)
        :unless j :return nil
        :do (setq pos j) (decf weight (aref cost-mx i j)) (incf len))
      (unless (= min weight) ; (< (/ (- min weight) (+ min weight)) eps)
        (warn "~S: rounding error ~5F * float-epsilon"
              'assignment (/ (- min weight) (+ (abs min) (abs weight))
                             (etypecase min
                               (short-float short-float-epsilon)
                               (single-float single-float-epsilon)
                               (double-float double-float-epsilon)
                               (long-float long-float-epsilon)))))
      (when out
        (format out "~&   => augmenting path: len=~:D weight=~S~%" len weight))
      (incf cost weight))))

(defun make-cost-matrix (a b change-cost-fn add-cost-fn)
  (let* ((al (length a)) (bl (length b)) (sl (+ al bl))
	 (m (make-array (list sl sl))))
    (dotimes (i al)
      (dotimes (j bl)
	(setf (aref m i j) (funcall change-cost-fn (svref a i) (svref b j))))
      (dotimes (j al)
	(setf 
	 (aref m i (+ bl j)) 
	 (if (= i j) (funcall add-cost-fn (svref a i)) most-positive-fixnum))))
    (dotimes (i bl)
      (dotimes (j bl)
	(setf 
	 (aref m (+ al i) j) 
	 (if (= i j) (funcall add-cost-fn (svref b i)) most-positive-fixnum)))
      (dotimes (j al)
	(setf (aref m (+ i al) (+ j bl)) 0)))
    m))

(defun regraph (graph)
  (let ((v (make-array (array-dimension graph 0))))
    (dotimes (i (array-dimension graph 0))
      (let (e)
	(dotimes (j (array-dimension graph 0))
	  (if (and (/= i j) (aref graph i j))
	      (push (aref graph i j) e)))
	(setf (svref v i) (cons (aref graph i i) (apply #'vector e)))))
    v))

(declaim (inline edge-cost))
(defun edge-cost (a b) (if (equal a b) 0 10))

(declaim (inline edge-price))
(defun edge-price (a) (declare (ignore a)) 10)

(defun vertex-cost (a b)
  (+ 
   (if (eq (unify-type (car a)) (unify-type (car b))) 0 10)
   (assignment (make-cost-matrix (cdr a) (cdr b) #'edge-cost #'edge-price)
		 :out nil)))

(declaim (inline vertex-price))
(defun vertex-price (a) (case (car a) (:dot 10) (t 100)))
		  
(defun graph-cost (a b)
  (assignment (make-cost-matrix a b #'vertex-cost #'vertex-price) :out nil))

(defun graph43 (g) (regraph (worsen-visibility g)))

(defun grafs->tree (lst)
  (grow 
   (make-instance 'burkhard-keller-tree)
   (mapcar #'(lambda (p) (cons (graph43 (cadr p)) (car p))) lst)
   #'graph-cost :key #'car))

(defun add23 (tree image char)
  (inoculate 
   tree (cons (graph43 (graphicalize image)) char) #'graph-cost :key #'car))

(defun file->tree (name) (grafs->tree (with-open-file (f name) (read f))))

(defun comb-char-list (lst)
  (sort
;  (let ((l
   (remove-duplicates 
    lst 
    :test 
    #'(lambda (a b) 
	(and 
	 (char= (car a) (car b)) 
	 (zerop (graph-cost (graph43 (cadr a)) (graph43 (cadr b))))))) ;))
;   (format t "~a~%" (length l))
;   (sort l 
   #'char< :key #'car))
;)

(defun grafs->file (lst name)
  (let ((*print-pretty* t))
    (with-open-file (f name :direction :output :if-exists :supersede :if-does-not-exist :create)
      (format f "~s" (comb-char-list lst)))))

(defun gather-fruit (tree graph &key (cost 150) char)
  (remove-duplicates
   (mapcar 
    #'(lambda (p)
	(let ((c (cdr (tree-fruit (second p)))))
	  (if char c (list (car p) c))))
    (sort (gather tree (graph43 graph) #'graph-cost cost :key #'car) 
	  #'< :key #'car))
   :from-end t :test #'equal :key (if char #'identity #'cadr)))

; спер из CLOCC
(defun subsets (set)
  "return a list of all subsets of the given set (represented as a list)"
  (let ((first (first set)) (rest (rest set)))
    (if rest
	(let ((others (subsets rest)))
	  (nconc others
		 (mapcar (lambda (subset)
			   (cons first subset))
			 others)))
	(list nil (list first)))))

; http://www.cs.cmu.edu/afs/cs/project/ai-repository/ai/lang/lisp/doc/notes/set_part.txt
(defun part (set)
  (if (null set) (list nil)
      (loop for partition in (part (cdr set))
	 nconc
	   (cons `((,(car set)) ,@partition)
		 (loop for piece in partition
		    for i from 0
		    collect
		      (nconc `(,(cons (car set) piece))
			     (subseq partition 0 i)
			     (subseq partition (1+ i))))))))

; http://www.informatik.uni-ulm.de/ni/Lehre/WS03/DMM/Software/partitions.pdf

(defun next-partition (n k m)
  (loop for i from (1- n) downto 1 do
       (when (<= (svref k i) (svref m (1- i)))
	 (setf (svref m i) (max (svref m i) (incf (svref k i))))
	 (loop for j from (1+ i) below n do
	      (setf (svref k j) (svref k 0) (svref m j) (svref m i)))
	 (return-from next-partition t)))
  nil)

(defun instantiate-partition (n k set)
  (let ((p (make-array n :initial-element nil)))
    (loop for e in set for i across k do (push e (svref p i)))
    (delete-if #'null (map 'list #'nreverse p))))

(defun do-partitions (set fn)
  (let* ((n (length set))
	 (k (make-array n :initial-element 0))
	 (m (make-array n :initial-element 0)))
    (funcall fn (instantiate-partition n k set))
    (loop while (next-partition n k m) do (funcall fn (instantiate-partition n k set)))))

(defun parth (set)
  (apply
   #'append
   (loop for i below (length set) collect 
	(if (zerop i)
	    (list (list set))
	    (mapcar
	     #'(lambda (p) (append (list (butlast set i)) p))
	     (parth (last set i)))))))
	   
(defun hor-extent (ss)
  (let ((x0 most-positive-fixnum) (x1 0))
    (dolist (b ss (list x0 x1))
      (setf x0 (min x0 (point-j (car (fourth b))))
	    x1 (max x1 (point-j (car (last (fourth b)))))))))

(defun range-intersection (r1 r2)
  (if (> (car r1) (car r2)) (rotatef r1 r2))
  (if (> (car r2) (cadr r1)) nil
      (- (min (cadr r1) (cadr r2)) (car r2))))

(defun range-length (r) (- (cadr r) (car r)))

(defun imagine-boxes (blocks)
  (parth blocks))

(defun image->sorted-blocks (image)
  (sort (image->blocks image) #'hor-point< 
	:key #'(lambda (b) (sb-average-centre (fourth b)))))

(defun choose-char-chains (tree img &key (cost 30))
  (let* ((blocks (image->sorted-blocks img)) (n (length blocks)))
    (mapcar
     #'(lambda (bs)
	 (let ((w (min 1 (abs (- 1 (/ (apply #'+ (mapcar #'length bs)) n))))))
	   (mapcar
	    #'(lambda (b)
		(mapcar 
		 #'(lambda (p) (list w (/ (car p) cost) (cadr p)))
		 (gather-fruit tree (blocks->graph b) :cost cost)))
	    bs)))
     (imagine-boxes blocks))))

(defun guess-char-chains (tree img &key (cost 30))
  (let ((h (array-dimension img 0)) (w (/ (array-dimension img 1) 6)))
    (remove-if
     #'(lambda (c) (some #'null c))
     (mapcar
      #'(lambda (ps)
	  (mapcar
	   #'(lambda (b)
	       (mapcar 
		#'(lambda (p) (list 0 (/ (car p) cost) (cadr p)))
		(gather-fruit 
		 tree 
		 (graphicalize 
		  (copy-rect 
		   (list 0 (floor (* w (1- (car b)))) 
			 h (floor (* w (car (last b)))))
		   img)) 
		 :cost cost)))
	   ps))
      (parth '(1 2 3 4 5 6))))))


(defun char-chains-p (tree img &key (cost 30))
  (let ((harvest (gather-fruit tree (graphicalize img) :cost cost)))
    (if harvest
	`((,(mapcar
	     #'(lambda (p) (list 0 (/ (car p) cost) (cadr p)))
	     harvest))))))

(defun link-chain (next chain)
  (second
   (reduce 
    #'(lambda (a b) 
	(list a (append (mapcar #'(lambda (e) (cons e (car b))) a) (second b))))
    chain :from-end t :initial-value next)))
  
(defun start-chains (&rest chains)
  (second
   (reduce
    #'(lambda (a b)
	(list
	 (list (apply #'append (mapcar #'car a)))
	 (append
	  (apply
	   #'append (mapcar #'(lambda (chain) (link-chain (first b) chain)) a))
	  (second b))))
    (append `(((( ,(list 0 0 #\Space) )))) chains) 
    :from-end t :initial-value (let ((e '(0 0 #\Space))) `((( ,e )) ((,e)))))))

(defun char-cost (freqs l k)
  (+ 0.2
     (* 0.3 
	(if (not (hash-table-p freqs)) 0
	    (let ((s (make-string 2)))
	      (setf (char s 0) (char-downcase (third l))
		    (char s 1) (char-downcase (third k)))
	      (gethash s freqs 1))))
     (* 0.5 (/ (+ (second l) (second k)) 2))))
;       (* 0.1 (/ (+ (first l) (first k)) 2)))))

(defun chains->graph (freqs links)
  (let ((ln (make-hash-table)) (len 0))
    (dolist (l links) (setf (gethash (car l) ln) len) (incf len))
    (let ((graph (make-array (list len len) :initial-element nil))
	  (vect (make-string len :initial-element #\Space)))
      (dolist (l links (values graph vect))
	(let* ((c (car l)) (i (gethash c ln)))
	  (setf (char vect i) (third c))
	  (dolist (k (cdr l))
	    (setf (aref graph i (gethash k ln)) (char-cost freqs c k))))))))

(defun shortest-chain (freqs &rest chains)
  (multiple-value-bind (g v) 
      (chains->graph freqs (apply #'start-chains chains))
    (string-trim " " (map 'string #'(lambda (i) (char v i))
			  (diter (dijkstra-shortest-path g))))))

; еще из deeplome

(defun blank-p (i max left getter)
  (if (= left 0) nil
    (if (= i max) t (blank-p (1+ i) max (- left (funcall getter i)) getter))))

(defun blank-line-p (image rect ln &optional (div 40))
  (blank-p 
   (second rect) (fourth rect) 
   (ceiling 
    (/ (min (- (third rect) (first rect)) (- (fourth rect) (second rect))) div))
   #'(lambda (i) (aref image ln i))))

(defun blank-column-p (image rect cn tol)
  (blank-p (first rect) (third rect) tol #'(lambda (i) (aref image i cn))))

(defun cut-into-lines (rect image &optional (div 40))
  (labels 
   ((linec
     (image rect i i1 lines)
     (if i
	 (if (= i (third rect))
	     (reverse (if i1 (push (list i1 (second rect)
					 i (fourth rect)) lines) lines))
	   (if (blank-line-p image rect i div)
	       (if i1 (linec image rect (1+ i) nil 
			     (push (list i1 (second rect)
					 i (fourth rect)) lines))
		 (linec image rect (1+ i) nil lines))
	     (if i1 (linec image rect (1+ i) i1 lines)
	       (linec image rect (1+ i) i lines))))
       (linec image rect (first rect) i1 lines))))
   (linec image (if rect rect (append '(0 0) (array-dimensions image)))
	   nil nil nil)))

(defun trim-box (rect image)
  (let* ((start (do ((i (first rect) (1+ i))) 
		    ((or (= i (third rect)) 
			 (not (blank-line-p image rect i))) i)))
	 (stop (do ((i (1- (third rect)) (1- i))) 
		   ((or (> start i) (not (blank-line-p image rect i))) i))))
    (if (> start stop) '(0 0 0 0) 
      `(,start ,(second rect) ,(1+ stop) ,(fourth rect)))))

(defun cut-into-boxes (rect image &optional (tol 1) (space 15))
  (labels
   ((boxec
     (image rect j j1 s1 boxes)
     (if j
	 (if (= j (fourth rect))
	     (mapcar
	      #'(lambda (r) (trim-box r image))
	       (reverse (if j1 (push (list (first rect) j1 
					   (third rect) j) boxes)
			  boxes)))
	   (if (blank-column-p image rect j tol)
	       (if j1 (boxec image rect (1+ j) nil j
			     (push (list (first rect) j1 
					 (third rect) j) boxes))
		 (if s1
		     (if (> (- j s1) space)
			 (boxec image rect (1+ j) nil j
				(push (list (first rect) s1 
					    (third rect) j) boxes))
			 (boxec image rect (1+ j) nil s1 boxes))
		   (boxec image rect (1+ j) nil j boxes)))
	     (if j1 (boxec image rect (1+ j) j1 nil boxes)
	       (boxec image rect (1+ j) j nil boxes))))
       (boxec image rect (second rect) j1 s1 boxes))))
   (boxec image (if rect rect (append '(0 0) (array-dimensions image)))
	  nil nil nil nil)))

(defun copy-rect (rect image)
  (let* ((height (- (third rect) (first rect))) 
	 (width (- (fourth rect) (second rect)))
	 (tmp (make-array `(,height ,width))))
    (do ((i (first rect) (1+ i)))
	((= i (third rect)) tmp)
	(do ((j (second rect) (1+ j))
	     (is (- i (first rect))))
	    ((= j (fourth rect)))
	    (setf (aref tmp is (- j (second rect))) (aref image i j))))))

(defmacro if-array (array dim var condition 1st 2nd)
  (let ((res (gensym)))
    `(do ((,var 0 (1+ ,var))
	  (,res nil ,condition))
	 ((or ,res (= ,var (array-dimension ,array ,dim)))
	  (setf ,var (1- ,var))
	  (if ,res ,1st ,2nd)))))

(defun myxor (image mask)
  (dotimes (i (array-dimension image 0))
    (dotimes (j (array-dimension image 1))
      (if (= (aref image i j) (aref mask i j))
	  (setf (aref image i j) 0)))))

(defun xtrim (image mask)
  ;(write-plain-pbm "i1.pnm" image)
  ;(write-plain-pbm "i2.pnm" mask)
  ;(bit-xor image mask t)
  (myxor image mask)
  ;(write-plain-pbm "i3.pnm" image)
  (let ((rect (append '(0 0) (array-dimensions image))))
    (if-array 
     image 1 j (not (blank-column-p image rect j 1))
     (progn
       (setf (second rect) j)
       (copy-rect (trim-box rect image) image))
     ())))

(defun look-around (i j image tmp)
  (remove-if
   #'(lambda (point)
       (let ((y (first point)) (x (second point)))
	 (or (not (array-in-bounds-p image y x))
	     (= 0 (aref image y x)) (= 1 (aref tmp y x)))))
  `((,(1- i) ,(1- j)) (,(1- i) ,j) (,(1- i) ,(1+ j))
    (,i ,(1- j)) (,i ,(1+ j)) 
    (,(1+ i) ,(1- j)) (,(1+ i) ,j) (,(1+ i) ,(1+ j)))))

(defun flood-fill (i j rect image tmp)
  (if (< i (nth 0 rect)) (setf (nth 0 rect) i))
  (if (> i (nth 2 rect)) (setf (nth 2 rect) i))
  (if (< j (nth 1 rect)) (setf (nth 1 rect) j))
  (if (> j (nth 3 rect)) (setf (nth 3 rect) j))
  (setf (aref tmp i j) 1)
  (setf (aref image i j) 0)
  (dolist (point (look-around i j image tmp))
    (let ((y (first point)) (x (second point))) 
      (flood-fill y x rect image tmp))))

(defun carve-into-boxes (image)
  (labels
   ((carvec
     (img boxes)
;     (format t "carvec ")
     (let ((tmp (make-array (array-dimensions img) ; :element-type 'bit
			    :initial-element 0))
	   (rct (append '(0 0) (array-dimensions img))))
       (if-array 
	img 0 i (= 1 (aref img i 0))
	(progn
;	  (write-plain-pbm "a1.pnm" img)
	  (flood-fill i 0 rct img tmp)
;	  (write-plain-pbm "a11.pnm" img)
;	  (write-plain-pbm "a12.pnm" tmp)
	  (if-array 
	   tmp 1 j (blank-column-p tmp rct j 1)
	   (progn
	     (setf (fourth rct) j)
;	     (format t "r: ~A~%" rct)
;	     (write-plain-pbm "a2.pnm" img)
;	     (format t "t: ~A~%" (trim-box rct tmp))
	     (let ((i4 (copy-rect (trim-box rct tmp) tmp)))
;	       (write-plain-pbm "aa.pnm" img)
;	       (write-plain-pbm "a.pnm" i4)
;	       (write-plain-pbm "aaa.pnm" tmp)
	       (push i4 boxes)
;	       (read)
	       )
	     (let ((nim (xtrim img tmp)))
	       (if nim (carvec nim boxes) (reverse boxes))))
	   (reverse (push tmp boxes))))
	(if (< 1 (array-dimension img 1))
	    (progn
	      (setf (second rct) 1)
	      (carvec (copy-rect rct img) boxes))
	  (reverse boxes))))))
   (let ((i5 (copy-image image)))
;     (write-plain-pbm "a0.pnm" i5)
     (carvec i5 nil))))

(defun boxify (image &key (tol 1) preserve-spaces
	       (rect (append '(0 0) (array-dimensions image))))
  (let ((boxes (mapcar #'(lambda (r) (copy-rect r image)) (cut-into-boxes rect image tol))))
    (if preserve-spaces
	boxes
	(remove-if 
	 #'(lambda (a) 
	     (or (zerop (array-dimension a 0)) (zerop (array-dimension a 1))))
	 boxes))))

(defun cut-into-words (image &key (tol 1)
		       (rect (append '(0 0) (array-dimensions image))))
  (labels
   ((space-p (box) (or (zerop (- (third box) (first box))) (zerop (- (fourth box) (second box)))))
    (wordec (boxes word words)
      (if (null boxes) (nreverse (if word (push word words) words))
	  (if (space-p (car boxes))
	      (wordec (cdr boxes) nil (if word (push word words) words))
	      (wordec (cdr boxes) (push (car boxes) word) words)))))
    (mapcar
     #'(lambda (w)
	 (loop for b in w
	      minimize (first b) into mini
	      maximize (third b) into maxi
	      minimize (second b) into minj
	      maximize (fourth b) into maxj
	      finally (return (list mini minj maxi maxj))))
    (wordec (cut-into-boxes rect image tol) nil nil))))
  
(defun observe-word (freqs tree image)
  (apply 
   #'shortest-chain freqs
   (apply
    #'append
    (mapcar
     #'(lambda (i0)
	 (let ((c0 (char-chains-p tree i0)))
;	   (format t "~%c0 ~a~%" c0)
	   (if (and c0 (zerop (second (caaar c0)))) (list c0)
	       (mapcar
		#'(lambda (i1)
		    (let ((c1 (char-chains-p tree i1)))
;		      (format t "~%c1 ~a~%" c1)
		      (if (and c1 (zerop (second (caaar c1)))) c1
			  (guess-char-chains tree i1))))
		(carve-into-boxes i0)))))
     (boxify image)))))


(defparameter english-alphabet "abcdefghijklmnopqrstuvwxyz")
(defparameter latin-alphabet "abcdefghijlmnopqrstuvxy")
(defparameter welsh-alphabet "abcdefghijlmnoprstuwyâôŵŷî")
(defparameter russian-alphabet "абвгдежзийклмнопрстуфхцчшщьъэюя")
(defparameter ukrainian-alphabet "абвгґдеєжзиiїйклмнопрстуфхцчшщьюя")
(defparameter persian-alphabet "ﺎآاﺍﺐﺒﺑﺏﭗﭙپﺖﺘﺗﺕﺚﺜﺛﺙﺞﺠﺟﺝﭻﭽﭼﭺﺢﺤﺣﺡﺦﺨﺧﺥﺪﺩﺬﺫﺮﺭﺰﺯﮋژﺲﺴﺳﺱﺶﺸﺷﺵﺺﺼﺻﺹﺾﻀﺿﺽﻂﻄﻃﻁﻆﻈﻇﻅﻊﻌﻋﻉﻎﻐﻏﻍﻒﻔﻓﻑﻖﻘﻗﻕﮏﮑﮐکﮓﮕﮔگﻞﻠﻟﻝﻢﻤﻣﻡﻦﻨﻧﻥﻮوﻪﻬﻫﻩﯽﻴﻳﻯ")

(defparameter persian-corpus " ﻰﺳﺭﺎﻓ ")

(defun calculate-frequencies (charset coprus-stream)
  (labels 
      ((get-char (stream)
	 (let ((c (read-char stream nil nil)))
	   (if (null c) 'eof
	       (if (position c " *\"'.,?!:;()-/") #\Space
		   (if (position (char-downcase c) charset :test #'char=) 
		       (char-downcase c)))))))
    (do ((hash (make-hash-table :test #'equal))
	 (max-freq 0)
	 (c1 #\Space c2)
	 (c2 (get-char coprus-stream) (get-char coprus-stream)))
	((eq c2 'eof)
	 (maphash
	  #'(lambda (k v) 
	      (setf (gethash k hash) (/ (- max-freq v) max-freq 1.0)))
	  hash)
	 hash)
      (if (and c1 c2 (not (char= #\Space c1 c2)))
	  (setf 
	   max-freq 
	   (max max-freq 
		(incf (gethash (coerce (vector c1 c2) 'string) hash 0))))))))

(defun hash->alist (hash)
  (loop for k being the hash-keys of hash using (hash-value v) 
     collect (cons k v)))
    
(defun alist->hash (alist &key (test #'equal))
  (let ((hash (make-hash-table :test test)))
    (dolist (p alist hash) (setf (gethash (car p) hash) (cdr p)))))

(defun file->freqs (name) (alist->hash (with-open-file (f name) (read f))))

(defun freqs->file (name freqs)
  (with-open-file (f name :direction :output :if-exists :supersede :if-does-not-exist :create)
    (format f "~s" (sort (hash->alist freqs) #'< :key #'cdr))))

(defun split-script (script charset0 charset1 
		     &key
		     (script0 (make-pathname :type "0" :defaults script))
		     (script1 (make-pathname :type "1" :defaults script)))
  (let (lst0 lst1)
    (dolist (c (with-open-file (f script) (read f)))
      (let ((p (position (car c) charset0)))
	(if (null p) (push c lst0)
	    (push (list (char charset1 p) (cadr c)) lst1))))
    (with-open-file (f script0 :direction :output :if-exists :supersede :if-does-not-exist :create)
      (format f "~s" lst0))
    (with-open-file (f script1 :direction :output :if-exists :supersede :if-does-not-exist :create)
      (format f "~s" lst1))))

(defparameter latin->cyrillic "ABCEHIJKMOPTXabceijmnouy")
(defparameter cyrillic->latin "АВСЕНIJКМОРТХаЬсеijтпоиУ")

;;
;; СДЕЛАТЬ МАКРОС MAP-IMAGE! или объединить с ch-image
;;

(defun bilinear (f x y x1 y1 x2 y2)
  (let ((dy21 (- y2 y1)) (dx21 (- x2 x1))
	(dy2 (- y2 y)) (dy1 (- y y1)) (dx2 (- x2 x)) (dx1 (- x x1)))
;    (format t "~a ~a ~a ~a ~a ~a ~a ~a ~a ~a " (funcall f x1 y1) (funcall f x2 y1) (funcall f x1 y2) (funcall f x2 y2) dx21 dy21 dx1 dx2 dy1 dy2)
    (cond
      ((and (= x x1 x2) (= y y1 y2)) (funcall f x1 y1))
      ((= x x1 x2) 
       (+ (* (/ (funcall f x1 y1) dy21) dy2) 
	  (* (/ (funcall f x1 y2) dy21) dy1)))
      ((= y y1 y2) 
       (+ (* (/ (funcall f x1 y1) dx21) dx2) 
	  (* (/ (funcall f x2 y1) dx21) dx1)))
      (t 
       (+ (* (/ (funcall f x1 y1) dx21 dy21) dx2 dy2)
	  (* (/ (funcall f x2 y1) dx21 dy21) dx1 dy2)
	  (* (/ (funcall f x1 y2) dx21 dy21) dx2 dy1)
	  (* (/ (funcall f x2 y2) dx21 dy21) dx1 dy1)))))) 

(defun poplar (image n)
  "return polar transform of cartesian IMAGE with N divisions along theta axis"
  (let ((rmax (ceiling (apply #'max (array-dimensions image)) 2))
	(i0 (floor (array-dimension image 0) 2))
	(j0 (floor (array-dimension image 1) 2)))
    (let ((p (make-array (list rmax n))))
      (dotimes (k n p)
	(let ((a (/ (* 2 pi k) n)))
	  (dotimes (r rmax)
	    (let ((x (+ j0 (* r (cos a)))) (y (+ i0 (* r (sin a)))))
	      (setf 
	       (aref p r k)
	       (round
		(bilinear
		 #'(lambda (i j) (cref image i j)) y x
		 (floor y) (floor x) (ceiling y) (ceiling x)))))))))))

; реализовать http://www.vuse.vanderbilt.edu/~rap2/papers/oncomplp.pdf?
(defun logpolar (image)
  "return the log-polar transform of cartesian IMAGE"
  (let* ((lmage (make-array (array-dimensions image)))
	 (h (array-dimension image 0))
	 (w (array-dimension image 1))
	 (dt (/ (* 2 pi) w))
	 (b (exp (/ (log h) h))))
    (dotimes (i h lmage)
      (dotimes (j w)
	(let* ((r (expt b i)) (th (* j dt))
	       (x (round (+ (* r (cos th)) (/ w 2))))
	       (y (round (+ (* r (sin th)) (/ h 2)))))
	  (if (and (< 0 x w) (< 0 y h)) 
	      (setf (aref lmage i j) (aref image y x))))))))

(defun poplar-log (m n image)
  (let* ((i0 (/ (array-dimension image 0) 2))
	 (j0 (/ (array-dimension image 1) 2))
	 (rmax (sqrt (+ (expt i0 2) (expt j0 2))))
	 (mm (/ m (log (1+ rmax)))) 
	 (lmage (make-array (list m n) :initial-element nil)))
    (dotimes (i (array-dimension image 0))
      (dotimes (j (array-dimension image 1))
	(let* ((di (- i j0)) (dj (- j j0))
	       (r (sqrt (+ (expt di 2) (expt dj 2))))
	       (rh (min (1- m) (max 0 (floor (+ 0.5 (* mm (log (1+ r))))))))
	       (th (min (1- n) (max 0 (floor (+ 0.5 (/ (* n (atan di dj)) 
						       ;(* 2 pi))))))))
						       pi)))))))
	  (push (aref image i j) (aref lmage rh th)))))
    (dotimes (i m lmage)
      (dotimes (j n)
	(setf (aref lmage i j) (round (apply #'average (aref lmage i j))))))))

(defun image-mean 
    (image &optional (i0 0) (j0 0) 
     (i1 (array-dimension image 0)) (j1 (array-dimension image 1)))
  (/ (loop for i from i0 below i1 sum 
	  (loop for j from j0 below j1 sum (aref image i j)))
     (array-total-size image)))

(defun image-deviation
    (image &optional mean (i0 0) (j0 0) 
     (i1 (array-dimension image 0)) (j1 (array-dimension image 1)))
  (if (null mean) (setf mean (image-mean image i0 j0 i1 j1)))
  (sqrt
   (/ (loop for i from i0 below i1 sum 
	   (loop for j from j0 below j1 sum (expt (- (aref image i j) mean) 2)))
      (array-total-size image))))

(defun cross-correlation 
    (mask image ic jc 
     &optional 
     (mask-av (image-mean mask))
     (mask-dv (image-deviation image mask-av)))
  (let ((i0 (- ic (floor (array-dimension mask 0) 2)))
	(i1 (+ ic (floor (array-dimension mask 0) 2)))
	(j0 (- jc (floor (array-dimension mask 1) 2)))
	(j1 (+ jc (floor (array-dimension mask 1) 2))))
    (let* ((image-av (image-mean image i0 j0 i1 j1))
	   (image-dv (image-deviation image image-av i0 j0 i1 j1)))
      (/
       (loop for i from i0 below i1 sum
	    (loop for j from j0 below j1 sum
		 (* (- (aref mask (- i i0) (- j j0)) mask-av)
		    (- (aref image i j) image-av))))
       (if (zerop mask-dv) 1 mask-dv) 
       (if (zerop image-dv) 1 image-dv)))))

(defun image-cross-correlation (mask image)
  (let* ((mask-av (image-mean mask))
	 (mask-dv (image-deviation mask mask-av))
	 (lmage (make-array (array-dimensions image)))
	 (i0 (+ 0 (floor (array-dimension mask 0) 2)))
	 (i1 (1+ (- (array-dimension image 0) (floor (array-dimension mask 0) 2))))
	 (j0 (+ 0 (floor (array-dimension mask 1) 2)))
	 (j1 (1+ (- (array-dimension image 1) (floor (array-dimension mask 1) 2)))))
    (loop for i from i0 below i1 do
	 (loop for j from j0 below j1 do
	      (setf (aref lmage i j) 
		    (cross-correlation mask image i j mask-av mask-dv))))
    lmage))

(defun slow-fourier-transform-2d (image)
  (let ((lmage (make-array (array-dimensions image)))
	(n (array-dimension image 0)) (n2 (apply #'* (array-dimensions image))))
    (dotimes (i n)
      (dotimes (j n)
	(setf
	 (aref lmage i j)
	 (/ (loop for k below n sum
		 (loop for l below n sum 
		      (* (aref image k l)
			 (exp (* #C(0 -1) 2 pi (/ (+ (* i k) (* j l)) n))))))
	    n2))))
    lmage))

(defun phase-difference (ft1 ft2)
  (declare (optimize (speed 3) (safety 0)))
  (let ((lmage (make-array (array-dimensions ft1)))
	(n0 (array-dimension ft1 0)) (n1 (array-dimension ft1 1)))
    (dotimes (i n0)
      (dotimes (j n1)
	(let ((p (* (aref ft1 i j) (conjugate (aref ft2 i j)))))
	  (setf (aref lmage i j) (if (zerop (abs p)) 0 (/ p (abs p)))))))
    lmage))

(defun phase-difference! (ft1 ft2)
  "maps into ft2"
  (declare (optimize (speed 3) (safety 0)))
  (let ((n0 (array-dimension ft1 0)) (n1 (array-dimension ft1 1)))
    (dotimes (i n0)
      (dotimes (j n1)
	(let ((p (* (aref ft1 i j) (conjugate (aref ft2 i j)))))
;	  (format t "~a ~a ~a~%" (aref ft1 i j) (aref ft2 i j)
		  (setf (aref ft2 i j) (if (zerop (abs p)) 0 (/ p (abs p))))
;		  )
	  )))
    ft2))

(defun slow-inverse-fourier-transform-2d (image)
  (let ((lmage (make-array (array-dimensions image)))
	(n (array-dimension image 0)) (n2 (apply #'* (array-dimensions image))))
    (dotimes (i n)
      (dotimes (j n)
	(setf
	 (aref lmage i j)
	 (/ (loop for k below n sum
		 (loop for l below n sum 
		      (* (aref image k l)
			 (exp (* #C(0 1) 2 pi (/ (+ (* i k) (* j l)) n))))))
	    n2))))
    lmage))

(defun fourier-magnitudes (image)
  (let ((lmage (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) (abs (aref image i j)))))))

(defun fourier-magnitudes! (image)
  (dotimes (i (array-dimension image 0) image)
    (dotimes (j (array-dimension image 1))
      (setf (aref image i j) (abs (aref image i j))))))
    
(defun image-max (image)
  (let ((y 0) (x 0) (m (aref image 0 0)))
    (dotimes (i (array-dimension image 0) (values y x m))
      (dotimes (j (array-dimension image 1))
	(if (> (aref image i j) m) (setf m (aref image i j) y i x j))))))

(defun rot-max (image)
  (let ((x 0) (m (aref image 0 0)))
    (dotimes (j (array-dimension image 1) (values 0 x m))
      (if (> (aref image 0 j) m) (setf m (aref image 0 j) x j)))))

(defun fourier-correlation (mask image)
  (fourier-magnitudes 
   (inverse-fft-2d! 
    (phase-difference 
     (fft-2d! (duplicate-image image))    ; почему тут наооборот?!
     (fft-2d! (duplicate-image mask)))))) ; внизу правильно!

(declaim (notinline fourier-correlation!))
(defun fourier-correlation! (f1 f2)
  "maps into ft1"
  (fourier-magnitudes! (inverse-fft-2d! (phase-difference! f2 f1))))

(declaim (notinline fft-n-mags!))
(defun fft-n-mags! (image)
  (let ((f (real-fft-2d! image)))
    (values (fourier-magnitudes f) f)))

(defun for-mellin (image)
;  (format t "4mellin~%")
  (multiple-value-bind (m f) (fft-n-mags! image)
;    (format t "info-mellin~%")
    (values
     (logpolar 
      (grayscale-image!
       (high-pass! 
	(fourier->image m)
	)
       )) 
     f)))

(defun logarithm-filter (image &optional (scale 0.0001))
  (let ((lmage (make-array (array-dimensions image)))
	(c (/ 255 (log (1+ (* scale (nth-value 2 (image-max image))))))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) 
	      (round (* c (log (1+ (* scale (aref image i j)))))))))))

(defun fourier->image (image)
  (let* ((lmage (make-array (array-dimensions image)))
;	 (m image) ;(fourier-magnitudes image))
;	 (c (/ 255 (max 1 (log (1+ (* 0.0001 (nth-value 2 (image-max m))))))))
	 (n (array-dimension image 0)))
    (dotimes (i n lmage)
      (dotimes (j n)
	(setf
	 (aref lmage 
	       (if (> (/ n 2) i) (+ (/ n 2) i) (- i (/ n 2)))
	       (if (> (/ n 2) j) (+ (/ n 2) j) (- j (/ n 2))))
;	 (round (* c (log (1+ (* 0.0001 (aref m i j)))))))))))
	 (aref image i j))))))

(defun fourier-for-rotation (image)
  (let ((lmage (make-array (array-dimensions image)))
	(n (array-dimension image 0)))
    (dotimes (i n lmage)
      (dotimes (j n)
	(setf
	 (aref lmage 
	       (if (> (/ n 2) i) (+ (/ n 2) i) (- i (/ n 2)))
	       (if (> (/ n 2) j) (+ (/ n 2) j) (- j (/ n 2))))
	 (aref image i j))))))

(defun grayscale-image (image)
  (let ((lmage (make-array (array-dimensions image)))
	(m (nth-value 2 (image-max image))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1)) 
	(setf (aref lmage i j) 
	      (floor (* 255 (/ (1+ (aref image i j)) (1+ m)))))))))

(defun grayscale-image! (image)
  (let ((m (nth-value 2 (image-max image))))
    (dotimes (i (array-dimension image 0) image)
      (dotimes (j (array-dimension image 1)) 
	(setf (aref image i j) 
	      (floor (* 255 (/ (1+ (aref image i j)) (1+ m)))))))))

(defun high-pass (image &optional (d (/ (array-dimension image 0) 8)))
  (let ((lmage (make-array (array-dimensions image)))
	(n/2 (/ (array-dimension image 0) 2)))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) 
	      (if (> d (sqrt (+ (expt (- i n/2) 2) (expt (- j n/2) 2)))) 0 
		  (aref image i j)))))))

(defun high-pass! (image &optional (d (/ (array-dimension image 0) 8)))
  (let ((n/2 (/ (array-dimension image 0) 2)))
    (dotimes (i (array-dimension image 0) image)
      (dotimes (j (array-dimension image 1))
	(setf (aref image i j) 
	      (if (> d (sqrt (+ (expt (- i n/2) 2) (expt (- j n/2) 2)))) 0 
		  (aref image i j)))))))

(defun low-pass (image &optional (d (/ (array-dimension image 0) 2)))
  (let ((lmage (make-array (array-dimensions image)))
	(n/2 (/ (array-dimension image 0) 2)))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) 
	      (if (<= d (sqrt (+ (expt (- i n/2) 2) (expt (- j n/2) 2)))) 0 
		  (aref image i j)))))))

(defun fourier-mellin-correlation (mask image)
  (multiple-value-bind (i1 f1) (for-mellin (duplicate-image mask))
    (multiple-value-bind (i2 f2) (for-mellin (duplicate-image image))
      (values 
       (fourier-correlation! 
	;(treshold-fourier! 0.01 i1)
	(real-fft-2d! i1) 
	(real-fft-2d! i2)) 
       f1 f2))))

(defun pre-fft
       (n
        &key
        (complex-exp-vector nil))
  (declare (optimize (speed 3) (safety 0)))
  (declare (type fixnum n))
  (let* ((the-vector (if complex-exp-vector complex-exp-vector
			 (make-array n)))
	 (s (/ (* 2 pi) n))
	 (c1 (complex (cos s) (- (sin s))))
	 (c2 (complex 1.0d0 0.0d0)))
    (declare (type (simple-array complex (*)) the-vector))
    (dotimes (i n the-vector)
      (setf (aref the-vector i) c2
	    c2 (* c2 c1)))))

(defun fft! (complex-vector)
  (let* ((n (length complex-vector))
	 (exponent (round (log n 2)))
         (W (pre-fft n))
	 (MM 1)
	 (LL n)
	 NN JJ
	 c1 c2)
    (declare (optimize (speed 3) (safety 0)))
    (dotimes (k exponent)
      (setf NN (/ LL 2)
            JJ MM)
      (do* ((i 0 (+ i LL))
            (kk NN (+ i NN)))
           ((>= i N))
        (setf c1 (+ (aref complex-vector i) (aref complex-vector kk))
              (aref complex-vector kk)
              (- (aref complex-vector i) (aref complex-vector kk))
              (aref complex-vector i) c1))
      (cond
       ((> NN 1)
        (do ((j 1 (1+ j)))
            ((>= j NN))
          (setf c2 (svref W JJ))
          (do* ((i j (+ i LL))
                (kk (+ j NN) (+ i NN)))
               ((>= i N))
            (setf c1 (+ (aref complex-vector i) (aref complex-vector kk))
                  (aref complex-vector kk)
                  (* (- (aref complex-vector i) (aref complex-vector kk))
                     c2)
                  (aref complex-vector i) c1))
          (incf jj MM))
        (setf LL NN)
        (setf MM (* MM 2)))))
    (let ((j 0)
          (nv2 (/ n 2))
          k)
      (dotimes (i (1- N))
        (if (< i j)
          (setf c1 (aref complex-vector j)
                (aref complex-vector j) (aref complex-vector i)
                (aref complex-vector i) c1))
        (setf k nv2)
        (loop
          (if (> k j) (return))
          (decf j k)
          (setf k (/ k 2)))
        (incf j k)))
    complex-vector))

(defun inverse-fft! (complex-vector)
  (let ((n (length complex-vector)))
    (dotimes (i N)
      (setf (aref complex-vector i)
	    (conjugate (aref complex-vector i))))
    (fft! complex-vector)
    (dotimes (i N complex-vector)
      (setf (aref complex-vector i)
	    (/ (conjugate (aref complex-vector i)) N)))))

(defun real-fft (real-vector)
  (declare (optimize (speed 3) (safety 0)))
  (let* ((n (length real-vector))
	 (n2 (/ n 2))
	 (tmp1 (make-array n2))
	 (tmp2 (make-array n)))
    (dotimes (i n2)
      (let ((j (* i 2)))
	(setf (aref tmp1 i) (complex (aref real-vector j) 
				     (aref real-vector (1+ j))))))
    (fft! tmp1)
    (setf (aref tmp2 0) (+ (realpart (aref tmp1 0)) (imagpart (aref tmp1 0)))
	  (aref tmp2 n2) (- (realpart (aref tmp1 0)) (imagpart (aref tmp1 0))))
    (loop for k from 1 to (/ n2 2) do
	 (let* ((fpk (aref tmp1 k)) (n2k (- n2 k))
		(fpnk (conjugate (aref tmp1 n2k)))
		(f1k (+ fpk fpnk)) (f2k (- fpk fpnk))
		(a (* (- pi) (+ 0.5 (/ k n2))))
		(tw (* f2k (complex (cos a) (sin a)))))
	   (setf (aref tmp2 k) (/ (+ f1k tw) 2) 
		 (aref tmp2 n2k) (conjugate (/ (- f1k tw) 2)))))
    (dotimes (k n2 tmp2) (setf (aref tmp2 (+ n2 k)) 
			       (conjugate (aref tmp2 (- n2 k)))))))

(declaim (notinline fft-2d!))
(defun fft-2d! (image)
  (let* ((n0 (array-dimension image 0)) (n1 (array-dimension image 1))
	 (tmp0 (make-array n1)) (tmp1 (make-array n0)))
    (dotimes (i n0)
      (dotimes (j n1) (setf (aref tmp0 j) (aref image i j)))
      (fft! tmp0)
      (dotimes (j n1) (setf (aref image i j) (aref tmp0 j))))
    (dotimes (j n1)
      (dotimes (i n0) (setf (aref tmp1 i) (aref image i j)))
      (fft! tmp1)
      (dotimes (i n0) (setf (aref image i j) (aref tmp1 i))))
    image))

(declaim (notinline real-fft-2d!))
(defun real-fft-2d! (image)
  (let* ((n0 (array-dimension image 0)) (n1 (array-dimension image 1))
	 (tmp0 (make-array n1)) (tmp1 (make-array n0)))
    (dotimes (i n0)
      (dotimes (j n1) (setf (aref tmp0 j) (aref image i j)))
      (setf tmp0 (real-fft tmp0))
      (dotimes (j n1) (setf (aref image i j) (aref tmp0 j))))
    (dotimes (j n1)
      (dotimes (i n0) (setf (aref tmp1 i) (aref image i j)))
      (fft! tmp1)
      (dotimes (i n0) (setf (aref image i j) (aref tmp1 i))))
    image))

(declaim (notinline inverse-fft-2d!))
(defun inverse-fft-2d! (image)
  (let* ((n0 (array-dimension image 0)) (n1 (array-dimension image 1))
	 (tmp0 (make-array n1)) (tmp1 (make-array n0)))
    (dotimes (i n0)
      (dotimes (j n1) (setf (aref tmp0 j) (aref image i j)))
      (inverse-fft! tmp0)
      (dotimes (j n1) (setf (aref image i j) (aref tmp0 j))))
    (dotimes (j n1)
      (dotimes (i n0) (setf (aref tmp1 i) (aref image i j)))
      (inverse-fft! tmp1)
      (dotimes (i n0) (setf (aref image i j) (aref tmp1 i))))
    image))
    
(defun %rotate-image (image a)
  (declare (optimize (speed 3) (safety 0)))
  (let* ((lmage (make-array (array-dimensions image)))
	    (h (array-dimension image 0))
	    (w (array-dimension image 1))
	    (c (cos (- a))) (s (sin (- a)))
	    (i0 (/ h 2)) (j0 (/ w 2)))
       (dotimes (i2 h lmage)
	 (dotimes (j2 w)
	   (let ((i1 (+ i0 (* s (- j2 j0)) (* c (- i2 i0))))
		 (j1 (+ j0 (* c (- j2 j0)) (* (- s) (- i2 i0)))))
	     (setf 
	      (aref lmage i2 j2)
	      (round
	       (bilinear
		#'(lambda (i j) (cref image i j)) i1 j1
		(floor i1) (floor j1) (ceiling i1) (ceiling j1))))
;	     (format t " ~a ~a ~a~%" i2 j2 (aref lmage i2 j2))
)))))

(declaim (notinline rotate-image))
(defun rotate-image (image a)
  (declare (optimize (speed 3) (safety 0)))
  (cond
    ((> 0.008 (abs a)) (duplicate-image image))
    ((> 0.008 (abs (- a pi)))
     (let ((lmage (make-array (array-dimensions image)))
	   (h (array-dimension image 0))
	   (w (array-dimension image 1)))
       (dotimes (i h lmage)
	 (dotimes (j w)
	   (setf (aref lmage i j) (aref image (- h i 1) (- w j 1)))))))
    (t (%rotate-image image a))))
     
(defun rotate-fourier (image a)
  (let ((lmage (make-array (array-dimensions image)))
	(1mage (fourier-for-rotation image))
	(c (cos (- a))) (s (sin (- a)))
	(i0 (/ (array-dimension image 0) 2))
	(j0 (/ (array-dimension image 1) 2)))
    (dotimes (i2 (array-dimension image 0) (fourier-for-rotation lmage))
      (dotimes (j2 (array-dimension image 1))
	(let ((i1 (+ i0 (* s (- j2 j0)) (* c (- i2 i0))))
	      (j1 (+ j0 (* c (- j2 j0)) (* (- s) (- i2 i0)))))
	  (setf 
	   (aref lmage i2 j2)
	    (bilinear
	     #'(lambda (i j) (cref 1mage i j)) i1 j1
	     (floor i1) (floor j1) (ceiling i1) (ceiling j1))))))))

(defun wolberg-zokai-fft-registration (mask image)
  (let ((lpm (fft-2d! (logpolar (low-pass mask))))
	(lmage (make-array (array-dimensions image)))
	(1mage (make-array (array-dimensions image) :initial-element nil))
	(i0 (+ 0 (floor (array-dimension mask 0) 2)))
	(i1 (1+ (- (array-dimension image 0) (floor (array-dimension mask 0) 2))))
	(j0 (+ 0 (floor (array-dimension mask 1) 2)))
	(j1 (1+ (- (array-dimension image 1) (floor (array-dimension mask 1) 2)))))
    (loop for i from i0 below i1 do
	 (loop for j from j0 below j1 do
	      (multiple-value-bind (y x m)
		  (image-max
		   (fourier-magnitudes 
		    (inverse-fft-2d! 
		     (phase-difference
		      (fft-2d!
		       (logpolar
			(low-pass
			 (copy-rect 
			  (list (- i (floor (array-dimension mask 0) 2))
				(- j (floor (array-dimension mask 1) 2))
				(+ i (floor (array-dimension mask 0) 2))
				(+ j (floor (array-dimension mask 1) 2)))
			  image))))
		      lpm))))
		(setf (aref lmage i j) m (aref 1mage i j) (list y x)))))
    (values lmage 1mage)))
    
(defun wolberg-zokai-registration (mask image)
  (let ((lpm (double-image (logpolar (low-pass mask))))
	(lmage (make-array (array-dimensions image)))
	(1mage (make-array (array-dimensions image) :initial-element nil))
	(i0 (+ 0 (floor (array-dimension mask 0) 2)))
	(i1 (1+ (- (array-dimension image 0) (floor (array-dimension mask 0) 2))))
	(j0 (+ 0 (floor (array-dimension mask 1) 2)))
	(j1 (1+ (- (array-dimension image 1) (floor (array-dimension mask 1) 2)))))
    (loop for i from i0 below i1 do
	 (loop for j from j0 below j1 do
	      (multiple-value-bind (y x m)
		  (image-max
		   (image-cross-correlation
		    (logpolar
		     (low-pass
		      (copy-rect 
		       (list (- i (floor (array-dimension mask 0) 2))
			     (- j (floor (array-dimension mask 1) 2))
			     (+ i (floor (array-dimension mask 0) 2))
			     (+ j (floor (array-dimension mask 1) 2)))
		       image)))
		    lpm))
		(setf (aref lmage i j) m (aref 1mage i j) (list y x)))))
    (values lmage 1mage)))

(defun double-image (image)
  (let ((lmage (make-array (list (* (array-dimension image 0) 2)
				 (* (array-dimension image 1) 2)))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(let ((i1 (+ i (/ (array-dimension image 0) 2)))
	      (j1 (+ j (/ (array-dimension image 1) 2)))
	      (i2 (if (< i (/ (array-dimension image 0) 2))
		      (+ i (* 3/2 (array-dimension image 0)))
		      (- i (/ (array-dimension image 0) 2))))
	      (j2 (if (< j (/ (array-dimension image 1) 2))
		      (+ j (* 3/2 (array-dimension image 1)))
		      (- j (/ (array-dimension image 1) 2)))))
	  (setf (aref lmage i1 j1) (aref image i j)
		(aref lmage i1 j2) (aref image i j)
		(aref lmage i2 j1) (aref image i j)
		(aref lmage i2 j2) (aref image i j)))))))

(defun pad-image (p image &optional (p1 p))
  (if (or (> (array-dimension image 0) p)
	  (> (array-dimension image 1) p1))
      (duplicate-image image)
      (let ((lmage (make-array (list p p1)))
	    (io (floor (- p (array-dimension image 0)) 2))
	    (jo (floor (- p1 (array-dimension image 1)) 2)))
	(dotimes (i (array-dimension image 0) lmage)
	  (dotimes (j (array-dimension image 1))
	    (setf (aref lmage (+ i io) (+ j jo)) (aref image i j)))))))

(defun get-angles (t0 t1 r)
  (let* ((mm (nth-value 2 (image-max r)))
	 (th0 (* mm t0)) (th1 (* mm t1))
         (h (array-dimension r 0)) (w (array-dimension r 1)) q
         (b (exp (/ (log h) h))))
    (dotimes (i h (group #'= q :key #'car))
      (let ((m (loop for j below w maximize (aref r i j))))
	(if (> m th0)
	    (dotimes (j w)
	      (if (> (aref r i j) th1)
		  (push (list (expt b (- i (/ h 2))) 
			      (* 2 pi (/ (- j (/ w 2)) w)))
			q))))))))
			  

(defun fourier-mellin-search (mask image &key (t0 0.3) (t1 0.5) (t2 0.5) (t3 0.2) (t4 1) debug-p dry)
  (let* ((n (power-of-2-ceiling (apply #'max (array-dimensions image))))
	 (io (floor (- n (array-dimension image 0)) 2)) (q ()) (ql 0) (qs 1)
	 (jo (floor (- n (array-dimension image 1)) 2)))
    (multiple-value-bind (r f1 f2)
	(fourier-mellin-correlation (pad-image n mask) (pad-image n image))
      (declare (ignore f1))
      ;(write-netpbm "~/Desktop/r.pnm" (grayscale-image! (fourier->image r)))
      (let ((gg (get-angles t2 t3 (fourier->image r))))
      (dolist (g gg
	       (progn
		 (if debug-p (format t "ql ~a ~a~%" ql qs))
		 (if dry qs
		     (remove-if 
		      #'(lambda (k) (/= qs (sixth k)))
		      (remove-if #'(lambda (k) (> (max 16 (* t1 ql)) k)) q 
				 :key #'third)))))
	(if (>= t4 (caar g) 1) ; пока только уменьшение... и то не работат :(
	    (let ((b (pad-image n (shrink-image (caar g) mask))))
	      (dolist (a (mapcar #'cadr g))
		(let* ((p (rotate-image b a))
		       (z (fourier-correlation! (real-fft-2d! p) f2))
		       (im (nth-value 2 (image-max z)))
		       (imm (image-mean z))
		       (m/m (round im imm))
		       (th1 (* t0 im)))
		  (if debug-p (format t "angle ~a ~a ~a~%" 
				      (caar g) (* 180 (/ a pi)) m/m))
		  (if (> m/m ql) (setf ql m/m qs (caar g)))
		  (if (not dry)
		      (dotimes (k n)
			(dotimes (l n)
			  (when (< th1 (aref z k l))
			    (let ((y (mod (+ k (/ n 2)) n)) 
				  (x (mod (+ l (/ n 2)) n)))
			      (push (list (- y io) (- x jo) m/m 
					  (/ (* a 180) pi)
					  (aref z k l) (caar g))
				    q)))))))))))))))

(defun fourier-search (mask image &key (t0 0.1) (t1 0.9) 
		       mask-cooked-p image-cooked-p)
  (let* ((n0 (power-of-2-ceiling (array-dimension image 0)))
	 (n1 (power-of-2-ceiling (array-dimension image 1)))
	 (io (floor (- n0 (array-dimension image 0)) 2)) (q) 
	 (jo (floor (- n1 (array-dimension image 1)) 2))
	 (m1 (if mask-cooked-p (duplicate-image mask)
		 (treshold-fourier! t0 (pad-image n0 mask n1))))
	 (f1 (if image-cooked-p (duplicate-image image-cooked-p)
		 (real-fft-2d! (pad-image n0 image n1))))
	 (z (fourier-correlation! m1 f1))
	 (th1 (* t1 (nth-value 2 (image-max z)))))
    (dotimes (k n0 q)
      (dotimes (l n1)
	(when (< th1 (aref z k l))
;	  (format t "~a ~a ~a~%" k l (aref z k l))
	  (let ((y (mod (+ k (/ n0 2)) n0)) 
		(x (mod (+ l (/ n1 2)) n1)))
	    (push (list (- y io) (- x jo) (aref z k l)) q)))))))

    
(defun hough-line-transform (image)
  "computes Hough transform for lines"
  (let* ((maxr (ceiling (sqrt (+ (expt (array-dimension image 0) 2) (expt (array-dimension image 1) 2)))))
	 (r0 (/ maxr 2))
	 (i0 (/ (array-dimension image 0) 2))
	 (j0 (/ (array-dimension image 1) 2))
	 (lmage (make-array (list maxr 360))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(if (not (zerop (aref image i j)))
	    (dotimes (a 360)
	      (let* ((aa (* (/ a 360) 2 pi)) 
		     (r (+ (* (- i i0) (sin aa)) (* (- j  j0) (cos aa)))))
		;(format t "~a ~a~%" a r)
		(incf (aref lmage (floor (+ r r0)) a) (aref image i j)))))))))

		
(defun sobel (image)
  (let ((gx (convolve sobel-gx image))
	(gy (convolve sobel-gy image))
	(g (make-array (array-dimensions image)))
	(q (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) (values g q))
      (dotimes (j (array-dimension image 1))
	(setf (aref g i j) 
	      (sqrt (+ (expt (aref gx i j) 2) (expt (aref gy i j) 2)))
	      (aref q i j) 
	      (atan (aref gy i j) 
		    (if (zerop (aref gx i j)) 1 (aref gx i j))))))))

(defun hough+sobel-line-transform (image)
  (multiple-value-bind (g q) (sobel image)
    (let* ((maxr (ceiling (sqrt (+ (expt (array-dimension image 0) 2) (expt (array-dimension image 1) 2)))))
	 (r0 (/ maxr 2))
	 (i0 (/ (array-dimension image 0) 2))
	 (j0 (/ (array-dimension image 1) 2))
	 (lmage (make-array (list maxr 360))))
      (dotimes (i (array-dimension image 0) lmage)
	(dotimes (j (array-dimension image 1))
	  (if (not (zerop (aref g i j)))
	      (let* ((u (aref q i j))
		     (a (floor (mod (+ (* (/ u pi) 180) 360) 360)))
		     (r (+ (* (- i i0) (sin u)) 
			   (* (- j j0) (cos u)))))
		(incf (aref lmage (floor (+ r r0)) a) (aref image i j)))))))))

(defun hough+sobel-circle-transform (r image)
  (multiple-value-bind (g q) (sobel image) ; порог на g
    (let ((lmage (make-array (array-dimensions image))))
      (dotimes (i (array-dimension image 0) lmage)
	(dotimes (j (array-dimension image 1))
	  (if (not (zerop (aref g i j)))
	      (let* ((c (cos (aref q i j))) 
		     (s (sin (aref q i j)))
		     (y (floor (- i (* r s)))) (x (floor (+ j (* r c)))))
		(if (array-in-bounds-p image y x)
		    (incf (aref lmage y x))))))))))

(defun hough+sobel-circles-transform (image)
  (multiple-value-bind (g q) (sobel image) ; порог на g
    (let ((lmage (make-array (array-dimensions image))))
      (dotimes (i (array-dimension image 0) lmage)
	(dotimes (j (array-dimension image 1))
	  (if (not (zerop (aref g i j)))
	      (let ((tg (tan (aref q i j))))
		(dotimes (x (array-dimension image 1))
		  (let ((y (floor (+ (- (* j tg) (* x tg)) i))))
		    (if (array-in-bounds-p image y x)
			(incf (aref lmage y x))))))))))))

(defun negate-grayscale-image (image)
  (let ((lmage (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) (abs (- (aref image i j) 255)))))))
	
(defun sobel-r-table 
    (image 
     &optional 
     (centre (mapcar #'(lambda (n) (/ n 2)) (array-dimensions image))))
  (multiple-value-bind (g q) (sobel image)
    (let ((lmage (make-array 360 :initial-element nil))
	  (i0 (car centre)) (j0 (cadr centre)))
      (dotimes (i (array-dimension image 0) lmage)
	(dotimes (j (array-dimension image 1))
	  (if (not (zerop (aref g i j)))
	      (let ((y (- i i0)) (x (- j j0)) (u (aref q i j)))
		(let ((a (floor (mod (+ (* (/ u pi) 180) 360) 360)))
		      (r (sqrt (+ (expt x 2) (expt y 2)))) ; может, вычитать
		      (th (if (zerop x) 0 (atan y x))))    ; градиент?
		  (push (list r th) (aref lmage a))))))))))

(defun generalized-hough+sobel-transform (r-table image)
  "computes the generalized Hough tranform, finding edges with Sobel operator"
  (multiple-value-bind (g q) (sobel image)
    (let ((lmage (make-array (array-dimensions image) :initial-element nil)))
      (dotimes (i (array-dimension image 0) lmage)
	(dotimes (j (array-dimension image 1))
	  (if (not (zerop (aref g i j)))
	      (let ((a (floor (mod (+ (* (/ (aref q i j) pi) 180) 360) 360))))
		(dotimes (u 360)
		  (dolist (p (aref r-table (mod (+ 360 (- a u)) 360)))
		    (let ((c (cos (cadr p))) (s (sin (cadr p))))
		      (dotimes (f 5)
			(let* ((r (* (expt 2 (- f 2)) (car p))) 
			       (y (floor (- i (* r s)))) 
			       (x (floor (+ j (* r c)))))
			  (when (array-in-bounds-p image y x)
			    (if (null (aref lmage y x))
				(setf (aref lmage y x) (make-array '(360 5))))
			    (incf (aref (aref lmage y x) u f)))))))))))))))
		      
(defun ght (treshold image)
  (let (r)
    (dotimes (i (array-dimension image 0) (sort r #'> :key #'fifth))
      (dotimes (j (array-dimension image 1))
	(if (not (null (aref image i j)))
	    (dotimes (k 360)
	      (dotimes (l 5)
		(if (< treshold (aref (aref image i j) k l))
		    (push (list i j k l (aref (aref image i j) k l)) r)))))))))

(defun mark-image (image lst)
  (let ((lmage (duplicate-image image)))
    (dolist (e lst lmage)
      (let ((y (car e)) (x (cadr e)))
	(loop for i from (- y 2) below (+ y 2) do
	     (loop for j from (- x 2) below (+ x 2) do
		  (if (array-in-bounds-p image i j)
		      (setf (aref lmage i j) 255))))))))

(defun binarize-image (image)
  (let ((lmage  (make-array (array-dimensions image))))
    (dotimes (i (array-dimension image 0) lmage)
      (dotimes (j (array-dimension image 1))
	(setf (aref lmage i j) (if (> 128 (aref image i j)) 0 255))))))

(defun treshold-fourier! (f image)
  (multiple-value-bind (a fft) (fft-n-mags! image)
    (let ((m (* f (nth-value 2 (image-max a)))))
      (dotimes (i (array-dimension image 0) image)
	(dotimes (j (array-dimension image 1))
	  (setf (aref image i j) 
		(if (> m (aref a i j)) 0 (aref fft i j))))))))

(defun join-centres (t0 t1 lst)
  (mapcar
   #'(lambda (g) ;(car (sort g #'> :key #'fifth)))
       (list (round (apply #'average (mapcar #'car g)))
	     (round (apply #'average (mapcar #'cadr g)))
	     1
	     (round (apply #'average (mapcar #'fourth g)))
	     (round (apply #'average (mapcar #'fifth g)))))
   (group 
    #'(lambda (a b)
	(and (> t0 (abs (- (fourth a) (fourth b))))
	     (> t1 (sqrt (+ (expt (- (first a) (first b)) 2)
			    (expt (- (second a) (second b)) 2))))))   
    lst :every t)))


(defun compare-complex (image image1)
  (dotimes (i (array-dimension image 0))
    (dotimes (j (array-dimension image 1))
      (let ((d (abs (- (aref image i j) (aref image1 i j))))
	    (m (min (abs (aref image i j)) (abs (aref image1 i j)))))
	(if (> (if (zerop m) d (/ d m)) 0.001)
	    (format t "~a ~a ~a ~a ~a ~a~%" (aref image i j) (aref image1 i j) i j d m))))))

(defun filter-centres (min max image lst)
  (remove-if-not  
   #'(lambda (c)
       ;(format t "~a~%" c)
       (and (array-in-bounds-p image (car c) (cadr c))
	    (< min (aref image (car c) (cadr c)) max))) 
   lst))

(defun shrink-image (f image)
  (let ((lmage (make-array (mapcar #'(lambda (d) (ceiling d f)) 
				   (array-dimensions image))
			   :initial-element nil))
	(lmage1 (make-array (mapcar #'(lambda (d) (ceiling d f)) 
				   (array-dimensions image))
			    :element-type '(unsigned-byte 8))))
    (dotimes (i (array-dimension image 0))
      (dotimes (j (array-dimension image 1)) 
	(push (aref image i j) (aref lmage (floor i f) (floor j f)))))
    (dotimes (i (array-dimension lmage 0) lmage1)
      (dotimes (j (array-dimension lmage 1))
	(setf (aref lmage1 i j) (round (apply #'average (aref lmage i j))))))))

(defun rotate-rect-points (a rect)
  (let* ((r (/ (* a pi) 180)) (c (cos r)) (s (sin r)))
    (values
     (%rotate-point c s (new-point (first rect) (second rect)))
     (%rotate-point c s (new-point (first rect) (fourth rect)))
     (%rotate-point c s (new-point (third rect) (fourth rect)))
     (%rotate-point c s (new-point (third rect) (second rect))))))

(defun rotate-roints (a rect)
  (let* ((c (cos a)) (s (sin a)))
    (values
     (%rotate-point c s (new-point (first rect) (second rect)))
     (%rotate-point c s (new-point (first rect) (fourth rect)))
     (%rotate-point c s (new-point (third rect) (fourth rect)))
     (%rotate-point c s (new-point (third rect) (second rect))))))

(defun restore-rect (rect a o image)
  (let ((p (multiple-value-list (rotate-rect-points a rect))))
    (let ((li (mapcar #'point-i p)) (lj (mapcar #'point-j p))
	  (h (array-dimension image 0))
	  (w (array-dimension image 1))
	  (hr (- (third rect) (first rect)))
	  (wr (- (fourth rect) (second rect))))
      (let ((maxi (min (1- h) (max 0 (+ (point-i o) (apply #'max li)))))
	    (mini (min (1- h) (max 0 (+ (point-i o) (apply #'min li)))))
	    (maxj (min (1- w) (max 0 (+ (point-j o) (apply #'max lj)))))
	    (minj (min (1- w) (max 0 (+ (point-j o) (apply #'min lj))))))
	(if (and (not (zerop (- maxi mini))) (not (zerop (- maxj minj))))
	    (let* ((ar (copy-rect (list mini minj maxi maxj) image))
		   (r (rotate-image 
		       (pad-image (max (- maxi mini) (- maxj minj)) ar)
		       (- (/ (* a pi) 180)))) ; хрень с углом
		   (r0 (array-dimension r 0))
		   (r1 (array-dimension r 1))
		   (i0 (max 0 (round (- (/ r0 2) (/ hr 2)))))
		   (i1 (min r0 (round (+ (/ r0 2) (/ hr 2)))))
		   (j0 (max 0 (round (- (/ r1 2) (/ wr 2)))))
		   (j1 (min r1 (round (+ (/ r1 2) (/ wr 2))))))
	      (pad-image hr (copy-rect (list i0 j0 i1 j1) r) wr)))))))

(defun dilate-points (pts)
  (remove-duplicates
   (apply
    #'append
    (mapcar
     #'(lambda (p)
	 (let ((i (point-i p)) (j (point-j p)))
	   (list 
	    p 
	    (new-point i (1+ j)) 
	    (new-point (1+ i) j)
	    (new-point (1+ i) (1+ j)))))
     pts))
   :test #'equal))

(defvar *this-file* (load-time-value
                     (or #.*compile-file-pathname* *load-pathname*)))


(defun pathify (prefix names)
  (map 
   'vector
   #'(lambda (name)
       (make-pathname 
	:name (concatenate 'string prefix name) 
	:type "lsp" :version nil :defaults *this-file*))
   names))

(defun slurp (proc paths)
  (map 
   'vector
   #'(lambda (path) (funcall proc (with-open-file (f path) (read f))))
   paths))

(defvar *omg* nil)
(defvar *scripts* 
  (pathify "script-" '("ogham" "runic" "latin" "cyrillic" "digits")))
(defvar *acorn* nil)
(defvar *oak* (slurp #'grafs->tree *scripts*))
(defvar *frg* 
  (slurp #'alist->hash 
	 (pathify "freq-" '("welsh" "english" "latin" "russian" "ukrainian"))))
(defvar *lng* #((2 . 0) (2 . 1) (2 . 2) (3 . 3) (3 . 4)))

(defun set-acorn! ()
  (setf *acorn* (slurp #'identity *scripts*))
  t)

(defun set-oak! ()
  (setf *oak* (slurp #'grafs->tree *scripts*))
  (when *acorn* (set-acorn!)))

(defun img<-path (path) 
  (instantiate-rectangle 
   (shrink-lst 10 (widen-lines 1 (render-svg-path path)))))

(defvar *plant* t)

(defun remember-blocks (i blocks &optional (stream t))
  (let* ((graph (blocks->graph blocks))
	 (c (read stream)))
    (push
     (list
      (etypecase c
	(character c)
	((integer 0 9) (digit-char c))
	(string (char c 0))
	(symbol (char (symbol-name c) 0)))
      graph) 
     (svref *acorn* i))
;    (format t "r ~a~%" (length *acorn*))
    (setf (svref *oak* i) (grafs->tree (svref *acorn* i)))
    (when *plant*  ;добавлять только новое?
      (grafs->file (svref *acorn* i) (svref *scripts* i)))
    graph))

(defun forget-last (i)
  (pop (svref *acorn* i))
  (setf (svref *oak* i) (grafs->tree (svref *acorn* i)))
  (grafs->file (svref *acorn* i) (svref *scripts* i)))

(defun remember-img (i img &optional (stream t))
  (remember-blocks i (image->blocks img) stream))

(defun stierlitz (image blocks)
  (let ((lmage (make-array (array-dimensions image) :initial-element 0)))
    (loop for b in blocks for c in '(1 2 3 4 5 6 7 8 9 a b c d e f g h i j k l m n o p q r s t u v w x y z) do
	 (dolist (p (fourth b))
	   (setf (aref lmage (point-i p) (point-j p)) c)))
    lmage))

(defun test-img (i image)
  (let ((blocks (image->blocks image)))
    (values 
     (stierlitz image blocks) 
     blocks (gather-fruit (svref *oak* i) (blocks->graph blocks)))))

(defun test-slice (i no &key box)
  (test-img i (nth no (funcall (if box #'boxify #'carve-into-boxes) *omg*))))

(defun prove-slice (i no &key box)
  (remember-img 
   i (nth no (funcall (if box #'boxify #'carve-into-boxes) *omg*))))

(defun approve-slice (i no &key box)
  (multiple-value-bind (k b c) (test-slice i no :box box)
    (format t "~%~a~%~a~%~a~%" k b c)
    (when (or (null c) (not (zerop (caar c)))) (remember-blocks i b))))

(defun grayscale->bitmap! (image)
  (dotimes (i (array-dimension image 0) image)
    (dotimes (j (array-dimension image 1))
      (setf (aref image i j) (round (aref image i j) 255)))))

;вэйвлеты
    
(defun haar-decomposition-step (h c)
  "переписал буквально"
  (if (<= h 1) c
    (let ((c/ (copy-seq c)) (half (/ h 2)))
      (haar-decomposition-step
       half
       (dotimes (i half c/)
	 (let ((2i (* i 2)) (2i+1 (1+ (* i 2))) (v2 (sqrt 2)))
	   (setf (svref c/ i) (/ (+ (svref c 2i) (svref c 2i+1)) v2))
	   (setf (svref c/ (+ half i)) 
		 (/ (- (svref c 2i) (svref c 2i+1)) v2))))))))

(defun haar-decomposition (c)
  (let* ((h (length c)) (s (sqrt h)))
    (haar-decomposition-step h (map 'vector #'(lambda (x) (/ x s)) c))))

(defun power-of-2-ceiling (n) (expt 2 (ceiling (log n 2))))

(defun row (image i)
  (let ((tmp (make-array (array-dimension image 1))))
    (dotimes (j (array-dimension image 1) tmp)
      (setf (svref tmp j) (aref image i j)))))

(defun col (image j)
  (let ((tmp (make-array (array-dimension image 0))))
    (dotimes (i (array-dimension image 0) tmp)
      (setf (svref tmp i) (aref image i j)))))

(defun (setf row) (r image i)
  (dotimes (j (length r) r)
    (setf (aref image i j) (svref r j))))

(defun (setf col) (c image j)
  (dotimes (i (length c) c)
    (setf (aref image i j) (svref c i))))

(defun haar-standard-decomposition! (image)
  (dotimes (i (array-dimension image 0))
    (setf (row image i) (haar-decomposition (row image i))))
  (dotimes (j (array-dimension image 1))
    (setf (col image j) (haar-decomposition (col image j))))
  image)

;(defun haar-standard-decomposition (image)
;  (haar-standard-decomposition! (copy-image-2 image)))

(defun l1-image-metric (i1 i2)
  (apply 
   #'+ 
   (mapcar 
    #'abs 
    (map 'list #'- 
	 (make-array (array-total-size i1) 
		     :element-type (array-element-type i1) 
		     :displaced-to i1)
	 (make-array (array-total-size i2) 
		     :element-type (array-element-type i2) 
		     :displaced-to i2)))))

;видимо, можно сделать вывод, что вэйвлеты L1 не слаще

; x-1, y+1 => i-1, j-1

(defparameter d8before '((-1 -1) (-1 0) (-1 1) (0 -1)))

(defparameter d8after '((0 0) (0 1) (1 -1) (1 0) (1 1)))

(defun coef-1 (p) (apply #'max (mapcar #'abs p)))

(defun coef-2 (p) (sqrt (apply #'+ (mapcar #'(lambda (x) (* x x)) p))))

(defun f (image cf half y x)
  (if (zerop (aref image y x)) 0
    (apply
     #'min
     (mapcar
      #'(lambda (p) 
	  (+ (funcall cf p) (cref image (+ y (car p)) (+ x (cadr p)) 100)))
      half))))  

(defun dransform! (image coef)
  (dotimes (y (array-dimension image 0))
      (dotimes (x (array-dimension image 1))
	(setf (aref image y x) (ceiling (f image coef d8before y x)))))
  (do ((y (1- (array-dimension image 0)) (1- y)))
      ((< y 0))
      (do ((x (1- (array-dimension image 1)) (1- x)))
	  ((< x 0))
	  (setf (aref image y x) (ceiling (f image coef d8after y x)))))  
  image)

(defun distance-transform! (image) (dransform! image #'coef-1)) 

(defun chamfer-transform! (image) (dransform! image #'coef-2))

(defun image-not (image)
  (let ((tmp (make-array (array-dimensions image) :element-type '(unsigned-byte 8))))
    (dotimes (i (array-dimension image 0) tmp)
      (dotimes (j (array-dimension image 1))
	(setf (aref tmp i j) (- 1 (aref image i j)))))))

(defun h (image dform)
;  (apply
;   #'max
;   (map 'list #'* 
;	(make-array 
;	 (array-total-size image) 
;	 :element-type (array-element-type image) 
;	 :displaced-to image)
;	(make-array 
;	 (array-total-size dform) 
;	 :element-type (array-element-type dform) 
;	 :displaced-to dform))))
  (let ((d 0))
    (dotimes (i (array-dimension image 0) d)
      (dotimes (j (array-dimension image 1))    
	(if (and (not (zerop (aref image i j)))
		 (< d (aref dform i j))) 
	    (setf d (aref dform i j)))))))

(defun hr (image r1 dform r2)
  (do ((d 0)
       (i1 (first r1) (1+ i1))
       (i2 (first r2) (1+ i2)))
      ((= i1 (third r1)) d)
      (do ((j1 (second r1) (1+ j1))
	   (j2 (second r2) (1+ j2)))
	  ((= j1 (fourth r1)))
	  (if (and (not (zerop (aref image i1 j1)))
		   (< d (aref dform i2 j2))) 
	      (setf d (aref dform i2 j2))))))

(defun ha (image r1 dform r2)
  (declare (optimize (speed 3) (safety 0)))
  (declare (type (simple-array (unsigned-byte 8) (* *)) image dform))
  (do ((d 0) (c 1)
       (i1 (the (unsigned-byte 24) (first r1)) 
	   (1+ (the (unsigned-byte 24) i1)))
       (i2 (the (unsigned-byte 24) (first r2)) 
	   (1+ (the (unsigned-byte 24) i2))))
      ((= i1 (the (unsigned-byte 24) (third r1))) (/ d c))
      (do ((j1 (the (unsigned-byte 24) (second r1)) 
	       (1+ (the (unsigned-byte 24) j1)))
	   (j2 (the (unsigned-byte 24) (second r2)) 
	       (1+ (the (unsigned-byte 24) j2))))
	  ((= j1 (the (unsigned-byte 24) (fourth r1))))
	  (if (not (zerop (aref image i1 j1)))
	      (setq d (+ (the (unsigned-byte 24) d) 
			 (the (unsigned-byte 8) (aref dform i2 j2))) 
		    c (1+ (the (unsigned-byte 24) c)))))))

(defun ht (image dform)
  (let ((d 0) (c 0))
    (dotimes (i (array-dimension image 0) (/ d c))
      (dotimes (j (array-dimension image 1))
	(if (not (zerop (aref image i j)))
	    (setq d (+ d (aref dform i j)) c (1+ c)))))))

(defun h-distance (i1 d1 i2 d2) (max (h i1 d2) (h i2 d1)))

(defun ht-distance (i1 d1 i2 d2) (max (ht i1 d2) (ht i2 d1)))

(defvar *treshold* 1)
(defvar *treshold2* 0.9)
(defvar *treshold3* 0.8)

(defun hr-distance (i1 d1 r1 i2 d2 r2) 
  (max (* *treshold3* (ha i1 r1 d2 r2)) (* *treshold2* (ha i2 r2 d1 r1))))

(defun hausdorff-distance (i1 i2)
  (h-distance i1 (chamfer-transform! (image-not i1))
	      i2 (chamfer-transform! (image-not i2))))

(defun hausdorth-distance (i1 i2)
  (ht-distance i1 (chamfer-transform! (image-not i1))
	       i2 (chamfer-transform! (image-not i2))))

(defun list-boxes (h w h1 w1)
  (let (result)
    (dotimes (i h (reverse result))
      (let ((ni (+ i h1)))
	(if (<= ni h)
	    (dotimes (j w)
	      (let ((nj (+ j w1)))
		(if (<= nj w) 
		    (push (list i j ni nj) result)))))))))