Index: /trunk/ppImage/src/Makefile.am =================================================================== --- /trunk/ppImage/src/Makefile.am (revision 25929) +++ /trunk/ppImage/src/Makefile.am (revision 25930) @@ -41,4 +41,5 @@ ppImageDetrendFringe.c \ ppImageDetrendFree.c \ + ppImageDetrendPattern.c \ ppImageRebinReadout.c \ ppImageMosaic.c \ Index: /trunk/ppImage/src/ppImage.h =================================================================== --- /trunk/ppImage/src/ppImage.h (revision 25929) +++ /trunk/ppImage/src/ppImage.h (revision 25930) @@ -159,4 +159,6 @@ bool ppImageBurntoolMask(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *mask); +bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options); + // Record which detrend file was used for the detrending bool ppImageDetrendRecord( Index: /trunk/ppImage/src/ppImageDetrendFringe.c =================================================================== --- /trunk/ppImage/src/ppImageDetrendFringe.c (revision 25929) +++ /trunk/ppImage/src/ppImageDetrendFringe.c (revision 25930) @@ -84,4 +84,6 @@ psArray *cells = chip->cells; // Component cells psArray *fringes = psArrayAlloc(cells->n); // Fringes, to return + int video_cell_zero = 0; + for (int i = 0; i < cells->n; i++) { fringes->data[i] = NULL; @@ -89,13 +91,42 @@ pmCell *cell = cells->data[i]; // Cell of interest + psTrace("psModules.detrend",7,"Readouts: Cell %d chip: %ld\n",i,cell->readouts->n); // XXX for now, skip the video cells (cell->readouts->n > 1) - if (cell->readouts->n > 1) { - psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i); - continue; - } - + // CZW: This mess creates a fake set of fringe stats by stealing the previous one. + // We let the fitting code know that this is all lies by scaling the weights by a crazy amount. + + if ( (cell->readouts->n > 1) ) { + psTrace("psModules.detrend",7,"Should be skipping scichip: %d\n",i); + psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i); + + if (i == 0) { + video_cell_zero = 1; + } + else { + pmFringeStats *prevFringe = fringes->data[i-1]; + pmFringeStats *fringe = pmFringeStatsAlloc(prevFringe->regions); + for (int j = 0; j < fringe->regions->nRequested; j++) { + fringe->f->data.F32[j] = prevFringe->f->data.F32[j]; + fringe->df->data.F32[j] = prevFringe->df->data.F32[j] / 1e6; + } + fringes->data[i] = fringe; + } + + continue; + } + fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis, source)); } + if (video_cell_zero == 1) { + pmFringeStats *prevFringe = fringes->data[1]; + pmFringeStats *fringe = pmFringeStatsAlloc(prevFringe->regions); + for (int j = 0; j < fringe->regions->nRequested; j++) { + fringe->f->data.F32[j] = NAN; + fringe->df->data.F32[j] = 1.0; + } + fringes->data[0] = fringe; + } + return fringes; } @@ -112,5 +143,5 @@ psArray *science = NULL; if (isResidual) { - science = getFringes(scienceChip, "FRINGE.RESIDUALS"); // Fringe residuals on science chip + science = getFringes(scienceChip, "FRINGE.RESIDUALS"); // Fringe residuals on science chip } else { science = getFringes(scienceChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on science chip @@ -124,12 +155,13 @@ psArray *references = getFringes(refChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on reference chip int numRefs = ((psArray*)references->data[0])->n; // Number of reference fringes - psArray *referencesCat = psArrayAlloc(numRefs); // Reference fringes - for (int i = 0; i < numRefs; i++) { // Iterate over fringes - psArray *refs = psArrayAlloc(references->n); // Array of fringes for each cell - for (int j = 0; j < references->n; j++) { // Iterate over cells - psArray *ref = references->data[j]; // Array of references for this cell + psArray *referencesCat = psArrayAlloc(numRefs); // Reference fringes + for (int i = 0; i < numRefs; i++) { // Iterate over fringes + psArray *refs = psArrayAlloc(references->n); // Array of fringes for each cell + for (int j = 0; j < references->n; j++) { // Iterate over cells + psArray *ref = references->data[j]; // Array of references for this cell + refs->data[j] = psMemIncrRefCounter(ref->data[i]); } - referencesCat->data[i] = pmFringeStatsConcatenate(refs, NULL, NULL); + referencesCat->data[i] = pmFringeStatsConcatenate(refs, NULL, NULL); psFree(refs); } Index: /trunk/ppImage/src/ppImageDetrendPattern.c =================================================================== --- /trunk/ppImage/src/ppImageDetrendPattern.c (revision 25930) +++ /trunk/ppImage/src/ppImageDetrendPattern.c (revision 25930) @@ -0,0 +1,69 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ppImage.h" + +#define ESCAPE(MESSAGE) { \ + psError(PS_ERR_UNKNOWN, false, MESSAGE); \ + psFree(view); \ + return false; \ +} + +bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options) { + + + pmCell *cell = NULL; + + assert (options->doPattern); // do not call if not needed + assert (inputView->chip != -1); + assert (inputView->cell == -1); + assert (inputView->readout == -1); + bool status; + pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.INPUT"); + + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + *view = *inputView; + + while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { + if (!cell->process || !cell->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + ESCAPE("load failure for Cell"); + } + + if (!cell->data_exists) { + continue; + } + + if (cell->readouts->n > 1) { + psWarning ("Skipping Video Cell for ppImageDetrendPatternApply"); + continue; + } + + // process each of the readouts + pmReadout *readout; // Readout from cell + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + ESCAPE("load failure for Readout"); + } + if (!readout->data_exists) { + continue; + } + + // perfore pattern correction + if (!pmPatternRow(readout, options->patternOrder, options->patternIter, options->patternRej, + options->patternThresh, options->patternMean, options->patternStdev, + options->maskValue, options->darkMask)) { + psFree(view); + return(false); + } + } + } + + psFree(view); + return(true); +} + + Index: /trunk/ppImage/src/ppImageDetrendReadout.c =================================================================== --- /trunk/ppImage/src/ppImageDetrendReadout.c (revision 25929) +++ /trunk/ppImage/src/ppImageDetrendReadout.c (revision 25930) @@ -128,13 +128,13 @@ } - // Pattern noise correction - if (options->doPattern) { - if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, - options->patternThresh, options->patternMean, options->patternStdev, - options->maskValue, options->darkMask)) { - psFree(detview); - return false; - } - } +/* // Pattern noise correction */ +/* if (options->doPattern) { */ +/* if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */ +/* options->patternThresh, options->patternMean, options->patternStdev, */ +/* options->maskValue, options->darkMask)) { */ +/* psFree(detview); */ +/* return false; */ +/* } */ +/* } */ // Normalization by a single (known) constant Index: /trunk/ppImage/src/ppImageLoop.c =================================================================== --- /trunk/ppImage/src/ppImageLoop.c (revision 25929) +++ /trunk/ppImage/src/ppImageLoop.c (revision 25930) @@ -163,4 +163,11 @@ } + // Apply the pattern noise correction + if (options->doPattern) { + if (!ppImageDetrendPatternApply(config,chip,view,options)) { + ESCAPE("Unable to apply pattern corrections"); + } + } + // measure various pixel-based statistics for this image if (!ppImagePixelStats(config, stats, options, view)) { Index: /trunk/ppImage/src/ppImageOptions.c =================================================================== --- /trunk/ppImage/src/ppImageOptions.c (revision 25929) +++ /trunk/ppImage/src/ppImageOptions.c (revision 25930) @@ -249,5 +249,6 @@ options->burntoolTrails = psMetadataLookupS32(&status, recipe, "BURNTOOL.TRAILS"); - fprintf(stderr,"TRAILS: %d %d %d\n",options->burntoolTrails,psMetadataLookupS32(&status,recipe,"BURNTOOL.TRAILS"),status); + psTrace("psModules.detrend", 7, "burntoolTrails: %d BURNTOOL.TRAILS: %d Status: %d\n", + options->burntoolTrails,psMetadataLookupS32(&status,recipe,"BURNTOOL.TRAILS"),status); if (!status) { psWarning("BURNTOOL.TRAILS not found in recipe: setting to default value.\n"); Index: /trunk/psModules/src/concepts/pmConcepts.c =================================================================== --- /trunk/psModules/src/concepts/pmConcepts.c (revision 25929) +++ /trunk/psModules/src/concepts/pmConcepts.c (revision 25930) @@ -249,4 +249,5 @@ CONCEPT_REGISTER_FUNCTION(S32, Enum, -1); // For enums: set default to -1 CONCEPT_REGISTER_FUNCTION(S32, S32, 0); // For values: set default to 0 +//CONCEPT_REGISTER_FUNCTION(Bool, Bool, NULL); // For values: set default to 0 static void conceptRegisterTime(const char *name, /* Name of concept */ \ @@ -333,5 +334,6 @@ conceptRegisterF32("FPA.TELTEMP.EXTRA", "Telescope Temperatures: extra", p_pmConceptParse_TELTEMPS, NULL, NULL, false, PM_FPA_LEVEL_FPA); conceptRegisterF32("FPA.PON.TIME", "Power On Time", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA); - conceptRegisterS32("FPA.BURNTOOL.APPLIED", "Status of burntool processing", p_pmConceptParse_BTOOLAPP,NULL,NULL,false,PM_FPA_LEVEL_FPA); + conceptRegisterS32("FPA.BURNTOOL.APPLIED", "Status of burntool processing", p_pmConceptParse_BTOOLAPP,p_pmConceptFormat_BTOOLAPP,NULL,false,PM_FPA_LEVEL_FPA); + // conceptRegisterBool("FPA.BURNTOOL.APPLIED", "Status of burntool processing", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA); conceptRegisterF32("FPA.EXPOSURE", "Exposure time (sec)", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA); } Index: /trunk/psModules/src/concepts/pmConceptsRead.c =================================================================== --- /trunk/psModules/src/concepts/pmConceptsRead.c (revision 25929) +++ /trunk/psModules/src/concepts/pmConceptsRead.c (revision 25930) @@ -47,4 +47,6 @@ case PS_DATA_F64: return psMetadataItemAllocF64(pattern->name, pattern->comment, psMetadataItemParseF64(concept)); + case PS_DATA_BOOL: + return psMetadataItemAllocBool(pattern->name, pattern->comment, psMetadataItemParseBool(concept)); default: psWarning("Concept %s (%s) is not of a standard type (%x)\n", Index: /trunk/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- /trunk/psModules/src/concepts/pmConceptsStandard.c (revision 25929) +++ /trunk/psModules/src/concepts/pmConceptsStandard.c (revision 25930) @@ -753,5 +753,17 @@ psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not BOOL\n", concept->name, concept->type); - return NULL; + if (concept->type != PS_DATA_S32) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Wasn't the type I'd guessed either.\n"); + return NULL; + } + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Looks like an S32 value? (%d)\n", + concept->data.S32); + + if (concept->data.S32 == 0) { + bt_status = 1; + } + else if (concept->data.S32 == 1) { + bt_status = -2; + } } @@ -765,4 +777,32 @@ return psMetadataItemAllocS32(concept->name, concept->comment, bt_status); } +psMetadataItem *p_pmConceptFormat_BTOOLAPP(const psMetadataItem *concept, + pmConceptSource source, + const psMetadata *cameraFormat, + const pmFPA *fpa, + const pmChip *chip, + const pmCell *cell) +{ + assert(concept); + + if (concept->type != PS_DATA_S32) { + return NULL; + } + + if (concept->data.S32 == 0) { + return NULL; + } + else if (concept->data.S32 == -2) { + return psMetadataItemAllocBool(concept->name,concept->comment,true); + } + else if (concept->data.S32 == 1) { + return psMetadataItemAllocBool(concept->name,concept->comment,false); + } + else { + return NULL; + } + +} + // Get the current value of a concept Index: /trunk/psModules/src/concepts/pmConceptsStandard.h =================================================================== --- /trunk/psModules/src/concepts/pmConceptsStandard.h (revision 25929) +++ /trunk/psModules/src/concepts/pmConceptsStandard.h (revision 25930) @@ -146,4 +146,12 @@ const pmCell *cell ///< Cell for concept, or NULL ); +psMetadataItem *p_pmConceptFormat_BTOOLAPP( + const psMetadataItem *concept, ///< Concept to format + pmConceptSource source, ///< Source for concept + const psMetadata *cameraFormat, ///< Camera format definition + const pmFPA *fpa, ///< FPA for concept, or NULL + const pmChip *chip, ///< Chip for concept, or NULL + const pmCell *cell ///< Cell for concept, or NULL + ); Index: /trunk/psModules/src/detrend/pmFringeStats.c =================================================================== --- /trunk/psModules/src/detrend/pmFringeStats.c (revision 25929) +++ /trunk/psModules/src/detrend/pmFringeStats.c (revision 25930) @@ -11,4 +11,9 @@ #include "pmFPA.h" #include "pmFringeStats.h" + +#include "psPolynomialMD.h" +#include "psMinimizePolyFit.h" +#include "psVector.h" + // Future optimisations for speed: @@ -331,6 +336,6 @@ dfPt[i] = 1.0 / medianSd->sampleStdev; - psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f\n", (int)region.x0, (int)region.x1, - (int)region.y0, (int)region.y1, fPt[i], dfPt[i]); + psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f : %s\n", (int)region.x0, (int)region.x1, + (int)region.y0, (int)region.y1, fPt[i], dfPt[i], readout->parent->hdu->extname); } psFree(sky); @@ -840,5 +845,5 @@ } } - B->data.F64[i] = vector; + B->data.F64[i] = vector; } @@ -966,4 +971,5 @@ } } + } @@ -988,7 +994,46 @@ psTrace("psModules.detrend", 9, "Masking region %d because not finite in fringe %d.\n", j, i); } - } - } - + else if (fabs(fringe->f->data.F32[j]) > 0.1) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 1; + psTrace("psModules.detrend", 9, "Masking region %d because too large fringe %d.\n", j, i); + } + // Mask bad points in the science data as well. + if ((i == 0) && (!isfinite(science->f->data.F32[j]))) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 1; + psTrace("psModules.detrend", 9, "Masking region %d because not finite in science fringe %d.\n", j, i); + } + psTrace("psModules.detrend", 7, "F %f %f %f %d\n", + fringe->f->data.F32[j], science->f->data.F32[j], + 1 / science->df->data.F32[j],(int) mask->data.PS_TYPE_VECTOR_MASK_DATA[j]); + } + } + +/* // Begin switch from old outlier removal and fitting code. */ + + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + pmFringeStats *fringe = fringes->data[0]; + psVector *errors = psVectorAlloc(science->df->n,PS_TYPE_F32); + for (int j = 0; j < errors->n; j++) { + errors->data.F32[j] = 1 / science->df->data.F32[j]; + } + psVectorFitPolynomial1D(poly,mask,0xff,science->f,errors,fringe->f); + + for (int i = 0; i <= poly->nX; i++) { + scale->coeff->data.F32[i] = poly->coeff[i]; + psTrace("psModules.detrend",7,"COEFFS: %d %g %g %g\n",i,scale->coeff->data.F32[i],poly->coeff[i],poly->coeffErr[i]); + } + + psFree(poly); + // psFree(fringe); + psFree(errors); + + psFree(median); + psFree(diff); + return scale; + // End switch from old code. + + + # if (0) // Write fringe data to file for a test @@ -1019,5 +1064,7 @@ } scale->coeff->data.F32[0] -= median->sampleMedian; - scale->coeff->data.F32[i] = 0.0; + if (i != 0) { + scale->coeff->data.F32[i] = 0.0; + } } psFree(median); @@ -1051,4 +1098,5 @@ return scale; + //# endif }