;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; $Id: morpho.lsh,v 1.3 2002/09/23 02:06:59 profshadoko Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? *** Morphological Operations
;; Erosion, dilation, distance transforms,
;; thresholding, and masking on on images of ubytes.
;;.AUTHOR Leon Bottou 09/96

(libload "libidx/idx-macros")


#? (ubim-internal-disttrans <im> <dist>)
;; Returns the internal distance transformation matrix of binary image <im>.
;; The elements of the output ubyte matrix <dist> contains the distance to the
;; closest background pixel in the image (background level is zero for this
;; purpose). Distances are clipped to 255. The dimensions of <dist> must be 
;; 2 larger than that of <im>. <dist> will be resized to the appropriate
;; size if necessary. 
(de ubim-internal-disttrans(im dist)
    ((-idx2- (-ubyte-)) im dist)
    (let* ((width (idx-dim im 1))
           (height (idx-dim im 0)))
      ((-int-) width height)
      (when (or (<> (idx-dim dist 0) (+ 2 height))
		(<> (idx-dim dist 1) (+ 2 width)))
	(idx-u2resize dist (+ 2 height) (+ 2 width)))
      #{ {
        unsigned char *p, *q, *qup, lastqup, minim;
        int i, j, qm0, qm1, pm0, pm1;
        pm0 = $im->mod[0];
        pm1 = $im->mod[1];
        qm0 = $dist->mod[0];
        qm1 = $dist->mod[1];
        for (j=1; j<=$height; j++) {
          p = IDX_PTR( $im, unsigned char) + (j-1)*pm0;
          q = IDX_PTR( $dist, unsigned char) + j*qm0 + qm1;
          qup = q - qm0 + qm1*((j-1)&1); 
          minim= *(q - qm1);
          lastqup= *(qup - qm1);
          for (i=1; i<=$width; i++, q+=qm1, p+=pm1, qup+=qm1) {
            if (*p == 0) {
              *q=minim=0;
              lastqup= *qup;
            } else {
              if ( lastqup<minim ) minim=lastqup;
              if ( (lastqup= *qup)<minim ) minim=lastqup;
              if (minim < 255) minim++;
              *q = minim;
            }
          }
        }
        for (j=$height; j>0; j--) {
          q  = IDX_PTR($dist, unsigned char) + j*qm0 + $width*qm1; 
          qup= q + qm0 + (((j-1)&1)-1)*qm1; 
          minim = *(q + qm1);
          lastqup= *(qup + qm1);
          for (i=$width; i>0; i--, q-=qm1, qup-=qm1) {
            if (*q == 0) {
              minim = 0; lastqup = *qup;
            } else {
              if (lastqup<minim ) minim=lastqup;
              if ( (lastqup= *qup)<minim ) minim=lastqup;
              if (minim < 255) minim++;
              if ( minim<*q ) *q=minim; else minim= *q;
            }
          }
        } 
      } #} 
      ()))


#? (ubim-external-disttrans <im> <dist>)
;; computes the external distance transformation matrix of binary image <im>.
;; The elements of the output ubyte matrix <dist> contains the distance to the
;; closest non-background pixel in the image (background level is zero for this
;; purpose). Distances are clipped to 255. The dimensions of <dist> must be 
;; 2 larger than that of <im>. <dist> will be resized to the appropriate
;; size if necessary. 
(de ubim-external-disttrans(im dist)
    ((-idx2- (-ubyte-)) im dist)
    (let* ((width (idx-dim im 1))
           (height (idx-dim im 0)))
      ((-int-) width height)
      (when (or (<> (idx-dim dist 0) (+ 2 height))
		(<> (idx-dim dist 1) (+ 2 width)))
	(idx-u2resize dist (+ 2 height) (+ 2 width)))
      #{ {
        unsigned char *p, *q, *qup, lastqup, minim;
        int i, j, qm0, qm1, pm0, pm1;
        pm0 = $im->mod[0];
        pm1 = $im->mod[1];
        qm0 = $dist->mod[0];
        qm1 = $dist->mod[1];
        q = IDX_PTR($dist, unsigned char);
        qup = q + $width*qm1 + qm1;
        for (i=0; i<$height+2; i++, q+=qm0, qup+=qm0)
           *q = *qup = 255;
        q = IDX_PTR($dist, unsigned char);
        qup = q + $height*qm0 + qm0;
        for (i=0; i<$width+2; i++, q+=qm1, qup+=qm1)
           *q = *qup = 255;
        for (j=1; j<=$height; j++) {
          p = IDX_PTR( $im, unsigned char) + (j-1)*pm0;
          q = IDX_PTR( $dist, unsigned char) + j*qm0 + qm1;
          qup = q - qm0 + qm1*((j-1)&1); 
          minim= *(q - qm1);
          lastqup= *(qup - qm1);
          for (i=1; i<=$width; i++, q+=qm1, p+=pm1, qup+=qm1) {
            if (*p > 0) {
              *q=minim=0;
              lastqup= *qup;
            } else {
              if ( lastqup<minim ) minim=lastqup;
              if ( (lastqup= *qup)<minim ) minim=lastqup;
              if (minim < 255) minim++;
              *q = minim;
            }
          }
        }
        for (j=$height; j>0; j--) {
          q  = IDX_PTR($dist, unsigned char) + j*qm0 + $width*qm1; 
          qup= q + qm0 + (((j-1)&1)-1)*qm1; 
          minim = *(q + qm1);
          lastqup= *(qup + qm1);
          for (i=$width; i>0; i--, q-=qm1, qup-=qm1) {
            if (*q == 0) {
              minim = 0; lastqup = *qup;
            } else {
              if (lastqup<minim ) minim=lastqup;
              if ( (lastqup= *qup)<minim ) minim=lastqup;
              if (minim < 255) minim++;
              if ( minim<*q ) *q=minim; else minim= *q;
            }
          }
        } 
      } #} 
      ()))



#? (ubim-positive-threshold <im> <thres>)
;; Applies a threshold on ubyte image <im>.
;; All values smaller than or equal to threshold are zeroed.
;; Other values are set to 255.

(de ubim-positive-threshold(im thres)
    ((-idx2- (-ubyte-)) im)
    ((-int-) thres)
    (idx-bloop ((imlin im))
      (cinline-idx1loop imlin "unsigned char" "mp" "i"
                        #{ if (*mp <= $thres) *mp = 0; else *mp = 255; #} ) )
    () )

#? (ubim-negative-threshold <im> <thres>)
;; Applies a threshold on ubyte image <im>.
;; All values smaller than or equal to threshold are set to 255.
;; Other values are set to 0

(de ubim-negative-threshold(im thres)
    ((-idx2- (-ubyte-)) im)
    ((-int-) thres)
    (idx-bloop ((imlin im))
      (cinline-idx1loop imlin "unsigned char" "mp" "i"
                        #{ if (*mp > $thres) *mp = 0; else *mp = 255; #} ) )
    () )




#? (ubim-erosion <im> <factor>)
;; Perform in place erosion of ubyte image <im> at distance <factor>.
(de ubim-erosion(im factor)
    ((-idx2- (-ubyte-)) im)
    ((-int-) factor)
    (let* ((w (idx-dim im 1))
	   (h (idx-dim im 0))
	   (dist (ubyte-matrix-nc (+ 2 h) (+ 2 w))))
      (ubim-internal-disttrans im dist)
      (ubim-positive-threshold dist factor)
      (idx-copy (narrow (narrow dist 0 h 1) 1 w 1) im)
      ()))


#? (ubim-dilation <im> <factor>)
;; Perform in place dilation of ubyte image <im> at distance <factor>.

(de ubim-dilation(im factor)
    ((-idx2- (-ubyte-)) im)
    ((-int-) factor)
    (let* ((w (idx-dim im 1))
	   (h (idx-dim im 0))
	   (dist (ubyte-matrix-nc (+ 2 h) (+ 2 w))))
      (ubim-external-disttrans im dist)
      (ubim-negative-threshold dist factor)
      (idx-copy (narrow (narrow dist 0 h 1) 1 w 1) im)
      ()))


#? (ubim-mask <im> <mask>)
;; Sets to zero all bytes of ubyte image <im> whose corresponding
;; pixel in <mask> is zero. This is a AND operation.
(de ubim-mask(im mask)
    ((-idx2- (-ubyte-)) im mask)
    (idx-bloop ((imlin im)(msklin mask))
      (cinline-idx1loop2 imlin "unsigned char" msklin "unsigned char" 
                        "m1" "m2" "i1" "i2"
                        #{ if (*m2==0) (*m1=0); #} ) ) ())






;;; ----------------------------------------
;;; COMPILE

(dhc-make () 
          ubim-internal-disttrans
          ubim-external-disttrans 
          ubim-positive-threshold 
          ubim-negative-threshold 
          ubim-erosion
          ubim-dilation
          ubim-mask
          )