include include "wpdef.h" .help maxrej .nf ---------------------------------------------------------------------------- COMBINING IMAGES: MAXIMUM REJECTION ALGORITHM For more than one input image they are combined by scaling and taking a weighted average excluding the maximum value and, if DQFs are used, flagged bad data. The exposure time of the output image is the scaled and weighted average of the input image exposure times. The average is computed in real arithmetic with trunction on output if the output image is an integer datatype. PROCEDURES: MAXREJ -- Combine image lines without weighting or scaling. DQMAXREJ -- Combine image lines using Data Quality flags, possibly with weighting or scaling. WTMAXREJ -- Combine image lines with weighting or scaling. .endhelp ----------------------------------------------------------------------- $for (silrd) ################################################################################# # MAXREJ -- Combine image lines after rejecting the maximum value at each # # pixel without weighting or scaling. This routine is based on # # the `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure maxrej$t (data, output, nimages, npts) # Passed arguments: pointer data[nimages] # IMIO data pointers $if (datatype == sil) real output[npts] # Output line (returned) $else PIXEL output[npts] # Output line (returned) $endif int nimages # Number of images to be combined int npts # Number of pixels per image line # Local variables: int i, j # Dummy loop counters int k # Index of maximum value int nims # Total of non-rejected images $if (datatype == sil) real maxval # Largest data value @pixel real sumval # Running sum of pixel values real val # Data value @pixel $else PIXEL maxval # Largest data value @pixel PIXEL sumval # Running sum of pixel values PIXEL val # Data value @pixel $endif begin nims = nimages - 1 do i = 1, npts { sumval = 0. maxval = Mem$t[data[1]+i-1] k = 1 do j = 2, nimages { val = Mem$t[data[j]+i-1] # $if (datatype == x) # if (abs(val) > abs(maxval)) { # $else if (val > maxval) { # $endif sumval = sumval + maxval maxval = val k = j } else sumval = sumval + val } # Store output value and set maximum value in working data array to INDEF for # future use (e.g., in SIGMA routine). output[i] = sumval / nims Mem$t[data[k]+i-1] = INDEF } end ################################################################################# # DQMAXREJ -- Combine image lines after rejecting the maximum value at each # # pixel, modulo the Data Quality flags, and possibly with # # weighting or scaling. This routine is based on the # # `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure dqmaxrej$t (data, dqfdata, output, nimages, npts) include "wpcom.h" # Passed arguments: pointer data[nimages] # IMIO data pointers int dqfdata[nimages] # Data Quality File pointer $if (datatype == sil) real output[npts] # Output line (returned) $else PIXEL output[npts] # Output line (returned) $endif int nimages # Number of images to be combined int npts # Number of pixels per image line # Local variables: int bflag[IMS_MAX] # int i, j # Dummy loop counters int k # Index of maximum value int ncount # Total of non-rejected pixels real netwt # Sum of weights minus maximum wt real sumwt # Sum of weights for each pixel $if (datatype == sil) real maxval # Largest data value @pixel real sumval # Running sum of non-flagged pixel values real val # Data value @pixel $else PIXEL maxval # Largest data value @pixel PIXEL sumval # Running sum of non-flagged pixel values PIXEL val # Data value @pixel $endif begin do i = 1, npts { # Select user-chosen Data Quality bits: do j = 1, nimages bflag[j] = Memi[dqfdata[j]+i-1] call aandki (bflag, BADBITS, bflag, nimages) # Initialize other variables maxval = -MAX_REAL sumval = 0. sumwt = 0. k = 1 ncount = 0 do j = 1, nimages { if (bflag[j] == 0) { ncount = ncount + 1 val = Mem$t[data[j]+i-1] / SCALES[j] - ZEROS[j] # $if (datatype == x) # if (abs(val) > abs(maxval)) { # $else if (val > maxval) { # $endif maxval = val k = j } sumval = sumval + WTS[j] * val sumwt = sumwt + WTS[j] } else Mem$t[data[j]+i-1] = INDEF } # Store output value and set maximum value in working data array to INDEF for # future use (e.g., in SIGMA routine). netwt = sumwt - WTS[k] if (netwt <= 0.) output[i] = BLANK else { output[i] = (sumval - WTS[k] * maxval) / netwt Mem$t[data[k]+i-1] = INDEF } } end ################################################################################# # WTMAXREJ -- Combine image lines after rejecting the maximum value at each # # pixel with weighting or scaling. This routine is based on # # the `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure wtmaxrej$t (data, output, nimages, npts) include "wpcom.h" # Passed arguments: pointer data[nimages] # IMIO data pointers $if (datatype == sil) real output[npts] # Output line (returned) $else PIXEL output[npts] # Output line (returned) $endif int nimages # Number of images to be combined int npts # Number of pixels per image line # Local variables: int i, j # Dummy loop counters int k # Index of maximum value $if (datatype == sil) real maxval # Largest data value @pixel real sumval # Running sum of non-flagged pixel values real val # Data value @pixel $else PIXEL maxval # Largest data value @pixel PIXEL sumval # Running sum of non-flagged pixel values PIXEL val # Data value @pixel $endif begin do i = 1, npts { # Initialize local variables maxval = Mem$t[data[1]+i-1] / SCALES[1] - ZEROS[1] sumval = 0. k = 1 do j = 2, nimages { val = Mem$t[data[j]+i-1] / SCALES[j] - ZEROS[j] # $if (datatype == x) # if (abs(val) > abs(maxval)) { # $else if (val > maxval) { # $endif sumval = sumval + WTS[k] * maxval maxval = val k = j } else sumval = sumval + WTS[j] * val } # Store output value and set maximum value in working data array to INDEF for # future use (e.g., in SIGMA routine). output[i] = sumval / (1. - WTS[k]) Mem$t[data[k]+i-1] = INDEF } end $endfor