Index: trunk/psModules/src/detrend/pmDark.c =================================================================== --- trunk/psModules/src/detrend/pmDark.c (revision 16948) +++ trunk/psModules/src/detrend/pmDark.c (revision 16949) @@ -1,4 +1,5 @@ #include #include +#include #include "psPolynomialMD.h" @@ -95,6 +96,6 @@ -bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, int iter, float rej, - psMaskType maskVal) +bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept, + int iter, float rej, psMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(output, false); @@ -109,9 +110,29 @@ int numBadInputs = 0; // Number of bad inputs psArray *values = psArrayAlloc(numInputs); + psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts + psVectorInit(roMask, 0); + psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each for (int i = 0; i < numInputs; i++) { values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32); - } - psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts - psVectorInit(roMask, 0); + if (norm) { + pmReadout *ro = inputs->data[i]; // Readout of interest + float normValue; // Normalisation value + if (!ordinateLookup(&normValue, normConcept, false, NAN, NAN, ro)) { + psWarning("Unable to find value of %s for readout %d", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + if (normValue == 0.0) { + psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + norm->data.F32[i] = 1.0 / normValue; + } + } psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept for (int i = 0; i < numOrdinates; i++) { @@ -135,7 +156,11 @@ for (int j = 0; j < numInputs; j++) { + psVector *val = values->data[j]; // Value vector for readout + if (roMask->data.U8[j]) { + val->data.F32[i] = NAN; + continue; + } + pmReadout *ro = inputs->data[j]; // Readout of interest - psVector *val = values->data[j]; // Value vector for readout - float value = NAN; // Value of ordinate if (!ordinateLookup(&value, ord->name, ord->scale, ord->min, ord->max, ro)) { @@ -148,5 +173,4 @@ } } - numInputs -= numBadInputs; if (psTraceGetLevel("psModules.detrend") > 9) { @@ -162,9 +186,10 @@ psFree(orders); int numTerms = poly->coeff->n; // Number of terms in polynomial - if (numTerms > numInputs) { + if (numTerms > numInputs - numBadInputs) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).", - numInputs, numTerms); + numInputs - numBadInputs, numTerms); psFree(values); psFree(roMask); + psFree(norm); return false; } @@ -176,4 +201,7 @@ inputs)) { psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); + psFree(values); + psFree(roMask); + psFree(norm); return false; } @@ -224,4 +252,7 @@ pixels->data.F32[r] = readout->image->data.F32[yIn][xIn]; + if (norm) { + pixels->data.F32[r] *= norm->data.F32[r]; + } if (readout->mask) { mask->data.PS_TYPE_MASK_DATA[r] = readout->mask->data.PS_TYPE_MASK_DATA[yIn][xIn]; @@ -240,4 +271,5 @@ } + psFree(norm); psFree(roMask); psFree(poly); @@ -248,4 +280,6 @@ psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); + psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, + PS_META_REPLACE, "Dark normalisation", normConcept); for (int i = 0; i < numTerms; i++) { @@ -286,4 +320,6 @@ return false; } + bool mdok; // Status of MD lookup + psString normConcept = psMetadataLookupStr(&mdok, dark->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation int numOrdinates = ordinates->n; // Number of ordinates @@ -299,4 +335,12 @@ values->data.F32[i] = value; } + float norm = NAN; // Normalisation value + if (normConcept) { + if (!ordinateLookup(&norm, normConcept, false, NAN, NAN, readout)) { + psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept); + psFree(values); + return false; + } + } psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Order for each polynomial @@ -322,4 +366,7 @@ } float value = psPolynomialMDEval(poly, values); // Value of dark current + if (normConcept) { + value *= norm; + } readout->image->data.F32[y][x] -= value; if (readout->mask && !isfinite(value)) { @@ -352,4 +399,5 @@ psArray *ordinates = NULL; // Dark ordinates, to write + const char *normConcept = NULL; // Normalisation concept psArray *chips = fpa->chips; // Component chips for (int i = 0; i < chips->n; i++) { @@ -359,15 +407,22 @@ pmCell *cell = cells->data[j]; // Cell of interest bool mdok; // Status of MD lookup - psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); + psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ords if (!mdok) { continue; } + psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Norm if (ordinates) { - if (new != ordinates) { + if (newOrd != ordinates) { psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells."); return false; } + if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) || + (!normConcept && newNorm)) { + psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells."); + return false; + } } else { - ordinates = new; + ordinates = newOrd; + normConcept = newNorm; } } @@ -379,5 +434,5 @@ } - return pmDarkWrite(fits, hdu->header, ordinates); + return pmDarkWrite(fits, hdu->header, ordinates, normConcept); } @@ -399,19 +454,27 @@ psArray *ordinates = NULL; // Dark ordinates, to write + const char *normConcept = NULL; // Normalisation concept psArray *cells = chip->cells; // Component cells for (int j = 0; j < cells->n; j++) { pmCell *cell = cells->data[j]; // Cell of interest bool mdok; // Status of MD lookup - psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); + psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ordinates if (!mdok) { continue; } + psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation if (ordinates) { - if (new != ordinates) { + if (newOrd != ordinates) { psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells."); return false; } + if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) || + (!normConcept && newNorm)) { + psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells."); + return false; + } } else { - ordinates = new; + ordinates = newOrd; + normConcept = newNorm; } } @@ -422,5 +485,5 @@ } - return pmDarkWrite(fits, hdu->header, ordinates); + return pmDarkWrite(fits, hdu->header, ordinates, normConcept); } @@ -447,9 +510,11 @@ return false; } - - return pmDarkWrite(fits, hdu->header, ordinates); -} - -bool pmDarkWrite(psFits *fits, const psMetadata *header, const psArray *ordinates) + bool mdok; // Status of MD lookup + psString normConcept = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation + + return pmDarkWrite(fits, hdu->header, ordinates, normConcept); +} + +bool pmDarkWrite(psFits *fits, psMetadata *header, const psArray *ordinates, const char *normConcept) { PS_ASSERT_FITS_NON_NULL(fits, false); @@ -472,24 +537,41 @@ } - // Format ordinates into FITS table - int numOrdinates = ordinates->n;// Number of ordinates - psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates - for (int i = 0; i < ordinates->n; i++) { - pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest - psMetadata *row = psMetadataAlloc(); // FITS table row - psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name); - psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order); - psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale); - psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min); - psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max); - table->data[i] = row; - } - - if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) { - psError(PS_ERR_IO, false, "Unable to write dark ordinates."); + if (!psMemIncrRefCounter((psMetadata*)header)) { + header = psMetadataAlloc(); + } + psMetadataAddStr(header, PS_LIST_TAIL, PM_DARK_HEADER_NORM, PS_META_REPLACE, + "Dark normalisation concept", normConcept); + + if (ordinates->n > 0) { + // Format ordinates into FITS table + int numOrdinates = ordinates->n;// Number of ordinates + psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates + for (int i = 0; i < ordinates->n; i++) { + pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest + psMetadata *row = psMetadataAlloc(); // FITS table row + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name); + psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order); + psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale); + psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min); + psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max); + table->data[i] = row; + } + + if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) { + psError(PS_ERR_IO, false, "Unable to write dark ordinates."); + psFree(table); + psFree(header); + return false; + } psFree(table); - return false; - } - psFree(table); + } else { + // No ordinates to write to a table, so write to a blank header. + if (!psFitsWriteBlank(fits, header, PM_DARK_FITS_EXTNAME)) { + psError(PS_ERR_IO, false, "Unable to write dark header."); + psFree(header); + return false; + } + } + psFree(header); return true; @@ -507,5 +589,6 @@ } - psArray *ordinates = pmDarkRead(fits); + psString normConcept = NULL; // Normalisation concept + psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates if (!ordinates) { psError(PS_ERR_IO, false, "Unable to read dark ordinates."); @@ -521,7 +604,10 @@ psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); + psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE, + "Dark normalisation", normConcept); } } psFree(ordinates); + psFree(normConcept); return true; @@ -538,5 +624,6 @@ } - psArray *ordinates = pmDarkRead(fits); + psString normConcept = NULL; // Normalisation concept + psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates if (!ordinates) { psError(PS_ERR_IO, false, "Unable to read dark ordinates."); @@ -549,6 +636,9 @@ psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); + psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE, + "Dark normalisation", normConcept); } psFree(ordinates); + psFree(normConcept); return true; @@ -567,5 +657,6 @@ } - psArray *ordinates = pmDarkRead(fits); + psString normConcept = NULL; // Normalisation concept + psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates if (!ordinates) { psError(PS_ERR_IO, false, "Unable to read dark ordinates."); @@ -574,12 +665,17 @@ psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); + psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE, + "Dark normalisation", normConcept); psFree(ordinates); + psFree(normConcept); return true; } -psArray *pmDarkRead(psFits *fits) -{ - PS_ASSERT_FITS_NON_NULL(fits, false); +psArray *pmDarkRead(psString *normConcept, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(normConcept, NULL); + PS_ASSERT_PTR_NULL(*normConcept, NULL); + PS_ASSERT_FITS_NON_NULL(fits, NULL); if (!psFitsMoveExtName(fits, PM_DARK_FITS_EXTNAME)) { @@ -589,26 +685,56 @@ } - psArray *table = psFitsReadTable(fits); // FITS Table with ordinates - int numOrdinates = table->n; // Number of ordinates - psArray *ordinates = psArrayAlloc(numOrdinates); // Parsed ordinates - - for (int i = 0; i < numOrdinates; i++) { - psMetadata *row = table->data[i]; // Row of interest - const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name - int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order - if (!name || order <= 0) { - psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table."); - psFree(table); - psFree(ordinates); - return false; - } - pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data - ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE); - ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN); - ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX); - - ordinates->data[i] = ord; - } - psFree(table); + psMetadata *header = psFitsReadHeader(NULL, fits); // Header + bool mdok; // Status of MD lookup + *normConcept = psMemIncrRefCounter(psMetadataLookupStr(&mdok, header, PM_DARK_HEADER_NORM)); + + psArray *ordinates = NULL; // Dark ordinates to return + + psFitsType type = psFitsGetExtType(fits); // Type of FITS extension + switch (type) { + case PS_FITS_TYPE_IMAGE: { + // Check that it's of zero size; otherwise we might have some conflict + int numCols = psMetadataLookupS32(&mdok, header, "NAXIS1"); + int numRows = psMetadataLookupS32(&mdok, header, "NAXIS2"); + if (numCols != 0 || numRows != 0) { + psError(PS_ERR_UNKNOWN, true, "Extension %s is not a DARK table.", PM_DARK_FITS_EXTNAME); + psFree(header); + return NULL; + } + // No ordinates fit --- only a constant term + ordinates = psArrayAlloc(0); + break; + } + case PS_FITS_TYPE_BINARY_TABLE: + case PS_FITS_TYPE_ASCII_TABLE: { + psArray *table = psFitsReadTable(fits); // FITS Table with ordinates + int numOrdinates = table->n; // Number of ordinates + ordinates = psArrayAlloc(numOrdinates); + + for (int i = 0; i < numOrdinates; i++) { + psMetadata *row = table->data[i]; // Row of interest + const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name + int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order + if (!name || order <= 0) { + psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table."); + psFree(table); + psFree(ordinates); + return false; + } + pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data + ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE); + ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN); + ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX); + + ordinates->data[i] = ord; + } + psFree(table); + break; + } + default: + psError(PS_ERR_UNKNOWN, true, "Unrecognised FITS extension type."); + return NULL; + } + psFree(header); return ordinates; Index: trunk/psModules/src/detrend/pmDark.h =================================================================== --- trunk/psModules/src/detrend/pmDark.h (revision 16948) +++ trunk/psModules/src/detrend/pmDark.h (revision 16949) @@ -7,4 +7,6 @@ #define PM_DARK_ANALYSIS_ORDINATES "DARK.ORDINATES" // Name for dark ordinates in the cell analysis metadata +#define PM_DARK_ANALYSIS_NORM "DARK.NORM" // Name for dark normalisation concept in cell analysis metadata +#define PM_DARK_HEADER_NORM "PSDRKNRM" // Header keyword for dark normalisation concept // An ordinate for fitting darks @@ -26,4 +28,5 @@ const psArray *inputs, // Input readouts for combination psArray *ordinates, // Ordinates for fitting + const char *normConcept, // Concept name to use to divide input pixel values int iter, // Number of rejection iterations float rej, // Rejection threshold (standard deviations) @@ -82,10 +85,12 @@ // Write dark table to FITS file bool pmDarkWrite(psFits *fits, // FITS file to which to write - const psMetadata *header, // Header to write - const psArray *ordinates // Dark ordinates to write + psMetadata *header, // Header to write + const psArray *ordinates, // Dark ordinates to write + const char *normConcept // Normalisation concept name ); // Read dark table from FITS file -psArray *pmDarkRead(psFits *fits // FITS file to read +psArray *pmDarkRead(psString *normConcept, // Normalisation concept name + psFits *fits // FITS file to read );