mrdfits.pro

来自「basic median filter simulation」· PRO 代码 · 共 1,852 行 · 第 1/5 页

    if error ne 0 then begin
     catch,/cancel
     status=-2
     return
    endif

    ; If necessary skip to beginning of desired data. 

    if range[0] gt 0 then mrd_skip, unit, range[0]*rsize

    status=-2
    if rsize eq 0 then return 

    on_ioerror,done
    readu, unit, table

    if N_elements(rows) GT 0 then begin
    row1 = rows- range[0]
    case size(table,/n_dimen) of 
    1: table = table[row1]
    2: table = table[*,row1]
    3: table = table[*,*,row1]
    4: table = table[*,*,*,row1]
    5: table = table[*,*,*,*,row1]
    6: table = table[*,*,*,*,*,row1]
    7: table = table[*,*,*,*,*,*,row1]
    8: table = table[*,*,*,*,*,*,*,row1]
    else: begin 
          print,'MRDFITS: Subscripted image must be between 1 and 8 dimensions'
          status = -1
          return
          end
    endcase
    endif

    ; Skip to the end of the data

    skipB = 2880 - (maxd*rsize) mod 2880
    if skipB eq 2880 then skipB = 0

    if range[1] lt maxd-1 then begin
        skipB = skipB + (maxd-range[1]-1)*rsize
    endif

    mrd_skip, unit, skipB
    if unixpipe then swap_endian_inplace, table,/swap_if_little

    ; Fix offset for unsigned data
    type = mrd_unsignedtype(table)
    if type gt 0 then begin
	table = table - mrd_unsigned_offset(type)
    endif
    
    status=0
    done:

;-- probably an EOF 

    if status ne 0 then begin 
          message,!ERROR_STATE.MSG,/CON
         free_lun,unit
    endif

    return
end 

; Truncate superfluous axes.

pro mrd_axes_trunc,naxis, dims, silent

    mysilent = silent
    for i=naxis-1,1,-1 do begin 

        if dims[i] eq 1 then begin
            if not mysilent then begin
                print, 'MRDFITS: Truncating unused dimensions'
                mysilent = 1
            endif
            dims = dims[0:i-1] 
            naxis = naxis - 1 
        
        endif else return
     
    endfor 
 
    return
end

; Define structure/array to hold a FITS image. 
pro mrd_image, header, range, maxd, rsize, table, scales, offsets, scaling, $
  status, silent=silent, unsigned=unsigned, rows = rows
 
    ; 
    ; Header                FITS header for table. 
    ; Range                 Range of data to be retrieved. 
    ; Rsize                 Size of a row or group. 
    ; Table                 Structure to be defined. 
    ; Status                Return status
    ; Silent=silent         Suppress info messages?
 
    table = 0

    ; type    0         1           2         3         4         5  6  7  8  9 10 11        12         13          14          15
    lens =  [ 0,        1,          2,        4,        4,        8, 0, 0, 0, 0, 0, 0,        2,         4,          8,          8] 
    typstrs=['',   'Byte',    'Int*2',  'Int*4', 'Real*4', 'Real*8','','','','','','', 'UInt*2',  'Uint*4',    'Int*8',    'Uint*8']
    typarr= ['', 'bytarr',   'intarr', 'lonarr', 'fltarr', 'dblarr','','','','','','','uintarr', 'ulonarr', 'lon64arr', 'ulon64arr'] 
 
    status = 0 
 
 
    naxis = fxpar(header, 'NAXIS') 
    bitpix= fxpar(header, 'BITPIX') 
    if naxis gt 0 then begin 
          dims = long64(fxpar(header, 'NAXIS*', Count = N_axis)) 
          if N_axis GT naxis then begin
; Check if extra NAXISn keywords are present (though this is not legal FITS)
                   nextra = N_axis - naxis
                   dim_extra = dims[naxis:N_axis-1]
                   if total(dim_extra) EQ nextra then $
                        dims = dims[0:naxis-1] else $
                   message,'ERROR - NAXIS = ' + strtrim(naxis,2) +  $
                          ' but NAXIS' + strtrim(N_axis,2) + ' keyword present'
          endif
   endif else dims = 0
    
    gcount = fxpar(header, 'GCOUNT') 
    pcount = fxpar(header, 'PCOUNT')
    isgroup = fxpar(header, 'GROUPS')
    gcount = long(gcount)

    xscale = fxpar(header, 'BSCALE', count=cnt)
    if cnt eq 0 then xscale = 1      ;Corrected 06/29/06
    
    xunsigned = mrd_chkunsigned(bitpix,  xscale, $
				fxpar(header, 'BZERO'), unsigned=unsigned)
    ; Note that type is one less than the type signifier returned in the size call.
    type = -1
    
    if not xunsigned then begin 
 
        if bitpix eq 8        then type = 1     $ 
        else if bitpix eq  16 then type = 2     $ 
        else if bitpix eq  32 then type = 3     $ 
        else if bitpix eq -32 then type = 4     $ 
        else if bitpix eq -64 then type = 5     $
        else if bitpix eq  64 then type = 14

    endif else begin

	if bitpix eq 16       then type = 12     $
	else if bitpix eq  32 then type = 13     $
	else if bitpix eq  64 then type = 15
	
    endelse

    if type eq -1 then begin
	print,'MRDFITS: Error: Invalid BITPIX: '+strtrim(bitpix)
	table = 0
	return
    endif

    ; Note that for random groups data we must ignore the first NAXISn keyword. 
    if isgroup GT 0  then begin 


        range[0] = range[0] > 0
        if (range[1] eq -1) then begin
            range[1] = gcount-1
        endif else begin
            range[1] = range[1] < gcount - 1
        endelse
	
	maxd = gcount
        
        if (n_elements(dims) gt 1) then begin
            dims = dims[1:*]
            naxis = naxis-1
        endif else begin
            print, 'MRDFITS: Warning: No data specified for group data.'
            dims = [0]
            naxis = 0
        endelse
        
        ; The last entry is the scaling for the sample data.
        
        if (pcount gt 0) then begin
            scales  = dblarr(pcount+1)
            offsets = dblarr(pcount+1)
        endif
        
        values = strarr(2)
        
        
        mrd_axes_trunc, naxis, dims, keyword_set(silent)
        
        values[0] = typarr[type] + "("+string(pcount)+")" 
        rsize = dims[0] 
        sarr = "(" + strcompress(string(dims[0]), /remo )
         
        for i=1, naxis-1 do begin
	    
            sarr = sarr + "," + strcompress(string(dims[i]),/remo)
            rsize = rsize*dims[i]
	    
        endfor 
         
        sarr = sarr + ")"

        if not keyword_set(silent) then print,'MRDFITS--Image with groups:', $
          ' Ngroup=',strcompress(string(gcount)),' Npar=',                   $
          strcompress(string(pcount),/remo), ' Group=', sarr, '  Type=',typstrs[type]

        sarr = typarr[type] + sarr
        values[1] = sarr 
        rsize = (rsize + pcount)*lens[type] 
         
        table = mrd_struct(['params','array'], values, range[1]-range[0]+1, $
                           silent=silent)

	if xunsigned then begin
	    fxaddpar,header, 'BZERO', 0, 'Reset by MRDFITS v'+mrd_version()
	endif
           

        for i=0, pcount-1 do begin
	    
            istr = strcompress(string(i+1),/remo)
	    
            scales[i] = fxpar(header, 'PSCAL'+istr)
            if scales[i] eq 0.0d0 then scales[i] =1.0d0
	    
            offsets[i] = fxpar(header, 'PZERO'+istr)
	    
            scales[pcount] = fxpar(header, 'BSCALE')
            if scales[pcount] eq 0.0d0 then scales[pcount] = 1.0d0
            offsets[pcount] = fxpar(header, 'BZERO')
	    
        endfor

        if scaling then begin
            w = where(scales ne 1.0d0 or offsets ne 0.0d0)
            if w[0] eq -1 then scaling = 0
        endif
        
    endif else begin 
 
        if naxis eq 0 then begin
	
            rsize = 0 
            table = 0
            if not keyword_set(silent) then begin
                print, 'MRDFITS: Null image, NAXIS=0'
            endif
            return
	    
        endif 
         
        if gcount gt 1 then begin 
            dims = [dims, gcount] 
            naxis = naxis + 1 
        endif 
         
        mrd_axes_trunc, naxis, dims, keyword_set(silent)

                
        maxd = dims[naxis-1] 
         
        if range[0] ne -1 then begin 
            range[0] = range[0]<(maxd-1) 
            range[1] = range[1]<(maxd-1) 
        endif else begin 
            range[0] = 0 
            range[1] = maxd - 1 
        endelse 

        Nlast = dims[naxis-1]   
        dims[naxis-1] = range[1]-range[0]+1
        pdims = dims
        if N_elements(rows) GT 0 then begin
             if max(rows) GE Nlast then begin 
               print, 'MRDFITS: Row numbers must be between 0 and ' + $
                      strtrim(Nlast-1,2)
               status = -1 & rsize = 0
               return
             endif
             pdims[naxis-1] = N_elements(rows)
        endif 
 
        if not keyword_set(silent) then begin
            str = '('
            for i=0, naxis-1 do begin
                if i ne 0 then str = str + ','
                str = str + strcompress(string(pdims[i]),/remo)
            endfor
            str = str+')'
            print, 'MRDFITS: Image array ',str, '  Type=', typstrs[type]
        endif
         
        rsize = 1
	
        if naxis gt 1 then for i=0, naxis - 2 do rsize=rsize*dims[i] 
        rsize = rsize*lens[type] 
        sz = lonarr(naxis+3) 
        sz[0] = naxis 
        sz[1:naxis] = dims 
        nele = 1l 
         
        for i=0, naxis-1 do begin 
            nele = nele*dims[i] 
        endfor 
         
        sz[naxis+1] = type   
        sz[naxis+2] = nele 
         
        if nele gt 0 then  begin
            table = make_array(size=sz) 
        endif else begin 
            table = 0 
        endelse 
        
        scales = dblarr(1)
        offsets = dblarr(1)

	if xunsigned then begin
	    fxaddpar,header, 'BZERO', 0, 'Updated by MRDFITS v'+mrd_version()
	endif
	
        scales[0] = fxpar(header, 'BSCALE')
        offsets[0] = fxpar(header, 'BZERO')
	
        if scales[0] eq 0.0d0 then scales[0] = 1.0d0
        if scaling and scales[0] eq 1.0d0 and offsets[0] eq 0.0d0 then scaling = 0
    endelse 
         
    status = 0 
    return 
 
end

; Scale an array of pointers
pro mrd_ptrscale, array, scale, offset

    for i=0, n_elements(array)-1 do begin
        if ptr_valid(array[i]) then begin
	    array[i] = ptr_new(*array[i] * scale + offset)
	endif
    endfor
end

; Scale a FITS array or table.
pro mrd_scale, type, scales, offsets, table, header,  $
               fnames, fvalues, nrec, dscale = dscale, structyp=structyp, silent=silent
    ;
    ; Type:         FITS file type, 0=image/primary array
    ;                               1=ASCII table
    ;                               2=Binary table
    ;
    ; scales:       An array of scaling info
    ; offsets:      An array of offset information
    ; table:        The FITS data.
    ; header:       The FITS header.
    ; dscale:       Should data be scaled to R*8?
    ; fnames:       Names of table columns.
    ; fvalues:      Values of table columns.
    ; nrec:         Number of records used.
    ; structyp:     Structure name.
 
    w = where( (scales ne 1.d0  or offsets ne 0.d0), Nw, $
                complement=ww, Ncomplement = Nww)
		
    if Nw EQ 0 then return    ;No tags require scaling?