fxbreadm.pro

basic median filter simulation

                    COLDIM[I,0:COLNDIM[I]-2] = COLDIM[I,1:COLNDIM[I]-1]
                    COLDIM[I,COLNDIM[I]-1]   = 0
                    COLNDIM[I] = COLNDIM[I] - 1
                ENDIF ELSE BEGIN  ;; Case of a single row
                    COLNDIM[I] = 1L
                    COLDIM[I,0] = NROWS
                ENDELSE
            ENDIF

            ;; Byte offsets
            BOFF1[I] = BYTOFF[ICOL[I],ILUN]
            IF ICOL[I] EQ TFIELDS[ILUN]-1 THEN $
              BOFF2[I] = NAXIS1[ILUN]-1 $
            ELSE $
              BOFF2[I] = BYTOFF[ICOL[I]+1,ILUN]-1

            LOOP_END_DIMS:

        ENDFOR

;
;  Construct any virtual columns first
;
        WC = WHERE(FOUND EQ 1 AND VIRTUAL EQ 1, WCCOUNT)
        FOR I = 0L, WCCOUNT-1 DO BEGIN
            ;; If it's virtual, then the value only needs to be
            ;; replicated
            EXTCMD = COLNAMES[WC[I]]+'= REPLICATE(D'+COLNAMES[WC[I]]+',NROWS)'
            ;; Run the command that selects the data
            RESULT = EXECUTE(EXTCMD)
            IF RESULT EQ 0 THEN BEGIN
                MESSAGE = 'ERROR: Could not extract data (column '+$
                  STRTRIM(MYCOL[WC[I]],2)+')'
                IF N_ELEMENTS(ERRMSG) NE 0 THEN BEGIN
                    ERRMSG = MESSAGE
                    RETURN
                ENDIF ELSE MESSAGE, MESSAGE
            ENDIF
            OUTSTATUS[I] = 1
        ENDFOR


;  Skip to processing variable-length columns if all other columns are virtual
        WC = WHERE(FOUND EQ 1 AND VIRTUAL EQ 0, WCCOUNT)
        IF WCCOUNT EQ 0 THEN GOTO, PROC_CLEANUP

;  Create NANVALUES, the template to use when a NAN is found
        IF N_ELEMENTS(NANVALUE) GE NUMCOLS THEN BEGIN
            NANVALUES = NANVALUE[0:NUMCOLS-1]
        ENDIF ELSE IF N_ELEMENTS(NANVALUE) GT 0 THEN BEGIN
            NANVALUES = REPLICATE(NANVALUE[0], NUMCOLS)
            NANVALUES[0] = NANVALUE
            I = N_ELEMENTS(NANVALUE)
            IF I LT NUMCOLS THEN $
              NANVALUES[I:*] = NANVALUE[0]
        ENDIF

;
;  Find the position of the first byte of the data array in the file.
;
        OFFSET0 = NHEADER[ILUN] + NAXIS1[ILUN]*(ROW1-1)
        POS = 0L
        NROWS0 = NROWS
        J = 0L
        FIRST = 1
        ;; Here, we constrain the buffer to be at least 16 rows long.
        ;; If we fill up 32 kB with fewer than 16 rows, then there
        ;; must be a lot of (big) columns in this table.  It's
        ;; probably a candidate for using FXBREAD instead.
        BUFFROWS = LONG((BUFFERSIZE/NAXIS1[ILUN]) > 16L)
        IF BUFFERSIZE LE 0 THEN BUFFROWS = NROWS0

;
;  Loop through the data in chunks
;
        WHILE NROWS GT 0 DO BEGIN
        J = J + 1
        NR  = NROWS < BUFFROWS
        OFFSET1 = NAXIS1[ILUN]*POS

;
;  Proceed by reading a byte array from the input data file
;  FXBREADM reads all columns from the specified rows, and
;  sorts out the details of which bytes belong to which columns
;  in the next FOR loop.
;
        BB = BYTARR(NAXIS1[ILUN], NR)
        POINT_LUN, UNIT, OFFSET0+OFFSET1
        READU, UNIT, BB

;
;  Now select out the desired columns
;
        FOR I = 0, NUMCOLS-1 DO BEGIN
           
            ;; Extract the proper rows and columns
            IF NOT FOUND[I] THEN GOTO, LOOP_END_STORE
            IF VIRTUAL[I]   THEN GOTO, LOOP_END_STORE

            ;; Extract the data from the byte array and convert it
            ;; The inner CALL_FUNCTION is to one of the coercion
            ;; functions, such as FIX(), DOUBLE(), STRING(), etc.,
            ;; which is called with an offset to force a conversion
            ;; from bytes to the data type.
            ;; The outer CALL_FUNCTION is to REFORM(), which makes
            ;; sure that the data structure is correct.
            ;;
            DIMS = COLDIM[I,0:COLNDIM[I]-1]
            PERROW = ROUND(PRODUCT(DIMS)/NROWS0)
            
            IF N_ELEMENTS(NANVALUES) GT 0 THEN $
              EXTRA={NANVALUE: NANVALUES[I]}

            FXBREADM_CONV, BB[BOFF1[I]:BOFF2[I], *], DD, COLTYPE[I], PERROW, NR,$
              NOIEEE=KEYWORD_SET(NOIEEE), NOSCALE=KEYWORD_SET(NOSCALE), $
              TZERO=TZERO[ICOL[I], ILUN], TSCAL=TSCAL[ICOL[I], ILUN], $
              VARICOL=VARICOL[I], DEFAULT_FLOAT=DEFAULT_FLOAT, _EXTRA=EXTRA

            ;; Initialize the output variable on the first chunk
            IF FIRST THEN BEGIN
                SZ = SIZE(DD)
                ;; NOTE: type could have changed if TSCAL/TZERO were used
                COLTYPEI = SZ(SZ[0]+1)  
                RESULT = EXECUTE(COLNAMES[I]+' = 0')
                RESULT = EXECUTE(COLNAMES[I]+' = '+$
                                 'MAKE_ARRAY(PERROW, NROWS0, TYPE=COLTYPEI)')
                RESULT = EXECUTE(COLNAMES[I]+' = '+$
                         'REFORM('+COLNAMES[I]+', PERROW, NROWS0,/OVERWRITE)')
            ENDIF

            ;; Finally, store this in the output variable
            RESULT = EXECUTE(COLNAMES[I]+'[0,POS] = DD')
            DD = 0
            IF RESULT EQ 0 THEN BEGIN
                MESSAGE = 'ERROR: Could not compose output data '+COLNAMES[I]
                IF N_ELEMENTS(ERRMSG) NE 0 THEN BEGIN
                    ERRMSG = MESSAGE
                    RETURN
                ENDIF ELSE MESSAGE, MESSAGE
            ENDIF

            OUTSTATUS[I] = 1

            LOOP_END_STORE:
        ENDFOR

        FIRST = 0
        NROWS = NROWS - NR
        POS   = POS + NR
        ENDWHILE

;
;  Read the variable-length columns from the heap.  Adjacent data are
;  coalesced into one read operation.  Note: this technique is thus
;  optimal for extensions with only one variable-length column.  If
;  there are more than one then coalescence will not occur.
;

        ;; Width of the various data types in bytes
        WIDARR = [0L, 1L, 2L, 4L, 4L, 8L, 8L, 1L, 0L,16L, 0L]
        WV = WHERE(OUTSTATUS EQ 1 AND VARICOL EQ 1, WVCOUNT)
        FOR J = 0, WVCOUNT-1 DO BEGIN
            I = WV[J]
            RESULT = EXECUTE('PDATA = '+COLNAMES[I])
            NVALS = PDATA[0,*]          ;; Number of values in each row
            NTOT  = ROUND(TOTAL(NVALS)) ;; Total number of values
            IF NTOT EQ 0 THEN BEGIN
                DD = {N_ELEMENTS: 0L, N_ROWS: NROWS0, $
                      INDICES: LONARR(NROWS0+1), DATA: 0L}
                GOTO, FILL_VARICOL
            ENDIF

            ;; Compute the width in bytes of the data value
            TYPE = IDLTYPE[ICOL[I], ILUN]
            WID = WIDARR[TYPE < 10]
            IF WID EQ 0 THEN BEGIN
                OUTSTATUS[I] = 0
                MESSAGE = 'ERROR: Column '+COLNAMES[I]+' has unknown data type'
                IF N_ELEMENTS(ERRMSG) NE 0 THEN BEGIN
                        ERRMSG = MESSAGE
                        RETURN
                END ELSE MESSAGE, MESSAGE
            ENDIF

            ;; Coalesce the data pointers
            BOFF1 = PDATA[1,*]
            BOFF2 = BOFF1 + NVALS*WID
            WH = WHERE(BOFF1[1:*] NE BOFF2, CT)
            IF CT GT 0 THEN BI = [-1L, WH, N_ELEMENTS(BOFF1)-1] $
            ELSE            BI = [-1L,     N_ELEMENTS(BOFF1)-1]
            CT = CT + 1

            ;; Create the output array
            BC = BOFF2[BI[1:*]] - BOFF1[BI[0:CT-1]+1] ;; Byte count
            NB = ROUND(TOTAL(BC))                     ;; Total # bytes
            BB = BYTARR(NB)                           ;; Byte array

            ;; Initialize the counter variables used in the read-loop
            CC = 0L & CC1 = 0L & K = 0L
            BUFFROWS = ROUND(BUFFERSIZE/WID) > 128L
            BASE = NHEADER[ILUN]+HEAP[ILUN]

            ;; Read data from file
            WHILE CC LT NB DO BEGIN
                NB1 = (BC[K]-CC1) < BUFFROWS
                BB1 = BYTARR(NB1)

                POINT_LUN, UNIT, BASE+BOFF1[BI[K]+1]+CC1
                READU, UNIT, BB1
                BB[CC] = TEMPORARY(BB1)

                CC  = CC  + NB1
                CC1 = CC1 + NB1
                IF CC1 EQ BC[K] THEN BEGIN
                    K = K + 1
                    CC1 = 0L
                ENDIF
            ENDWHILE

            ;; Convert the data
            IF N_ELEMENTS(NANVALUES) GT 0 THEN $
              EXTRA={NANVALUE: NANVALUES[I]}

            FXBREADM_CONV, BB, DD, TYPE, NTOT, 1L, $
              NOIEEE=KEYWORD_SET(NOIEEE), NOSCALE=KEYWORD_SET(NOSCALE), $
              TZERO=TZERO[ICOL[I], ILUN], TSCAL=TSCAL[ICOL[I], ILUN], $
              DEFAULT_FLOAT=DEFAULT_FLOAT, _EXTRA=EXTRA
            
            ;; Ensure the correct dimensions, now that we know them
            COLNDIM[I] = 1
            COLDIM[I,0] = NTOT
            
            ;; Construct the indices; unfortunately we need to make an
            ;; accumulant with a FOR loop
            INDICES = LONARR(NROWS0+1)
            FOR K = 1L, NROWS0 DO $
              INDICES[K] = INDICES[K-1] + NVALS[K-1]

            ;; Construct a structure with additional data
            DD = {N_ELEMENTS: NTOT, N_ROWS: NROWS0, TYPE: TYPE, $
                  INDICES: INDICES, DATA: TEMPORARY(DD)}

            FILL_VARICOL:
            RESULT = EXECUTE(COLNAMES[I] +' = TEMPORARY(DD)')
        ENDFOR

;
;  Compose the output columns, which might need reforming
;
        FOR I = 0, NUMCOLS-1 DO BEGIN
            IF OUTSTATUS[I] NE 1 THEN GOTO, LOOP_END_FINAL

            ;; Extract the dimensions and name of the column data
            DIMS = COLDIM[I,0:COLNDIM[I]-1]
            NEL  = PRODUCT(DIMS)
            CNAME = COLNAMES[I]
            IF VARICOL[I] THEN CNAME = CNAME + '.DATA'

            ;; Compose the reforming part
            IF NEL EQ 1 THEN $
              CMD = CNAME+'[0]' $
            ELSE $
              CMD = 'REFORM(TEMPORARY('+CNAME+'),DIMS,/OVERWRITE)'

            ;; Variable-length columns return extra information
            IF VARICOL[I] THEN BEGIN
                CMD = ('{VARICOL:    1,'+$
                       ' N_ELEMENTS: '+COLNAMES[I]+'.N_ELEMENTS, '+$
                       ' TYPE:       '+COLNAMES[I]+'.TYPE, '+$
                       ' N_ROWS:     '+COLNAMES[I]+'.N_ROWS, '+$
                       ' INDICES:    '+COLNAMES[I]+'.INDICES, '+$
                       ' DATA:       '+CMD+'}')
            ENDIF

            ;; Assign to pointer, or re-assign to column
            IF PASS EQ 'ARGUMENT' THEN $
              CMD = COLNAMES[I]+' = ' + CMD $
            ELSE IF PASS EQ 'POINTER' THEN $
              CMD = '*(POINTERS[I]) = ' + CMD 

            RESULT = EXECUTE(CMD)
            LOOP_END_FINAL:
        ENDFOR

        PROC_CLEANUP:
;
        IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = ''
        RETURN
        
        END