Changeset 1880

Timestamp:

Sep 24, 2004, 11:05:00 AM (22 years ago)

Author:

gusciora

Message:

...

Location:

trunk/psModules/src

Files:

• Makefile (added)
• Makefile.am (modified) (2 diffs)
• pmSubtractBias.c (modified) (9 diffs)

trunk/psModules/src/Makefile.am

@@ -2 +2 @@
 bin_PROGRAMS = tst_pmFlatField tst_pmMaskBadPixels
 
-libpsmodule_a_SOURCES = pmFlatField.c pmMaskBadPixels.c psTest.c
-libpsmodule_a_HEADERS = pmFlatField.h pmMaskBadPixels.h psTest.h
+libpsmodule_a_SOURCES = pmFlatField.c pmMaskBadPixels.c pmSubtractBias.c psTest.c
+libpsmodule_a_HEADERS = pmFlatField.h pmMaskBadPixels.h pmSubtractBias.h psTest.h
 libpsmodule_adir = .
 
@@ -11 +11 @@
 tst_pmMaskBadPixels_SOURCES = tst_pmMaskBadPixels.c
 tst_pmMaskBadPixels_LDFLAGS = libpsmodule.a
+
+tst_pmSB_SOURCES = tst_pmMaskBadPixels.c
+tst_pmSB_LDFLAGS = libpsmodule.a

trunk/psModules/src/pmSubtractBias.c

@@ -8 +8 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-24 18:31:02 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-24 21:05:00 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -35 +35 @@
 //    PM_FIT_SPLINE                             ///< Fit cubic splines
 
+#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
+#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
 
 psReadout *psSubtractFrame(psReadout *in,
@@ -84 +86 @@
 }
 
+/******************************************************************************
+GenNewStatOptions(): this routine will take as input the options member of
+the stat data structure, determine if multiple options have been specified,
+issue a warning message if so, and return the highest priority option
+(according to the order of the if-statements in this code).
+ *****************************************************************************/
+psStatsOptions GenNewStatOptions(const psStats *stat)
+{
+    int numOptions = 0;
+    psStatsOptions opt;
+
+    if (stat->options & PS_STAT_SAMPLE_MEAN) {
+        numOptions++;
+        opt = PS_STAT_SAMPLE_MEAN;
+    }
+
+    if (stat->options & PS_STAT_SAMPLE_MEDIAN) {
+        if (numOptions == 0) {
+            opt = PS_STAT_SAMPLE_MEDIAN;
+        }
+        numOptions++;
+    }
+
+    if (stat->options & PS_STAT_CLIPPED_MEAN) {
+        if (numOptions == 0) {
+            opt = PS_STAT_CLIPPED_MEAN;
+        }
+        numOptions++;
+    }
+
+    if (stat->options & PS_STAT_ROBUST_MEAN) {
+        if (numOptions == 0) {
+            opt = PS_STAT_ROBUST_MEAN;
+        }
+        numOptions++;
+    }
+
+    if (stat->options & PS_STAT_ROBUST_MEDIAN) {
+        if (numOptions == 0) {
+            opt = PS_STAT_ROBUST_MEDIAN;
+        }
+        numOptions++;
+    }
+
+    if (stat->options & PS_STAT_ROBUST_MODE) {
+        if (numOptions == 0) {
+            opt = PS_STAT_ROBUST_MODE;
+        }
+        numOptions++;
+    }
+
+    if (numOptions != 1) {
+\\ XXX:
+        Generate warning message.
+    }
+
+    return(opt);
+}
+
+// XXX: Must scale the overscan vector
+psVector *ScaleOverscanVector(overscanVector,
+                              (in->image)->numCols,
+                              pmFit fit)
+{
+    return(overscanVector);
+}
+
 
 psReadout *pmSubtractBias(psReadout *in,
@@ -98 +167 @@
     int nBin;
     int numBins;
-    static psVector *tmp = NULL;
-    static psVector *tmpRow = NULL;
-    static psVector *tmpCol = NULL;
-    static psVector *tmpMask = NULL;
+    static psVector *overscanVector = NULL;
+    psVector *tmpRow = NULL;
+    psVector *tmpCol = NULL;
+    psVector *tmpMask = NULL;
+    psVector *myBin = NULL;
+    psVector *binVec = NULL;
+    psList *tmpOverscan = NULL;
 
     if (in == NULL) {
@@ -107 +179 @@
     }
 
+    //
+    //  Determine if multiple options are specified in stat.
+    //
+    if (stat != NULL) {
+        stat->options = GenNewStatOptions(stat);
+    }
+
     if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
         //XXX: psErrorMsg()
-    }
-    if ((overscans != NULL) && (overScanAxis == PM_OVERSCAN_NONE)) {
-        //XXX: psWarningMsg()
-        return(psSubtractFrame(in, bias));
-    }
-
-    if (overScanAxis == PM_OVERSCAN_ALL) {
-        psImageStats(stat, in->image, in->mask, 0xffffffff);
-        p_psImageSubtractScalar(in->image, stat->value);
         return(in);
     }
@@ -129 +199 @@
     }
 
+    if (overScanAxis == PM_OVERSCAN_NONE) {
+        if (overscans != NULL) {
+            //XXX: psWarningMsg()
+        }
+        return(psSubtractFrame(in, bias));
+    }
+
+    if (overScanAxis == PM_OVERSCAN_ALL) {
+        psImageStats(stat, in->image, in->mask, 0xffffffff);
+        p_psImageSubtractScalar(in->image, stat->value);
+        return(in);
+    }
+
     // XXX: Is there a better way to extract a psVector from a psImage without
     // having to copy every element in that vector?
@@ -135 +218 @@
     // or PM_OVERSCAN_COL?
 
+    // XXX: How to use multiple overscans?  Currently, only the first is used.
     if ((overScanAxis == PM_OVERSCAN_ROW) || (overScanAxis == PM_OVERSCAN_COL)) {
+        myOverscanImage = (psImage *) overscans->data;
+
         if (overScanAxis == PM_OVERSCAN_ROW) {
-            tmp = psVectorAlloc((in->image)->numCols, PS_TYPE_F32);
+            overscanVector = psVectorAlloc((myOverscanImage->image)->numCols, PS_TYPE_F32);
+            for (i=0;i<overscanVector->n;i++) {
+                overscanVector->data.F32[i] = 0.0;
+            }
             tmpRow = psVectorAlloc((in->image)->numCols, PS_TYPE_F32);
             tmpMask = psVectorAlloc((in->image)->numCols, PS_TYPE_U8);
-            for (i=0;i<(in->image)->numCols;i++) {
-                for (j=0;j<(in->image)->numRows;j++) {
-                    tmpRow->data.F32[j] = (in->image)->data.F32[i][j];
-                    tmpMask->data.U8[j] = (in->mask)->data.U8[i][j];
-                }
-                stat = psVectorStats(stat, tmpRow, tmpMask, 0xffffffff);
-                tmp->data.F32[i] = stat->value;
-            }
+
+            tmpOverscan = (psList *) overscans;
+            while (tmpOverscan) {
+                myOverscanImage = (psImage *) tmpOverscan->data;
+                for (i=0;i<(myOverscanImage->image)->numCols;i++) {
+                    for (j=0;j<(myOverscanImage->image)->numRows;j++) {
+                        tmpRow->data.F32[j] = (myOverscanImage->image)->data.F32[i][j];
+                        tmpMask->data.U8[j] = (myOverscanImage->mask)->data.U8[i][j];
+                    }
+                    stat = psVectorStats(stat, tmpRow, tmpMask, 0xffffffff);
+                    overscanVector->data.F32[i] = stat->value;
+                }
+                //                tmpOverscan = tmpOverscan->next;
+tmpOverscan = NULL:
+                          }
+
+                          psFree(tmpRow);
+            psFree(tmpMask);
+            ScaleOverscanVector(overscanVector, (in->image)->numCols, fit);
         }
         if (overScanAxis == PM_OVERSCAN_COL) {
-            tmp = psVectorAlloc((in->image)->numRows, PS_TYPE_F32);
+            overscanVector = psVectorAlloc((myOverscanImage->image)->numRows, PS_TYPE_F32);
+            for (i=0;i<overscanVector->n;i++) {
+                overscanVector->data.F32[i] = 0.0;
+            }
             tmpcol = psVectorAlloc((in->image)->numRows, PS_TYPE_F32);
             tmpMask = psVectorAlloc((in->image)->nuRows, PS_TYPE_U8);
-            for (i=0;i<(in->image)->numRows;i++) {
-                for (j=0;j<(in->image)->numCols;j++) {
-                    tmpRow->data.F32[j] = (in->image)->data.F32[i][j];
-                    tmpMask->data.U8[j] = (in->mask)->data.U8[i][j];
-                }
-                stat = psVectorStats(stat, tmpRow, tmpMask, 0xffffffff);
-                tmp->data.F32[i] = stat->value;
-            }
-        }
-
-        if (nBin > 0) {
-            numBins = 1+((tmp->n)/nBin);
+
+            tmpOverscan = (psList *) overscans;
+            while (tmpOverscan) {
+                for (i=0;i<(myOverscanImage->image)->numRows;i++) {
+                    for (j=0;j<(myOverscanImage->image)->numCols;j++) {
+                        tmpRow->data.F32[j] = (myOverscanImage->image)->data.F32[i][j];
+                        tmpMask->data.U8[j] = (myOverscanImage->mask)->data.U8[i][j];
+                    }
+                    stat = psVectorStats(stat, tmpRow, tmpMask, 0xffffffff);
+                    overscanVector->data.F32[i] = stat->value;
+                }
+                //                tmpOverscan = tmpOverscan->next;
+tmpOverscan = NULL:
+                          }
+                          psFree(tmpCol);
+            psFree(tmpMask);
+            ScaleOverscanVector(overscanVector, (in->image)->numRows, fit);
+        }
+
+        if ((nBin > 0) && (nBin < overscanVector->n)) {
+            numBins = 1+((overscanVector->n)/nBin);
             psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
+            binVec = psVectorAlloc(nBin, PS_TYPE_F32);
 
             for (i=0;i<numBins;i++) {
-                // XXX: do this
-            }
-            // XXX: resize tmp to myBin, and store myBin elements there.
+                for(j=0;j<nBin;j++) {
+                    if (overscanVector->n > ((i*nBin)+j)) {
+                        binVec->data.F32[j] = overscanVector->data.F32[(i*nBin)+j];
+                    } else {
+                        // XXX: we get here if nBin 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, 0);
+                myBin->data.F32[i] = stat->value;
+            }
+            // 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);
         } else {
             nBin = 1;
         }
 
-        if (!((fitSpec == NULL) || (fit ==PM_FIT_NONE))) {
-            fit the vector to a linear or cubic spline.
+        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((float) i, myPoly);
+                }
+
+            } else if (fit == PM_FIT_SPLINE) {
+                mySpline = (psSPline1D *) fitSpec;
+                // XXX: What is the point of doing this?
+            }
         }
 
@@ -182 +322 @@
             for (i=0;i<(in->image)->numCols;i++) {
                 for (j=0;j<(in->image)->numRows;j++) {
-                    (in->image)->data.F32[i][j]-= tmp->data.F32[j/nBin];
+                    (in->image)->data.F32[i][j]-= overscanVector->data.F32[j/nBin];
                 }
             }
@@ -190 +330 @@
             for (i=0;i<(in->image)->numRows;i++) {
                 for (j=0;j<(in->image)->numCols;j++) {
-                    (in->image)->data.F32[i][j]-= tmp->data.F32[j/nBin];
-                }
-            }
-        }
+                    (in->image)->data.F32[i][j]-= overscanVector->data.F32[j/nBin];
+                }
+            }
+        }
+        psFree(overscanVector);
         return(in);
     }