Changeset 7060 for trunk/psModules/src/detrend/pmFlatNormalize.c

Timestamp:

May 3, 2006, 5:57:32 PM (20 years ago)

Author:

Paul Price

Message:

Updates for ppMerge

File:

• trunk/psModules/src/detrend/pmFlatNormalize.c (modified) (6 diffs)

trunk/psModules/src/detrend/pmFlatNormalize.c

@@ -7 +7 @@
 
 
-// Estimate the flat-field normalisation
-bool pmFlatNormalize(psVector *sourceFlux, // The source flux in each image; modified for return
-                     psVector *chipGains, // Initial guess of the chip gains; modified for return
-                     psImage *fluxLevels, // Fluxes for each integration (row) and chip (col); modified
-                     unsigned int maxIter, // Maximum number of iterations
-                     double tolerance   // Tolerance level before dying
-                    )
+// Estimate the flat-field normalisation; return the source flux in each integration
+psVector *pmFlatNormalize(bool *converge, // Did we converge?
+                          psVector *chipGains, // Initial guess of the chip gains; modified for return
+                          const psImage *fluxLevels, // Fluxes for each integration (row) and chip (col)
+                          unsigned int maxIter, // Maximum number of iterations
+                          double tolerance   // Tolerance level before dying
+                         )
 {
-    int numSources = sourceFlux->n;     // Number of measurements made
-    int numChips = chipGains->n;        // Number of chips with which each measurement is made
+    int numSources = fluxLevels->numRows; // Number of integrations
+    int numChips = fluxLevels->numCols; // Number of chips with which each integration is made
 
     // Sanity checks
-    assert(fluxLevels->numCols == numSources);
-    assert(fluxLevels->numRows == numChips);
-    assert(sourceFlux->type.type == PS_TYPE_F64);
+    assert(chipGains->n == numChips);
     assert(chipGains->type.type == PS_TYPE_F64);
-    assert(fluxLevels->type.type == PS_TYPE_F64);
     assert(maxIter >= 1);
     assert(tolerance > 0);
 
     // Take the logarithms
+    psImage *flux = psImageCopy(NULL, fluxLevels, PS_TYPE_F64); // Copy of the input flux levels matrix
     psImage *fluxMask = psImageAlloc(numSources, numChips, PS_TYPE_U8); // Mask for bad measurements
     psImageInit(fluxMask, 0);
@@ -44 +42 @@
 
         for (int j = 0; j < numSources; j++) {
-            if (isfinite(fluxLevels->data.F64[j][i]) && fluxLevels->data.F64[j][i] > 0) {
-                fluxLevels->data.F64[j][i] = log(fluxLevels->data.F64[j][i]);
+            if (isfinite(flux->data.F64[j][i]) && flux->data.F64[j][i] > 0) {
+                flux->data.F64[j][i] = log(flux->data.F64[j][i]);
             } else {
                 // Blank out this measurement
                 fluxMask->data.U8[j][i] = 1;
-                fluxLevels->data.F64[j][i] = NAN;
+                flux->data.F64[j][i] = NAN;
             }
         }
     }
-    // Don't need to initialise sourceFlux, since that is changed immediately
 
 
     double diff = INFINITY;             // Difference from previous iteration
+    psVector *sourceFlux = psVectorAlloc(numSources, PS_TYPE_F64); // The flux in each integration
+    // Don't need to initialise sourceFlux, since that is changed immediately
     for (int iter = 0; iter < maxIter && diff > tolerance; iter++) {
         diff = 0.0;
@@ -70 +69 @@
             for (int j = 0; j < numChips; j++) {
                 if (!gainMask->data.U8[j] && !fluxMask->data.U8[i][j]) {
-                    sum += fluxLevels->data.F64[i][j] - chipGains->data.F64[j];
+                    sum += flux->data.F64[i][j] - chipGains->data.F64[j];
                     number++;
                 }
@@ -93 +92 @@
             for (int j = 0; j < numSources; j++) {
                 if (!fluxMask->data.U8[j][i]) {
-                    sum += fluxLevels->data.F64[j][i] - sourceFlux->data.F64[j];
+                    sum += flux->data.F64[j][i] - sourceFlux->data.F64[j];
                     number++;
                 }
@@ -106 +105 @@
         }
     }
+    psFree(flux);
     psFree(fluxMask);
 
-    // Un-log everything
+    // Un-log the vectors
     for (int i = 0; i < numChips; i++) {
         if (!gainMask->data.U8[i]) {
@@ -122 +122 @@
     psFree(sourceMask);
 
-    return (diff <= tolerance);         // Did we converge?
+    if (converge) {
+        *converge = (diff < tolerance); // Did we converge?
+    }
+
+    return sourceFlux;
 }