Changeset 6921 for trunk/psModules/src/detrend/pmFringeStats.c

Timestamp:

Apr 19, 2006, 6:28:00 PM (20 years ago)

Author:

Paul Price

Message:

Fringe subtraction is working; needs clipping and refinement.

File:

• trunk/psModules/src/detrend/pmFringeStats.c (modified) (11 diffs)

trunk/psModules/src/detrend/pmFringeStats.c

@@ -3 +3 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-20 04:28:00 $
  *
  *  Copyright 2004 IfA
@@ -12 +12 @@
 #include "pmFPA.h"
 #include "pmFringeStats.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
 
 static void fringeRegionsFree(pmFringeRegions *fringe)
@@ -56 +59 @@
     stats->regions = psMemIncrRefCounter(regions);
     stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32);
+    stats->f->n = regions->x->n;
     stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32);
+    stats->df->n = regions->x->n;
 
     return stats;
@@ -77 +82 @@
     scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
     scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
+    scale->coeff->n = nFringeFrames + 1;
+    scale->coeff->n = nFringeFrames + 1;
 
     return scale;
 }
 
-bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image)
+bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, psImage *image)
 {
 
@@ -92 +99 @@
     fringe->y = psVectorRecycle(fringe->y, fringe->nRequested, PS_TYPE_F32);
     fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8);
+    fringe->x->n = fringe->y->n = fringe->mask->n = fringe->nRequested;
 
     int nX = image->numCols;
@@ -116 +124 @@
     if (!fringe->x || !fringe->y) {
         // create the fringe vectors for this image
-        pmFringeStatsCreatePoints(fringe, readout->image);
+        pmFringeRegionsCreatePoints(fringe, readout->image);
     }
 
@@ -138 +146 @@
 
     for (int i = 0; i < fringe->x->n; i++) {
-        psRegion region = psRegionSet(xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, yPt[i] + dY + 1);
+        psRegion region = psRegionSet(xPt[i] - dX, xPt[i] + dX + 1, yPt[i] - dY, yPt[i] + dY + 1);
         region = psRegionForImage(image, region);
         psImage *subImage = psImageSubset(image, region);
@@ -144 +152 @@
         psImageStats(stats, subImage, subMask, maskVal);
 
-        fPt[i] = stats->sampleMedian;
-        dfPt[i] = stats->sampleStdev;
+        fPt[i] = stats->robustMedian;
+        dfPt[i] = stats->robustStdev;
+
+        psTrace(__func__, 7, "[%d:%d,%d:%d]: %f %f\n", region.x0, region.x1, region.y0, region.y1,
+                fPt[i], dfPt[i]);
     }
 
@@ -154 +165 @@
 pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes)
 {
-    double sum;
-    psVector *v1;
-    psVector *v2;
-
-    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n);
-
-    int nCof = fringes->n + 1;
-    int nPts = science->f->n;
-
-    psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64);
-    psVector *B = psVectorAlloc(nCof, PS_TYPE_F64);
-
-    for (int i = 0; i < nCof; i++) {
-        psF32 *p1 = NULL;
+    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n); // The resultant fringe scales
+
+    int numCoeffs = fringes->n + 1;
+    int numPoints = science->f->n;
+
+    psImage *A = psImageAlloc(numCoeffs, numCoeffs, PS_TYPE_F64);
+    psVector *B = psVectorAlloc(numCoeffs, PS_TYPE_F64);
+    B->n = numCoeffs;
+
+    for (int i = 0; i < numCoeffs; i++) {
+        psVector *fringe1 = NULL;
         if (i != 0) {
-            v1 = fringes->data[i - 1];
-            p1 = v1->data.F32;
+            pmFringeStats *fringe = fringes->data[i - 1];
+            fringe1 = fringe->f;
         }
-        for (int j = 0; j < nCof; j++) {
-            psF32 *p2 = NULL;
+        for (int j = 0; j < numCoeffs; j++) {
+            psVector *fringe2 = NULL;
             if (j != 0) {
-                v2 = fringes->data[j - 1];
-                p2 = v2->data.F32;
+                pmFringeStats *fringe = fringes->data[j - 1];
+                fringe2 = fringe->f;
             }
 
-            sum = 0;
-            for (int k = 0; k < nPts; k++) {
-                psF32 f1 = (p1 == NULL) ? 1.0 : p1[k];
-                psF32 f2 = (p2 == NULL) ? 1.0 : p2[k];
+            double sum = 0.0;
+            for (int k = 0; k < numPoints; k++) {
+                psF32 f1 = (fringe1 == NULL) ? 1.0 : fringe1->data.F32[k];
+                psF32 f2 = (fringe2 == NULL) ? 1.0 : fringe2->data.F32[k];
                 sum += f1*f2;
             }
-            A->data.F32[i][j] = sum;
+            A->data.F64[i][j] = sum;
         }
 
-        sum = 0;
-        psF32 *p2 = science->f->data.F32;
-        for (int j = 0; j < nPts; j++) {
-            psF32 f1 = (p1 == NULL) ? 1.0 : p1[j];
-            psF32 f2 = p2[j];
+        double sum = 0.0;
+        for (int j = 0; j < numPoints; j++) {
+            psF32 f1 = (fringe1 == NULL) ? 1.0 : fringe1->data.F32[j];
+            psF32 f2 = science->f->data.F32[j];
             sum += f1*f2;
         }
-        B->data.F32[i] = sum;
+        B->data.F64[i] = sum;
     }
 
@@ -202 +209 @@
     }
 
-
-    for (int i = 0; i < nCof; i++) {
-        scale->coeff->data.F32[i] = B->data.F32[i];
-        scale->coeffErr->data.F32[i] = sqrt(A->data.F32[i][i]);
+    for (int i = 0; i < numCoeffs; i++) {
+        scale->coeff->data.F32[i] = B->data.F64[i];
+        scale->coeffErr->data.F32[i] = sqrt(A->data.F64[i][i]);
     }
 
@@ -219 +225 @@
     psFree(scienceStats);
 
+    psTrace(__func__, 7, "Fringe solution:\n");
+    for (int i = 0; i < fringeImages->n; i++) {
+        psTrace(__func__, 7, "%d: %f %f\n", i, scale->coeff->data.F32[i],
+                scale->coeffErr->data.F32[i]);
+    }
+
     // build the fringe correction image
     // XXX we could save data space by making the first image the output image
-    psImage *sumFringe = NULL;
+    psImage *sumFringe = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
     for (int i = 0; i < fringeImages->n; i++) {