Changeset 5516 for trunk/psModules/src/imsubtract

Timestamp:

Nov 15, 2005, 10:09:03 AM (20 years ago)

Author:

gusciora

Message:

SubtractBias was recoded. Significant mods to removeBadPixels.
Additional mods to other files.

Location:

trunk/psModules/src/imsubtract

Files:

• pmSubtractBias.c (modified) (11 diffs)
• pmSubtractBias.h (modified) (2 diffs)
• pmSubtractSky.c (modified) (11 diffs)

trunk/psModules/src/imsubtract/pmSubtractBias.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-20 23:06:24 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-15 20:09:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
  */
-
+/*****************************************************************************/
+/* INCLUDE FILES                                                             */
+/*****************************************************************************/
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include "pslib.h"
 #if HAVE_CONFIG_H
 #include <config.h>
 #endif
-
 #include "pmSubtractBias.h"
 
-#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
-#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
-
+/*****************************************************************************/
+/* DEFINE STATEMENTS                                                         */
+/*****************************************************************************/
 // XXX: put these in psConstants.h
-void PS_POLY1D_PRINT(psPolynomial1D *poly)
+void PS_POLY1D_PRINT(
+    psPolynomial1D *poly)
 {
     printf("-------------- PS_POLY1D_PRINT() --------------\n");
@@ -51 +57 @@
 }\
 
-/******************************************************************************
-psSubtractFrame(): this routine will take as input a readout for the input
-image and a readout for the bias image.  The bias image is subtracted in
-place from the input image.
+/*****************************************************************************/
+/* TYPE DEFINITIONS                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* GLOBAL VARIABLES                                                          */
+/*****************************************************************************/
+psS32 currentId = 0;                // XXX: remove
+psS32 memLeaks = 0;                 // XXX: remove
+//PRINT_MEMLEAKS(8); XXX
+/*****************************************************************************/
+/* FILE STATIC VARIABLES                                                     */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - LOCAL                                           */
+/*****************************************************************************/
+
+/******************************************************************************
+psSubtractFrame(): this routine will take as input the pmReadout for the input
+image and a pmReadout for the bias image.  The bias image is subtracted in
+place from the input image.  We assume that sizes and types are checked
+elsewhere.
+ 
+XXX: Verify that the image and readout offsets are being used the right way.
+ 
+XXX: Ensure that it does the correct thing with image size.
 *****************************************************************************/
-static pmReadout *SubtractFrame(pmReadout *in,
-                                const pmReadout *bias)
-{
-    psS32 i;
-    psS32 j;
-
-    if (bias == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL.  Returning original image.\n");
-        return(in);
-    }
-
-
-    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
-        return(in);
-    }
-    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
-        return(in);
-    }
-
-    for (i=0;i<in->image->numRows;i++) {
-        for (j=0;j<in->image->numCols;j++) {
+static pmReadout *SubtractFrame(
+    pmReadout *in,
+    const pmReadout *bias)
+{
+    // XXX: When did the ->row0 and ->col0 offsets get coded?
+    for (psS32 i=0;i<in->image->numRows;i++) {
+        for (psS32 j=0;j<in->image->numCols;j++) {
             in->image->data.F32[i][j]-=
                 bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
+
             if ((in->mask != NULL) && (bias->mask != NULL)) {
                 (in->mask->data.U8[i][j])|=
@@ -92 +105 @@
 }
 
+
+/******************************************************************************
+psSubtractDarkFrame(): this routine will take as input the pmReadout for the
+input image and a pmReadout for the dark image.  The dark image is scaled and
+subtracted in place from the input image.
+ 
+XXX: Verify that the image and readout offsets are being used the right way.
+ 
+XXX: Ensure that it does the correct thing with image size.
+*****************************************************************************/
+static pmReadout *SubtractDarkFrame(
+    pmReadout *in,
+    const pmReadout *dark,
+    psF32 scale)
+{
+    // XXX: When did the ->row0 and ->col0 offsets get coded?
+    if (fabs(scale) > FLT_EPSILON) {
+        for (psS32 i=0;i<in->image->numRows;i++) {
+            for (psS32 j=0;j<in->image->numCols;j++) {
+                in->image->data.F32[i][j]-=
+                    (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]);
+
+                if ((in->mask != NULL) && (dark->mask != NULL)) {
+                    (in->mask->data.U8[i][j])|=
+                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
+                }
+            }
+        }
+    } else {
+        for (psS32 i=0;i<in->image->numRows;i++) {
+            for (psS32 j=0;j<in->image->numCols;j++) {
+                in->image->data.F32[i][j]-=
+                    dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0];
+
+                if ((in->mask != NULL) && (dark->mask != NULL)) {
+                    (in->mask->data.U8[i][j])|=
+                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
+                }
+            }
+        }
+    }
+
+    return(in);
+}
+
 /******************************************************************************
 ImageSubtractScalar(): subtract a scalar from the input image.
  
-XXX: Use a psLib function for this.
- 
-XXX: This should
- *****************************************************************************/
-static psImage *ImageSubtractScalar(psImage *image,
-                                    psF32 scalar)
+XXX: Is there a psLib function for this?
+ *****************************************************************************/
+static psImage *ImageSubtractScalar(
+    psImage *image,
+    psF32 scalar)
 {
     for (psS32 i=0;i<image->numRows;i++) {
@@ -164 +221 @@
 
     if (numOptions == 0) {
-        psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n");
+        psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
     }
     if (numOptions != 1) {
@@ -173 +230 @@
 }
 
-
+/******************************************************************************
+Polynomial1DCopy(): This private function copies the members of the existing
+psPolynomial1D "in" into the existing psPolynomial1D "out".  The previous
+members of the existing psPolynomial1D "out" are psFree'ed.
+ *****************************************************************************/
+static psBool Polynomial1DCopy(
+    psPolynomial1D *out,
+    psPolynomial1D *in)
+{
+    psFree(out->coeff);
+    psFree(out->coeffErr);
+    psFree(out->mask);
+
+    out->type = in->type;
+    out->nX = in->nX;
+
+    out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
+    // XXX: use memcpy
+    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
+        out->coeff[i] = in->coeff[i];
+    }
+
+    out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
+    // XXX: use memcpy
+    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
+        out->coeffErr[i] = in->coeffErr[i];
+    }
+
+    out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType));
+    // XXX: use memcpy
+    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
+        out->mask[i] = in->mask[i];
+    }
+
+    return(true);
+}
+
+/******************************************************************************
+Polynomial1DDup(): This private function duplicates and then returns the input
+psPolynomial1D "in".
+ *****************************************************************************/
+static psPolynomial1D *Polynomial1DDup(
+    psPolynomial1D *in)
+{
+    psPolynomial1D *out = psPolynomial1DAlloc(in->nX, in->type);
+    Polynomial1DCopy(out, in);
+    return(out);
+}
+
+
+/******************************************************************************
+SplineCopy(): This private function copies the members of the existing
+psSpline in into the existing psSpline out.
+ *****************************************************************************/
+static psBool SplineCopy(
+    psSpline1D *out,
+    psSpline1D *in)
+{
+    PS_ASSERT_PTR_NON_NULL(out, false);
+    PS_ASSERT_PTR_NON_NULL(in, false);
+
+    for (psS32 i = 0 ; i < out->n ; i++) {
+        psFree(out->spline[i]);
+    }
+    psFree(out->spline);
+    psFree(out->knots);
+    psFree(out->p_psDeriv2);
+
+    out->n = in->n;
+    out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *));
+    for (psS32 i = 0 ; i < in->n ; i++) {
+        out->spline[i] = Polynomial1DDup(in->spline[i]);
+    }
+
+    out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
+
+    out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
+    // XXX: use memcpy
+    for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
+        out->p_psDeriv2[i] = in->p_psDeriv2[i];
+    }
+
+    return(true);
+}
 
 /******************************************************************************
@@ -181 +321 @@
     PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
     PM_FIT_SPLINE: fit splines to the input vector data.
-XXX: Doesn't it make more sense to do polynomial interpolation on a few
-elements of the input vector, rather than fit a polynomial to the entire
-vector?
- *****************************************************************************/
-static psVector *ScaleOverscanVector(psVector *overscanVector,
-                                     psS32 n,
-                                     void *fitSpec,
-                                     pmFit fit)
+The resulting spline or polynomial is set in the fitSpec argument.
+ *****************************************************************************/
+static psVector *ScaleOverscanVector(
+    psVector *overscanVector,
+    psS32 n,
+    void *fitSpec,
+    pmFit fit)
 {
     psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
             "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
-    //    PS_VECTOR_PRINT_F32(overscanVector);
 
     if (NULL == overscanVector) {
@@ -205 +343 @@
     //
     if (n == overscanVector->n) {
-        for (psS32 i = 0 ; i < n ; i++) {
-            newVec->data.F32[i] = overscanVector->data.F32[i];
-        }
-        return(newVec);
-    }
-    psPolynomial1D *myPoly;
-    psSpline1D *mySpline;
+        return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
+    }
     psF32 x;
-    psS32 i;
 
     if (fit == PM_FIT_POLYNOMIAL) {
         // Fit a polynomial to the old overscan vector.
-        myPoly = (psPolynomial1D *) fitSpec;
+        psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
         PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+        PS_ASSERT_POLY1D(myPoly, NULL);
         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
         if (myPoly == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n");
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n");
             return(NULL);
         }
@@ -228 +361 @@
         // of the old vector, use the fitted polynomial to determine the
         // interpolated value at that point, and set the new vector.
-        for (i=0;i<n;i++) {
+        for (psS32 i=0;i<n;i++) {
             x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
             newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
         }
     } else if (fit == PM_FIT_SPLINE) {
-        psS32 mustFreeSpline = 0;
-        // Fit a spline to the old overscan vector.
-        mySpline = (psSpline1D *) fitSpec;
-        // XXX: Does it make any sense to have a psSpline argument?
-        if (mySpline == NULL) {
-            mustFreeSpline = 1;
-        }
-
         //
         // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
@@ -246 +371 @@
         // properly when doing the spline eval.
         //
-        //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
-        mySpline = psVectorFitSpline1D(NULL, overscanVector);
+        psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
         if (mySpline == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n");
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
             return(NULL);
         }
-        //        PS_PRINT_SPLINE(mySpline);
 
         // For each element of the new vector, convert the x-ordinate to that
-        // of the old vector, use the fitted polynomial to determine the
+        // of the old vector, use the fitted spline to determine the
         // interpolated value at that point, and set the new vector.
-        for (i=0;i<n;i++) {
+        for (psS32 i=0;i<n;i++) {
             // Scale to [0 : overscanVector->n - 1]
             x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
             newVec->data.F32[i] = psSpline1DEval(mySpline, x);
         }
-        if (mustFreeSpline ==1) {
-            psFree(mySpline);
-        }
-        //        PS_VECTOR_PRINT_F32(newVec);
-
-
+
+        psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
+        if (ptrSpline != NULL) {
+            // Copy the resulting spline fit into ptrSpline.
+            PS_ASSERT_SPLINE(ptrSpline, NULL);
+            SplineCopy(ptrSpline, mySpline);
+        }
+        psFree(mySpline);
     } else {
         psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
@@ -280 +405 @@
 
 /******************************************************************************
-XXX: The SDRS does not specify type support.  F32 is implemented here.
+ *****************************************************************************/
+static psS32 GetOverscanSize(
+    psImage *inImg,
+    pmOverscanAxis overScanAxis)
+{
+    if (overScanAxis == PM_OVERSCAN_ROWS) {
+        return(inImg->numCols);
+    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
+        return(inImg->numRows);
+    } else if (overScanAxis == PM_OVERSCAN_ALL) {
+        return(1);
+    }
+    return(0);
+}
+
+/******************************************************************************
+GetOverscanAxis(in) this private routine determines the appropiate overscan
+axis from the parent cell metadata.
+ 
+XXX: Verify the READDIR corresponds with my overscan axis.
+ *****************************************************************************/
+static pmOverscanAxis GetOverscanAxis(pmReadout *in)
+{
+    psBool rc;
+    if ((in->parent != NULL) && (in->parent->concepts)) {
+        psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
+        if (rc == true) {
+            if (dir == 1) {
+                return(PM_OVERSCAN_ROWS);
+            } else if (dir == 2) {
+                return(PM_OVERSCAN_COLUMNS);
+            } else if (dir == 3) {
+                return(PM_OVERSCAN_ALL);
+            }
+        }
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN,
+             "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
+    return(PM_OVERSCAN_NONE);
+}
+
+/******************************************************************************
+psListLength(list): determine the length of a psList.
+ 
+XXX: Put this elsewhere.
+ *****************************************************************************/
+static psS32 psListLength(
+    psList *list)
+{
+    psS32 length = 0;
+    psListElem *tmpElem = (psListElem *) list->head;
+    while (NULL != tmpElem) {
+        tmpElem = tmpElem->next;
+        length++;
+    }
+    return(length);
+}
+
+/******************************************************************************
+Note: this isn't needed anymore as of psModule SDRS 12-09.
+ *****************************************************************************/
+static psBool OverscanReducePixel(
+    psImage *in,
+    psList *bias,
+    psStats *myStats)
+{
+    PS_ASSERT_PTR_NON_NULL(in, NULL);
+    PS_ASSERT_PTR_NON_NULL(bias, NULL);
+    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
+    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
+
+    // Allocate a psVector with one element per overscan image.
+    psS32 numOverscanImages = psListLength(bias);
+    psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
+    psListElem *tmpOverscan = (psListElem *) bias->head;
+    psS32 i = 0;
+    psF64 statValue;
+    //
+    // We loop through each overscan image, calculating the specified
+    // statistic on that image.
+    //
+    while (NULL != tmpOverscan) {
+        psImage *myOverscanImage = (psImage *) tmpOverscan->data;
+
+        PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
+        myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
+        if (myStats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
+            psFree(statsAll);
+            return(false);
+        }
+        if (false == p_psGetStatValue(myStats, &statValue)) {
+            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
+            psFree(statsAll);
+            return(false);
+        }
+        statsAll->data.F32[i] = statValue;
+        i++;
+        tmpOverscan = tmpOverscan->next;
+    }
+
+    //
+    // We reduce the individual stats for each overscan image to
+    // a single psF32.
+    //
+    myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
+    if (myStats == NULL) {
+        psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
+        psFree(statsAll);
+        return(false);
+    }
+    if (false == p_psGetStatValue(myStats, &statValue)) {
+        psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
+        psFree(statsAll);
+        return(false);
+    }
+
+    //
+    // Subtract the result and return.
+    //
+    ImageSubtractScalar(in, statValue);
+    psFree(statsAll);
+    return(in);
+}
+
+/******************************************************************************
+ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
+a single psImage to a column by combining all pixels from each row into a
+single pixel via requested statistic in myStats.
+ *****************************************************************************/
+static psVector *ReduceOverscanImageToCol(
+    psImage *overscanImage,
+    psStats *myStats)
+{
+    psF64 statValue;
+    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
+    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
+
+    //
+    // For each row, we store all pixels in that row into a temporary psVector,
+    // then we run psVectorStats() on that vector.
+    //
+    for (psS32 i=0;i<overscanImage->numRows;i++) {
+        for (psS32 j=0;j<overscanImage->numCols;j++) {
+            tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
+        }
+
+        psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
+        if (rc == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
+            return(NULL);
+        }
+
+        if (false ==  p_psGetStatValue(rc, &statValue)) {
+            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
+            return(NULL);
+        }
+
+        tmpCol->data.F32[i] = (psF32) statValue;
+    }
+    psFree(tmpRow);
+
+    return(tmpCol);
+}
+
+/******************************************************************************
+ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
+a single psImage to a row by combining all pixels from each column into a
+single pixel via requested statistic in myStats.
+ *****************************************************************************/
+static psVector *ReduceOverscanImageToRow(
+    psImage *overscanImage,
+    psStats *myStats)
+{
+    psF64 statValue;
+    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
+    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
+
+    //
+    // For each column, we store all pixels in that column into a temporary psVector,
+    // then we run psVectorStats() on that vector.
+    //
+    for (psS32 i=0;i<overscanImage->numCols;i++) {
+        for (psS32 j=0;j<overscanImage->numRows;j++) {
+            tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
+        }
+
+        psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
+        if (rc == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
+            return(NULL);
+        }
+
+        if (false ==  p_psGetStatValue(rc, &statValue)) {
+            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
+            return(NULL);
+        }
+
+        tmpRow->data.F32[i] = (psF32) statValue;
+    }
+    psFree(tmpCol);
+
+    return(tmpRow);
+}
+
+/******************************************************************************
+OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
+overscan images (in bias) and reduces them to a single psVector via the
+specified psStats struct.  The vector is then scaled to the length or the
+row/column in inImg.
+ *****************************************************************************/
+static psVector* OverscanReduce(
+    psImage *inImg,
+    pmOverscanAxis overScanAxis,
+    psList *bias,
+    void *fitSpec,
+    pmFit fit,
+    psStats *myStats)
+{
+    if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
+        psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
+        return(NULL);
+    }
+    PS_ASSERT_PTR_NON_NULL(inImg, NULL);
+    PS_ASSERT_PTR_NON_NULL(bias, NULL);
+    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
+    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
+    //
+    // Allocate a psVector for the output of this routine.
+    //
+    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
+    psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
+
+    //
+    // Allocate an array of psVectors with one psVector per element of the
+    // final oversan column vector.  These psVectors will be used with
+    // psStats to reduce the multiple elements from each overscan column
+    // vector to a single final column vector.
+    //
+    psS32 numOverscanImages = psListLength(bias);
+    psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
+    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+        overscanVectors[i] = NULL;
+    }
+
+    //
+    // We iterate through the list of overscan images.  For each image,
+    // we reduce it to a single column or row.  Save the overscan vector
+    // in overscanVectors[].
+    //
+    psListElem *tmpOverscan = (psListElem *) bias->head;
+    psS32 overscanID = 0;
+    while (tmpOverscan != NULL) {
+        psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
+        if (overScanAxis == PM_OVERSCAN_ROWS) {
+            overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
+        } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
+            overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
+        }
+
+        tmpOverscan = tmpOverscan->next;
+        overscanID++;
+    }
+
+    //
+    // For each overscan vector, if necessary, we scale that column or
+    // row to vecSize.  Note: we should have already ensured that the
+    // fit is poly or spline.
+    //
+    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+        psVector *tmpOverscanVector = overscanVectors[i];
+
+        if (tmpOverscanVector->n != vecSize) {
+            overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
+            if (overscanVectors[i] == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
+                for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+                    psFree(overscanVectors[i]);
+                }
+                psFree(overscanVectors);
+                psFree(tmpOverscanVector);
+                return(NULL);
+            }
+            psFree(tmpOverscanVector);
+        }
+    }
+
+    //
+    // We collect all elements in the overscan vectors for the various
+    // overscan images into a single psVector (tmpVec).  Then we call
+    // psStats on that vector to determine the final values for the
+    // overscan vector.
+    //
+    psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
+    psF64 statValue;
+    for (psS32 i = 0 ; i < vecSize ; i++) {
+        // Collect the i-th elements from each overscan vector into a single vector.
+        for (psS32 j = 0 ; j < numOverscanImages ; j++) {
+            tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
+        }
+
+        if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
+            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+                psFree(overscanVectors[i]);
+            }
+            psFree(overscanVectors);
+            psFree(tmpVec);
+            return(NULL);
+        }
+        if (false == p_psGetStatValue(myStats, &statValue)) {
+            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
+            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+                psFree(overscanVectors[i]);
+            }
+            psFree(overscanVectors);
+            psFree(tmpVec);
+            return(NULL);
+        }
+
+        overscanVector->data.F32[i] = (psF32) statValue;
+    }
+
+    //
+    // We're done.  Free the intermediate overscan vectors.
+    //
+    psFree(tmpVec);
+    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
+        psFree(overscanVectors[i]);
+    }
+    psFree(overscanVectors);
+
+    //
+    // Return the computed overscanVector
+    //
+    return(overscanVector);
+}
+
+/******************************************************************************
+RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
+takes groups of nBinOrig elements in the input vector, combines them into a
+single pixel via myStats and psVectorStats(), and then outputs a vector of
+those pixels.
+ *****************************************************************************/
+static psS32 RebinOverscanVector(
+    psVector *overscanVector,
+    psS32 nBinOrig,
+    psStats *myStats)
+{
+    psF64 statValue;
+    psS32 nBin;
+    if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
+        psS32 numBins = 1+((overscanVector->n)/nBinOrig);
+        psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
+        psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
+
+        for (psS32 i=0;i<numBins;i++) {
+            for(psS32 j=0;j<nBinOrig;j++) {
+                if (overscanVector->n > ((i*nBinOrig)+j)) {
+                    binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
+                } else {
+                    // XXX: we get here if nBinOrig does not evenly divide
+                    // the overscanVector vector.  This is the last bin.  Should
+                    // we change the binVec->n to acknowledge that?
+                    binVec->n = j;
+                }
+            }
+            psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
+            if (rc == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
+                return(-1);
+            }
+            if (false ==  p_psGetStatValue(rc, &statValue)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
+                return(-1);
+            }
+            myBin->data.F32[i] = statValue;
+        }
+
+        // Change the effective size of overscanVector.
+        overscanVector->n = numBins;
+        for (psS32 i=0;i<numBins;i++) {
+            overscanVector->data.F32[i] = myBin->data.F32[i];
+        }
+        psFree(binVec);
+        psFree(myBin);
+        nBin = nBinOrig;
+    } else {
+        nBin = 1;
+    }
+
+    return(nBin);
+}
+
+/******************************************************************************
+FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
+nBin):  this private routine fits a psPolynomial or psSpline to the overscan
+vector.  It then creates a new vector, with a size determined by the input
+image, evaluates the psPolynomial or psSpline at each element in that vector,
+then returns that vector.
+ *****************************************************************************/
+static psVector *FitOverscanVectorAndUnbin(
+    psImage *inImg,
+    psVector *overscanVector,
+    pmOverscanAxis overScanAxis,
+    void *fitSpec,
+    pmFit fit,
+    psS32 nBin)
+{
+    psPolynomial1D* myPoly = NULL;
+    psSpline1D *mySpline = NULL;
+    //
+    // Fit a polynomial or spline to the overscan vector.
+    //
+    if (fit == PM_FIT_POLYNOMIAL) {
+        myPoly = (psPolynomial1D *) fitSpec;
+        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+        PS_ASSERT_POLY1D(myPoly, NULL);
+        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
+        if (myPoly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
+            return(NULL);
+        }
+    } else if (fit == PM_FIT_SPLINE) {
+        mySpline = psVectorFitSpline1D(NULL, overscanVector);
+        if (mySpline == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
+            return(NULL);
+        }
+        if (fitSpec != NULL) {
+            // Copy the resulting spline fit into fitSpec.
+            psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
+            PS_ASSERT_SPLINE(ptrSpline, NULL);
+            SplineCopy(ptrSpline, mySpline);
+        }
+    }
+
+    //
+    // Evaluate the poly/spline at each pixel in the overscan row/column.
+    //
+    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
+    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
+    if ((nBin > 1) && (nBin < overscanVector->n)) {
+        for (psS32 i = 0 ; i < vecSize ; i++) {
+            if (fit == PM_FIT_POLYNOMIAL) {
+                newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
+            } else if (fit == PM_FIT_SPLINE) {
+                newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
+            }
+        }
+    } else {
+        for (psS32 i = 0 ; i < vecSize ; i++) {
+            if (fit == PM_FIT_POLYNOMIAL) {
+                newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
+            } else if (fit == PM_FIT_SPLINE) {
+                newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
+            }
+        }
+    }
+
+    psFree(mySpline);
+    psFree(overscanVector);
+    return(newVec);
+}
+
+
+
+/******************************************************************************
+UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
+routine takes a psVector overscanVector that was previously binned by a factor
+of nBin, and then expands it to its original size, duplicated elements nBin
+times for each element in the input vector overscanVector.
+ *****************************************************************************/
+static psVector *UnbinOverscanVector(
+    psImage *inImg,
+    psVector *overscanVector,
+    pmOverscanAxis overScanAxis,
+    psS32 nBin)
+{
+    psS32 vecSize;
+
+    if (overScanAxis == PM_OVERSCAN_ROWS) {
+        vecSize = inImg->numCols;
+    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
+        vecSize = inImg->numRows;
+    }
+
+    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
+    for (psS32 i = 0 ; i < vecSize ; i++) {
+        newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
+    }
+
+    psFree(overscanVector);
+    return(newVec);
+}
+
+
+/******************************************************************************
+SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
+routine subtracts the overscanVector column-wise or row-wise from inImg.
+ *****************************************************************************/
+static psImage *SubtractVectorFromImage(
+    psImage *inImg,
+    psVector *overscanVector,
+    pmOverscanAxis overScanAxis)
+{
+    //
+    // Subtract overscan vector row-wise from the image.
+    //
+    if (overScanAxis == PM_OVERSCAN_ROWS) {
+        for (psS32 i=0;i<inImg->numCols;i++) {
+            for (psS32 j=0;j<inImg->numRows;j++) {
+                inImg->data.F32[j][i]-= overscanVector->data.F32[i];
+            }
+        }
+    }
+
+    //
+    // Subtract overscan vector column-wise from the image.
+    //
+    if (overScanAxis == PM_OVERSCAN_COLUMNS) {
+        for (psS32 i=0;i<inImg->numRows;i++) {
+            for (psS32 j=0;j<inImg->numCols;j++) {
+                inImg->data.F32[i][j]-= overscanVector->data.F32[i];
+            }
+        }
+    }
+
+    return(inImg);
+}
+
+
+
+typedef enum {
+    PM_ERROR_NO_SUBTRACTION,
+    PM_WARNING_NO_SUBTRACTION,
+    PM_ERROR_NO_BIAS_SUBTRACT,
+    PM_WARNING_NO_BIAS_SUBTRACT,
+    PM_ERROR_NO_DARK_SUBTRACT,
+    PM_WARNING_NO_DARK_SUBTRACT,
+    PM_OKAY
+} pmSubtractBiasAssertStatus;
+/******************************************************************************
+AssertCodeOverscan(....) this private routine verifies that the various input
+parameters to pmSubtractBias() are correct for overscan subtraction.
+ *****************************************************************************/
+pmSubtractBiasAssertStatus AssertCodeOverscan(
+    pmReadout *in,
+    void *fitSpec,
+    pmFit fit,
+    bool overscan,
+    psStats *stat,
+    int nBinOrig,
+    const pmReadout *bias,
+    const pmReadout *dark)
+{
+
+    PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
+    PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
+    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
+    PS_WARN_PTR_NON_NULL(in->parent);
+    if (in->parent != NULL) {
+        PS_WARN_PTR_NON_NULL(in->parent->concepts);
+    }
+
+    if (overscan == true) {
+        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
+        PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
+        PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
+        PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
+        //
+        // Check the type, size of each bias image.
+        //
+        psListElem *tmpOverscan = (psListElem *) in->bias->head;
+        psS32 numOverscans = 0;
+        while (NULL != tmpOverscan) {
+            numOverscans++;
+            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
+            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
+            // XXX: Get this right with the rows and columns.
+            if (overScanAxis == PM_OVERSCAN_ROWS) {
+                if (myOverscanImage->numRows != in->image->numRows) {
+                    psLogMsg(__func__, PS_LOG_WARN,
+                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
+                             numOverscans, myOverscanImage->numCols, in->image->numRows);
+                    if (fit == PM_FIT_NONE) {
+                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
+                        return(PM_ERROR_NO_SUBTRACTION);
+                    }
+                }
+            } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
+                if (myOverscanImage->numCols != in->image->numCols) {
+                    psLogMsg(__func__, PS_LOG_WARN,
+                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
+                             numOverscans, myOverscanImage->numCols, in->image->numCols);
+                    if (fit == PM_FIT_NONE) {
+                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
+                        return(PM_ERROR_NO_SUBTRACTION);
+                    }
+                }
+            } else if (overScanAxis != PM_OVERSCAN_ALL) {
+                psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
+                return(PM_ERROR_NO_SUBTRACTION);
+            }
+            tmpOverscan = tmpOverscan->next;
+        }
+    } else {
+        if (fit != PM_FIT_NONE) {
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
+            return(PM_WARNING_NO_SUBTRACTION);
+        }
+    }
+
+    // XXX: I do not like the following spec since it's useless to specify
+    // a psSpline as the fitSpec.
+    if (0) {
+        if ((fitSpec == NULL) &&
+                ((fit != PM_FIT_NONE) || (overscan == true))) {
+            psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
+            return(PM_ERROR_NO_SUBTRACTION);
+        }
+    }
+
+    return(PM_OKAY);
+}
+
+/******************************************************************************
+AssertCodeBias(....) this private routine verifies that the various input
+parameters to pmSubtractBias() are correct for bias subtraction.
+ *****************************************************************************/
+static pmSubtractBiasAssertStatus AssertCodeBias(
+    pmReadout *in,
+    void *fitSpec,
+    pmFit fit,
+    bool overscan,
+    psStats *stat,
+    int nBinOrig,
+    const pmReadout *bias,
+    const pmReadout *dark)
+{
+    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
+        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
+        return(PM_ERROR_NO_BIAS_SUBTRACT);
+    }
+    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
+        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
+        return(PM_ERROR_NO_BIAS_SUBTRACT);
+    }
+
+    if (bias != NULL) {
+        PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
+        PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
+    }
+    return(PM_OKAY);
+}
+
+/******************************************************************************
+AssertCodeDark(....) this private routine verifies that the various input
+parameters to pmSubtractBias() are correct for dark subtraction.
+ *****************************************************************************/
+pmSubtractBiasAssertStatus AssertCodeDark(
+    pmReadout *in,
+    void *fitSpec,
+    pmFit fit,
+    bool overscan,
+    psStats *stat,
+    int nBinOrig,
+    const pmReadout *bias,
+    const pmReadout *dark)
+{
+    if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
+        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
+        return(PM_ERROR_NO_DARK_SUBTRACT);
+    }
+    if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
+        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
+        return(PM_ERROR_NO_DARK_SUBTRACT);
+    }
+
+    if (dark != NULL) {
+        PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
+        PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
+    }
+    return(PM_OKAY);
+}
+
+/******************************************************************************
+p_psDetermineTrimmedImage(): global routine: determines the region of the
+input pmReadout which will be operated on by the various detrend modules.  It
+does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
+ 
+Use it this way:
+    PS_WARN_PTR_NON_NULL(in->parent);
+    if (in->parent != NULL) {
+        PS_WARN_PTR_NON_NULL(in->parent->concepts);
+    }
+    //
+    // Determine trimmed image from metadata.
+    //
+    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
+ 
+XXX: Create a pmUtils.c file and put this routine there.
+ *****************************************************************************/
+psImage *p_psDetermineTrimmedImage(pmReadout *in)
+{
+    if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
+        return(in->image);
+    }
+
+    psBool rc = false;
+    psImage *trimmedImg = NULL;
+    psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
+                           "CELL.TRIMSEC");
+    if (rc == false) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
+        trimmedImg = in->image;
+    } else {
+        trimmedImg = psImageSubset(in->image, *trimRegion);
+    }
+
+    return(trimmedImg);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/******************************************************************************
+pmSubtractBias(....): see SDRS for complete specification.
+ 
+XXX: Code and assert type support: U16, S32, F32.
+XXX: Add trace messages.
  *****************************************************************************/
 pmReadout *pmSubtractBias(
@@ -290 +1164 @@
     int nBin,
     const pmReadout *bias,
-    const pmReadout *dark
-)
+    const pmReadout *dark)
 {
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
             "---- pmSubtractBias() begin ----\n");
-    PS_ASSERT_READOUT_NON_NULL(in, NULL);
-    PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
-    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
-
-    //
-    // If the overscans != NULL, then check the type of each image.
-    //
-    if (overscans != NULL) {
-        psListElem *tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            tmpOverscan = tmpOverscan->next;
-        }
-    }
-
-    if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
-        psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE).  Returning in image\n");
+    //
+    // Check input parameters, generate warnings and errors.
+    //
+    if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
         return(in);
     }
-
-    // Check for an unallowable pmFit.
-    if ((fit != PM_OVERSCAN_NONE) &&
-            (fit != PM_OVERSCAN_ROWS) &&
-            (fit != PM_OVERSCAN_COLUMNS) &&
-            (fit != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d).  Returning in image.\n", fit);
-        return(in);
-    }
-    // Check for an unallowable pmOverscanAxis.
-    if ((overScanAxis != PM_OVERSCAN_NONE) &&
-            (overScanAxis != PM_OVERSCAN_ROWS) &&
-            (overScanAxis != PM_OVERSCAN_COLUMNS) &&
-            (overScanAxis != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).  Returning in image.\n", overScanAxis);
-        return(in);
-    }
-    psS32 i;
-    psS32 j;
-    psS32 numBins = 0;
-    static psVector *overscanVector = NULL;
-    psVector *tmpRow = NULL;
-    psVector *tmpCol = NULL;
-    psVector *myBin = NULL;
-    psVector *binVec = NULL;
-    psListElem *tmpOverscan = NULL;
-    double statValue;
-    psImage *myOverscanImage = NULL;
-    psPolynomial1D *myPoly = NULL;
-    psSpline1D *mySpline = NULL;
-    psS32 nBin;
-
-    //
-    //  Create a static stats data structure and determine the highest
-    //  priority stats option.
-    //
-    static psStats *myStats = NULL;
-    if (myStats == NULL) {
-        myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        p_psMemSetPersistent(myStats, true);
-    }
-    if (stat != NULL) {
-        myStats->options = GenNewStatOptions(stat);
-    }
-
-
-    if (overScanAxis == PM_OVERSCAN_NONE) {
-        if (fit != PM_FIT_NONE) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE.  Proceeding to full fram subtraction.\n");
-        }
-
-        if (overscans != NULL) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL.  Proceeding to full fram subtraction.\n");
-        }
-        return(SubtractFrame(in, bias));
-    }
-
-    if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE.  Proceeding with the rest of the module.\n");
-    }
-
-
-    //
-    // We subtract each overscan region from the image data.
-    // If we get here we know that overscans != NULL.
-    //
-
-    if (overScanAxis == PM_OVERSCAN_ALL) {
-        tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            myOverscanImage = (psImage *) tmpOverscan->data;
-
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
-            if (rc == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
+    //
+    // Determine trimmed image from metadata.
+    //
+    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
+
+    //
+    // Subtract overscan frames if necessary.
+    //
+    if (overscan == true) {
+        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
+
+        //
+        //  Create a psStats data structure and determine the highest
+        //  priority stats option.
+        //
+        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+        if (stat != NULL) {
+            myStats->options = GenNewStatOptions(stat);
+        }
+
+        //
+        // Reduce overscan images to a single pixel, then subtract.
+        // This code is no longer required as of SDRS 12-09.
+        //
+        if (overScanAxis == PM_OVERSCAN_ALL) {
+            if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) {
                 return(in);
             }
-            if (false == p_psGetStatValue(myStats, &statValue)) {
-                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
+            psFree(myStats);
+        } else {
+            //
+            // Reduce the overscan images to a single overscan vector.
+            //
+            psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
+                                       in->bias, fitSpec,
+                                       fit, myStats);
+            if (overscanVector == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
+                psFree(myStats);
                 return(in);
             }
-            ImageSubtractScalar(in->image, statValue);
-
-            tmpOverscan = tmpOverscan->next;
-        }
-        return(in);
-    }
-
-    // This check is redundant with above code.
-    if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis);
-        return(in);
-    }
-
-    tmpOverscan = (psListElem *) overscans->head;
-    while (NULL != tmpOverscan) {
-        //        PS_IMAGE_PRINT_F32_HIDEF(in->image);
-        myOverscanImage = (psImage *) tmpOverscan->data;
-
-        if (overScanAxis == PM_OVERSCAN_ROWS) {
-            if (myOverscanImage->numCols != (in->image)->numCols) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n",
-                         myOverscanImage->numCols, in->image->numCols);
-            }
-
-            // We create a row vector and subtract this vector from image.
-            // XXX: Is there a better way to extract a psVector from a psImage without
-            // having to copy every element in that vector?
-            overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-
-            // For each column of the input image, loop through every row,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numCols;i++) {
-                for (j=0;j<myOverscanImage->numRows;j++) {
-                    tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
-                }
-                psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
+
+            //
+            // Rebin the overscan vector if necessary.
+            //
+            psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
+            if (newBin < 0) {
+                psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
+                psFree(myStats);
+                return(in);
+            }
+
+            //
+            // If necessary, fit a psPolynomial or psSpline to the overscan vector.
+            // Then, unbin the overscan vector to appropriate length for the in image.
+            //
+            if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
+                overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
+                if (overscanVector == NULL) {
+                    psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
+                    psFree(myStats);
                     return(in);
                 }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
+            } else {
+                overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
+            }
+
+            //
+            // Subtract the overscan vector from the input image.
+            //
+            SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
+            psFree(myStats);
+            psFree(overscanVector);
+        }
+    }
+
+    //
+    // Perform bias subtraction if necessary.
+    //
+    if (bias != NULL) {
+        if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
+            SubtractFrame(in, bias);
+        }
+    }
+
+    //
+    // Perform dark subtraction if necessary.
+    //
+    if (dark != NULL) {
+        if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
+            psBool rc;
+            psF32 scale = 0.0;
+            if (in->parent != NULL) {
+                scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
+                if (rc == false) {
+                    psLogMsg(__func__, PS_LOG_WARN,
+                             "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
                 }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpRow);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numCols) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numCols,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-            if (myOverscanImage->numRows != (in->image)->numRows) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n",
-                         myOverscanImage->numRows, in->image->numRows);
-            }
-
-            // We create a column vector and subtract this vector from image.
-            overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-
-            // For each row of the input image, loop through every column,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numRows;i++) {
-                for (j=0;j<myOverscanImage->numCols;j++) {
-                    tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
-                }
-                psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpCol);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numRows) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numRows,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        //
-        // Re-bin the overscan vector (change its length).
-        //
-        // Only if nBinOrig > 1.
-        if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
-            numBins = 1+((overscanVector->n)/nBinOrig);
-            myBin = psVectorAlloc(numBins, PS_TYPE_F32);
-            binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
-
-            for (i=0;i<numBins;i++) {
-                for(j=0;j<nBinOrig;j++) {
-                    if (overscanVector->n > ((i*nBinOrig)+j)) {
-                        binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
-                    } else {
-                        // XXX: we get here if nBinOrig does not evenly divide
-                        // the overscanVector vector.  This is the last bin.  Should
-                        // we change the binVec->n to acknowledge that?
-                        binVec->n = j;
-                    }
-                }
-                psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                myBin->data.F32[i] = statValue;
-            }
-
-            // Change the effective size of overscanVector.
-            overscanVector->n = numBins;
-            for (i=0;i<numBins;i++) {
-                overscanVector->data.F32[i] = myBin->data.F32[i];
-            }
-            psFree(binVec);
-            psFree(myBin);
-            nBin = nBinOrig;
-        } else {
-            nBin = 1;
-        }
-
-        // At this point the number of data points in overscanVector should be
-        // equal to the number of rows/columns (whatever is appropriate) in the
-        // image divided by numBins.
-        //
-
-
-        //
-        // This doesn't seem right.  The only way to do a spline fit is if,
-        // by SDRS requirements, fitSpec is not-NULL>  But in order for it
-        // to be non-NULL, someone must have called psSpline1DAlloc() with
-        // the min, max, and number of splines.
-        //
-        if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
-            //
-            // Fit a polynomial or spline to the overscan vector.
-            //
-            if (fit == PM_FIT_POLYNOMIAL) {
-                myPoly = (psPolynomial1D *) fitSpec;
-                myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
-                if (myPoly == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            } else if (fit == PM_FIT_SPLINE) {
-                // XXX: This makes no sense
-                // XXX: must free mySpline?
-                mySpline = (psSpline1D *) fitSpec;
-                mySpline = psVectorFitSpline1D(NULL, overscanVector);
-                if (mySpline == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= tmpF32;
-                    }
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= tmpF32;
-                    }
-                }
-            }
-        } else {
-            //
-            // If we get here, then no polynomials were fit to the overscan
-            // vector.  We simply subtract it, taking into account binning,
-            // from the image.
-            //
-
-            //
-            // Subtract overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-
-            //
-            // Subtract overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-        }
-
-        psFree(overscanVector);
-
-        tmpOverscan = tmpOverscan->next;
-    }
-
+            }
+            SubtractDarkFrame(in, dark, scale);
+        }
+    }
+
+    //
+    // All done.
+    //
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
             "---- pmSubtractBias() exit ----\n");
-
-    if (bias != NULL) {
-        return(SubtractFrame(in, bias));
-    }
     return(in);
 }

trunk/psModules/src/imsubtract/pmSubtractBias.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-20 23:06:24 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-15 20:09:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -47 +47 @@
     const pmReadout *bias,              ///< A possibly NULL bias pmReadout which is to be subtracted
     const pmReadout *dark               ///< A possibly NULL bias pmReadout which is to be subtracted
-)
+);
 // OLD: remove this
 //    const psList *overscans,      ///< A psList of overscan images
 //    pmOverscanAxis overScanAxis,  ///< Defines how overscans are applied
 
+/******************************************************************************
+DetermineTrimmedImageg() This private routine determines the region of the
+input pmReadout which will be operated on by the various detrend modules.  It
+does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
+ 
+XXX: Consider making a pmUtils.c file and put this routine there.
+ *****************************************************************************/
+psImage *p_psDetermineTrimmedImage(
+    pmReadout *in
+);
+
 #endif

trunk/psModules/src/imsubtract/pmSubtractSky.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-12 21:02:04 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-15 20:09:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18 +18 @@
 #include<math.h>
 #include "pslib.h"
-#include "psConstants.h"
 #include "pmSubtractSky.h"
+
+// XXX: Get rid of the.  Create pmUtils.h
+psImage *p_psDetermineTrimmedImage(
+    pmReadout *in
+);
 
 /******************************************************************************
@@ -152 +156 @@
  
 XXX: Use your brain and figure out the analytical expression.
+ 
+XXX: Why isn't it simply (xOrder+1) * (yOrder+1)?
  *****************************************************************************/
 static psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
@@ -170 +176 @@
         }
     }
-
     psTrace("SubtractSky.CalculatePolyTerms", 4,
             "Exiting CalculatePolyTerms(%d, %d) -> %d\n", xOrder, yOrder, localPolyTerms);
     return(localPolyTerms);
+
+    //    return((xOrder+1) * (yOrder+1));
 }
 
@@ -283 +290 @@
 XXX: Different trace message facilities in use here.
  *****************************************************************************/
-static psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
-        psImage *dataImage,
-        psImage *maskImage)
+static psPolynomial2D *ImageFitPolynomial(
+    psPolynomial2D *myPoly,
+    psImage *dataImage,
+    psImage *maskImage)
 {
     psTrace("SubtractSky.ImageFitPolynomial", 4,
@@ -471 +479 @@
     PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
+    PS_WARN_PTR_NON_NULL(in->parent);
+    if (in->parent != NULL) {
+        PS_WARN_PTR_NON_NULL(in->parent->concepts);
+    }
     psTrace(".psModule.pmSubtractSky", 4,
             "---- pmSubtractSky() begin ----\n");
@@ -492 +504 @@
         return(in);
     }
-    psImage *origImage = in->image;
+
+    //
+    // Determine trimmed image from metadata.
+    //
+
+    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
     psImage *binnedImage = NULL;
     psPolynomial2D *myPoly = NULL;
@@ -539 +556 @@
         // No binning is required here.  Simply create a copy of the image
         // and a mask.
-        binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
+        binnedImage = psImageCopy(binnedImage, trimmedImg, PS_TYPE_F32);
         if (binnedImage == NULL) {
             psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL.  Returning in image.\n");
@@ -559 +576 @@
         }
     } else {
-        binnedImage = psImageRebin(NULL, origImage, in->mask, 0, binFactor, stats);
+        binnedImage = psImageRebin(NULL, trimmedImg, in->mask, 0, binFactor, stats);
         if (binnedImage == NULL) {
             psError(PS_ERR_UNKNOWN, false, "psImageRebin() returned NULL.  Returning in image.\n");
@@ -677 +694 @@
     if (binFactor <= 1) {
         // The binned image is the same size as the original image.
-        for (psS32 row = 0; row < origImage->numRows ; row++) {
-            for (psS32 col = 0; col < origImage->numCols ; col++) {
-                origImage->data.F32[row][col]-= binnedImage->data.F32[row][col];
+        for (psS32 row = 0; row < trimmedImg->numRows ; row++) {
+            for (psS32 col = 0; col < trimmedImg->numCols ; col++) {
+                trimmedImg->data.F32[row][col]-= binnedImage->data.F32[row][col];
             }
         }
     } else {
-        for (psS32 row = 0; row < origImage->numRows ; row++) {
-            for (psS32 col = 0; col < origImage->numCols ; col++) {
+        for (psS32 row = 0; row < trimmedImg->numRows ; row++) {
+            for (psS32 col = 0; col < trimmedImg->numCols ; col++) {
                 // We calculate the F32 value of the pixel coordinates in the
                 // binned image and then use a pixel interpolation routine to
@@ -700 +717 @@
                                      binnedImage, binColF64, binRowF64,
                                      NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
-                origImage->data.F32[row][col]-= binPixel;
+                trimmedImg->data.F32[row][col]-= binPixel;
 
                 psTrace(".psModule.pmSubtractSky", 8,