Index: trunk/psModules/src/detrend/pmFlatNormalize.c =================================================================== --- trunk/psModules/src/detrend/pmFlatNormalize.c (revision 7025) +++ trunk/psModules/src/detrend/pmFlatNormalize.c (revision 7060) @@ -7,25 +7,23 @@ -// 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,17 +42,18 @@ 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,5 +69,5 @@ 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,5 +92,5 @@ 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,7 +105,8 @@ } } + psFree(flux); psFree(fluxMask); - // Un-log everything + // Un-log the vectors for (int i = 0; i < numChips; i++) { if (!gainMask->data.U8[i]) { @@ -122,4 +122,8 @@ psFree(sourceMask); - return (diff <= tolerance); // Did we converge? + if (converge) { + *converge = (diff < tolerance); // Did we converge? + } + + return sourceFlux; } Index: trunk/psModules/src/detrend/pmFlatNormalize.h =================================================================== --- trunk/psModules/src/detrend/pmFlatNormalize.h (revision 7025) +++ trunk/psModules/src/detrend/pmFlatNormalize.h (revision 7060) @@ -5,11 +5,11 @@ // Normalise the flat-field measurements (f_ij = g_i s_j where f_ij is the flux recorded for chip i and // 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). -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 - ); - +// 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 + ); #endif Index: trunk/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 7025) +++ trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 7060) @@ -24,6 +24,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-19 20:37:35 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-04 03:57:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -39,5 +39,5 @@ /** Mask values */ typedef enum { - PM_MASK_CLEAR = 0x00, ///< The pixel is a charge trap. + PM_MASK_CLEAR = 0x00, ///< The pixel is not masked PM_MASK_TRAP = 0x01, ///< The pixel is a charge trap. PM_MASK_BADCOL = 0x02, ///< The pixel is a bad column.