IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8397


Ignore:
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.

Location:
trunk/psModules/src/detrend
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/detrend/pmFlatNormalize.c

    r8246 r8397  
    66#include "pmFlatNormalize.h"
    77
     8// I'm not sure that many many iterations are required, but rather suspect that the system converges within a
     9// few with absolutely no trouble (it *is* over-constrained).  For this reason, I'm putting the maximum number
     10// of iterations and tolerance as preset values.
     11#define MAXITER 10
     12#define TOLERANCE 1e-3
    813
    914// Estimate the flat-field normalisation; return the source flux in each integration
    10 psVector *pmFlatNormalize(bool *converge, // Did we converge?
    11                           psVector *chipGains, // Initial guess of the chip gains; modified for return
    12                           const psImage *fluxLevels, // Fluxes for each integration (row) and chip (col)
    13                           unsigned int maxIter, // Maximum number of iterations
    14                           float tolerance   // Tolerance level before dying
    15                          )
     15bool pmFlatNormalize(psVector **expFluxesPtr, // Flux in each exposure; modified for return
     16                     psVector **chipGainsPtr, // Initial guess of the chip gains; modified for return
     17                     const psImage *bgMatrix
     18                    )
    1619{
    17     PS_ASSERT_PTR_NON_NULL(chipGains, NULL);
    18     PS_ASSERT_PTR_NON_NULL(fluxLevels, NULL);
     20    PS_ASSERT_PTR_NON_NULL(bgMatrix, false);
     21    PS_ASSERT_IMAGE_NON_NULL(bgMatrix, false);
    1922
    20     int numSources = fluxLevels->numRows; // Number of integrations
    21     int numChips = fluxLevels->numCols; // Number of chips with which each integration is made
     23    int numExps = bgMatrix->numRows; // Number of exposures
     24    int numChips = bgMatrix->numCols; // Number of chips with which each exposure is made
    2225
    23     // Sanity checks
    24     assert(chipGains->n == numChips);
    25     assert(chipGains->type.type == PS_TYPE_F32);
    26     assert(maxIter >= 1);
    27     assert(tolerance > 0);
     26    psVector *expFluxes;                // Dereferenced version of expFluxesPtr
     27    if (expFluxesPtr) {
     28        if (*expFluxesPtr) {
     29            PS_ASSERT_VECTOR_TYPE(*expFluxesPtr, PS_TYPE_F32, false);
     30            PS_ASSERT_VECTOR_SIZE(*expFluxesPtr, numExps, false);
     31        } else {
     32            *expFluxesPtr = psVectorAlloc(numExps, PS_TYPE_F32);
     33            (*expFluxesPtr)->n = numExps;
     34        }
     35        expFluxes = psMemIncrRefCounter(*expFluxesPtr);
     36    } else {
     37        expFluxes = psVectorAlloc(numExps, PS_TYPE_F32);
     38        expFluxes->n = numExps;
     39    }
     40
     41    psVector *chipGains;                // Dereferenced version of chipGainsPtr
     42    if (chipGainsPtr) {
     43        if (*chipGainsPtr) {
     44            PS_ASSERT_VECTOR_TYPE(*chipGainsPtr, PS_TYPE_F32, false);
     45            PS_ASSERT_VECTOR_SIZE(*chipGainsPtr, numChips, false);
     46        } else {
     47            *chipGainsPtr = psVectorAlloc(numChips, PS_TYPE_F32);
     48            (*chipGainsPtr)->n = numChips;
     49            psVectorInit(*chipGainsPtr, 1.0);
     50        }
     51        chipGains = psMemIncrRefCounter(*chipGainsPtr);
     52    } else {
     53        chipGains = psVectorAlloc(numChips, PS_TYPE_F32);
     54        chipGains->n = numChips;
     55        psVectorInit(chipGains, 1.0);
     56    }
    2857
    2958    // Take the logarithms
    30     psImage *flux = psImageCopy(NULL, fluxLevels, PS_TYPE_F32); // Copy of the input flux levels matrix
    31     psImage *fluxMask = psImageAlloc(numChips, numSources, PS_TYPE_U8); // Mask for bad measurements
     59    psImage *flux = psImageCopy(NULL, bgMatrix, PS_TYPE_F32); // Copy of the input flux levels matrix
     60    psImage *fluxMask = psImageAlloc(numChips, numExps, PS_TYPE_U8); // Mask for bad measurements
    3261    psImageInit(fluxMask, 0);
    3362    psVector *gainMask = psVectorAlloc(numChips, PS_TYPE_U8); // Mask for bad gains
    3463    gainMask->n = numChips;
    3564    psVectorInit(gainMask, 0);
    36     psVector *sourceMask = psVectorAlloc(numSources, PS_TYPE_U8); // Mask for bad integrations
    37     sourceMask->n = numSources;
    38     psVectorInit(sourceMask, 0);
     65    psVector *expMask = psVectorAlloc(numExps, PS_TYPE_U8); // Mask for bad exposures
     66    expMask->n = numExps;
     67    psVectorInit(expMask, 0);
    3968    for (int i = 0; i < numChips; i++) {
    4069        // Note: the input gains are in e/ADU; we want to work with ADU/e (bg [ADU] = g [ADU/e] * f [e])
     70        // Hence the minus sign
    4171        if (isfinite(chipGains->data.F32[i]) && chipGains->data.F32[i] > 0) {
    42             chipGains->data.F32[i] = -log(chipGains->data.F32[i]);
     72            chipGains->data.F32[i] = -logf(chipGains->data.F32[i]);
    4373        } else {
    44             chipGains->data.F32[i] = 0.0; // Take a wild guess
     74            chipGains->data.F32[i] = 0.0; // Take a wild guess, gain ~ 1 e/ADU
    4575        }
    4676
    47         for (int j = 0; j < numSources; j++) {
     77        for (int j = 0; j < numExps; j++) {
    4878            if (isfinite(flux->data.F32[j][i]) && flux->data.F32[j][i] > 0) {
    49                 flux->data.F32[j][i] = log(flux->data.F32[j][i]);
     79                flux->data.F32[j][i] = logf(flux->data.F32[j][i]);
    5080            } else {
    5181                // Blank out this measurement
    … …  
    5888    // Not really sure that we need to iterate, but here we go anyway...
    5989
    60     float diff = INFINITY;             // Difference from previous iteration
    61     psVector *sourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The flux in each integration
    62     sourceFlux->n = numSources;
    63     psVector *oldSourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The fluxes in the previous iteration
    64     oldSourceFlux->n = numSources;
    65     psVectorInit(oldSourceFlux, 0.0);
    66     for (int iter = 0; iter < maxIter && diff > tolerance; iter++) {
    67         diff = 0.0;
    68 
    69         // Improve on the fluxes
    70         long numFluxes = 0;             // Number of fluxes
    71         for (int i = 0; i < numSources; i++) {
    72             if (sourceMask->data.U8[i]) {
     90    float diff = INFINITY;              // Difference from previous iteration
     91    psVector *oldExpFluxes = NULL;      // The fluxes in the previous iteration
     92    psVector *oldChipGains = NULL;      // Chip gains in the previous iteration
     93    for (int iter = 0; iter < MAXITER && diff > TOLERANCE; iter++) {
     94        // Improve on the exposure fluxes
     95        int numFluxes = 0;              // Number of fluxes
     96        for (int i = 0; i < numExps; i++) {
     97            if (expMask->data.U8[i]) {
    7398                psTrace("psModules.detrend", 7, "Flux for exposure %d is masked.\n", i);
    7499                continue;
    75100            }
    76101            numFluxes++;
    77             float sum = 0.0;           // Sum of F_ij - G_j
     102            float sum = 0.0;            // Sum of F_ij - G_j
    78103            int number = 0;             // Number of chips contributing
    79104            for (int j = 0; j < numChips; j++) {
    … …  
    84109            }
    85110            if (number > 0) {
    86                 sourceFlux->data.F32[i] = sum / (float)number;
     111                expFluxes->data.F32[i] = sum / (float)number;
    87112            } else {
    88                 sourceMask->data.U8[i] = 1;
    89                 sourceFlux->data.F32[i] = NAN;
     113                expMask->data.U8[i] = 1;
     114                expFluxes->data.F32[i] = NAN;
    90115            }
    91             psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, exp(sourceFlux->data.F32[i]));
    92         }
    93 
    94         for (int i = 0; i < numSources; i++) {
    95             if (!sourceMask->data.U8[i]) {
    96                 diff += abs((sourceFlux->data.F32[i] - oldSourceFlux->data.F32[i]) / sourceFlux->data.F32[i]);
    97             }
     116            psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, expf(expFluxes->data.F32[i]));
    98117        }
    99118
    100119        // Improve on the gains
     120        float meanGain = 0.0;           // Mean gain
     121        int numGains = 0;               // Number of gains
    101122        for (int i = 0; i < numChips; i++) {
    102123            if (gainMask->data.U8[i]) {
    … …  
    105126            float sum = 0.0;           // Sum of F_ji - S_j
    106127            int number = 0;             // Numer of sources contributing
    107             for (int j = 0; j < numSources; j++) {
     128            for (int j = 0; j < numExps; j++) {
    108129                if (!fluxMask->data.U8[j][i]) {
    109                     sum += flux->data.F32[j][i] - sourceFlux->data.F32[j];
     130                    sum += flux->data.F32[j][i] - expFluxes->data.F32[j];
    110131                    number++;
    111132                }
    … …  
    117138                chipGains->data.F32[i] = NAN;
    118139            }
    119             psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, exp(-chipGains->data.F32[i]));
     140            psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, expf(-chipGains->data.F32[i]));
     141            meanGain += expf(chipGains->data.F32[i]);
     142            numGains++;
     143        }
     144
     145        // Normalise the mean gain to unity, and measure the difference
     146        meanGain /= (float)numGains;
     147        meanGain = logf(meanGain);
     148        if (iter > 0) {
     149            diff = 0.0;
     150            for (int i = 0; i < numChips; i++) {
     151                if (gainMask->data.U8[i]) {
     152                    continue;
     153                }
     154                chipGains->data.F32[i] -= meanGain;
     155                diff += abs((chipGains->data.F32[i] - oldChipGains->data.F32[i]) / chipGains->data.F32[i]);
     156            }
     157            for (int i = 0; i < numExps; i++) {
     158                if (expMask->data.U8[i]) {
     159                    continue;
     160                }
     161                diff += abs((expFluxes->data.F32[i] - oldExpFluxes->data.F32[i]) / expFluxes->data.F32[i]);
     162            }
    120163        }
    121164
    122165        psTrace("psModules.detrend", 2, "Iteration %d: difference is %e\n", iter, diff);
    123166
    124         // Switch the old and new
    125         psVector *temp = sourceFlux;
    126         sourceFlux = oldSourceFlux;
    127         oldSourceFlux = temp;
     167        // Copy the new to the old
     168        oldChipGains = psVectorCopy(oldChipGains, chipGains, PS_TYPE_F32);
     169        oldExpFluxes = psVectorCopy(oldExpFluxes, expFluxes, PS_TYPE_F32);
    128170    }
    129171    psFree(flux);
    130172    psFree(fluxMask);
    131     psFree(oldSourceFlux);
     173    psFree(oldChipGains);
     174    psFree(oldExpFluxes);
    132175
    133176    // Un-log the vectors
    134177    for (int i = 0; i < numChips; i++) {
    135178        if (!gainMask->data.U8[i]) {
    136             chipGains->data.F32[i] = exp(chipGains->data.F32[i]);
     179            chipGains->data.F32[i] = expf(chipGains->data.F32[i]);
    137180        }
    138181    }
    139     for (int i = 0; i < numSources; i++) {
    140         if (!sourceMask->data.U8[i]) {
    141             sourceFlux->data.F32[i] = exp(sourceFlux->data.F32[i]);
     182    for (int i = 0; i < numExps; i++) {
     183        if (!expMask->data.U8[i]) {
     184            expFluxes->data.F32[i] = expf(expFluxes->data.F32[i]);
    142185        }
    143186    }
    144187    psFree(gainMask);
    145     psFree(sourceMask);
     188    psFree(expMask);
    146189
    147     if (converge) {
    148         *converge = (diff < tolerance); // Did we converge?
    149     }
     190    psFree(chipGains);
     191    psFree(expFluxes);
    150192
    151     return sourceFlux;
     193    return (diff < TOLERANCE); // Did we converge?
    152194}

  • trunk/psModules/src/detrend/pmFlatNormalize.h

    r7179 r8397  
    66// Normalise the flat-field measurements (f_ij = g_i s_j where f_ij is the flux recorded for chip i and
    77// integration j, g_i is the gain for the i-th chip, s_j is the flux of the source in the j-th integration).
    8 // Return the source flux in each integration.
    9 psVector *pmFlatNormalize(bool *converge, // Did we converge?
    10                           psVector *chipGains, // Initial guess of the chip gains; modified for return
    11                           const psImage *fluxLevels, // Fluxes for each integration (row) and chip (col)
    12                           unsigned int maxIter, // Maximum number of iterations
    13                           float tolerance // Tolerance level before dying
    14                          );
     8bool pmFlatNormalize(psVector **expFluxesPtr, // Flux in each exposure; modified for return
     9                     psVector **chipGainsPtr, // Initial guess of the chip gains; modified for return
     10                     const psImage *bgMatrix
     11                    );
    1512
    1613#endif
Note: See TracChangeset for help on using the changeset viewer.