; $Id: warp.pro,v 1.1 1993/04/02 18:38:01 idl Exp $

; Warp demo for Sun display with mouse.  Modifications
; will have to be made for other displays.
;
; Compile subroutines that we use:
pro get_cp,a,b,x0,y0,x1,y1 ;Obtain control points for two images
		;in frame buffers 0 and 1.
; It is assumed that the images have been loaded into the two frame
; buffers.

; Each control point is marked with a 3 by 3 rectangle of 0's.
; On exit:
; x0,y0 = vectors containing coordinates of control points in buffer 0.
; x1,y1 = vectors containing coordinates of control points in buffer 1.
;

common colors, r_orig, g_orig, b_orig, r_curr, g_curr, b_curr

device,set_w=255		;Use all colors
if n_elements(r_orig) eq 0 then begin	;colors defined?
	r_orig=indgen(256) & g_orig = r_orig & b_orig = r_orig
	tvlct,r_orig, g_orig, b_orig
	endif

p1 = 16 * [ lindgen(256)/16, lindgen(256) and 15 ] ;512 array of 2 tables

;	** Load both images by writing 4 bits from each **
device,set_write=15	;load bottom 4 bits from b
tv,b/16b
device,set_write=240
tv,a		;Show A in top 4 bits
empty		;Flush out display

mx = 200		;Maximum # of points to accept
xx = intarr(mx,2)	;Accept up to mx points
yy = xx
;
n = 0			;count of points
print,'Enter each control point in the two images by pressing the left mouse button.'
print,'Press the right mouse button when entering the last control point to terminate.'
repeat begin		;point loop
	for i=0,1 do begin	;get each point
		p = p1(i*256 : i*256 + 255)
		tvlct,r_orig(p), g_orig(p), b_orig(p)
		empty
		print,'Enter point',n+1,', Image ',i
		tvrdc,x,y,/dev	;Get coordinate, dev coords
		xx(n,i) = x
		yy(n,i) = y
		device,set_write=([240, 15])(i)
		tv,bytarr(5,5),x-1,y-1	 ;Show rectangle
		empty
		wait,.5		;Try to eliminate contact bounce
		endfor
	n=n+1
	endrep until !err eq 4
x0 = xx(0:n-1,0) & y0 = yy(0:n-1,0)
x1 = xx(0:n-1,1) & y1 = yy(0:n-1,1)
device,set_write=255		;Restore write mask
tvlct,r_orig, g_orig, b_orig
return
end	
pro rms2d, xi, yi, xo, yo, kx, ky
deg = fix(sqrt(n_elements(kx)))-1 ;Get degree
n = n_elements(xi)
xp = fltarr(n)	;X,Y predicted
yp = xp
xii = replicate(1.,n)	;X^i
for i=0,deg do begin
	yii = xii	;X^i * y^j
	for j=0,deg do begin
		xp = xp + kx(j,i) * yii
		yp = yp + ky(j,i) * yii
		yii = yii * yi
		endfor
	xii = xii * xi
	endfor
xrms = sqrt(total((xo - xp)^2)) / n
yrms = sqrt(total((yo - yp)^2)) / n
print,'X RMS error = ',xrms
print,'Y RMS error = ',yrms
return
end


;
; This IDL program illustrates polynomial image warping.
;
; It reads and displays two aerial images of 
; the same landscape taken from two different
; flights into the two frame buffers.
;
; Then the user enters a bunch of control points in each frame and
; the program warps them and displays the result.
; The user should enter at least four control points for a linear fit,
;	and 9 points for a quadratic fit.
;

;
; This reads in two images, A (512,392), and B (512,355).
;
if !d.window lt 0 then window	;Why this is necessary, I don't know
openr, 1, filepath('warp_data.dat', subdirectory='images')
aa=assoc(1,bytarr(512,392))
a=aa(0)
b=aa(1)
close,1
get_cp,a,b,x0,y0,x1,y1	;Get control points
y0=y0 & y1 = y1		;Remove bias from Y because of display offset
c = b			;Make image to mess with
a0=a			;First image
for i=0,n_elements(x1)-1 do begin
	c(x1(i)-1, y1(i)-1) = bytarr(5,5) ;image for play
	a0(x0(i)-1, y0(i)-1) = bytarr(5,5)
	endfor
for i=0,511,100 do c(i,*)=255	;Make grid lines.
for i=0,300,100 do c(*,i)=255
ndeg = fix(sqrt(n_elements(x0)))-1	;Degree of polynomial obtained from
			;the number of control points
print,'Using polynomial of order',ndeg
polywarp,x1,y1,x0,y0, ndeg, kx,ky ;Least square fit for coeff.
rms2d,x0,y0,x1,y1,kx,ky	;Print error statistics
			;Warp it & also make result same size as A
c = poly_2d(c,kx,ky,0,512,392,miss=0) 

flick,a0,c,8			;And flicker result, 8 frames/sec
end