Index: /trunk/psModules/src/pmReadoutCombine.c
===================================================================
--- /trunk/psModules/src/pmReadoutCombine.c	(revision 2828)
+++ /trunk/psModules/src/pmReadoutCombine.c	(revision 2829)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-24 02:39:17 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-27 20:14:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20,4 +20,6 @@
 p_psDetermineNumBits(data): This routine takes an enum psStatsOptions as an
 argument and returns the number of non-zero bits.
+ 
+XXX: Use static vectors.
  *****************************************************************************/
 psStatsOptions p_psDetermineNumBits(psStatsOptions data)
@@ -36,6 +38,4 @@
 }
 
-// XXX: Whats this?
-#define MAX_INT 10000
 psImage *pmReadoutCombine(psImage *output,
                           const psList *inputs,
@@ -44,6 +44,6 @@
                           const psVector *scale,
                           bool applyZeroScale,
-                          float gain,
-                          float readnoise)
+                          psF32 gain,
+                          psF32 readnoise)
 {
     PS_PTR_CHECK_NULL(inputs, NULL);
@@ -63,14 +63,14 @@
 
     psStats *stats = params->stats;
-    int i;
-    int j;
-    int maxInputCols = 0;
-    int maxInputRows = 0;
-    int minInputCols = MAX_INT;
-    int minInputRows = MAX_INT;
+    psS32 i;
+    psS32 j;
+    psS32 maxInputCols = 0;
+    psS32 maxInputRows = 0;
+    psS32 minInputCols = PS_MAX_S32;
+    psS32 minInputRows = PS_MAX_S32;
     psListElem *tmpInput = NULL;
     psReadout *tmpReadout = NULL;
-    int numInputs = 0;
-    int tmpI;
+    psS32 numInputs = 0;
+    psS32 tmpI;
     psElemType outputType = PS_TYPE_F32;
 
@@ -153,6 +153,6 @@
         output = psImageAlloc(maxInputCols-minInputCols,
                               maxInputRows-minInputRows, outputType);
-        *(int *) &(output->col0) = minInputCols;
-        *(int *) &(output->row0) = minInputRows;
+        *(psS32 *) &(output->col0) = minInputCols;
+        *(psS32 *) &(output->row0) = minInputRows;
     } else {
         if (((output->col0 + output->numCols) < maxInputCols) ||
@@ -210,5 +210,5 @@
     for (i = output->row0; i < (output->row0 + output->numRows) ; i++) {
         for (j = output->col0; j < (output->col0 + output->numCols) ; j++) {
-            for (int r = 0; r < numInputs ; r++) {
+            for (psS32 r = 0; r < numInputs ; r++) {
                 //  printf("[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outRowLower->data.U32[r], outColLower->data.U32[r], outRowUpper->data.U32[r], outColUpper->data.U32[r]);
                 if ((outRowLower->data.U32[r] <= i) &&
@@ -217,8 +217,8 @@
                         (outColUpper->data.U32[r] > j)) {
 
-                    int imageRow = i - (tmpReadouts[r]->row0 +
-                                        tmpReadouts[r]->image->row0);
-                    int imageCol = j - (tmpReadouts[r]->col0 +
-                                        tmpReadouts[r]->image->col0);
+                    psS32 imageRow = i - (tmpReadouts[r]->row0 +
+                                          tmpReadouts[r]->image->row0);
+                    psS32 imageCol = j - (tmpReadouts[r]->col0 +
+                                          tmpReadouts[r]->image->col0);
 
                     if ((NULL == tmpReadouts[r]->mask) ||
@@ -240,6 +240,6 @@
 
             // Determine how many pixels lie between fracLow and fracHigh.
-            int pixelCount = 0;
-            for (int r = 0; r < numInputs ; r++) {
+            psS32 pixelCount = 0;
+            for (psS32 r = 0; r < numInputs ; r++) {
                 if (tmpPixelMask->data.U8[r] == 0) {
                     if ((params->fracLow <= tmpPixels->data.F32[r]) &&
@@ -254,5 +254,5 @@
             // that range.
             if (pixelCount >= params->nKeep) {
-                for (int r = 0; r < numInputs ; r++) {
+                for (psS32 r = 0; r < numInputs ; r++) {
                     if (tmpPixelMask->data.U8[r] == 0) {
                         if ((params->fracLow <= tmpPixels->data.F32[r]) &&
@@ -266,31 +266,57 @@
             }
 
-            // XXX: Is this correct?
-            // We scale the pixels by the scale vector, if not-NULL.
-            if (scale != NULL) {
-                for (int r = 0; r < numInputs ; r++) {
-                    tmpPixels->data.F32[r]*= scale->data.F32[r];
-                }
-            }
-
-            // XXX: Is this correct?
-            // We add the zero vector, if non-NULL.
-            if (zero != NULL) {
-                for (int r = 0; r < numInputs ; r++) {
-                    tmpPixels->data.F32[r]+= zero->data.F32[r];
-                }
-            }
-
-            // XXX: Is this correct?
-            if (gain > 0.0) {
-                for (int r = 0; r < numInputs ; r++) {
-                    tmpPixels->data.F32[r]*= gain;
-                }
-            }
-
-            // XXX: Is this correct?
-            if (!(readnoise < 0.0)) {
-                for (int r = 0; r < numInputs ; r++) {
-                    tmpPixels->data.F32[r]+= readnoise;
+            // XXX: Still waiting on clarification of the algorithm.
+            #define SOME_FUNC(X, SIGMA) ((X))
+            if ((gain > 0.0) && (readnoise >= 0.0))
+            {
+                psF32 x;
+                psF32 sigma;
+                if (applyZeroScale == true) {
+                    for (psS32 r = 0; r < numInputs ; r++) {
+                        if (zero != NULL) {
+                            x = zero->data.F32[r];
+                        } else {
+                            x = 0.0;
+                        }
+                        if (scale != NULL) {
+                            x+= tmpPixels->data.F32[r] * scale->data.F32[r];
+                        } else {
+                            x+= tmpPixels->data.F32[r];
+                        }
+                        sigma = PS_SQRT_F32((readnoise*readnoise) + gain * x) / gain;
+
+                        tmpPixels->data.F32[r]= SOME_FUNC(x, sigma);
+                    }
+                } else {
+                    for (psS32 r = 0; r < numInputs ; r++) {
+                        x= tmpPixels->data.F32[r];
+
+                        if (zero != NULL) {
+                            sigma = zero->data.F32[r];
+                        } else {
+                            sigma = 0.0;
+                        }
+                        if (scale != NULL) {
+                            sigma+= tmpPixels->data.F32[r] * scale->data.F32[r];
+                        } else {
+                            sigma+= tmpPixels->data.F32[r];
+                        }
+                        sigma = PS_SQRT_F32((readnoise*readnoise) + (gain * sigma)) / gain;
+
+                        tmpPixels->data.F32[r]= SOME_FUNC(x, sigma);
+                    }
+                }
+            } else {
+                if (scale != NULL) {
+                    for (psS32 r = 0; r < numInputs ; r++) {
+                        tmpPixels->data.F32[r]*= scale->data.F32[r];
+                    }
+                }
+
+                // We add the zero vector, if non-NULL.
+                if (zero != NULL) {
+                    for (psS32 r = 0; r < numInputs ; r++) {
+                        tmpPixels->data.F32[r]+= zero->data.F32[r];
+                    }
                 }
             }
@@ -309,5 +335,5 @@
                 return(NULL);
             } else {
-                output->data.F32[i-output->row0][j-output->col0] = (float) statValue;
+                output->data.F32[i-output->row0][j-output->col0] = (psF32) statValue;
             }
         }
Index: /trunk/psModules/test/tst_pmReadoutCombine.c
===================================================================
--- /trunk/psModules/test/tst_pmReadoutCombine.c	(revision 2828)
+++ /trunk/psModules/test/tst_pmReadoutCombine.c	(revision 2829)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-24 00:58:08 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-27 20:14:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -118,7 +118,9 @@
            minOutRow, minOutCol, maxOutRow, maxOutCol);
 
-    output = pmReadoutCombine(output, list, params, zero, scale, true, 1.0, 0.0);
+    output = pmReadoutCombine(output, list, params, zero, scale, true, 0.0, 0.0);
     psF32 NR = (psF32) NUM_READOUTS;
-    psF32 expectedPixel = ((((NR-1.0) * (NR/2.0)) / NR) + VEC_ZERO) * VEC_SCALE;
+    //    psF32 expectedPixel = ((((NR-1.0) * (NR/2.0)) / NR) + VEC_ZERO) * VEC_SCALE;
+    psF32 expectedPixel = ((NR/2.0) * (VEC_ZERO + (VEC_ZERO + VEC_SCALE * (NR - 1)))) / NR;
+
     VerifyTheOutput(output, expectedPixel);
 
@@ -179,8 +181,12 @@
         zero->data.F32[i] = 3.0;
         zeroBig->data.F32[i] = 3.0;
+        zero->data.F32[i] = VEC_ZERO;
+        zeroBig->data.F32[i] = VEC_ZERO;
     }
     for (i=0;i<NUM_READOUTS;i++) {
         scale->data.F32[i] = 6.0;
         scaleBig->data.F32[i] = 6.0;
+        scale->data.F32[i] = VEC_SCALE;
+        scaleBig->data.F32[i] = VEC_SCALE;
     }
 
@@ -220,19 +226,26 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with NULL zero vector.\n");
-    rc = pmReadoutCombine(NULL, list, params, NULL, scale, true, 1.0, 0.0);
-    if (rc == NULL) {
-        //XXX: We should verify the output image here.
+    rc = pmReadoutCombine(NULL, list, params, NULL, scale, true, 0.0, 0.0);
+    if (rc != NULL) {
+        psF32 NR = (psF32) NUM_READOUTS;
+        psF32 expectedPixel = ((NR/2.0) * (0.0 + (0.0 + VEC_SCALE * (NR - 1)))) / NR;
+        if (false == VerifyTheOutput(rc, expectedPixel)) {
+            testStatus = false;
+        }
+
+        if (rc->type.type != scale->type.type) {
+            printf("ERROR: output readout->image has incorrect type.\n");
+            testStatus = false;
+        }
+        psFree(rc);
+    } else {
         printf("ERROR: pmReadoutCombine() returned NULL\n");
         testStatus = false;
-    }
-    if (rc->type.type != scale->type.type) {
-        printf("ERROR: output readout->image has incorrect type.\n");
-        testStatus = false;
-    }
-    psFree(rc);
+
+    }
 
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect length zero vector (too small).  Should generate error.\n");
-    rc = pmReadoutCombine(NULL, list, params, zeroHalf, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(NULL, list, params, zeroHalf, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -242,5 +255,5 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect type zero vector.  Should generate error.\n");
-    rc = pmReadoutCombine(NULL, list, params, zeroF64, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(NULL, list, params, zeroF64, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -250,7 +263,10 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect length zero vector (too big).  Should generate warning.\n");
-    rc = pmReadoutCombine(output, list, params, zeroBig, scale, true, 1.0, 0.0);
-    if (rc != NULL) {
-        if (false == VerifyTheOutput(rc, 45.0)) {
+    rc = pmReadoutCombine(output, list, params, zeroBig, scale, true, 0.0, 0.0);
+    if (rc != NULL) {
+        psF32 NR = (psF32) NUM_READOUTS;
+        psF32 expectedPixel = ((NR/2.0) * (VEC_ZERO + (VEC_ZERO + VEC_SCALE * (NR - 1)))) / NR;
+
+        if (false == VerifyTheOutput(rc, expectedPixel)) {
             testStatus = false;
         }
@@ -264,15 +280,27 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with NULL scale vector.\n");
-    rc = pmReadoutCombine(output, list, params, zero, NULL, true, 1.0, 0.0);
-    if (rc == NULL) {
-        //XXX: We should verify the output image here.
+    rc = pmReadoutCombine(output, list, params, zero, NULL, true, 0.0, 0.0);
+
+    if (rc != NULL) {
+        psF32 NR = (psF32) NUM_READOUTS;
+        psF32 expectedPixel = ((NR/2.0) * (VEC_ZERO + (VEC_ZERO + 1.0 * (NR - 1)))) / NR;
+        if (false == VerifyTheOutput(rc, expectedPixel)) {
+            testStatus = false;
+        }
+
+        if (rc->type.type != scale->type.type) {
+            printf("ERROR: output readout->image has incorrect type.\n");
+            testStatus = false;
+        }
+        psFree(rc);
+    } else {
         printf("ERROR: pmReadoutCombine() returned NULL\n");
         testStatus = false;
-    }
-    psFree(rc);
+
+    }
 
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect length scale vector (too small).  Should generate error.\n");
-    rc = pmReadoutCombine(output, list, params, zero, scaleHalf, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, params, zero, scaleHalf, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -282,5 +310,5 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect type scale vector.  Should generate error.\n");
-    rc = pmReadoutCombine(output, list, params, zero, scaleF64, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, params, zero, scaleF64, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -290,7 +318,10 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with incorrect length scale vector (too big).  Should generate warning.\n");
-    rc = pmReadoutCombine(output, list, params, zero, scaleBig, true, 1.0, 0.0);
-    if (rc != NULL) {
-        if (false == VerifyTheOutput(rc, 45.0)) {
+    rc = pmReadoutCombine(output, list, params, zero, scaleBig, true, 0.0, 0.0);
+    if (rc != NULL) {
+        psF32 NR = (psF32) NUM_READOUTS;
+        psF32 expectedPixel = ((NR/2.0) * (VEC_ZERO + (VEC_ZERO + VEC_SCALE * (NR - 1)))) / NR;
+
+        if (false == VerifyTheOutput(rc, expectedPixel)) {
             testStatus = false;
         }
@@ -305,5 +336,5 @@
     printf("Calling pmReadoutCombine() insufficient size output image.  Should generate error, return NULL.\n");
     output = psImageAlloc(1, 1, PS_TYPE_F32);
-    rc = pmReadoutCombine(output, list, params, zero, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, params, zero, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -315,5 +346,5 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with NULL input list.  Should generate error, return NULL.\n");
-    rc = pmReadoutCombine(output, NULL, params, zero, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, NULL, params, zero, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -323,5 +354,5 @@
     printf("----------------------------------------------------------------------------\n");
     printf("Calling pmReadoutCombine() with NULL params.  Should generate error, return NULL.\n");
-    rc = pmReadoutCombine(output, list, NULL, zero, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, NULL, zero, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -332,5 +363,5 @@
     psStatsOptions oldStatsOpts = params->stats->options |= PS_STAT_MIN;
     printf("Calling pmReadoutCombine() with multiple stats->options.  Should generate error, return NULL.\n");
-    rc = pmReadoutCombine(output, list, params, zero, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, params, zero, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");
@@ -342,5 +373,5 @@
     psStats *oldStats = params->stats;
     printf("Calling pmReadoutCombine() with NULL param->stats.  Should generate error, return NULL.\n");
-    rc = pmReadoutCombine(output, list, params, zero, scale, true, 1.0, 0.0);
+    rc = pmReadoutCombine(output, list, params, zero, scale, true, 0.0, 0.0);
     if (rc != NULL) {
         printf("ERROR: pmReadoutCombine() did not return NULL\n");