Changeset 2848 for trunk/psModules/src/pmSubtractBias.c

Timestamp:

Dec 29, 2004, 12:44:33 PM (22 years ago)

Author:

gusciora

Message:

Code.

File:

• trunk/psModules/src/pmSubtractBias.c (modified) (8 diffs)

trunk/psModules/src/pmSubtractBias.c

@@ -6 +6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-29 21:18:33 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-29 22:44:33 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -205 +205 @@
         for (i=0;i<n;i++) {
             x = ((float) i) * ((float) overscanVector->n) / ((float) n);
-            newVec->data.F32[i] = psSpline1DEval(x, mySpline);
+            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
         }
     } else {
@@ -234 +234 @@
         return(in);
     }
+
+    // Check for an unallowable pmFit.
     if ((fit != PM_OVERSCAN_NONE) &&
             (fit != PM_OVERSCAN_ROWS) &&
@@ -241 +243 @@
         return(in);
     }
+    // Check for an unallowable pmOverscanAxis.
     if ((overScanAxis != PM_OVERSCAN_NONE) &&
             (overScanAxis != PM_OVERSCAN_ROWS) &&
@@ -259 +262 @@
     double statValue;
     psImage *myOverscanImage = NULL;
-    psPolynomial1D *myPoly;
-    psSpline1D *mySpline;
+    psPolynomial1D *myPoly = NULL;
+    psSpline1D *mySpline = NULL;
     int nBin;
     //  XXX: This comment isn't right?
@@ -308 +311 @@
     }
 
+    // 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)\n", overScanAxis);
+        return(in);
+    }
 
     tmpOverscan = (psListElem *) overscans->head;
@@ -313 +321 @@
         myOverscanImage = (psImage *) tmpOverscan->data;
 
-        // XXX: Is there a better way to extract a psVector from a psImage without
-        // having to copy every element in that vector?
-        // XXX: If we get here, do we know that overScanAxis == PM_OVERSCAN_ROWS
-        // or PM_OVERSCAN_COL?
-
-        if ((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS)) {
-            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.
-                overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-                for (i=0;i<overscanVector->n;i++) {
-                    overscanVector->data.F32[i] = 0.0;
-                }
-                tmpRow = psVectorAlloc((in->image)->numRows, PS_TYPE_F32);
-
-                for (i=0;i<myOverscanImage->numCols;i++) {
-                    for (j=0;j<myOverscanImage->numRows;j++) {
-                        tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
-                    }
-                    stat = psVectorStats(stat, tmpRow, NULL, NULL, 0);
-                    p_psGetStatValue(stat, &statValue);
-                    overscanVector->data.F32[i] = statValue;
-                }
-                psFree(tmpRow);
+        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((in->image)->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];
+                }
+                stat = psVectorStats(stat, tmpRow, NULL, NULL, 0);
+                p_psGetStatValue(stat, &statValue);
+                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,
+                                                           in->image->numCols,
                                                            fitSpec, fit);
                     psFree(overscanVector);
                     overscanVector = newVec;
-                }
-            }
-
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                // We create a column vector and subtract this vector from image.
-                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);
-                }
-                overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-                for (i=0;i<overscanVector->n;i++) {
-                    overscanVector->data.F32[i] = 0.0;
-                }
-                tmpCol = psVectorAlloc((in->image)->numCols, PS_TYPE_F32);
-
-                for (i=0;i<myOverscanImage->numRows;i++) {
-                    for (j=0;j<myOverscanImage->numCols;j++) {
-                        tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
-                    }
-                    stat = psVectorStats(stat, tmpCol, NULL, NULL, 0);
-                    p_psGetStatValue(stat, &statValue);
-                    overscanVector->data.F32[i] = statValue;
-                }
-
-                psFree(tmpCol);
+                } 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.\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((in->image)->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];
+                }
+                stat = psVectorStats(stat, tmpCol, NULL, NULL, 0);
+                p_psGetStatValue(stat, &statValue);
+                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,
+                                                           in->image->numRows,
                                                            fitSpec, fit);
                     psFree(overscanVector);
                     overscanVector = newVec;
-                }
-            }
-
-            //
-            // Re-bin the overscan vector (change its length).
-            //
-            if ((nBinOrig > 0) && (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;
-                        }
-                    }
-                    stat = psVectorStats(stat, binVec, NULL, NULL, 0);
-                    p_psGetStatValue(stat, &statValue);
-                    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;
-            }
-
-            //
-            // Fit a polynomial to the overscan vector.
-            //
-            if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
-                if (fit == PM_FIT_POLYNOMIAL) {
-                    myPoly = (psPolynomial1D *) fitSpec;
-                    myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL);
-                    for (i=0;i<numBins;i++) {
-                        overscanVector->data.F32[i] = psPolynomial1DEval(myPoly, (float) i);
-                    }
-
-                } else if (fit == PM_FIT_SPLINE) {
-                    mySpline = (psSpline1D *) fitSpec;
-                    // XXX: The spline functions do not work if the number of
-                    // data points and the number of splines are not equal.
-                }
-            }
+                } 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.\n");
+                    psFree(overscanVector);
+                    return(in);
+                }
+            }
+        }
+
+        //
+        // Re-bin the overscan vector (change its length).
+        //
+        if ((nBinOrig > 0) && (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;
+                    }
+                }
+                stat = psVectorStats(stat, binVec, NULL, NULL, 0);
+                p_psGetStatValue(stat, &statValue);
+                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.
+        //
+        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, overscanVector, NULL);
+                if (myPoly == NULL) {
+                    psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.\n");
+                    psFree(overscanVector);
+                    return(in);
+                }
+            } else if (fit == PM_FIT_SPLINE) {
+                mySpline = (psSpline1D *) fitSpec;
+                mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
+                if (mySpline == NULL) {
+                    psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.\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, ((float) i) / ((float) nBin));
+                    } else if (fit == PM_FIT_SPLINE) {
+                        tmpF32 = psSpline1DEval(mySpline, ((float) i) / ((float) 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, ((float) i) / ((float) nBin));
+                    } else if (fit == PM_FIT_SPLINE) {
+                        tmpF32 = psSpline1DEval(mySpline, ((float) i) / ((float) 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.
 
             //
@@ -457 +533 @@
                 }
             }
-            psFree(overscanVector);
-        }
+        }
+
+        psFree(overscanVector);
+
         tmpOverscan = tmpOverscan->next;
     }
-
 
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,