mrdfits.pro

basic median filter simulation

; First do ASCII and Binary tables.    We need to create a new structure 
; because scaling will change the tag data types.

    if type ne 0 then begin
        
        if type eq 1 then begin
            if keyword_set(dscale) then begin
                fvalues[w] = '0.0d0'
            endif else begin
                fvalues[w] = '0.0'
            endelse
        endif else if type eq 2 then begin

            if keyword_set(dscale) then begin
                sclr = '0.d0'
                vc = 'dblarr'
            endif else begin
                sclr = '0.0'
                vc = 'fltarr'
            endelse
                
           for i=0, Nw-1 do begin
                col = w[i]
                sz = size(table[0].(col),/str)

		; Handle pointer columns
		if sz.type eq 10 then begin
		    fvalues[col] = 'ptr_new()'

		; Scalar columns
		endif else if sz.N_dimensions eq 0 then begin
                    fvalues[col] = sclr

		; Vectors
                endif else begin
		    dim = sz.dimensions[0:sz.N_dimensions-1]
                    fvalues[col] = vc + $
		      '(' + strjoin(strtrim(dim,2),',') + ')'
		    
                endelse
            endfor
        endif

        tabx = mrd_struct(fnames, fvalues, nrec, structyp=structyp, silent=silent )

; First copy the unscaled columns indexed by ww.     This is actually more 
; efficient than using STRUCT_ASSIGN since the tag names are all identical,
; so STRUCT_ASSIGN would copy everything (scaled and unscaled).
 	    
       for i=0, Nww - 1 do tabx.(ww[i]) = table.(ww[i])
       
; Now copy the scaled items indexed by w after applying the scaling.        
       
        for i=0, Nw - 1 do begin	    
 		
		dtype = size(tabx.(w[i]),/type)
		if dtype eq 10 then $
		    mrd_ptrscale, table.(w[i]), scales[w[i]], offsets[w[i]]
		
                tabx.(w[i]) = table.(w[i])*scales[w[i]] + offsets[w[i]]
		
            istr = strtrim(w[i]+1,2)
            fxaddpar, header, 'TSCAL'+istr, 1.0, ' Set by MRD_SCALE'
            fxaddpar, header, 'TZERO'+istr, 0.0, ' Set by MRD_SCALE'
	    
        endfor

        table = temporary(tabx)   ;Remove original structure from memory
    endif else begin
    ; Now process images and random groups.

        sz = size(table[0])
        if sz[sz[0]+1] ne 8 then begin
            ; Not a structure so we just have an array of data.
            if keyword_set(dscale) then begin
                table = temporary(table)*scales[0]+offsets[0]
            endif else begin
                table = tempoary(table)*float(scales[0]) + float(offsets[0])
            endelse
            fxaddpar, header, 'BSCALE', 1.0, 'Set by MRD_SCALE'
            fxaddpar, header, 'BZERO', 0.0, 'Set by MRD_SCALE'

        endif else begin
            ; Random groups.  Get the number of parameters by looking
            ; at the first element in the table.
            nparam = n_elements(table[0].(0))
            if keyword_set(dscale) then typ = 'dbl' else typ='flt'
            s1 = typ+'arr('+string(nparam)+')'
            ngr = n_elements(table)
            sz = size(table[0].(1))
            if sz[0] eq 0 then dims = [1] else dims=sz[1:sz[0]]
            s2 = typ + 'arr('
            for i=0, n_elements(dims)-1 do begin 
                if i ne 0 then s2 = s2+ ','
                s2 = s2+string(dims[i])
            endfor
            s2 = s2+')'
            tabx = mrd_struct(['params', 'array'],[s1,s2],ngr, silent=silent)

            for i=0, nparam-1 do begin
                istr = strcompress(string(i+1),/remo)
                fxaddpar, header, 'PSCAL'+istr, 1.0, 'Added by MRD_SCALE'
                fxaddpar, header, 'PZERO'+istr, 0.0, 'Added by MRD_SCALE'
                tabx.(0)[i] = table.(0)[i]*scales[i]+offsets[i]
            endfor
	    
            tabx.(1) = table.(1)*scales[nparam] + offsets[nparam]
            fxaddpar, header, 'BSCALE', 1.0, 'Added by MRD_SCALE'
            fxaddpar, header, 'BZERO', 0.0, 'Added by MRD_SCALE'
            table = temporary(tabx)
        endelse
    endelse

end

; Read a variable length column into a pointer array.
pro mrd_varcolumn, vtype, array, heap, off, siz

    ; Guaranteed to have at least one non-zero length column
    w   = where(siz gt 0)
    nw  = n_elements(w)
    
    if vtype eq 'X' then siz = 1 + (siz-1)/8
    
    siz = siz[w]
    off = off[w]

    unsigned = 0
    if vtype eq '1' then begin
	unsigned = 12
    endif else if vtype eq '2' then begin
	unsigned = 13
    endif else if vtype eq '3' then begin
	unsigned = 15;
    endif
    unsigned = mrd_unsigned_offset(unsigned)
    

    for j=0, nw-1 do begin

        case vtype of

            'L': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
            'X': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
            'B': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
	
            'I': array[w[j]] = ptr_new(  fix(heap, off[j], siz[j]) )
            'J': array[w[j]] = ptr_new(  long(heap, off[j], siz[j]) )
            'K': array[w[j]] = ptr_new(  long64(heap, off[j], siz[j]) )
			  
            'E': array[w[j]] = ptr_new(  float(heap, off[j], siz[j]) )
            'D': array[w[j]] = ptr_new(  double(heap, off[j], siz[j]) )
			  
            'C': array[w[j]] = ptr_new(  complex(heap, off[j], siz[j]) )
            'M': array[w[j]] = ptr_new(  dcomplex(heap, off[j], siz[j]) )
			   
            '1': array[w[j]] = ptr_new(  uint(heap, off[j], siz[j]) )
            '2': array[w[j]] = ptr_new(  ulong(heap, off[j], siz[j]) )
            '3': array[w[j]] = ptr_new(  ulong64(heap, off[j], siz[j]) )
      
        endcase

	; Fix endianness.
        if vtype ne 'B' and vtype ne 'X' and vtype ne 'L' then begin
	    swap_endian_inplace, *array[w[j]],/swap_if_little
        endif

	; Scale unsigneds.
	if unsigned gt 0 then *array[w[j]] = *array[w[j]] - unsigned
	
    endfor
end

; Read a variable length column into a fixed length array.
pro mrd_fixcolumn, vtype, array, heap, off, siz

    w   = where(siz gt 0)
    if w[0] eq -1 then return
    
    nw  = n_elements(w)

    
    if vtype eq 'X' then siz = 1 + (siz-1)/8
    
    siz = siz[w]
    off = off[w]

    for j=0, nw-1 do begin
        case vtype of
            'L': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j])  
            'X': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j])
            'B': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j]) 
	
            'I': array[0:siz[j]-1,w[j]] = fix(heap, off[j], siz[j]) 
            'J': array[0:siz[j]-1,w[j]] = long(heap, off[j], siz[j]) 
            'K': array[0:siz[j]-1,w[j]] = long64(heap, off[j], siz[j]) 
			  
            'E': begin                  ;Delay conversion until after byteswapping to avoid possible math overflow   Feb 2005
	         temp = heap[off[j]: off[j] + 4*siz[j]-1 ]
		 byteorder, temp, /LSWAP, /SWAP_IF_LITTLE	 
	         array[0:siz[j]-1,w[j]] = float(temp,0,siz[j]) 
                 end			  
           'D': begin 
	         temp = heap[off[j]: off[j] + 8*siz[j]-1 ]
		 byteorder, temp, /L64SWAP, /SWAP_IF_LITTLE		 
	         array[0:siz[j]-1,w[j]] = double(temp,0,siz[j]) 
                 end			  
            'C': array[0:siz[j]-1,w[j]] = complex(heap, off[j], siz[j]) 
            'M': array[0:siz[j]-1,w[j]] = dcomplex(heap, off[j], siz[j]) 
			   
            'A': array[w[j]] = string(byte(heap,off[j],siz[j])) 

            '1': array[0:siz[j]-1,w[j]] = uint(heap, off[j], siz[j]) 
            '2': array[0:siz[j]-1,w[j]] = ulong(heap, off[j], siz[j])
            '3': array[0:siz[j]-1,w[j]] = ulong64(heap, off[j], siz[j])
      
        endcase

    endfor

    ; Fix endianness
    if (vtype ne 'A') and (vtype ne 'B') and (vtype ne 'X') and (vtype ne 'L')  and $
       (vtype NE 'D')  and (vtype NE 'E') then begin
	swap_endian_inplace, array, /swap_if_little
    endif

    ; Scale unsigned data
    unsigned = 0
    if vtype eq '1' then begin
	unsigned = 12
    endif else if vtype eq '2' then begin
	unsigned = 13
    endif else if vtype eq '3' then begin
	unsigned = 15;
    endif
    
    if unsigned gt 0 then begin
        unsigned = mrd_unsigned_offset(unsigned)
    endif
    
    if unsigned gt 0 then begin
        for j=0, nw-1 do begin
            array[0:siz[j]-1,w[j]] = array[0:siz[j]-1,w[j]] - unsigned
	endfor
    endif


end
		
; Read the heap area to get the actual values of variable 
; length arrays. 
pro mrd_read_heap, unit, header, range, fnames, fvalues, vcls, vtpes, table, $ 
   structyp, scaling, scales, offsets, status, silent=silent,                $
   columns=columns, rows = rows, pointer_var=pointer_var, fixed_var=fixed_var

    ; 
    ; Unit:         FITS unit number. 
    ; header:       FITS header. 
    ; fnames:       Column names. 
    ; fvalues:      Column values. 
    ; vcols:        Column numbers of variable length columns. 
    ; vtypes:       Actual types of variable length columns 
    ; table:        Table of data from standard data area, on output 
    ;               contains the variable length data. 
    ; structyp:     Structure name. 
    ; scaling:      Is there going to be scaling of the data?
    ; status:       Set to -1 if an error occurs.
    ;
    typstr = 'LXBIJKAEDCM123' 
    prefix = ['bytarr(', 'bytarr(', 'bytarr(', 'intarr(',     $ 
              'lonarr(', 'lon64arr(', 'string(bytarr(', 'fltarr(',         $ 
              'dblarr(', 'cmplxarr(', 'dblarr(2,',            $
	      'uintarr(', 'ulonarr(', 'ulon64arr(']
    
    status = 0 

    ; Convert from a list of indicators of whether a column is variable
    ; length to pointers to only the variable columns.

    vcols = where(vcls eq 1)
    vtypes = vtpes[vcols]

    nv = n_elements(vcols) 
 
    ; Find the beginning of the heap area. 
 
    heapoff = long64(fxpar(header, 'THEAP')) 
    sz = fxpar(header, 'NAXIS1')*fxpar(header, 'NAXIS2')
    
    if heapoff ne 0 and heapoff lt sz then begin 
        print, 'MRDFITS: ERROR Heap begins within data area' 
        status = -1 
        return 
    endif

    ; Skip to beginning.
    if (heapoff > sz) then begin
        mrd_skip, unit, heapoff-sz
    endif
 
    ; Get the size of the heap. 
    pc = long64(fxpar(header, 'PCOUNT')) 
    if heapoff eq 0 then heapoff = sz 
    hpsiz = pc - (heapoff-sz) 
 
    if (hpsiz gt 0) then heap = bytarr(hpsiz) 
 
 
    ; Read in the heap 
    readu, unit, heap

    ; Skip to the end of the data area.
    skipB = 2880 - (sz+pc) mod 2880
    if skipB ne 2880 then begin
        mrd_skip, unit, skipB
    endif
 
    ; Find the maximum dimensions of the arrays. 
    ; 
    ; Note that the variable length column currently has fields which 
    ; are I*4 2-element arrays where the first element is the 
    ; length of the field on the current row and the second is the 
    ; offset into the heap. 

    vdims = lonarr(nv)
    for i=0, nv-1 do begin 
        col = vcols[i]
        curr_col = table.(col)
        vdims[i] = max(curr_col[0,*])
        w = where(curr_col[0,*] ne vdims[i])
        if w[0] ne -1 then begin
            if n_elements(lencols) eq 0 then begin
                lencols = [col]
            endif else begin
                lencols=[lencols,col]
            endelse
        endif

        if vtypes[i] eq 'X' then vdims[i]=(vdims[i]+7)/8 
        ind = strpos(typstr, vtypes[i])
    
        ; Note in the following that we ensure that the array is
        ; at least one element long.

        fvalues[col] = prefix[ind] + string((vdims[i] > 1)) + ')'
        if vtypes[i] eq 'A' then fvalues[col] = fvalues[col] + ')' 
    
    endfor 
 
    nfld = n_elements(fnames) 

    ; Get rid of columns which have no actual data. 
    w= intarr(nfld) 
    w[*] = 1
    corres = indgen(nfld)

    
    ; Should we get rid of empty columns?
    delete = 1
    if keyword_set(pointer_var) then delete = pointer_var eq 1

    if delete then begin
	
        ww = where(vdims eq 0) 
        if ww[0] ne -1 then  begin
            w[vcols[ww]] = 0
            if not keyword_s