;
;+
;
;  Routine matches a pair images to have the same
;	magnification, rotation, and XY position.

PRO image_match, image1, image2, outimg1, outimg2, scale1, scale2, angle1, angle2, xyref1, xyref2, xyout, scaleout

;
; Input Parameters:
;	image1, image2 = 2-d input images.  May be of any type.
;
;	scale1, scale2 = image scales in units per pixel.  
;			Assumes that both x and y axes have same scale.
;
;	angle1, angle2 = Rotation angle in degrees, measured from the
;			+Y direction toward -X,  to a common reference
;			direction.  
;
;	xyref1, xyref2 = 2-element vectors giving the location of a common
;			reference location in original pixels.  The
;			2 elements are the (x,y) distance from the (0,0)
;			pixel.  The reference location can be outside
;			both images.
;
; Output Parameters:
;	outimg1, outimg2 = output images, magnified, shifted, rotated, to
;			a common coordinate system.  Output images
;			are always larger than the input images to
;			avoid losing data.  The reference direction is
;			toward +Y.
;
;	xyout   =  2-element vector giving the location of the reference
;			location in output pixels.  Same convention as
;			input reference vector.
;
;	scaleout  =  output image scale, units per pixel.  This is
;			always the smaller of the 2 input scales.
;
;
; Routines Called:
;	ROT_TRIM
;
; Procedure:
;	Rotates and magnifies images to common system.  Computes
;	vector offset of image centers, inserts images into a common
;	large format that completely includes both fields.
;
; History:
;	Written Jun 8, 1993  Barry LaBonte
;
;-
;

size1 = SIZE(image1)
size2 = SIZE(image2)

; Final scale is smallest of all
fscale = scale1 <  scale2

; First magnify and rotate images
ROT_TRIM, image1, angle1, scale1/fscale, out1, x1off, y1off
ROT_TRIM, image2, angle2, scale2/fscale, out2, x2off, y2off
osiz1 = SIZE(out1)
osiz2 = SIZE(out2)

; Find offset of image centers
; Vectors in inital frames
x1c = (size1(1)/2. - xyref1(0)) * scale1
y1c = (size1(2)/2. - xyref1(1)) * scale1
x2c = (size2(1)/2. - xyref2(0)) * scale2
y2c = (size2(2)/2. - xyref2(1)) * scale2

; Rotate into reference frame
c1 = COS(!DTOR * angle1)
s1 = SIN(!DTOR * angle1)
c2 = COS(!DTOR * angle2)
s2 = SIN(!DTOR * angle2)

x1r = x1c * c1 - y1c * s1
y1r = y1c * c1 + x1c * s1
x2r = x2c * c2 - y2c * s2
y2r = y2c * c2 + x2c * s2

xoff = FIX( (x1r - x2r)/fscale )
yoff = FIX( (y1r - y2r)/fscale )
ix1 = xoff > 0
iy1 = yoff > 0
ix2 = (-xoff) > 0
iy2 = (-yoff) > 0



; Now make final arrays and insert with offsets
dimx =  ( ix1 + osiz1(1)) > ( ix2 + osiz2(1) )
dimy =  ( iy1 + osiz1(2)) > ( iy2 + osiz2(2) )
outimg1 = MAKE_ARRAY( DIMENSION=[dimx, dimy], TYPE=size1(size1(0)+1) )
outimg2 = MAKE_ARRAY( DIMENSION=[dimx, dimy], TYPE=size2(size2(0)+1) )

outimg1( ix1:ix1+osiz1(1)-1, iy1:iy1+osiz1(2)-1 ) = out1
outimg2( ix2:ix2+osiz2(1)-1, iy2:iy2+osiz2(2)-1 ) = out2

; Figure out where (0,0) is
scaleout = fscale
xyout = FLTARR(2)
xyout(0) =  (ix2 + osiz2(1)/2) - x2r/fscale
xyout(1) =  (iy2 + osiz2(2)/2) - y2r/fscale

RETURN
END