;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; 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: rgbafimage.lsh,v 1.6 2006/02/27 15:43:47 laseray Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? *** RGBA Images of floats
;; basic image processing on RGBA images
;; where each pixel component is a float.


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; warping functions for RGBA images of floats
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(libload "libidx/idx-macros")
(libload "libidx/idx-sort")
(libload "libidx/idx-convol")
(libload "libimage/image-transform")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


#? rgbafimage
;; a rgbafimage is an idx3 of flts intended to store images.
;; The last dimension is 4, and contains the R,G,B, and Alpha
;; channels respectively.

;; transformation matrices between color spaces

(defvar rgb-uvw 
 [f[ 0.405 0.116 0.133]
   [ 0.299 0.587 0.114]
   [ 0.145 0.827 0.627]])

(defvar rgb-yuv
  [f[  0.299  0.587  0.114 ] 
    [ -0.147 -0.289  0.437 ] 
    [  0.615 -0.515 -0.100 ]])

(defvar rgb-yiq
 [f[  0.299  0.587  0.114 ]
   [  0.596 -0.274 -0.322 ]
   [  0.212 -0.523  0.311 ]])

(defvar rgb-ycrcb 
 [f[ 0.257   0.504   0.098  16]
   [ 0.439  -0.368  -0.071  128]
   [-0.148  -0.291   0.439  128]
   [ 0.0     0.0     0.0    1 ]])

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; subsampling
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (rgbafim-subsample <in> <nlin> <ncol>)
;; subsamples image <in> with integer ratios <nlin> (vertical) <ncol>
;; (horizontal). the subsampled image is returned.
;; The horizontal (resp. vertical) size of the returned image is 
;; equal to the integer part of the horizontal (resp vertical)
;; size of <in> divided by <ncol> (resp <nlin>).
(de rgbafim-subsample (in nlin ncol)
    ((-idx3- (-flt-)) in)
    ((-int-) nlin ncol)
    (let* ((h (idx-dim in 0))
           (w (idx-dim in 1))
           (nh (int (/ h nlin)))
           (nw (int (/ w ncol)))
           (inp (narrow (narrow in 0 (* nlin nh) 0) 1 (* ncol nw) 0))
           (out (float-matrix-nc nh nw 4)))
      #{{
      flt *_idx2loopc1, *pin;
      flt *_idx2loopc2, *pout;
      int i, _imax = ($out)->dim[0];
      int j, _jmax = ($out)->dim[1];
      int _imat1_m0 = ($inp)->mod[0];
      int _imat1_m1 = ($inp)->mod[1];
      int _imat2_m0 = ($out)->mod[0];
      int _imat2_m1 = ($out)->mod[1];
      int pin_incr = $ncol * _imat1_m1;
      int norm = 1.0/($ncol * $nlin);
      register int acc0, acc1, acc2, acc3;
      register int k,l;
      register flt *pinptr;
      register int pinptr_incr = _imat1_m0 - $ncol* _imat1_m1;

      _idx2loopc1 = IDX_PTR (($inp), flt);
      _idx2loopc2 = IDX_PTR (($out), flt);
      for (i = 0; i < _imax; i++) {
        pin = _idx2loopc1;
        pout = _idx2loopc2;
        for (j = 0; j < _jmax; j++) {
          acc0 = acc1 = acc2 = acc3 =0;
          pinptr = pin;
          for (k=0; k<$nlin; k++) {
            for (l=0; l<$ncol; l++) { 
              acc0 += pinptr[0];
              acc1 += pinptr[1];
              acc2 += pinptr[2];
              acc3 += pinptr[3];
              pinptr += _imat1_m1;
              }
            pinptr += pinptr_incr;
            }
          pout[0] = acc0 * norm;
          pout[1] = acc1 * norm;
          pout[2] = acc2 * norm;
          pout[3] = acc3 * norm;
          pin += pin_incr;
          pout += _imat2_m1;
	  }
	  _idx2loopc1 += _imat1_m0 * $nlin;
	  _idx2loopc2 += _imat2_m0;
	}
        } #}
  out))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; color processing
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; look in Digital Image Processing by W. Pratt, page 66.

#? (rgbafim-luminance <in> <out>)
;; compute luminance of each pixel and put result in image <out>
;; the formula used is lum = 0.299 * R + 0.587 * G + 0.114 * B
(de rgbafim-luminance (in out)
    ((-idx3- (-flt-)) in)
    ((-idx2- (-flt-)) out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{ *pout = 0.299 * pin[0] + 0.587 * pin[1] + 0.114 * pin[2]; #})
      ()))

#? (rgbafim-lum2rgba <in> <out>)
;; transform a grey-level image <in> (idx2 of flts) into
;; an RGBA image. <in> and <out> must have identical sizes,
;; though no check of that is done
(de rgbafim-lum2rgba (in out)
    ((-idx2- (-flt-)) in)
    ((-idx3- (-flt-)) out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{ pout[0] = pout[1] = pout[2] = *pin; pout[3] = 0; #})
      ()))

#? (rgbafim-chromin-broken <in> <z> <k> <out>)
;; attempts to compute the chrominance info in <in>
;; (normalized colors with almost identical luminance)
;; and put the result in <out>. Both <in> and <out>
;; are idx3 of flts. <z> is a multiplier for the
;; average pixel component value. <k> is a fudge factor
;; that prevents the normalisation from blowing up
;; if the color is black or almost black.
(de rgbafim-chromin-broken (in z k out)
    ((-flt-) z k)
    ((-idx3- (-flt-)) in out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{{ flt x, n, r, g, b;
         r = pin[0]; g = pin[1]; b = pin[2];
         x = 0.0841*r*r + 0.3481*g*g + 0.0144*b*b;
         n = x+1.0;
         n = $z*sqrt(n/(x*n+$k));
         pout[0] = r*n; pout[1] = g*n; pout[2] = b*n; pout[3] = pin[3];
         }#}
      ()))

#? (rgbafim-chrominance <in> <z> <k> <out>)
;; attempts to compute the chrominance info in <in>
;; (normalized colors with almost identical luminance)
;; and put the result in <out>. Both <in> and <out>
;; are idx3 of flts. <z> is a multiplier for the
;; average pixel component value. <k> is a fudge factor
;; that prevents the normalisation from blowing up
;; if the color is black or almost black.
(de rgbafim-chrominance (in z k out)
    ((-flt-) z k)
    ((-idx3- (-flt-)) in out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{{ flt x, n, r, g, b;
         r = pin[0]; g = pin[1]; b = pin[2];
         x = 0.299*r + 0.587*g + 0.114*b;
         n = $z/(x+$k);
         pout[0] = r*n; pout[1] = g*n; pout[2] = b*n; pout[3] = pin[3];
         }#}
      ()))

#? (rgbafim-uvw <in> <out>)
;; transform an RGBA image into an UVWA image, where the V component
;; contains the luminance information, and the U and W components 
;; collectively contain the chrominance information.
;; In the UVW color coding scheme the euclidean distance is meant
;; to reflect the subjective "perceptual" distance.
;; See Digital Image Processing by W. Pratt, page 66.
(de rgbafim-uvw (in out)
    ((-idx3- (-flt-)) in)
    ((-idx3- (-flt-)) out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{{ register float p0=pin[0],p1=pin[1],p2=pin[2];
         pout[0] = 0.405 * p0 + 0.116 * p1 + 0.133 * p2;
         pout[1] = 0.299 * p0 + 0.587 * p1 + 0.114 * p2;
         pout[2] = 0.145 * p0 + 0.827 * p1 + 0.627 * p2; 
         pout[3] = pin[3]; }#})
      ()))

#? (rgbafim-fromuvw <in> <out>)
;; converts back from UVWA to RGBA.
;; It is advasable to saturate the resulting image, as 
;; numerical errors may cause the pixel values to exceed the
;; range 0-255.
(de rgbafim-fromuvw (in out)
    ((-idx3- (-flt-)) in)
    ((-idx3- (-flt-)) out)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{{ register float p0=pin[0],p1=pin[1],p2=pin[2];
	 pout[0] =  2.432 * p0 + 0.332 * p1 - 0.576 * p2;
         pout[1] = -1.519 * p0 + 2.083 * p1 - 0.057 * p2;
         pout[2] =  1.440 * p0 - 2.823 * p1 + 1.803 * p2; 
         pout[3] = pin[3]; }#})
      ()))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; pixel processing (saturation, transformation...)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (rgbafim-saturate <in> <min> <max> <out>)
;; Saturate pixel component values in image <in> to the range
;; [<min>,<max>]. Result is put in <out>.
(de rgbafim-saturate (in min max out)
    ((-idx3- (-flt-)) in out)
    ((-flt-) min max)
    (cinline-idx2loop2 in "flt" out "flt" "pin" "pout" "i" "j"
      #{{ register flt r;
         r = pin[0]; pout[0] = (r>$max)? $max : (r<$min) ? $min : r;
         r = pin[1]; pout[1] = (r>$max)? $max : (r<$min) ? $min : r;
         r = pin[2]; pout[2] = (r>$max)? $max : (r<$min) ? $min : r;
         r = pin[3]; pout[3] = (r>$max)? $max : (r<$min) ? $min : r;
         } #})
      ()))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; blitting
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (rgbafim-blit <x> <y> <in> <out>)
;; blit <in> into <out> at position <x>,<y>
;; with alpha blending.
;; an alpha value of 0 in <in> means opaque, 
;; 255 means totally transparent.
(de rgbafim-blit (x y in out)
    ((-idx3- (-flt-)) in out)
    ((-int-) x y)
    (let* ((hin (idx-dim in 0))
           (win (idx-dim in 1))
           (hout (idx-dim out 0))
           (wout (idx-dim out 1))
           (lx (if (< y 0) 0 y))
           (ly (if (< x 0) 0 x))
           (ux (if (> (+ x win) wout) wout (+ x win)))
           (uy (if (> (+ y hin) hout) hout (+ y hin)))
           (dest (narrow (narrow out 0 (- uy ly) ly) 1 (- ux lx) lx))
           (src  (narrow (narrow in  0 (- uy ly) (- ly y)) 
                         1 (- ux lx) (- lx x))))
      ((-int-) hin hout win wout lx ly ux uy)
      (cinline-idx2loop2 src "flt" dest "flt" "pin" "pout" "i" "j"
        #{{ register flt r, r1;
            r = 0.00392157 * pin[3]; r1 = 1-r;
            pout[0] = r * pout[0] + r1 * pin[0];
            pout[1] = r * pout[1] + r1 * pin[1];
            pout[2] = r * pout[2] + r1 * pin[2];
            } #})
      ()))

#? (rgbafim-blitcolor <x> <y> <in> <r> <g> <b> <out>)
;; blit color defined by <r>,<g>,<b> using idx2 <in> as a stencil
;; into RGBA image <out> at position <x> <y>
;; The values in <in> must be between 0 and 255. 
;; if the value is 255, no paint is put in; if the value is 0,
;; the pixel takes the value specified by <r>,<g>,<b>.
(de rgbafim-blitcolor (x y in r g b out)
    ((-idx2- (-flt-)) in)
    ((-idx3- (-flt-)) out)
    ((-int-) x y)
    ((-flt-) r g b)
    (let* ((hin (idx-dim in 0))
           (win (idx-dim in 1))
           (hout (idx-dim out 0))
           (wout (idx-dim out 1))
           (lx (if (< y 0) 0 y))
           (ly (if (< x 0) 0 x))
           (ux (if (> (+ x win) wout) wout (+ x win)))
           (uy (if (> (+ y hin) hout) hout (+ y hin)))
           (dest (narrow (narrow out 0 (- uy ly) ly) 1 (- ux lx) lx))
           (src  (narrow (narrow in  0 (- uy ly) (- ly y)) 
                         1 (- ux lx) (- lx x))))
      ((-int-) hin hout win wout lx ly ux uy)
      (cinline-idx2loop2 src "flt" dest "flt" "pin" "pout" "i" "j"
        #{{ register flt v, v1;
            v = 0.00392157*(*pin); v1 = 1-v;
            pout[0] = v*pout[0] + v1*$r;
            pout[1] = v*pout[1] + v1*$g;
            pout[2] = v*pout[2] + v1*$b;
            } #})
      ()))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; color histogram in 32x32x32
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (rgbafim-init-lum2rgb )
;; build a table of 32x32x32 elements containing RGB cell indices
;; in ascending luminosity.
(de rgbafim-init-lum2rgb ()
    (cheader "#define NCOLCELL 32768")
    (cheader "static int lum2rgb[NCOLCELL], lum2rgbinit = 0;")
    (let ((lum (float-matrix (* 32 32 32)))
	  (ind (int-matrix (* 32 32 32)))
	  (n 0))
      ((-int-) n)
      (for (r 0 31) (for (g 0 31) (for (b 0 31)
	   (ind n n)
	   (lum n (+ (* .299 r) (* .587 g) (* .114 b)))
	   (incr n))))
      (idx-f1i1sortup lum ind)
      #{{
        int *pind, i;
	pind = IDX_PTR($ind, int);
	for (i=0; i<NCOLCELL; i++){ lum2rgb[i] = pind[i]; }
	lum2rgbinit = 1;
      }#}
      ()))

#? (rgbafim-lum2rgb <n>)
;; converts a luminosity index to a 5-bit RGB cell number.
(de rgbafim-lum2rgb (n) ((-int-) n) (int #{ lum2rgb[(int)$n] #} ))

#? (rgbafim-histo32 <rgbafim> <ppal> <count>)
;; compute histogram of image <rgbafim>.
;; <ppal> must be a 32768 by 3 matrix which will contain the list of pixel
;; values found. <count> must be a 32768 matrix which will contain
;; the corresponding pixel counts.
;; pixel colors are sorted by luminance.
;; This function really computes a 32768 color palette by assigning each pixel
;; to the cube it belongs to in the 32x32x32 RBG color cube.
;; The color prototype assigned to each cube is the mean of the pixels 
;; belonging to the cube.
(de rgbafim-histo32 (rgbafim ppal count)
    ((-idx3- (-flt-)) rgbafim)
    ((-idx2- (-flt-)) ppal)
    ((-idx1- (-flt-)) count)
    (let* ((ncolcell (* 32 32 32))
	   (histo (float-matrix ncolcell 3))
	   (hcount (float-matrix ncolcell))
	   (phisto (gptr #{ IDX_PTR($histo,flt) #}))
	   (phcount (gptr #{ IDX_PTR($hcount,flt) #})))
      ;; loop over image pixels and compute histogram
      ;; for each histogram cell in the 32x32x32 table we compute the
      ;; average color value for the pixels (it's like a 32k color palette).
      (cinline-idx2loop rgbafim "flt" "p" "i" "j"
			#{{ register flt *ph, *pc ;
			int offset ;
			offset = ((((int)(p[0]))&0xf8)<<7) | 
                                 ((((int)(p[1]))&0xf8)<<2) | 
                                  (((int)(p[2]))>>3) ;
			ph = ((flt*)$phisto) + offset*3 ;
			pc = ((flt*)$phcount) + offset ;
			/* printf("i=%d, j=%d, offset=%d, *pc=%g, ph[]=%g, %g, %g\n",
				  i,j,offset,*pc,ph[0],ph[1],ph[2]); */
			ph[0] += (flt)p[0];
			ph[1] += (flt)p[1];
			ph[2] += (flt)p[2];
			*pc += 1 ; 
			} #} )
      ;; now compact the list of color cells and sort them by luminosity
      ;; (using a precomputed table of cell index -> luminosity)
      (when (= 0 (int #{ lum2rgbinit #})) (rgbafim-init-lum2rgb))
      (let ((ccount 0))
	((-int-) ccount)
	(for (i 0 (1- ncolcell))
	  ((-int-) i)
	  (let* ((c (rgbafim-lum2rgb i)))
	    ((-int-) c)
	    (let ((cnt (hcount c)))
	    ((-flt-) cnt)
	    (when (> cnt 0)
	      ;; (printf "i=%d, c=%d, ccount=%d\n" i c ccount)
	      (ppal ccount 0 (/ (histo c 0) cnt))
	      (ppal ccount 1 (/ (histo c 1) cnt))
	      (ppal ccount 2 (/ (histo c 2) cnt))
	      (count ccount cnt)
	      (incr ccount)))))
	(idx-f2resize ppal ccount 3)
	(idx-f1resize count ccount)) ()))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; non-linear warping functions with and without bilinear interpolation.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (rgbafim-interpolate-bilin <background> <pin> <indimi> <indimj> <inmodi> <inmodj> <ppi> <ppj> <rez>)
;;.VP
;; (-gptr- pin rez background):
;; (-int-  indimi indimj inmodi inmodj ppi ppj):
;;
;; SIDE EFFECTS: <rez>
;; AUTHOR: Y. LeCun
;; COMPILABLE: Yes
;; DESCRIPTION:
;; returns a bilinearly interpolated RGBA pixel value for coordinate
;; <ppi> <ppj>. The image data is pointed to by <pin>, with 
;; <indimi> <indimj> <inmodi> <inmodj> being the dimensions and modulos.
;; This function clips automatically if <ppi> <ppj> are outside of the
;; bounds by assuming the outside values are equal to <background>.
;; pixel values are flts, while coordinates are 32 bit fixed point
;; with 16 bit integer part and 16 bit fractional part.
(de rgbafim-interpolate-bilin 
    (background pin indimi indimj inmodi inmodj ppi ppj rez)
    ((-gptr- "float *") pin background rez)
    ((-int-)  indimi indimj inmodi inmodj ppi ppj)
    #{{
         int li0, lj0;
	 register int li1, lj1;
	 flt deltai, ndeltai;
	 flt deltaj, ndeltaj;
	 register flt *pin00, *back;
	 register flt *v00, *v01, *v10, *v11;
         back = (flt *)$background;
	 li0 = $ppi >> 16;
	 li1 = li0+1;
	 deltai = 0.0000152587890625*(flt)($ppi & 0x0000ffff);
	 ndeltai = 1.0 - deltai;
	 lj0 = $ppj  >> 16;
	 lj1 = lj0+1;
	 deltaj = 0.0000152587890625*(flt)($ppj & 0x0000ffff);
	 ndeltaj = 1.0 - deltaj;
	 pin00 = (flt*)($pin)+$inmodi*li0+$inmodj*lj0;
         if ((li1>0)&&(li1<$indimi)) {
	     if ((lj1>0)&&(lj1<$indimj)) { 
		 v00 = (pin00); 
		 v01 = (pin00+$inmodj); 
		 v11 = (pin00+$inmodi+$inmodj); 
		 v10 = (pin00+$inmodi);
		 v00 = (pin00); 
	     } else if (lj1==0) { 
		 v00 = back;
		 v01 = (pin00+$inmodj); 
		 v11 = (pin00+$inmodi+$inmodj);
		 v10 = back;
	     } else if (lj1==$indimj) {
		 v00 = (pin00); 
		 v01 = back;
		 v11 = back;
		 v10 = (pin00+$inmodi); 
	     } else {
		 v00 = back;
		 v01 = back;
		 v11 = back;
		 v10 = back;
	     }
	 } else if (li1==0) {
	     if ((lj1>0)&&(lj1<$indimj)) { 
		 v00 = back;
		 v01 = back;
		 v11 = (pin00+$inmodi+$inmodj); 
		 v10 = (pin00+$inmodi);
	     } else if (lj1==0) { 
		 v00 = back;
		 v01 = back;
		 v11 = (pin00+$inmodi+$inmodj);
		 v10 = back;
	     } else if (lj1==$indimj) {
		 v00 = back;
		 v01 = back;
		 v11 = back;
		 v10 = (pin00+$inmodi); 
	     } else {
		 v00 = back;
		 v01 = back;
		 v11 = back;
		 v10 = back;
	     }
	 } else if (li1==$indimi) {
	     if ((lj1>0)&&(lj1<$indimj)) { 
		 v00 = (pin00); 
		 v01 = (pin00+$inmodj); 
		 v11 = back;
		 v10 = back;
	     } else if (lj1==0) { 
		 v00 = back;
		 v01 = (pin00+$inmodj); 
		 v11 = back;
		 v10 = back;
	     } else if (lj1==$indimj) {
		 v00 = (pin00); 
		 v01 = back;
		 v11 = back;
		 v10 = back;
	     } else {
		 v00 = back;
		 v01 = back;
		 v11 = back;
		 v10 = back;
	     }
         } else {
	     v00 = back;
	     v01 = back;
	     v11 = back;
	     v10 = back;
	 }
         /* printf("deltai=%d, deltaj=%d, r00=%d, r10=%d, r01=%d, r11=%d\n",
                deltai, deltaj, v00[0], v10[0], v01[0], v11[0]); */
         { register flt *ubout;
           ubout = (flt *)($rez);
	   ubout[0] = (ndeltaj * (( v10[0]*deltai + v00[0]*ndeltai )) +
                       deltaj  * (( v11[0]*deltai + v01[0]*ndeltai )));
	   ubout[1] = (ndeltaj * (( v10[1]*deltai + v00[1]*ndeltai )) +
                       deltaj  * (( v11[1]*deltai + v01[1]*ndeltai )));
	   ubout[2] = (ndeltaj * (( v10[2]*deltai + v00[2]*ndeltai )) +
                       deltaj  * (( v11[2]*deltai + v01[2]*ndeltai )));
	   ubout[3] = (ndeltaj * (( v10[3]*deltai + v00[3]*ndeltai )) +
                       deltaj  * (( v11[3]*deltai + v01[3]*ndeltai )));
         }
     }

  #} ())


#? (rgbafim-warp <in> <out> <background> <pi> <pj>)
;;.VP
;; ((-idx1- (-flt)) background):
;; ((-idx3- (-flt-)) in out):
;; ((-idx2- (-int-)) pi pj):
;;
;; RETURNS: Null
;; SIDE EFFECTS: <out> 
;; AUTHOR: Y. LeCun
;; COMPILABLE: Yes
;; DESCRIPTION:
;; Warps an image using <pi> and <pj> as tabulated coordinate transforms.
;; <in> and <out> are idx2 of flts. <background> is the value assumed outside
;; of the input image. <pi> and <pj> are tabulated coordinates which can
;; be filled up using compute-bilin-transform or similar functions.
;; Pixel values are antialiased using bilinear interpolation.
(de rgbafim-warp (in out background pi pj)
    ((-idx1- (-flt-)) background)
    ((-idx3- (-flt-)) in out)
    ((-idx2- (-int-)) pi pj)
    (when (or (<> 4 (idx-dim in 2)) 
              (<> 4 (idx-dim out 2)) 
              (<> 4 (idx-dim background 0)))
      (error "last dimension of in, out and background must be 4"))
    (let ((ni (idx-dim in 0))
          (nj (idx-dim in 1))
          (pin (to-gptr (cinline "IDX_PTR(%s, flt)" in)))
          (pback (to-gptr (cinline "IDX_PTR(%s, flt)" background)))
          (indimi (idx-dim in 0))
          (indimj (idx-dim in 1))
          (inmodi (idx-modulo in 0))
          (inmodj (idx-modulo in 1)))
      ((-gptr- "float *") pin pback)
      ((-int-)  indimi indimj inmodi inmodj)
      (idx-bloop ((lout out) (lpi pi) (lpj pj))
        (idx-bloop ((llout lout) (llpi lpi) (llpj lpj))
          (let ((ppi (llpi)) (ppj (llpj))
                (pout (to-gptr (cinline "IDX_PTR(%s, flt)" llout))))
            ((-int-) ppi ppj) ((-gptr- "float *") pout)
            (rgbafim-interpolate-bilin 
             pback pin indimi indimj inmodi inmodj ppi ppj pout)))))
    ())

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; fast version, not interpolation.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


#? (rgbafim-warp-fast <in> <out> <background> <pi> <pj>)
;;.VP
;; ((-int-) background):
;; ((-idx2- (-flt-)) in out):
;; ((-idx2- (-int-)) pi pj):
;;
;; RETURNS: Null
;; SIDE EFFECTS: <out>
;; AUTHOR: Y. LeCun
;; COMPILABLE: Yes
;; DESCRIPTION:
;; Warps an image using <pi> and <pj> as tabulated coordinate transforms.
;; <in> and <out> are idx2 of flts. <background> is the value assumed outside
;; of the input image. <pi> and <pj> are tabulated coordinates which can
;; be filled up using compute-bilin-transform or similar functions.
;; This is essentially identical to warp-rgbafimage, except no antialiasing
;; is performed (it goes a lot faster, but is not nearly as nice).
(de rgbafim-warp-fast (in out background pi pj)
    ((-idx1- (-flt-)) background)
    ((-idx3- (-flt-)) in out)
    ((-idx2- (-int-)) pi pj)
    (unprotect
      (let ((ni (idx-dim in 0))
            (nj (idx-dim in 1))
            (pin (to-gptr (cinline "IDX_PTR(%s, flt)" in)))
            (back (to-gptr (cinline "IDX_PTR(%s, flt)" background)))
            (indimi (idx-dim in 0))
            (indimj (idx-dim in 1))
            (inmodi (idx-modulo in 0))
            (inmodj (idx-modulo in 1)))
        ((-gptr- "float *") pin back)
        ((-int-)  indimi indimj inmodi inmodj)
        (idx-bloop ((lout out) (lpi pi) (lpj pj))
          (idx-bloop ((llout lout) (llpi lpi) (llpj lpj))
            (let ((ppi (llpi)) (ppj (llpj)))
              ((-int-) ppi ppj) 
              #{{
  		     register int li, lj;
                     register flt *out, *in;
                     out = IDX_PTR($llout, flt);
  		     li = ($ppi+0x7f) >> 16;
  		     lj = ($ppj+0x7f) >> 16;
  		     if ((li>=0)&&(li<$indimi)&&(lj>=0)&&(lj<$indimj)) {
                       in = (flt*)($pin)+$inmodi*li+$inmodj*lj;
                       out[0] = in[0]; out[1] = in[1]; 
                       out[2] = in[2]; out[3] = in[3]; 
		     } else {
                       out[0] = *((flt *)($back));
                       out[1] = *((flt *)($back)+1);
                       out[2] = *((flt *)($back)+2);
                       out[3] = *((flt *)($back)+3);
		     }
  		   }
                 #} 
  	         )))))
      ())

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; top-level function


#? (rgbafim-warp-quad <in> <out> <background> <mode> <x1> <y1> <x2> <y2> <x3> <y3> <x4> <y4> <p1> <q1> <p3> <q3>)
;;.VP
;; ((-flt-) x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3):
;; ((-int-) background mode):
;; ((-idx2- (-flt-)) in out):
;;
;; RETURNS: Null
;; SIDE EFFECTS: <out>
;; AUTHOR: Y. LeCun
;; COMPILABLE: Yes
;; DESCRIPTION:
;; transforms rgbafimage <in> (idx2 of flt) mapping quadrilateral <x1> <y1> <x2>
;; <y2> <x3> <y3> <x4> (points are numbered clockwise starting from upper left)
;; to rectangle whose upper left and lower right points are <p1><q1>, <p3><q3>.
;; result is put in <out> (idx2 of flt).  Clipping is automatically performed.
;; pixels outside of the bounds of <in> are assumed to have the value
;; <background>.  <mode>=0: no antialiasing, <mode>=1: antialiasing with
;; bilinear interpolation (2 times slower).  execution time on sparc 10 is about
;; 5 ms in mode 0 and 10 ms in mode 1 for a 32x32 target image.
(de rgbafim-warp-quad (in out background mode x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3)
    ((-flt-) x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3)
    ((-idx1- (-flt-)) background)
    ((-int-) mode)
    ((-idx3- (-flt-)) in out)

    (let* ((outi (idx-dim out 0))
	   (outj (idx-dim out 1))
	   (dispi (int-matrix-nc outi outj))
	   (dispj (int-matrix-nc outi outj)))
      (compute-bilin-transform dispi dispj x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3)
      (if (= 0 mode)
	  (rgbafim-warp-fast in out background dispi dispj)
	(rgbafim-warp in out background dispi dispj)))
    ())

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(dhc-make "rgbafimage" 
          rgbafim-subsample rgbafim-luminance
          rgbafim-lum2rgba rgbafim-chromin-broken
          rgbafim-chrominance
          rgbafim-uvw rgbafim-fromuvw
          rgbafim-saturate
          rgbafim-blit
          rgbafim-blitcolor
          rgbafim-init-lum2rgb 
          rgbafim-lum2rgb
          rgbafim-histo32
          rgbafim-interpolate-bilin 
          rgbafim-warp rgbafim-warp-fast
          rgbafim-warp-quad )
