include include "wpdef.h" .help minrej .nf ---------------------------------------------------------------------------- COMBINING IMAGES: MINIMUM REJECTION ALGORITHM For more than two input images they are combined by scaling and taking a weighted average excluding the minimum 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: MINREJ -- Combine image lines after rejecting the minimum value at each pixel, without weighting or scaling. DQMINREJ -- Combine image lines using Data Quality flags, after rejecting the minimum value at each pixel, possibly with weighting and/or scaling. WTMINREJ -- Combine image lines after rejecting the minimum value at each pixel, with weighting and/or scaling. .endhelp ----------------------------------------------------------------------- #$for (silrdx) $for (silrd) ################################################################################# # MINREJ -- Combine image lines after rejecting the minimum value at each # # pixel without weighting or scaling. This routine is based # # upon the `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure minrej$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 minimum value int nims # Total of non-rejected images $if (datatype == sil) real val # Data value @pixel real sumval # Running sum of pixel values real minval # Smallest data value @pixel $else PIXEL val # Data value @pixel PIXEL sumval # Running sum of pixel values PIXEL minval # Smallest data value @pixel $endif begin nims = nimages - 1 do i = 1, npts { sumval = 0. minval = 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(minval)) { # $else if (val < minval) { # $endif sumval = sumval + minval minval = val k = j } else sumval = sumval + val } # Save output value and set minimum 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 ################################################################################# # DQMINREJ -- Combine image lines, modulo the Data Quality flags, after # # rejecting the minimum value at each pixel and possibly with # # weighting or scaling. This routine is based upon the # # `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure dqminrej$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 minimum value int ncount # Total of non-rejected pixels real netwt # Sum of weights minus minimum wt real sumwt # Sum of weights for each pixel $if (datatype == sil) real val # Data value @pixel real sumval # Running sum of non-flagged pixel values real minval # Smallest data value @line $else PIXEL val # Data value @pixel PIXEL sumval # Running sum of non-flagged pixel values PIXEL minval # Smallest data value @line $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 minval = +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(minval)) { # $else if (val < minval) { # $endif minval = val k = j } sumval = sumval + WTS[j] * val sumwt = sumwt + WTS[j] } else Mem$t[data[j]+i-1] = INDEF } # Save output value and set minimum 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] * minval) / netwt Mem$t[data[k]+i-1] = INDEF } } end ################################################################################# # WTMINREJ -- Combine image lines after rejecting the minimum value at each # # pixel with weighting and/or scaling. This routine is based # # upon the `images.imcombine' package. # # # # Development version: 1/91 RAShaw # procedure wtminrej$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 minimum value $if (datatype == sil) real val # Data value @pixel real minval # Smallest data value @pixel real sumval # Running sum of non-flagged pixel values $else PIXEL val # Data value @pixel PIXEL sumval # Running sum of non-flagged pixel values PIXEL minval # Smallest data value @pixel $endif begin do i = 1, npts { # Initialize local variables. minval = 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(minval)) { # $else if (val < minval) { # $endif sumval = sumval + WTS[k] * minval minval = val k = j } else sumval = sumval + WTS[j] * val } # Save output value and set minimum 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