;+
;
; $Id: hi_align_image.pro,v 1.4 2011/11/17 21:03:34 nathan Exp $
;
; Project: STEREO-SECCHI
;
;
; Name: hi_align_image
;
; Purpose: Align image from one header to be consistent with another 
;
; Category: STEREO, SECCHI, HI
;
; Explanation: The wcs coordinates vary from one image to another
;              this routine allows a shift in the image to be applied
;              to aling the images.  If you are aligning images using
;              the radec option you should have removed the f-corona first
;
; Syntax: shifted=hi_align_image(hdr1,image1,hdr2[,/radec])
;
; Example: shifted=hi_align_image(hdr1,image1,hdr2,/radec)
;
; Inputs: hdr1 - the secchi header for the input image
;         image1 - the input image
;         hdr2 - the secchi header for the target image
;
; Keywords: radec - use ra/dec coordinates for alignment
;         : shift - use nearest neighbour rather than interpolate
;   	/ZEROMISSING	Set missing values=0 not NaN
;   	/VERBOSE    Print shift info
;
; Calls: fitshead2wcs
;
; Common: None
;
; Restrictions: None 
;
; Side effects: None
;
; Prev. Hist.: None
;
; History: Written 10 August 2008 by Steve Crothers 
;
; $Log: hi_align_image.pro,v $
; Revision 1.4  2011/11/17 21:03:34  nathan
; add /VERBOSE and clarify inputs description
;
; Revision 1.3  2011/11/15 23:05:22  nathan
; add /ZEROMISSING
;
; Revision 1.2  2009/10/07 13:45:43  crothers
; added nearest neighbour option for interpolation
;
; Revision 1.1  2008-11-19 12:29:52  crothers
; initial release
;
;
;-

function hi_align_image,src_hdr,image,dest_hdr,radec=radec,shift=shift,ZEROMISSING=zeromissing, $
    VERBOSE=verbose, _extra=_extra

if keyword_set(radec) then begin

    src_wcs=fitshead2wcs(src_hdr,system='A')

    dest_wcs=fitshead2wcs(dest_hdr,system='A')
endif else begin

    src_wcs=fitshead2wcs(src_hdr)

    dest_wcs=fitshead2wcs(dest_hdr)

endelse

IF keyword_set(VERBOSE) THEN print,'Shift x,y pixels to match '+dest_hdr.filename+': ' $
    +trim((dest_hdr.crval1-src_hdr.crval1)/src_hdr.cdelt1)+ ', ' $
    +trim((dest_hdr.crval2-src_hdr.crval2)/src_hdr.cdelt2)

; get location of destination points in source image coordinates
; pixels is [2,nx,ny] where first coordinate is x/y

coords=wcs_get_coord(dest_wcs)
pixels=wcs_get_pixel(src_wcs,coords)

IF keyword_set(ZEROMISSING) THEN missval=0 ELSE missval=!values.F_NAN
if keyword_set(shift) then begin
    dest=image*missval

    pixels=round(pixels)

    s=size(image)

    w=where(pixels[0,*,*] ge 0 and pixels[0,*,*] lt s[1] and pixels[1,*,*] ge 0 and pixels[1,*,*] le s[2])

    posn=pixels[1,*,*]*s[1]+pixels[0,*,*]

    dest[w]=image[posn[w]]
endif else begin

    dest=interpolate(image,pixels[0,*,*],pixels[1,*,*],missing=missval,_extra=_extra)
endelse

; drop leading single dimension

dest=reform(dest)

return,dest

end