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

Timestamp:

Aug 16, 2006, 2:49:01 PM (20 years ago)

Author:

Paul Price

Message:

pmFlatNormalize now normalizes the mean of the gains to unity. Removed external iteration control, since I'm somewhat confident that the system is stable. Fixed bug that produced memory errors if the number of iterations was odd (due to swapping pointers around). Output both flux and gain solutions.

File:

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

trunk/psModules/src/detrend/pmFlatNormalize.c

@@ -6 +6 @@
 #include "pmFlatNormalize.h"
 
+// I'm not sure that many many iterations are required, but rather suspect that the system converges within a
+// few with absolutely no trouble (it *is* over-constrained).  For this reason, I'm putting the maximum number
+// of iterations and tolerance as preset values.
+#define MAXITER 10
+#define TOLERANCE 1e-3
 
 // 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
-                          float tolerance   // Tolerance level before dying
-                         )
+bool pmFlatNormalize(psVector **expFluxesPtr, // Flux in each exposure; modified for return
+                     psVector **chipGainsPtr, // Initial guess of the chip gains; modified for return
+                     const psImage *bgMatrix
+                    )
 {
-    PS_ASSERT_PTR_NON_NULL(chipGains, NULL);
-    PS_ASSERT_PTR_NON_NULL(fluxLevels, NULL);
+    PS_ASSERT_PTR_NON_NULL(bgMatrix, false);
+    PS_ASSERT_IMAGE_NON_NULL(bgMatrix, false);
 
-    int numSources = fluxLevels->numRows; // Number of integrations
-    int numChips = fluxLevels->numCols; // Number of chips with which each integration is made
+    int numExps = bgMatrix->numRows; // Number of exposures
+    int numChips = bgMatrix->numCols; // Number of chips with which each exposure is made
 
-    // Sanity checks
-    assert(chipGains->n == numChips);
-    assert(chipGains->type.type == PS_TYPE_F32);
-    assert(maxIter >= 1);
-    assert(tolerance > 0);
+    psVector *expFluxes;                // Dereferenced version of expFluxesPtr
+    if (expFluxesPtr) {
+        if (*expFluxesPtr) {
+            PS_ASSERT_VECTOR_TYPE(*expFluxesPtr, PS_TYPE_F32, false);
+            PS_ASSERT_VECTOR_SIZE(*expFluxesPtr, numExps, false);
+        } else {
+            *expFluxesPtr = psVectorAlloc(numExps, PS_TYPE_F32);
+            (*expFluxesPtr)->n = numExps;
+        }
+        expFluxes = psMemIncrRefCounter(*expFluxesPtr);
+    } else {
+        expFluxes = psVectorAlloc(numExps, PS_TYPE_F32);
+        expFluxes->n = numExps;
+    }
+
+    psVector *chipGains;                // Dereferenced version of chipGainsPtr
+    if (chipGainsPtr) {
+        if (*chipGainsPtr) {
+            PS_ASSERT_VECTOR_TYPE(*chipGainsPtr, PS_TYPE_F32, false);
+            PS_ASSERT_VECTOR_SIZE(*chipGainsPtr, numChips, false);
+        } else {
+            *chipGainsPtr = psVectorAlloc(numChips, PS_TYPE_F32);
+            (*chipGainsPtr)->n = numChips;
+            psVectorInit(*chipGainsPtr, 1.0);
+        }
+        chipGains = psMemIncrRefCounter(*chipGainsPtr);
+    } else {
+        chipGains = psVectorAlloc(numChips, PS_TYPE_F32);
+        chipGains->n = numChips;
+        psVectorInit(chipGains, 1.0);
+    }
 
     // Take the logarithms
-    psImage *flux = psImageCopy(NULL, fluxLevels, PS_TYPE_F32); // Copy of the input flux levels matrix
-    psImage *fluxMask = psImageAlloc(numChips, numSources, PS_TYPE_U8); // Mask for bad measurements
+    psImage *flux = psImageCopy(NULL, bgMatrix, PS_TYPE_F32); // Copy of the input flux levels matrix
+    psImage *fluxMask = psImageAlloc(numChips, numExps, PS_TYPE_U8); // Mask for bad measurements
     psImageInit(fluxMask, 0);
     psVector *gainMask = psVectorAlloc(numChips, PS_TYPE_U8); // Mask for bad gains
     gainMask->n = numChips;
     psVectorInit(gainMask, 0);
-    psVector *sourceMask = psVectorAlloc(numSources, PS_TYPE_U8); // Mask for bad integrations
-    sourceMask->n = numSources;
-    psVectorInit(sourceMask, 0);
+    psVector *expMask = psVectorAlloc(numExps, PS_TYPE_U8); // Mask for bad exposures
+    expMask->n = numExps;
+    psVectorInit(expMask, 0);
     for (int i = 0; i < numChips; i++) {
         // Note: the input gains are in e/ADU; we want to work with ADU/e (bg [ADU] = g [ADU/e] * f [e])
+        // Hence the minus sign
         if (isfinite(chipGains->data.F32[i]) && chipGains->data.F32[i] > 0) {
-            chipGains->data.F32[i] = -log(chipGains->data.F32[i]);
+            chipGains->data.F32[i] = -logf(chipGains->data.F32[i]);
         } else {
-            chipGains->data.F32[i] = 0.0; // Take a wild guess
+            chipGains->data.F32[i] = 0.0; // Take a wild guess, gain ~ 1 e/ADU
         }
 
-        for (int j = 0; j < numSources; j++) {
+        for (int j = 0; j < numExps; j++) {
             if (isfinite(flux->data.F32[j][i]) && flux->data.F32[j][i] > 0) {
-                flux->data.F32[j][i] = log(flux->data.F32[j][i]);
+                flux->data.F32[j][i] = logf(flux->data.F32[j][i]);
             } else {
                 // Blank out this measurement
@@ -58 +88 @@
     // Not really sure that we need to iterate, but here we go anyway...
 
-    float diff = INFINITY;             // Difference from previous iteration
-    psVector *sourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The flux in each integration
-    sourceFlux->n = numSources;
-    psVector *oldSourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The fluxes in the previous iteration
-    oldSourceFlux->n = numSources;
-    psVectorInit(oldSourceFlux, 0.0);
-    for (int iter = 0; iter < maxIter && diff > tolerance; iter++) {
-        diff = 0.0;
-
-        // Improve on the fluxes
-        long numFluxes = 0;             // Number of fluxes
-        for (int i = 0; i < numSources; i++) {
-            if (sourceMask->data.U8[i]) {
+    float diff = INFINITY;              // Difference from previous iteration
+    psVector *oldExpFluxes = NULL;      // The fluxes in the previous iteration
+    psVector *oldChipGains = NULL;      // Chip gains in the previous iteration
+    for (int iter = 0; iter < MAXITER && diff > TOLERANCE; iter++) {
+        // Improve on the exposure fluxes
+        int numFluxes = 0;              // Number of fluxes
+        for (int i = 0; i < numExps; i++) {
+            if (expMask->data.U8[i]) {
                 psTrace("psModules.detrend", 7, "Flux for exposure %d is masked.\n", i);
                 continue;
             }
             numFluxes++;
-            float sum = 0.0;           // Sum of F_ij - G_j
+            float sum = 0.0;            // Sum of F_ij - G_j
             int number = 0;             // Number of chips contributing
             for (int j = 0; j < numChips; j++) {
@@ -84 +109 @@
             }
             if (number > 0) {
-                sourceFlux->data.F32[i] = sum / (float)number;
+                expFluxes->data.F32[i] = sum / (float)number;
             } else {
-                sourceMask->data.U8[i] = 1;
-                sourceFlux->data.F32[i] = NAN;
+                expMask->data.U8[i] = 1;
+                expFluxes->data.F32[i] = NAN;
             }
-            psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, exp(sourceFlux->data.F32[i]));
-        }
-
-        for (int i = 0; i < numSources; i++) {
-            if (!sourceMask->data.U8[i]) {
-                diff += abs((sourceFlux->data.F32[i] - oldSourceFlux->data.F32[i]) / sourceFlux->data.F32[i]);
-            }
+            psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, expf(expFluxes->data.F32[i]));
         }
 
         // Improve on the gains
+        float meanGain = 0.0;           // Mean gain
+        int numGains = 0;               // Number of gains
         for (int i = 0; i < numChips; i++) {
             if (gainMask->data.U8[i]) {
@@ -105 +126 @@
             float sum = 0.0;           // Sum of F_ji - S_j
             int number = 0;             // Numer of sources contributing
-            for (int j = 0; j < numSources; j++) {
+            for (int j = 0; j < numExps; j++) {
                 if (!fluxMask->data.U8[j][i]) {
-                    sum += flux->data.F32[j][i] - sourceFlux->data.F32[j];
+                    sum += flux->data.F32[j][i] - expFluxes->data.F32[j];
                     number++;
                 }
@@ -117 +138 @@
                 chipGains->data.F32[i] = NAN;
             }
-            psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, exp(-chipGains->data.F32[i]));
+            psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, expf(-chipGains->data.F32[i]));
+            meanGain += expf(chipGains->data.F32[i]);
+            numGains++;
+        }
+
+        // Normalise the mean gain to unity, and measure the difference
+        meanGain /= (float)numGains;
+        meanGain = logf(meanGain);
+        if (iter > 0) {
+            diff = 0.0;
+            for (int i = 0; i < numChips; i++) {
+                if (gainMask->data.U8[i]) {
+                    continue;
+                }
+                chipGains->data.F32[i] -= meanGain;
+                diff += abs((chipGains->data.F32[i] - oldChipGains->data.F32[i]) / chipGains->data.F32[i]);
+            }
+            for (int i = 0; i < numExps; i++) {
+                if (expMask->data.U8[i]) {
+                    continue;
+                }
+                diff += abs((expFluxes->data.F32[i] - oldExpFluxes->data.F32[i]) / expFluxes->data.F32[i]);
+            }
         }
 
         psTrace("psModules.detrend", 2, "Iteration %d: difference is %e\n", iter, diff);
 
-        // Switch the old and new
-        psVector *temp = sourceFlux;
-        sourceFlux = oldSourceFlux;
-        oldSourceFlux = temp;
+        // Copy the new to the old
+        oldChipGains = psVectorCopy(oldChipGains, chipGains, PS_TYPE_F32);
+        oldExpFluxes = psVectorCopy(oldExpFluxes, expFluxes, PS_TYPE_F32);
     }
     psFree(flux);
     psFree(fluxMask);
-    psFree(oldSourceFlux);
+    psFree(oldChipGains);
+    psFree(oldExpFluxes);
 
     // Un-log the vectors
     for (int i = 0; i < numChips; i++) {
         if (!gainMask->data.U8[i]) {
-            chipGains->data.F32[i] = exp(chipGains->data.F32[i]);
+            chipGains->data.F32[i] = expf(chipGains->data.F32[i]);
         }
     }
-    for (int i = 0; i < numSources; i++) {
-        if (!sourceMask->data.U8[i]) {
-            sourceFlux->data.F32[i] = exp(sourceFlux->data.F32[i]);
+    for (int i = 0; i < numExps; i++) {
+        if (!expMask->data.U8[i]) {
+            expFluxes->data.F32[i] = expf(expFluxes->data.F32[i]);
         }
     }
     psFree(gainMask);
-    psFree(sourceMask);
+    psFree(expMask);
 
-    if (converge) {
-        *converge = (diff < tolerance); // Did we converge?
-    }
+    psFree(chipGains);
+    psFree(expFluxes);
 
-    return sourceFlux;
+    return (diff < TOLERANCE); // Did we converge?
 }