fxbreadm.pro

basic median filter simulation

;+
; NAME:
;       FXBREADM
; PURPOSE: 
;       Read multiple columns/rows from a disk FITS binary table file.
; EXPLANATION : 
;       A call to FXBREADM will read data from multiple rows and
;       multiple columns in a single procedure call.  Up to forty-nine
;       columns may be read in a single pass; the number of rows is
;       limited essentially by available memory.  The file should have
;       already been opened with FXBOPEN.  FXBREADM optimizes reading
;       multiple columns by first reading a large chunk of data from
;       the FITS file directly, and then slicing the data into columns
;       within memory.  FXBREADM can read variable-length arrays (see
;       below).
;
;       The number of columns is limited to 49 if data are passed by
;       positional argument.  However, this limitation can be overcome
;       by having FXBREADM return the data in an array of pointers.
;       The user should set the PASS_METHOD keyword to 'POINTER', and an 
;       array of pointers to the data will be returned in the POINTERS keyword.
;       The  user is responsible for freeing the pointers; however,
;       FXBREADM will reuse any pointers  passed into the procedure, and 
;       hence any pointed-to data will be destroyed.
;
;       FXBREADM can also read variable-length columns from FITS
;       binary tables.  Since such data is not of a fixed size, it is
;       returned as a structure.  The structure has the following
;       elements:
;
;              VARICOL:    ;; Flag: variable length column (= 1)
;              N_ELEMENTS: ;; Total number of elements returned
;              TYPE:       ;; IDL data type code (integer)
;              N_ROWS:     ;; Number of rows read from table (integer)
;              INDICES:    ;; Indices of each row's data (integer array)
;              DATA:       ;; Raw data elements (variable type array)
;
;       In order to gain access to the Ith row's data, one should
;       examine DATA(INDICES(I):INDICES(I+1)-1), which is similar in
;       construct to the REVERSE_INDICES keyword of the HISTOGRAM
;       function.
;
; CALLING SEQUENCE: 
;       FXBREADM, UNIT, COL, DATA1, [ DATA2, ... DATA48, ROW=, BUFFERSIZE = ]
;           /NOIEEE, /NOSCALE, /VIRTUAL, NANVALUE=, PASS_METHOD = POINTERS=, 
;           ERRMSG = , WARNMSG = , STATUS = , /DEFAULT_FLOAT]
;
; INPUT PARAMETERS : 
;       UNIT    = Logical unit number corresponding to the file containing the
;                 binary table.
;       COL     = An array of columns in the binary table to read data
;                 from, either as character strings containing column
;                 labels (TTYPE), or as numerical column indices
;                 starting from column one.
; Outputs     : 
;       DATA1, DATA2...DATA48 = A named variable to accept the data values, one
;                 for each column.  The columns are stored in order of the
;                 list in COL.  If the read operation fails for a
;                 particular column, then the corresponding output Dn
;                 variable is not altered.  See the STATUS keyword.
;                 Ignored if PASS_METHOD is 'POINTER'.
;
; OPTIONAL INPUT KEYWORDS: 
;       DEFAULT_FLOAT = If set, then scaling with TSCAL/TZERO is done with
;                   floating point rather than double precision.
;       ROW     = Either row number in the binary table to read data from,
;                 starting from row one, or a two element array containing a
;                 range of row numbers to read.  If not passed, then the entire
;                 column is read in.
;       /NOIEEE = If set, then then IEEE floating point data will not
;                be converted to the host floating point format (and
;                this by definition implies NOSCALE).  The user is
;                responsible for their own floating point conversion.
;       /NOSCALE = If set, then the output data will not be scaled using the
;                 optional TSCAL and TZERO keywords in the FITS header.
;                 Default is to scale.
;       VIRTUAL = If set, and COL is passed as a name rather than a number,
;                 then if the program can't find a column with that name, it
;                 will then look for a keyword with that name in the header.
;                 Such a keyword would then act as a "virtual column", with the
;                 same value for every row.
;       DIMENSIONS = FXBREADM ignores this keyword.  It is here for
;	          compatibility only.
;       NANVALUE= Value signalling data dropout.  All points corresponding to
;                 IEEE NaN (not-a-number) are converted to this number.
;                 Ignored unless DATA is of type float, double-precision or
;                 complex.
;       PASS_METHOD = A scalar string indicating method of passing
;                 data from FXBREADM.  Either 'ARGUMENT' (indicating
;                 pass by positional argument), or 'POINTER' (indicating
;                 passing an array of pointers by the POINTERS
;                 keyword).
;                 Default: 'ARGUMENT'
;       POINTERS = If PASS_METHOD is 'POINTER' then an array of IDL
;                 pointers is returned in this keyword, one for each
;                 requested column.    Any pointers passed into FXBREADM will 
;                 have their pointed-to data destroyed.  Ultimately the
;                 user is responsible for deallocating pointers. 
;       BUFFERSIZE = Raw data are transferred from the file in chunks
;                 to conserve memory.  This is the size in bytes of
;                 each chunk.  If a value of zero is given, then all
;                 of the data are transferred in one pass.  Default is
;                 32768 (32 kB).
; OPTIONAL OUTPUT KEYWORDS:
;       ERRMSG  = If defined and passed, then any error messages will be
;                 returned to the user in this parameter rather than
;                 depending on the MESSAGE routine in IDL.  If no errors are
;                 encountered, then a null string is returned.  In order to
;                 use this feature, ERRMSG must be defined first, e.g.
;
;                       ERRMSG = ''
;                       FXBREAD, ERRMSG=ERRMSG, ...
;                       IF ERRMSG NE '' THEN ...
;       WARNMSG = Messages which are considered to be non-fatal
;                 "warnings" are returned in this output string.
;                 Note that if some but not all columns are
;                 unreadable, this is considered to be non-fatal.
;       STATUS  = An output array containing the status for each
;                 column read, 1 meaning success and 0 meaning failure.
;
; Calls       : 
;       IEEE_TO_HOST, FXPAR(), WHERENAN()
; Common      : 
;       Uses common block FXBINTABLE--see "fxbintable.pro" for more
;       information.
; Restrictions: 
;       The binary table file must have been opened with FXBOPEN.
;
;       The data must be consistent with the column definition in the binary
;       table header.
;
;       The row number must be consistent with the number of rows stored in the
;       binary table header.
;
;       Generaly speaking, FXBREADM will be faster than iterative
;       calls to FXBREAD when (a) a large number of columns is to be
;       read or (b) the size in bytes of each cell is small, so that
;       the overhead of the FOR loop in FXBREAD becomes significant.
;
; SIDE EFFECTS: 
;       If there are no elements to read in (the number of elements is zero),
;       then the program sets !ERR to -1, and DATA is unmodified.
;
; Category    : 
;       Data Handling, I/O, FITS, Generic.
; Prev. Hist. : 
;       C. Markwardt, based in concept on FXBREAD version 12 from
;                              IDLASTRO, but with significant and
;                              major changes to accomodate the
;                              multiple row/column technique.  Mostly
;                              the parameter checking and general data
;                              flow remain.
;       C. Markwardt, updated to read variable length arrays, and to
;                              pass columns by handle or pointer.
;                              20 Jun 2001
;       C. Markwardt, try to conserve memory when creating the arrays
;                              13 Oct 2001
;   Handle case of GE 50 columns, C. Markwardt, 18 Apr 2002
;   Handle case where TSCAL/TZERO changes type of column,
;       C. Markwardt, 23 Feb 2003
;   Fix bug in handling of FOUND and numeric columns, 
;       C. Markwardt 12 May 2003
;   Removed pre-V5.0 HANDLE options  W. Landsman July 2004
;   Fix bug when HANDLE options were removed, July 2004
;   Handle special cases of TSCAL/TZERO which emulate unsigned
;      integers, Oct 2003
;   Add DEFAULT_FLOAT keyword to select float values instead of double
;      for TSCAL'ed, June 2004
;  read 64bit integer columns, E. Hivon, Mar 2008
;
;
;-
;


;; This is a utility routine which converts the data from raw bytes to
;; IDL variables.
PRO FXBREADM_CONV, BB, DD, CTYPE, PERROW, NROWS, $
                   NOIEEE=NOIEEE, NOSCALE=NOSCALE, VARICOL=VARICOL, $
                   NANVALUE=NANVALUE, TZERO=TZERO, TSCAL=TSCAL, $
                   DEFAULT_FLOAT=DF

  COMMON FXBREADM_CONV_COMMON, DTYPENAMES
  IF N_ELEMENTS(DTYPENAMES) EQ 0 THEN $
    DTYPENAMES = [ '__BAD', 'BYTE', 'FIX', 'LONG', $
                   'FLOAT', 'DOUBLE', 'COMPLEX', 'STRING', $
                   '__BAD', 'DCOMPLEX', '__BAD', '__BAD', '__BAD', '__BAD', 'LONG64' ]
  
  TYPENAME = DTYPENAMES[CTYPE]

  IF CTYPE EQ 7 THEN BEGIN
      DD = STRING(TEMPORARY(BB))
  ENDIF ELSE BEGIN
      DD = CALL_FUNCTION(TYPENAME, TEMPORARY(BB), 0, PERROW*NROWS)
  ENDELSE
  IF N_ELEMENTS(DD) EQ 1 THEN DD = [DD]
  DD = REFORM(DD, PERROW, NROWS, /OVERWRITE)

  ;; Now perform any type-specific conversions, etc.
  COUNT = 0L
  CASE 1 OF
      ;; Integer types
      (CTYPE EQ 2 OR CTYPE EQ 3 or ctype eq 14): BEGIN
          IF NOT KEYWORD_SET(NOIEEE) OR KEYWORD_SET(VARICOL) THEN $
            IEEE_TO_HOST, DD 
      END

      ;; Floating and complex types
      (CTYPE GE 4 OR CTYPE LE 6 OR CTYPE EQ 9): BEGIN
          IF NOT KEYWORD_SET(NOIEEE) THEN BEGIN
              IF N_ELEMENTS(NANVALUE) GT 0 THEN W=WHERENAN(DD,COUNT)
              IEEE_TO_HOST, DD
          ENDIF
      END

      ;; String types (CTYPE EQ 7) have already been converted
      ;; in the above CALL_FUNCTION.  No further conversion
      ;; is necessary here.
  ENDCASE

;
;  If the parameters TZERO and TSCAL are non-trivial, then adjust the array by
;  these values.
;
  IF ((NOT KEYWORD_SET(NOIEEE) AND NOT KEYWORD_SET(NOSCALE)) AND $
      (NOT KEYWORD_SET(VARICOL)) AND $
      (N_ELEMENTS(TZERO) EQ 1 AND N_ELEMENTS(TSCAL) EQ 1)) THEN BEGIN

      IF KEYWORD_SET(DF) THEN BEGIN
          ;; Default to float
          TSCAL = FLOAT(TSCAL)
          TZERO = FLOAT(TZERO)
      ENDIF

      FORWARD_FUNCTION UINT, ULONG
      IF CTYPE EQ 2 AND TSCAL[0] EQ 1 AND TZERO[0] EQ 32768 THEN BEGIN
          ;; SPECIAL CASE: Unsigned 16-bit integer
          DD = UINT(DD) - UINT(32768)
      ENDIF ELSE IF CTYPE EQ 3 AND TSCAL[0] EQ 1 AND $
        TZERO[0] EQ 2147483648D THEN BEGIN
          ;; SPECIAL CASE: Unsigned 32-bit integer
          DD = ULONG(DD) - ULONG(2147483648)
      ENDIF ELSE BEGIN
          IF (TSCAL[0] NE 0) AND (TSCAL[0] NE 1) THEN DD = TSCAL[0]*DD
          IF TZERO[0] NE 0 THEN DD = DD + TZERO[0]
      ENDELSE
  ENDIF

;
;  Store NANVALUE everywhere where the data corresponded to IEEE NaN.
;
  IF COUNT GT 0 AND N_ELEMENTS(NANVALUE) GT 0 THEN DD[W] = NANVALUE
  
END

PRO FXBREADM, UNIT, COL, $
              D0,  D1,  D2,  D3,  D4,  D5,  D6,  D7,  D8,  D9, $
              D10, D11, D12, D13, D14, D15, D16, D17, D18, D19, $
              D20, D21, D22, D23, D24, D25, D26, D27, D28, D29, $
              D30, D31, D32, D33, D34, D35, D36, D37, D38, D39, $
              D40, D41, D42, D43, D44, D45, D46, D47, $
              ROW=ROW, VIRTUAL=VIR, DIMENSIONS=DIM, $
              NOSCALE=NOSCALE, NOIEEE=NOIEEE, DEFAULT_FLOAT=DEFAULT_FLOAT, $
              PASS_METHOD=PASS_METHOD, POINTERS=POINTERS, $
              NANVALUE=NANVALUE, BUFFERSIZE=BUFFERSIZE, $
              ERRMSG=ERRMSG, WARNMSG=WARNMSG, STATUS=OUTSTATUS

@fxbintable
        ON_ERROR, 2
;
;  Check the number of parameters.
;
        IF N_PARAMS() LT 2 THEN BEGIN
                MESSAGE = 'Syntax:  FXBREADM, UNIT, COL, D0, D1, ... [, ROW= ]'
                IF N_ELEMENTS(ERRMSG) NE 0 THEN BEGIN
                        ERRMSG = MESSAGE
                        RETURN
                END ELSE MESSAGE, MESSAGE
        ENDIF
        IF N_ELEMENTS(BUFFERSIZE) EQ 0 THEN BUFFERSIZE = 32768L

;
;  COL may be one of several descriptors:
;     * a list of column numbers, beginning with 1
;     * a list of column names
;
        MYCOL = [ COL ]    ; Make sure it is an array

        SC = SIZE(MYCOL)
        NUMCOLS = N_ELEMENTS(MYCOL)
        OUTSTATUS = LONARR(NUMCOLS)
        COLNAMES = 'D'+STRTRIM(LINDGEN(NUMCOLS),2)