;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; 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: image-transform.lsh,v 1.1 2002/08/19 07:01:28 profshadoko Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; bilinear linear coordinate transform
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(libload "libidx/idx-macros")

#? *** Computing Image Warping Maps
;; compute displacement tables for 2D geometric 
;; transformations of images, such as bilinear
;; warps that transform any quadrilateral into 
;; any rectangle.

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

#? (compute-bilin-transform <dispi> <dispj> <x1> <y1> <x2> <y2> <x3> <y3> <x4> <y4> <p1> <q1> <p3> <q3>)
;;.VP
;; ((-idx2- (-int-)) dispi):
;; ((-idx2- (-int-)) dispj):
;; ((-flt-) x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3):
;;
;; RETURNS: Null
;; SIDE EFFECTS: <dispi> <dispj>
;; AUTHOR: Y. LeCun
;; COMPILABLE: Yes
;; DESCRIPTION:
;; Tabulates a bilinear transformation that maps the quadrilateral defined
;; by the <xi> <yi> to a rectangle whose upper left point is <p1> <q1> and
;; lower right point is <p3> <q3>.
;; <x1> <y1> is the upper left point on the quadrilateral, and the points
;; are numbered clockwise.
;; <dispi> and <dispj> must be idx2 of int. On output, element (i,j)
;; of <dispi> and <dispj> will be filled respectively with the abscissa
;; and ordinate of the point in the quadrilateral that maps to point (i,j).
;; Values in these idx2 are interpreted as 32bit fixed point numbers
;; with 16 bit integer part and 16 bit fractional part.
;; <dispi> and <dispj> can subsequently be used as input to 
;; warp-shimage, or warp-shimage-fast.
;; This function makes minimal use of floating point arithmetics.
(de compute-bilin-transform (dispi dispj x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3)
    ((-idx2- (-int-)) dispi)
    ((-idx2- (-int-)) dispj)
    ((-flt-) x1 y1 x2 y2 x3 y3 x4 y4 p1 q1 p3 q3)

    ;; compute transformation matrix from coordinates
    ;; in target (rectangular) space to coordinates
    ;; in original (irregular quadrilateral) image
    ;; transformation matrix is in 16.16 fixed point format.
    (let* ((k (/ 65536 (* (- p3 p1) (- q3 q1))))
	   (x41 (- x4 x1))
	   (x21 (- x2 x1))
	   (x43 (- x4 x3))
	   (x23 (- x2 x3))
	   (mx0 (* k (+ (* q3 x21) (* q1 x43))))
	   (mx1 (* k (+ (* p3 x41) (* p1 x23))))
	   (mx2 (* (- k) (+ x41 x23)))
	   (mx3 (* k (- (+ (* p3 q3 x1) (* p1 q1 x3)) 
			(+ (* p1 q3 x2) (* p3 q1 x4)))))
	   (y41 (- y4 y1))
	   (y21 (- y2 y1))
	   (y43 (- y4 y3))
	   (y23 (- y2 y3))
	   (my0 (* k (+ (* q3 y21) (* q1 y43))))
	   (my1 (* k (+ (* p3 y41) (* p1 y23))))
	   (my2 (* (- k) (+ y41 y23)))
	   (my3 (* k (- (+ (* p3 q3 y1) (* p1 q1 y3)) 
			 (+ (* p1 q3 y2) (* p3 q1 y4))))))
      ((-flt-) k)
      ((-flt-) x41 x21 x43 x23)
      ((-int-) mx0 mx1 mx2 mx3)
      ((-flt-) y41 y21 y43 y23)
      ((-int-) my0 my1 my2 my3)
      (cinline-idx2loop2 dispi "int" dispj "int" "di" "dj" "q" "p"
	#{ 
	*di = $my0*(p) + $my1*(q) + $my2*(p)*(q) + $my3;
	*dj = $mx0*(p) + $mx1*(q) + $mx2*(p)*(q) + $mx3;
	#}))
    ())

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

(dhc-make () compute-bilin-transform)


