Changeset 6873 for trunk/psModules/src/imcombine/pmReadoutCombine.c

Timestamp:

Apr 17, 2006, 8:10:08 AM (20 years ago)

Author:

magnier

Message:

fixed up conflicts

File:

• trunk/psModules/src/imcombine/pmReadoutCombine.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  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);
@@ -418 +690 @@
     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);
@@ -454 +726 @@
         }
 
-        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;
     }
@@ -464 +736 @@
     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);
     }
     return(diffStats);
 }
-
-
 
 /**