fits_read.pro

basic median filter simulation

              end else begin
;
;       case 3: read next extension
;
                enum = fcb.last_extension + 1
                if (enum eq 0) and (fcb.naxis[0] eq 0) then enum = 1
            end
        end
;
; check to see if it is a valid extension
;
        if enum gt fcb.nextend then begin
                message='EOF encountered'
                goto,error_exit
        end
;
; extract information from FCB for the extension
;
        bitpix = fcb.bitpix[enum]
        naxis = fcb.naxis[enum]
        if naxis gt 0 then axis = fcb.axis[0:naxis-1,enum]
        gcount = fcb.gcount[enum]
        pcount = fcb.pcount[enum]
        xtension = fcb.xtension[enum]
;
; read header ================================================================
;
        if data_only then goto,read_data
        h = bytarr(80,36,/nozero)
        nbytes_in_file = fcb.nbytes
        position = fcb.start_header[enum]
        point_lun,fcb.unit,position
        first_block = 1         ; first block in header flag
        repeat begin
             if position ge nbytes_in_file then begin
                 message = 'EOF encountered while reading header'
                 goto,error_exit
             endif

             readu,fcb.unit,h
             position = position + 2880
             hdr = string(h>32b)
             endline = where(strmid(hdr,0,8) eq 'END     ',nend)
             if nend gt 0 then hdr = hdr[0:endline[0]]
             if first_block then header = hdr else header = [header,hdr]
             first_block = 0
        end until (nend gt 0)
;
; extract some header information
;
        bscale = sxpar(header,'bscale', Count = N_bscale)
        bzero = sxpar(header,'bzero', Count = N_bzero)
        if bscale eq 0.0 then bscale = 1.0
        unsgn_int = (bitpix EQ 16) and (Bzero EQ 32768) and (bscale EQ 1)
        unsgn_lng = (bitpix EQ 32) and (Bzero EQ 2147483648) and (bscale EQ 1)
        if (unsgn_int or unsgn_lng) then $
	        if not keyword_set(no_unsigned) then noscale = 1
        if (N_bscale gt 0) and (noscale eq 0) and (data_only eq 0) and $
           (last eq 0) and (header_only eq 0) then sxaddpar,header,'bscale',1.0
        if (N_bzero gt 0) and (noscale eq 0) and (data_only eq 0) and $
           (last eq 0) and (header_only eq 0) then sxaddpar,header,'bzero',0.0
        groups = sxpar(header,'groups')
;
; create header with form:
;       ! Required Keywords
;       ! BEGIN MAIN HEADER ------------------------------------------
;       ! Primary data unit header keywords
;       ! BEGIN EXTENSION HEADER -------------------------------------
;       ! Extension header keywords
;       ! END           
;
;
; add Primary Data Unit header to it portion of the header to it, unless the
; NO_PDU keyword is set, or the INHERIT keyword is not found or set to false
;
	       
	if no_pdu EQ 0 then no_pdu = 1 - (sxpar(header,'INHERIT') > 0)
        if pdu then no_pdu = 0
        if (no_pdu eq 0) and (enum gt 0) then begin
	
;
; delete required keywords
;
        sxdelpar,header,['SIMPLE','BITPIX','NAXIS','NAXIS1', $
                         'NAXIS2','NAXIS3','NAXIS4','NAXIS5', $
                         'NAXIS6','NAXIS7','NAXIS8','EXTEND', $
                         'PCOUNT','GCOUNT','GROUPS', $
                         'XTENSION']
	

; create required keywords
;
        hreq = strarr(20)
        hreq[0] = 'END     '

        if enum eq 0 then $
                sxaddpar,hreq,'SIMPLE','T','image conforms to FITS standard' $
           else sxaddpar,hreq,'XTENSION',xtension,'extension type'

        sxaddpar,hreq,'bitpix',bitpix,'bits per data value'
        sxaddpar,hreq,'naxis',naxis,'number of axes'
        if naxis gt 0 then for i=1,naxis do $
                sxaddpar,hreq,'naxis'+strtrim(i,2),axis[i-1]
        if (enum eq 0)and (fcb.nextend GE 1) then $
                sxaddpar,hreq,'EXTEND','T','file may contain extensions'
        if groups then sxaddpar,hreq,'GROUPS','T','Group format'
        if (enum gt 0) or (pcount gt 0) then $
                     sxaddpar,hreq,'PCOUNT',pcount,'Number of group parameters'
        if (enum gt 0) or (gcount gt 0) then $
                    sxaddpar,hreq,'GCOUNT',gcount,'Number of groups'
       n0 = where(strmid(hreq,0,8) eq 'END     ') & n0=n0[0]
            hpdu = fcb.hmain
            n1 = n_elements(hpdu)
            if n1 gt 1 then begin               
                hreq = [hreq[0:n0-1], $
                        'BEGIN MAIN HEADER ---------------------------------', $
                        hpdu[0:n1-2], $
                        'BEGIN EXTENSION HEADER ----------------------------', $
                        'END     ']
                n0 = n0 + n1 + 1
            end
;
; add extension header
;
        header = [hreq[0:n0-1],header]
        end
        if header_only then begin
                data = 0
                goto,done
        endif
;
; Read Data ===================================================================
;
read_data:
        if naxis eq 0 then begin        ;null image?
                data = 0
;
; check for implicit data specified by NPIX1, NPIX2, and PIXVALUE (provided
; the header was red, i.e. data_only was not specified)
;
                if data_only eq 0 then begin
                        NPIX1 = sxpar(header,'NPIX1')
                        NPIX2 = sxpar(header,'NPIX2')
                        PIXVALUE = sxpar(header,'PIXVALUE')
                        if (NPIX1*NPIX2) gt 0 then $
                                data = replicate(pixvalue,npix1,npix2)
                end
                goto,done
        endif

        case BITPIX of
           8:   IDL_type = 1          ; Byte
          16:   IDL_type = 2          ; Integer*2
          32:   IDL_type = 3          ; Integer*4
         -32:   IDL_type = 4          ; Real*4
         -64:   IDL_type = 5          ; Real*8
        else:   begin
                message = 'ERROR - Illegal value of BITPIX (= ' +  $
                               strtrim(bitpix,2) + ') in FITS header'
                goto,error_exit
                end
        endcase

        ndata = long64( axis[0] )
        bytes_per_word = (abs(bitpix)/8)
        if naxis gt 1 then $
           for i=2,naxis do ndata = ndata*axis[i-1]
        nbytes_per_group = bytes_per_word * (pcount + ndata)
        nbytes = (gcount>1) * nbytes_per_group
        nwords = nbytes / bytes_per_word
;
; starting data position
;
        position = fcb.start_data[enum]
;
; find correct group
;
        if last eq 0 then begin
                if group ge (gcount>1) then begin
                        message='INVALID group number specified'
                        goto,error_exit
                end
                position = position + group * nbytes_per_group 
        end
;
; read group parameters
;
        if (enum eq 0) and (fcb.random_groups eq 1) and (pcount gt 0) and $
           (last eq 0) then begin
            if n_params(0) gt 3 then begin
                group_par = make_array( dim = [pcount], type = idl_type, /nozero)
                point_lun,fcb.unit,position
                readu,fcb.unit,group_par
            endif
            position = position + pcount * bytes_per_word
        endif
;
; create data array
;
        if last gt 0 then begin
;
; user specified first and last
;
                if (first lt 0) or (last le 1) or (first gt last) or $
                   (last gt nwords-1) then begin
                        message = 'INVALID value for parameters FIRST & LAST'
                        goto,error_exit
                endif
                data = make_array(dim = [last-first+1], type=idl_type, /nozero)
                position = position + first * bytes_per_word
            endif else begin
;
; full array
;
                if ndata eq 0 then begin
                        data = 0
                        goto,done
                endif 
                if naxis gt 8 then begin
                        message = 'Maximum value of NAXIS allowed is 8'
                        goto,error_exit
                endif
                data = make_array(dim = axis, type = idl_type, /nozero)
        endelse
;
; read array
;
        point_lun,fcb.unit,position
        readu,fcb.unit,data
        if not keyword_set(No_Unsigned) and (not data_only) then begin
        if unsgn_int then begin 
                data =  uint(data) - uint(32768) 
        endif else if unsgn_lng then begin 
                data = ulong(data) - ulong(2147483648)
        endif
	endif
;
; scale data if header was read and first and last not used.   Do a special
; check of an unsigned integer (BZERO = 2^15) or unsigned long (BZERO = 2^31) 
;
        if (data_only eq 0) and (last eq 0) and (noscale eq 0) then begin

                if bitpix lt 32 then begin      ;use real*4 for bitpix<32
                        bscale = float(bscale)
                        bzero = float(bzero)
                endif
                if bscale ne 1.0 then data = temporary(data)*bscale
                if bzero ne 0.0 then data = temporary(data) + bzero
 	endif
;
; done
;
done:   
        if fcbtype eq 7 then fits_close,fcb else file_or_fcb.last_extension=enum
        !err = 1
        return

;
; error exit
;
ioerror:
        message = !err_string
error_exit:
        if (fcbtype eq 7) and (N_elements(fcb) GT 0) then  $
                   fits_close,fcb, no_abort=no_abort
        !err = -1
        if keyword_set(no_abort) then return
        print,'FITS_READ ERROR: '+message
        retall
end