hextract.pro

basic median filter simulation

pro hextract, oldim, oldhd, newim, newhd, x0, x1, y0, y1, SILENT = silent, $
    ERRMSG = errmsg,ALT = alt  
;+
; NAME:
;       HEXTRACT
; PURPOSE:
;       Extract a subimage from an array and update astrometry in FITS header
; EXPLANATION:
;       Extract a subimage from an array and create a new FITS header with
;       updated astrometry for the subarray
; CALLING SEQUENCE:
;       HEXTRACT, Oldim, Oldhd, [ Newim, Newhd, x0, x1, y0, y1, /SILENT ]
;               or
;       HEXTRACT, Oldim, Oldhd, [x0, x1, y0, y1, /SILENT, ERRMSG =  ]    
;
; INPUTS:
;       Oldim - the original image array
;       Oldhd - the original image header
;
; OPTIONAL INPUTS:
;       x0, x1, y0, y1 - respectively, first and last X pixel, and first and
;       last Y pixel to be extracted from the original image, integer scalars.
;       HEXTRACT will convert these values to long integers. 
;       If omitted,  HEXTRACT will prompt for these parameters
;
; OPTIONAL OUTPUTS:
;       Newim - the new subarray extracted from the original image 
;       Newhd - header for newim containing updated astrometry info
;               If output parameters are not supplied or set equal to
;               -1, then the HEXTRACT will modify the input parameters 
;               OLDIM and OLDHD to contain the subarray and updated header.
;
; OPTIONAL INPUT KEYWORD:
;      ALT - Single character 'A' through 'Z' or ' ' specifying which astrometry
;          system to modify in the FITS header.    The default is to use the
;          primary astrometry or ALT = ' '.    See Greisen and Calabretta (2002)
;          for information about alternate astrometry keywords.
;      /SILENT - If set and non-zero, then a message describing the extraction
;               is not printed at the terminal.   This message can also be 
;               suppressed by setting !QUIET.
; OPTIONAL KEYWORD OUTPUT:
;       ERRMSG - If this keyword is supplied, then any error mesasges will be
;               returned to the user in this parameter rather than depending on
;               on the MESSAGE routine in IDL.   If no errors are encountered
;               then a null string is returned.               
;
; PROCEDURE:
;       The FITS header parameters NAXIS1, NAXIS2, CRPIX1, and CRPIX2 are
;       updated for the extracted image.
;
; EXAMPLE:  
;       Read an image from a FITS file 'IMAGE', extract a 512 x 512 subimage 
;       with the same origin, and write to a new FITS file 'IMAGENEW'
;
;       IDL> im = READFITS( 'IMAGE', hdr )      ;Read FITS files into IDL arrays
;       IDL> hextract, im, h, 0, 511, 0, 511    ;Extract 512 x 512 subimage
;       IDL> writefits, 'IMAGENEW', im ,h       ;Write subimage to a FITS file
;
; PROCEDURES CALLED
;       CHECK_FITS, STRN(), SXPAR(), SXADDPAR, SXADDHIST
; MODIFICATION HISTORY:
;       Written, Aug. 1986 W. Landsman, STX Corp.
;       Use astrometry structure,   W. Landsman      Jan, 1994
;       Minor fix if bad Y range supplied   W. Landsman    Feb, 1996
;       Added /SILENT keyword              W. Landsman     March, 1997
;       Added ERRMSG keyword    W. Landsman   May 2000
;       Work for dimensions larger than 32767   W.L., M.Symeonidis Mar 2007
;       Added ALT keyword  W.L. April 2007
;- 
 On_error, 2
 compile_opt idl2
 npar = N_params()

 if (npar EQ 3) or (npar LT 2) then begin       ;Check # of parameters
    print,'Syntax - HEXTRACT, oldim, oldhd, [ newim, newhd, x0, x1, y0, y1]'
    print,'   or    HEXTRACT, oldim, oldhd, x0, x1, y0, y1, [/SILENT, ERRMSG=]'
    return
 endif
 
 save_err = arg_present(errmsg)      ;Does user want to return error messages?
;                                    Check for valid 2-D image & header
  check_FITS, oldim, oldhd, dimen, /NOTYPE, ERRMSG = errmsg
  if errmsg NE '' then begin
        if not save_err then message,'ERROR - ' + errmsg,/CON
        return
  endif

  if N_elements(dimen) NE 2 then begin 
           errmsg = 'Input image array must be 2-dimensional'
           if not save_err then message,'ERROR - ' + errmsg,/CON
           return
  endif

  xsize = dimen[0]  &  ysize = dimen[1]


 if ( npar LT 4 ) then Update = 1 else Update = 0     ;Update old array?

 if ( npar EQ 6 ) then begin                 ;Alternative calling sequence ?

     if ( N_elements(newim) EQ 1 ) and ( N_elements(newhd) EQ 1 ) and $
        ( N_elements(x0) EQ 1 ) and ( N_elements(x1) EQ 1 ) then begin
              y0 = x0   &  y1 = x1
              x0 = newim   &   x1 = newhd
              Update = 1
      endif 

 endif

 RDX: 
 if ( npar LE 5 )  then begin

      message, /INF, $ 
           'Original array size is ' + strn(xsize) + ' by ' + strn(ysize) 
      x0 = 0l & x1 = 0l
      read,'% HEXTRACT: Enter first and last X pixel to be extracted: ',x0,x1

 endif

 if ( x1 LT x0 ) or ( x0 LT 0 ) or ( x1 GE xsize ) then begin

     message,'ERROR - Illegal pixel range: X direction',  /CON
     print, ' '
     message, /INF,   $
     ' Legal Range is 0 < First Pixel < Last Pixel < ' + strn(xsize-1)
     if update then npar = npar < 2 else npar = npar < 4
     goto, RDX 

 endif

 RDY: if ((not update) and ( npar LE 7 )) or (update and (npar LT 6) ) then $ 
    read,'% HEXTRACT: Enter first and last Y pixel to be extracted: ',y0,y1

 if ( y1 LT y0 ) or ( y0 LT 0 ) or ( y1 GE ysize ) then begin

     message,'ERROR - Illegal pixel range: Y direction', /CON
     message, /INF,     $ 
      'Legal Range is 0 < First Pixel < Last Pixel < ' + strn(ysize-1)
     if update then npar = npar < 4 else npar = npar < 6 
     goto, RDY

 endif

 x0 = long(x0) & x1 = long(x1)
 y0 = long(y0) & y1 = long(y1)                                          

 naxis1 = x1 - x0 + 1 
 naxis2 = y1 - y0 + 1   ;New dimensions

 if not keyword_set(SILENT) then message, /INF,        $
      'Now extracting a '+ strn(naxis1) + ' by ' + strn(naxis2) + ' subarray'

  if Update then oldim = oldim[ x0:x1,y0:y1 ]        $
            else newim = oldim[ x0:x1,y0:y1 ]

 newhd = oldhd
 sxaddpar, newhd, 'NAXIS1', naxis1                                   
 sxaddpar, newhd, 'NAXIS2', naxis2
 label = 'HEXTRACT: ' + systime(0)

 hist = [label,'Original image size was '+ strn(xsize) + ' by ' + strn(ysize), $
         'Extracted Image: [' + strn(x0) + ':'+ strn(x1) +  $
         ',' + strn(y0) + ':'+ strn(y1) + ']'  ]

 sxaddhist, hist, newhd


;GSSS image uses CNPIX instead of CRPIX
   cnpix1 = sxpar( oldhd, 'CNPIX1', COUNT = Ncnpix1)
         if ( Ncnpix1 EQ 1 ) then begin   ;Shift position of reference pixel

                sxaddpar, newhd, 'CNPIX1', cnpix1+x0
                cnpix2 = sxpar( oldhd, 'CNPIX2' )
                sxaddpar, newhd, 'CNPIX2', cnpix2+y0
        endif

; Update astrometry info if it exists

  if N_elements(alt) EQ 0 then alt = ''
  extast, newhd, astr, noparams, ALT = alt
  if noparams GE 0 then begin

  sxaddpar, newhd, 'CRPIX1'+alt, astr.crpix[0]-x0
  sxaddpar, newhd, 'CRPIX2'+alt, astr.crpix[1]-y0


 endif 
 
 if Update then begin

      oldhd = newhd
      newim = x0 & newhd = x1
      x0 = y0 & x1 = y1

 endif

 return
 end