Changeset 6448 for branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.c

Timestamp:

Feb 17, 2006, 7:13:42 AM (20 years ago)

Author:

magnier

Message:

bulk merge of eam_rel9_p0 onto this branch

File:

• branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.c (modified) (6 diffs)

branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.c

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-06 21:03:25 $
+ *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 17:13:41 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -36 +36 @@
 }
 
+static void pmCombineParamsFree (pmCombineParams *params)
+{
+
+    if (params == NULL)
+        return;
+
+    psFree (params->stats);
+    return;
+}
+
+pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions)
+{
+
+    pmCombineParams *params = psAlloc (sizeof(pmCombineParams));
+    psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree);
+
+    params->stats = psStatsAlloc (statsOptions);
+    params->maskVal = 0;
+    params->fracHigh = 0.25;
+    params->fracHigh = 0.25;
+    params->nKeep = 3;
+
+    return (params);
+}
+
 /******************************************************************************
 XXX: Must add support for S16 and S32 types.  F32 currently supported.
  *****************************************************************************/
 psImage *pmReadoutCombine(psImage *output,
-                          const psList *inputs,
-                          psCombineParams *params,
+                          const psArray *inputs,
                           const psVector *zero,
                           const psVector *scale,
+                          pmCombineParams *params,
                           bool applyZeroScale,
                           psF32 gain,
                           psF32 readnoise)
+{
+    PS_ASSERT_PTR_NON_NULL(inputs, NULL);
+    PS_ASSERT_PTR_NON_NULL(params, NULL);
+    PS_ASSERT_PTR_NON_NULL(params->stats, NULL);
+    if (zero != NULL) {
+        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
+        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL);
+    }
+    if (scale != NULL) {
+        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
+        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
+    }
+    if ((zero != NULL) && (scale != NULL)) {
+        PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL);
+        // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
+    }
+
+    psStats *stats = params->stats;
+    psS32 maxInputCols = 0;
+    psS32 maxInputRows = 0;
+    psS32 minInputCols = PS_MAX_S32;
+    psS32 minInputRows = PS_MAX_S32;
+    pmReadout *tmpReadout = NULL;
+    psS32 tmpI;
+    psElemType outputType = PS_TYPE_F32;
+
+    if (DetermineNumBits(stats->options) != 1) {
+        psError(PS_ERR_UNKNOWN, true,
+                "Multiple statistical options have been requested.  Returning NULL.\n");
+        return(NULL);
+    }
+
+    // We step through each readout in the input image list.  If any readout is
+    // NULL, empty, or has the wrong type, we generate an error and return
+    // NULL.  We determine how big of an output image is needed to combine
+    // these input images.  We do this by taking the
+    //     max(readout->col0 + readout->numCols + image->col0 + image->numCols)
+    //     max(readout->row0 + readout->numRows + image->row0 + image->numRows)
+    //
+    for (int i = 0; i < inputs->n; i++) {
+        tmpReadout = inputs->data[i];
+        PS_ASSERT_READOUT_NON_NULL(tmpReadout, output);
+        PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output);
+        PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output);
+
+        minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0));
+        tmpI = tmpReadout->row0 +
+               tmpReadout->image->row0 +
+               tmpReadout->image->numRows;
+        maxInputRows = PS_MAX(maxInputRows, tmpI);
+
+        minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0));
+        tmpI = tmpReadout->col0 +
+               tmpReadout->image->col0 +
+               tmpReadout->image->numCols;
+        maxInputCols = PS_MAX(maxInputCols, tmpI);
+    }
+
+    // We ensure that the zero vector is of the proper size.
+    if (zero != NULL) {
+        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
+        if (zero->n < inputs->n) {
+            psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d).  Returning NULL.\n", zero->n);
+            return(NULL);
+        } else if (zero->n > inputs->n) {
+            // XXX EAM : abort on this condition? is probably an error
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: the zero vector too many elements (%d)\n", zero->n);
+        }
+    }
+
+    // We ensure that the scale vector is of the proper size.
+    if (scale != NULL) {
+        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
+        if (scale->n < inputs->n) {
+            psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d).  Returning NULL.\n", scale->n);
+            return(NULL);
+        } else if (scale->n > inputs->n) {
+            // XXX EAM : abort on this condition? is probably an error
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: the scale vector has too many elements (%d)\n", scale->n);
+        }
+    }
+
+    // At this point, the following variables have been computed:
+    // maxInputRows: the largest input row value, in output image space.
+    // maxInputCols: the largest input column value, in output image space.
+    // minInputRows: the smallest input row value, in output image space.
+    // minInputCols: the smallest input column value, in output image space.
+    //
+    if (output == NULL) {
+        output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType);
+        *(psS32 *) &(output->col0) = minInputCols;
+        *(psS32 *) &(output->row0) = minInputRows;
+    } else {
+
+        // XXX EAM : recycle the existing output image?  why would we care about the existing pixels?
+        PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL);
+        if (((output->col0 + output->numCols) < maxInputCols) ||
+                ((output->row0 + output->numRows) < maxInputRows)) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "Output image (%d, %d) is too small to hold combined images.  Returning NULL.\n",
+                    output->row0 + output->numRows,
+                    output->col0 + output->numCols);
+            return(NULL);
+        }
+
+        // reset output origin using logic of above
+        *(psS32 *) &(output->col0) = minInputCols;
+        *(psS32 *) &(output->row0) = minInputRows;
+    }
+
+    psVector *tmpPixels     = psVectorAlloc(inputs->n, PS_TYPE_F32);
+    psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32);
+    psVector *outRowLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outRowUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outColLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outColUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+
+    // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates,
+    // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper).
+    for (int i = 0; i < inputs->n; i++) {
+        tmpReadout = (pmReadout *) inputs->data[i];
+        outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0;
+        outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0;
+        outRowUpper->data.U32[i] = tmpReadout->row0 +
+                                   tmpReadout->image->row0 +
+                                   tmpReadout->image->numRows;
+        outColUpper->data.U32[i] = tmpReadout->col0 +
+                                   tmpReadout->image->col0 +
+                                   tmpReadout->image->numCols;
+    }
+
+    // We loop through each pixel in the output image.  We loop through each
+    // input readout.  We determine if that output pixel is contained in the
+    // image from that readout.  If so, we save it in psVector tmpPixels.
+    // If not, we set a mask for that element in tmpPixels.  Then, we mask off
+    // pixels not between fracLow and fracHigh.  Then we call the vector
+    // stats routine on those pixels/mask.  Then we set the output pixel value
+    // to the result of the stats call.
+
+    int nx, ny;
+    int nKeep, nMin;
+    float keepFrac = 1.0 - params->fracLow - params->fracHigh;
+    float value = 0;
+    psF32 *saveVector = tmpPixelsKeep->data.F32;
+
+    for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) {
+        if (j % 10 == 0)
+            fprintf (stderr, ".");
+        for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) {
+            int nPix = 0;
+            for (int r = 0; r < inputs->n; r++) {
+                tmpReadout = (pmReadout *) inputs->data[r];
+
+                // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]);
+                if (i <  outColLower->data.U32[r])
+                    continue;
+                if (i >= outColUpper->data.U32[r])
+                    continue;
+                if (j <  outRowLower->data.U32[r])
+                    continue;
+                if (j >= outRowUpper->data.U32[r])
+                    continue;
+
+                nx = i - (tmpReadout->col0 + tmpReadout->image->col0);
+                ny = j - (tmpReadout->row0 + tmpReadout->image->row0);
+
+                if (tmpReadout->mask != NULL) {
+                    if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal)
+                        continue;
+                }
+
+                tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx];
+                // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]);
+                nPix ++;
+            }
+            tmpPixels->n = nPix;
+
+            // are there enough valid pixels to apply fracLow,fracHigh?
+            nKeep = nPix * keepFrac;
+            if (nKeep >= params->nKeep) {
+                psVectorSort (tmpPixels, tmpPixels);
+                nMin = nPix * params->fracLow;
+                tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin];
+                tmpPixelsKeep->n = nKeep;
+            } else {
+                tmpPixelsKeep->data.F32 = tmpPixels->data.F32;
+                tmpPixelsKeep->n = nPix;
+            }
+
+            // tmpPixelsKeep is already sorted.  sample mean and median are very easy
+            if (stats->options & PS_STAT_SAMPLE_MEAN) {
+                value = 0;
+                for (int r = 0; r < tmpPixelsKeep->n; r++) {
+                    value += tmpPixelsKeep->data.F32[r];
+                }
+                if (tmpPixelsKeep->n == 0) {
+                    value = 0;
+                } else {
+                    value = value / tmpPixelsKeep->n;
+                }
+            }
+            if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
+                int r = tmpPixelsKeep->n / 2;
+                if (tmpPixelsKeep->n == 0) {
+                    value = 0;
+                    goto got_value;
+                }
+                if (tmpPixelsKeep->n % 2 == 1) {
+                    int r = 0.5*tmpPixelsKeep->n;
+                    value = tmpPixelsKeep->data.F32[r];
+                    goto got_value;
+                }
+                if (tmpPixelsKeep->n % 2 == 0) {
+                    value = 0.5*(tmpPixelsKeep->data.F32[r] +
+                                 tmpPixelsKeep->data.F32[r-1]);
+                    goto got_value;
+                }
+            }
+got_value:
+            output->data.F32[j-output->row0][i-output->col0] = value;
+        }
+    }
+    tmpPixelsKeep->data.F32 = saveVector;
+
+    psFree(tmpPixels);
+    psFree(tmpPixelsKeep);
+    psFree(outRowLower);
+    psFree(outRowUpper);
+    psFree(outColLower);
+    psFree(outColUpper);
+
+    return(output);
+}
+
+/******************************************************************************
+XXX: Must add support for S16 and S32 types.  F32 currently supported.
+ *****************************************************************************/
+psImage *pmReadoutCombine_OLD(psImage *output,
+                              const psList *inputs,
+                              pmCombineParams *params,
+                              const psVector *zero,
+                              const psVector *scale,
+                              bool applyZeroScale,
+                              psF32 gain,
+                              psF32 readnoise)
 {
     PS_ASSERT_PTR_NON_NULL(inputs, NULL);
@@ -410 +682 @@
     psRegion minRegion;
     psRegion maxRegion;
-    psStats *minStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
-    psStats *maxStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
+    psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32);
@@ -445 +717 @@
         }
 
-        fp->midValue = 0.5 * (maxStats->fittedMean + minStats->fittedMean);
-        fp->delta = maxStats->fittedMean - minStats->fittedMean;
+        fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian);
+        fp->delta = maxStats->robustMedian - minStats->robustMedian;
         diffs->data.F32[i] = fp->delta;
     }
@@ -455 +727 @@
     psFree(diffs);
     if (diffStats == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not determine fitted median of the differences.\n");
+        psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n");
         return(NULL);
     }
@@ -461 +733 @@
 }
 
-
-
-/**
- *
- * The input array fluxLevels consists of Ni vectors, one per mosaic image.
- * Each vector consists of Nj elements, each a measurement of the input
- * flat-field image flux levels. All of these vectors must be constructed with
- * the same number of elements, or the function will return an error. If a chip
- * is missing from a particular image, that element should be set to NaN. The
- * vector chipGains supplies initial guesses for the chip gains. If the vector
- * contains the values 0.0 or NaN for any of the elements, the gain is set to the
- * mean of the valid values. If the vector length does not match the number of
- * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the
- * vector length modified to match the number of chips defined by the supplied
- * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the
- * routine will set the vector length to the number of source images regardless
- * of the initial state of the vector. All vectors used by this function must be
- * of type PS_DATA_F64.
- *
- 
-fluxLevels(i, j): for each flat field image i, this psArray contains a vector
-with an elemenmt for each chip j.  So, fluxLevels(i, j) corresponds to the
-measured flux M_(i, j) for flat image i, chip j.
- 
-chipGains[]: has j elements, one for each chip.
- 
- 
-They have the observed flux levels for each chip of each image.  They want to
-solve for the actual flux levels and the gain of each chip.
- 
-Okay, they want to solve for source fluxes and chip gains.
- 
- *
- */
-bool pmFlatNormalization(
-    psVector *sourceFlux,
-    psVector *chipGains,
-    psArray *fluxLevels)
-{
-    PS_ASSERT_PTR_NON_NULL(fluxLevels, false);
-    psS32 numImages = fluxLevels->n;
-    psS32 numChips = ((psVector *) fluxLevels->data[0])->n;
-    for (psS32 i = 0 ; i < numImages ; i++) {
-        psVector *tmpVec = (psVector *) fluxLevels->data[i];
-        PS_ASSERT_VECTOR_NON_NULL(tmpVec, false);
-        PS_ASSERT_VECTOR_TYPE(tmpVec, PS_TYPE_F64, false);
-        PS_ASSERT_VECTOR_SIZE(tmpVec, numChips, false);
-    }
-
-    //
-    // Ensure that *localChipGains points to a vector of the same length as numImages.
-    //
-    PS_ASSERT_PTR_NON_NULL(chipGains, false);
-    PS_ASSERT_VECTOR_TYPE(chipGains, PS_TYPE_F64, false);
-    psVector *localChipGains = chipGains;
-    if (numChips != chipGains->n) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: the chipGains vector length does not match the number of chips.\n");
-        localChipGains = psVectorAlloc(numChips, PS_TYPE_F64);
-        psBool rc = psVectorInit(localChipGains, 1.0);
-        if (rc == false) {
-            printf("XXX: gen error\n");
-        }
-    }
-
-    //
-    // If the chipGains vector contains the values 0.0 or NaN for any of the elements,
-    // the gain is set to the mean of the valid values.
-    //
-    psBool meanFlag = false;
-    psVector *chipGainsMask = psVectorAlloc(chipGains->n, PS_TYPE_U8);
-    for (psS32 i = 0 ; i < chipGains->n ; i++) {
-        if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||
-                (isnan(chipGains->data.F64[i]))) {
-            chipGainsMask->data.U8[i] = 1;
-            meanFlag = true;
-        }
-    }
-    // Must calculate the mean.
-    if (meanFlag == true) {
-        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        stats = psVectorStats(stats, chipGains, NULL, chipGainsMask, 1);
-        if (stats == NULL) {
-            printf("XXX: gen error\n");
-        }
-        psF64 mean;
-        psBool rc = p_psGetStatValue(stats, &mean);
-        if (rc == false) {
-            printf("XXX: gen error\n");
-        }
-        // Set the gain to this mean for chips with a gain of 0.0 or NAN
-
-        for (psS32 i = 0 ; i < chipGains->n ; i++) {
-            if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||
-                    (isnan(chipGains->data.F64[i]))) {
-                chipGains->data.F64[i] = mean;
-            }
-        }
-    }
-
-    //
-    // Assert that sourceFlux is non-NULL, correct type, correct size.
-    //
-    PS_ASSERT_PTR_NON_NULL(sourceFlux, false);
-    PS_ASSERT_VECTOR_TYPE(sourceFlux, PS_TYPE_F64, false);
-    psVectorRealloc(sourceFlux, numImages);
-
-    //    psFree(psVector);
-    if (numImages != chipGains->n) {
-        psFree(localChipGains);
-    }
-
-    return(true);
-}