postagestamp.pro

;+
; NAME:
;  POSTAGESTAMP
;
; PURPOSE:
;  This function extracts a subregion from a 2 or 3-dimensional FITS
;  image. It uses a nearest-neighbor regridding scheme. 
;
; CATEGORY:
;  image processing
;
; CALLING SEQUENCE:
;  result=POSTAGESTAMP(im, head, cen, wid, scale, [/NAN])
;
; INPUTS:
;  IM: A 2- or 3-dimensional image array
;  HEAD: The FITS header corresponding to IM (as is generated by
;        MRDFITS, for example)
;  CEN: A  2 or 3 element array giving the center of the desired
;       subfield in 'real sky units' (i.e. whatever real sky units are
;       assumed within the CRVAL keywords of HEAD)
;  WID: A 2 or 3 element array giving the size of the extracted
;       subfield in 'real sky units' (see note for CEN)
;  SCALE: A 2 or 3 element array giving the size of the output pixels
;         in real sky units (whatever real sky units are assumed in
;         the CD or CDELT keywords of HEAD) 
;
; KEYWORD PARAMETERS:
;  NAN: If set and nonzero, regions of the subimage which lie outside
;  the boundaries of the input image are output as NaN. Otherwise,
;  they are output as zero.
;
;  NPIX: If set, this is a 2 or 3 element array giving the pixel size
;  of the output image. It replaces the role of the WID parameter.
;
;  OUTHEAD: A named variable to hold the postage stamp header
;
; OUTPUTS:
;  The postage stamp described by cen, wid, and scale
;
; EXAMPLE:
;  Given im and head, extract a .25 x .25 degree field centered at
;  (a,d)=(10.,20.) where each pixel is .002 degrees
;
;    result=postagestamp(in,head,[10.,20.],[.25,.25],.[.002,.002])
;  
;  Repeat, but this time preserve the convention that RA increases to
;  the left in sky images:
;
;    result=postagestamp(in,head,[10.,20.],[.25,.25],[-.002,.002])
;
; RESTRICTIONS:
;  The coordinates of CEN will always be the exact center of the
;  output image, and the pixels will always be exactly the size given
;  by scale. However, the output images have an odd number of pixels
;  along each axis, and thus may be a pixel or two different from the
;  size requested by WID.
;
;  This function uses SKY2PIX, which provides rather limited
;  astrometry capabilities (e.g. no corrections for field curvature or
;  other higher order distortion terms). If you need more robust
;  astrometry, consider using the IDL astro library routine
;  HEXTRACT. However, that routine will not handle 3D images.
;
; MODIFICATION HISTORY:
;  Written by:       Chris Beaumont August 9, 2008
;  December 17 2008: Added NPIX keyword. cnb
;  Feb 16, 2009: Added a warning if the requested postage stamp lies
;                outside the bounds of the input image. cnb.
;  Feb 16, 2009: Added the OUTHEAD keyword. cnb.
;-
FUNCTION postagestamp, im, head, cen, wid, scale, $
                       nan=nan, npix = npix, outhead = outhead
compile_opt idl2
on_error, 2

;-check inputs
if n_params() ne 5 then begin
    print,'Calling Squence: '
    print,'  subim = postagestamp(im, head, cen, wid, scale)'
    print,'     im: Image to sample'
    print,'   head: Header corresponding to image'
    print,'    cen: Center of output image in real sky coordinates'
    print,'    wid: Width of window in real sky coordinates'
    print,'  scale: Pixel scale. Same units as cdelt in head (usually degrees/pix)'
    return,0
endif

;-get header info
ast=nextast(head)

;-check parameters
if (n_elements(cen) ne ast.naxis) or $
  (n_elements(wid) ne ast.naxis) or $
  (n_elements(scale) ne ast.naxis) then $
  message,'Image described by cen, wid and scale does not have a consistent number '+$
          'of dimensions, or has a dimensionality different from the input image'

if (ast.naxis lt 2) or (ast.naxis gt 3) then $
  message, 'Input image restricted to be 2 or three dimensional'

;-determine pixel length of putput images. Force it to be an odd
;number, so that the location of cen is the exact center of the output image
if ~keyword_set(npix) then npix=ceil(abs(wid/scale/2))*2+1

;-output reference values.
lout = ( findgen(npix[0])-(npix[0]-1) / 2. ) * scale[0] + cen[0]
mout = ( findgen(npix[1])-(npix[1]-1) / 2. ) * scale[1] + cen[1]
if (ast.naxis eq 3) then $
  nout = ( findgen(npix[2])-(npix[2]-1) / 2. ) * scale[2] + cen[2]

;-convert to pixel values
lout = reform( rebin(  lout, npix[0], npix[1] ), npix[0] * npix[1], /over)
mout = reform( rebin(1#mout, npix[0], npix[1] ), npix[0] * npix[1], /over)

xy = round(sky2pix(head, transpose([[lout],[mout]])))
x= -1 > reform(xy[0,*],npix[0],npix[1]) < ast.sz[0]
y= -1 > reform(xy[1,*],npix[0],npix[1]) < ast.sz[1]

;- determine which pixels lie outside the bounds of the input image
off = (x eq -1) or (x eq ast.sz[0]) or (y eq -1) or (y eq ast.sz[1])
bad=where(off, badct, ncomplement = goodct)

if badct ne 0 then begin
    x[bad]= 0
    y[bad] = 0
endif

if goodct eq 0 then begin
   message,/continue, 'WARNING: Requested postage stamp and input'$
           +' do not overlap.'
endif

if (ast.naxis eq 2) then begin
    result=im[x,y] 
    if badct ne 0 then result[bad] = keyword_set(nan) ? !values.f_nan : 0
endif else begin
    mask=bytarr(npix[0],npix[1])+1
    if badct ne 0 then mask[bad]=(keyword_set(nan)?!values.f_nan:0)
    z=0 > round((nout-ast.crval[2])/ast.cd[2,2]+ast.crpix[2]-1) < (ast.sz[2]-1)
    zplane=fltarr(npix[0],npix[1])
    result=fltarr(npix[0],npix[1],npix[2])
    for i=0,npix[2]-1, 1 do $
      result[0,0,i]=im[x,y,zplane+z[i]]*mask
endelse

;- create new header, if requested
if arg_present(outhead) then begin
   mkhdr, outhead, result
   sxaddpar, outhead, 'crval1', cen[0]
   sxaddpar, outhead, 'crval2', cen[1]
   sxaddpar, outhead, 'crpix1', (npix[0] + 1) / 2.
   sxaddpar, outhead, 'crpix2', (npix[1] + 1) / 2.
   sxaddpar, outhead, 'cdelt1', scale[0]
   sxaddpar, outhead, 'cdelt2', scale[1]
   
   if (ast.naxis eq 3) then begin
      sxaddpar, outhead, 'crval3', cen[2]
      sxaddpar, outhead, 'crpix3', (npix[2] + 1) / 2.
      sxaddpar, outhead, 'cdelt3', scale[2]
   endif

   sxaddhist, 'Created with postagestamp.pro on ' + systime(), outhead
endif

return, result
end