Index: branches/eam_branches/stackphot.20100406/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/Makefile.am (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/Makefile.am (revision 27649) @@ -87,13 +87,14 @@ # a psphot-variant for stack photometry psphotStack_SOURCES = \ - psphotStack.c \ - psphotStackArguments.c \ - psphotStackParseCamera.c \ - psphotStackImageLoop.c \ - psphotStackReadout.c \ - psphotStackChisqImage.c \ - psphotCleanup.c - -# psphotFitSourceLinearStack.c + psphotStack.c \ + psphotStackArguments.c \ + psphotStackParseCamera.c \ + psphotStackImageLoop.c \ + psphotStackReadout.c \ + psphotStackChisqImage.c \ + psphotFitSourcesLinearStack.c \ + psphotSourceMatch.c \ + psphotCleanup.c + Index: branches/eam_branches/stackphot.20100406/psphot/src/psphot.h =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphot.h (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphot.h (revision 27649) @@ -350,8 +350,18 @@ bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration const pmFPAview *view, - pmReadout *chiReadout, + pmReadout **chiReadout, char *filename, int index); +bool psphotStackRemoveChisqFromInputs (pmConfig *config); +bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num); + +psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view); +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index); +bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS); + +bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final); +int pmPhotObjSortBySN (const void **a, const void **b); +int pmPhotObjSortByX (const void **a, const void **b); #endif Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotApResid.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotApResid.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotApResid.c (revision 27649) @@ -16,6 +16,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotChoosePSF.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotChoosePSF.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotChoosePSF.c (revision 27649) @@ -13,6 +13,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotEfficiency.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotEfficiency.c (revision 27649) @@ -160,6 +160,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c (revision 27649) @@ -3,9 +3,9 @@ // XXX define the 'good' / 'bad' flags? +# define COVAR_FACTOR 1.0 + bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) { bool status; - float x; - float y; float f; @@ -29,5 +29,5 @@ // analysis is done in spatial order (to speed up overlap search) // sort by first element in each source list - objects = psArraySort (objects, pmPhotObjSortByY); + objects = psArraySort (objects, pmPhotObjSortByX); // storage array for fitSources @@ -66,5 +66,5 @@ } } - psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n); + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); if (fitSources->n == 0) { @@ -86,10 +86,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor); + float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -98,5 +98,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseVectorElement (sparse, i, f); @@ -109,11 +109,11 @@ // skip over disjoint source images, break after last possible overlap - if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break; - if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; - if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue; - if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; + if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; // source(i) is above source(j) + if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) continue; // source(i) is below source(j) + if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; // source(i) is right of source(j) + if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) break; // source(i) is left of source(j) [no other source(j) can overlap source(i)] // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseMatrixElement (sparse, j, i, f); } @@ -157,5 +157,5 @@ if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor); + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -172,24 +172,12 @@ } -// sort by Y (ascending) -int pmPhotObjSortBySN (const void **a, const void **b) +// sort by X (ascending) +int pmPhotObjSortByX (const void **a, const void **b) { pmPhotObj *objA = *(pmPhotObj **)a; pmPhotObj *objB = *(pmPhotObj **)b; - psAssert (objA->sources, "missing sources?"); - psAssert (objB->sources, "missing sources?"); - - psAssert (objA->sources->n, "missing sources"); - psAssert (objB->sources->n, "missing sources"); - - psAssert (objA->sources->data[0], "missing sources"); - psAssert (objB->sources->data[0], "missing sources"); - - pmSource *A = objA->sources->data[0]; - pmSource *B = objB->sources->data[0]; - - psF32 fA = (A->peak == NULL) ? 0 : A->peak->y; - psF32 fB = (B->peak == NULL) ? 0 : B->peak->y; + psF32 fA = objA->x; + psF32 fB = objB->x; psF32 diff = fA - fB; @@ -198,3 +186,2 @@ return (0); } - Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotGuessModels.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotGuessModels.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotGuessModels.c (revision 27649) @@ -15,6 +15,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotImageLoop.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotImageLoop.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotImageLoop.c (revision 27649) @@ -91,15 +91,15 @@ } - status = true; - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); - if (!status) { - psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); - psFree (view); - return false; - } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } } - // save output which is saved at the chip level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotImageQuality.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotImageQuality.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotImageQuality.c (revision 27649) @@ -16,6 +16,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotMagnitudes.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotMagnitudes.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotMagnitudes.c (revision 27649) @@ -12,6 +12,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image // find the currently selected readout Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotModelBackground.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotModelBackground.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotModelBackground.c (revision 27649) @@ -371,6 +371,6 @@ psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters - pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning); - pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning); + pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index); + pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotReadoutCleanup.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotReadoutCleanup.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotReadoutCleanup.c (revision 27649) @@ -53,5 +53,5 @@ pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psArray *sources = detections->allSources; + psArray *sources = detections ? detections->allSources : NULL; // XXX where do we free these, in here (psMetadataRemove?) @@ -73,5 +73,5 @@ // Check to see if the image quality was measured // XXX not sure we want / need this test - if (!psf) { + if (0 && !psf) { bool mdok; // Status of MD lookup int nIQ = psMetadataLookupS32(&mdok, recipe, "IQ_NSTAR"); // Number of stars for IQ measurement Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotRoughClass.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotRoughClass.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotRoughClass.c (revision 27649) @@ -19,6 +19,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSkyReplace.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotSkyReplace.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSkyReplace.c (revision 27649) @@ -8,6 +8,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c (revision 27649) @@ -14,56 +14,128 @@ if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i); - return false; + psFree (objects); + return NULL; } } + + // psphotMatchSourcesGenerate (objects, config, view, "PSPHOT.INPUT", i); + + return objects; +} + +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int RADIUS = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.MATCH.RADIUS"); + psAssert (status, "programming error: must define PSPHOT.STACK.MATCH.RADIUS"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + psAssert (detections->newSources, "new sources not defined?"); + psAssert (!detections->allSources, "all sources already defined?"); + + // XXX TEST: + if (detections->newSources) { + psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS); + } + return true; } -// we need a couple of functions to distinguish coincident sources: -// XXX these are similar (identical?) to the goals of pmSourceMatch.c - # define NEXT1 { i++; continue; } # define NEXT2 { j++; continue; } - -bool psphotSourceMerge (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { +bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { float dx, dy; + float RADIUS2 = RADIUS*RADIUS; + // sort the source list by X - pmSourceSortByX (sources1); - pmSourceSortByX (sources2); + sources = psArraySort (sources, pmSourceSortByX); + objects = psArraySort (objects, pmPhotObjSortByX); + psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8); + psVectorInit (found, 0); + + // match sources to existing objects + + psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld)", sources->n, objects->n); + int i, j; - for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { + for (i = j = 0; (i < sources->n) && (j < objects->n); ) { - pmSource *src1 = sources1->data[i]; - pmSource *src2 = sources2->data[j]; + pmSource *src = sources->data[i]; + pmPhotObj *obj = objects->data[j]; - if (!src1) NEXT1; - if (!src1->peak) NEXT1; - if (!finite(src1->peak->xf)) NEXT1; - if (!finite(src1->peak->yf)) NEXT1; + if (!src) NEXT1; + if (!src->peak) NEXT1; + if (!finite(src->peak->xf)) NEXT1; + if (!finite(src->peak->yf)) NEXT1; - if (!src2) NEXT2; - if (!src2->peak) NEXT2; - if (!finite(src2->peak->xf)) NEXT2; - if (!finite(src2->peak->yf)) NEXT2; + if (!obj) NEXT2; + if (!finite(obj->x)) NEXT2; + if (!finite(obj->y)) NEXT2; - dx = src1->peak->xf - src2->peak->xf; + dx = src->peak->xf - obj->x; if (dx < -1.02*RADIUS) NEXT1; if (dx > +1.02*RADIUS) NEXT2; // we are within match range, look for matches: - for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { + int Jmin = -1; + float Rmin = RADIUS2; + for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { - dx = src1->peak->xf - src2->peak->xf; - dy = src1->peak->yf - src2->peak->yf; + obj = objects->data[J]; + + dx = src->peak->xf - obj->x; + dy = src->peak->yf - obj->y; - dr = dx*dx + dy*dy; + float dr = dx*dx + dy*dy; if (dr > RADIUS2) continue; - - // add to group? - } + if (dr < Rmin) { + Rmin = dr; + Jmin = J; + } + } + + // no match, try next source + if (Jmin == -1) { + i++; + continue; + } + obj = objects->data[Jmin]; + + // add to object + pmPhotObjAddSource (obj, src); + found->data.U8[i] = 1; i++; } + + // add missed sources to new objects + + for (i = 0; i < sources->n; i++) { + + if (found->data.U8[i]) continue; + + pmSource *src = sources->data[i]; + + pmPhotObj *obj = pmPhotObjAlloc(); + pmPhotObjAddSource(obj, src); + psArrayAdd (objects, 100, obj); + } + psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n); + + return true; } Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceSize.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceSize.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceSize.c (revision 27649) @@ -44,6 +44,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) { psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotStackChisqImage.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotStackChisqImage.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotStackChisqImage.c (revision 27649) @@ -12,16 +12,8 @@ psTimerStart ("psphot.chisq.image"); - pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.OUTPUT", 0); + pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.IMAGE", 0); psAssert (chisqFile, "missing chisq image FPA?"); - pmCell *chisqCell = pmFPAviewThisCell(view, chisqFile->fpa); - - // create a holder for the image - pmReadout *chiReadout = pmFPAviewThisReadout(view, chisqFile->fpa); - if (!chiReadout) { - chiReadout = pmReadoutAlloc(chisqCell); - } else { - psMemIncrRefCounter(chiReadout); - } + pmReadout *chiReadout = NULL; int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); @@ -30,5 +22,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotStackChisqImageAddReadout(config, view, chiReadout, "PSPHOT.INPUT", i)) { + if (!psphotStackChisqImageAddReadout(config, view, &chiReadout, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); return false; @@ -51,5 +43,4 @@ psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image")); - psFree (chiReadout); return true; } @@ -57,5 +48,5 @@ bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration const pmFPAview *view, - pmReadout *chiReadout, + pmReadout **chiReadout, char *filename, int index) @@ -71,10 +62,24 @@ psImage *inImage = inReadout->image; + psAssert (inImage, "missing image?"); + psImage *inVariance = inReadout->variance; + psAssert (inVariance, "missing variance?"); + psImage *inMask = inReadout->mask; + psAssert (inMask, "missing mask?"); - psImage *chiImage = chiReadout->image; - psImage *chiVariance = chiReadout->variance; - psImage *chiMask = chiReadout->mask; + if (*chiReadout == NULL) { + *chiReadout = pmFPAGenerateReadout(config, view, "PSPHOT.CHISQ.IMAGE", input->fpa, NULL, 0); + } + + psImage *chiImage = (*chiReadout)->image; + psAssert (chiImage, "missing chi image"); + + psImage *chiVariance = (*chiReadout)->variance; + psAssert (chiVariance, "missing chi variance"); + + psImage *chiMask = (*chiReadout)->mask; + psAssert (chiMask, "missing chi mask"); // select the appropriate recipe information @@ -105,2 +110,72 @@ return true; } + +bool psphotStackRemoveChisqFromInputs (pmConfig *config) { + + bool status = false; + + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM"); + + int inputNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + pmFPAfileRemoveSingle (config->files, "PSPHOT.INPUT", chisqNum); + + inputNum --; + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum); + + return true; +} + +bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num) +{ + PS_ASSERT_PTR_NON_NULL(files, NULL); + PS_ASSERT_INT_NONNEGATIVE(num, NULL); + + psList* mdList = files->list; + psHash* mdHash = files->hash; + + // Generate a REGEX to select only items that match 'name' + psString regex = NULL; // Regular expression + if (name) { + if (!psMetadataLookup(files, name)) { + // No files match the requested name + return false; + } + psStringAppend(®ex, "^%s$", name); + } + + psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + + bool found = false; + int i = 0; // Counter + for (i = 0; !found && (item = psMetadataGetAndIncrement(iter)); i++) { + if (i == num) found = true; + } + psFree(iter); + if (!found) { + return false; + } + + char *key = item->name; + + // look up the name via hash to see if we have a multi or not + psMetadataItem* hashItem = psHashLookup(mdHash, name); + if (hashItem == NULL) { + psError(PS_ERR_UNKNOWN, false, _("Failed to remove metadata item, %s, from metadata table."), key); + return false; + } + + if (hashItem->type == PS_DATA_METADATA_MULTI) { + // multiple entries with same key, remove just the specified one + psListRemoveData(hashItem->data.list, item); + } else { + psHashRemove(mdHash, key); + } + psListRemoveData(mdList, item); + + return true; +} Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotStackImageLoop.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotStackImageLoop.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotStackImageLoop.c (revision 27649) @@ -9,4 +9,15 @@ bool psphotStackImageLoop (pmConfig *config) { + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + pmFPAview *view = pmFPAviewAlloc (0); @@ -14,6 +25,44 @@ if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); - // XXX for now, we assume there is only a single chip in the PHU: - psphotStackReadout (config, view); + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack."); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack."); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 6, "Readout %d: %x %x\n", view->readout, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + // XXX for now, we assume there is only a single chip in the PHU: + if (!psphotStackReadout (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + + } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } + } + // save output which is saved at the chip level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + } + // save output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); // fail if we failed to handle an error @@ -27,4 +76,4 @@ the easiest way to implement this is to assume we can pre-load the full set of images up front. with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad. - */ +*/ Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotStackParseCamera.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotStackParseCamera.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotStackParseCamera.c (revision 27649) @@ -15,7 +15,7 @@ } - psMetadataItem *item = NULL; - int nInputs = 0; - while ((item = psMetadataGet(inputs, nInputs)) != NULL) { + int nInputs = inputs->list->n; + for (int i = 0; i < nInputs; i++) { + psMetadataItem *item = psMetadataGet(inputs, i); if (item->type != PS_DATA_METADATA) { psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name); @@ -34,5 +34,5 @@ pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image if (!imageFile) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", nInputs, image); + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image); return false; } @@ -40,6 +40,6 @@ psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask if (mask && strlen(mask) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.INPUT.MASK", mask, PM_FPA_FILE_MASK)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", nInputs, mask); + if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); return false; } @@ -48,13 +48,25 @@ psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map if (variance && strlen(variance) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.INPUT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", nInputs, variance); + if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); return false; } } - nInputs ++; + // the output sources are carried on the input->fpa structures + pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); + return false; + } + outsources->save = true; + outsources->fileID = i; // this is used to generate output names } psMetadataRemoveKey(config->arguments, "FILENAMES"); psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs); + + if (!psphotSetMaskBits (config)) { + psError (PS_ERR_UNKNOWN, false, "failed to set mask bit values"); + return NULL; + } // generate an pmFPAimage for the chisqImage @@ -64,22 +76,27 @@ return false; } - pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.MASK"); + chisqImage->save = true; + + pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); if (!chisqMask) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPIMAGE.CHISQ.MASK")); + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); return NULL; } if (chisqMask->type != PM_FPA_FILE_MASK) { - psError(PS_ERR_IO, true, "PPIMAGE.CHISQ.MASK is not of type MASK"); + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); return NULL; } - pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.VARIANCE"); + chisqMask->save = true; + + pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); if (!chisqVariance) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPIMAGE.CHISQ.VARIANCE")); + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); return NULL; } if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { - psError(PS_ERR_IO, true, "PPIMAGE.CHISQ.VARIANCE is not of type VARIANCE"); + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); return NULL; } + chisqVariance->save = true; # if (0) Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotStackReadout.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotStackReadout.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotStackReadout.c (revision 27649) @@ -59,4 +59,5 @@ // find the detections (by peak and/or footprint) in the image. + // This finds the detections on Chisq image as well as the individuals if (!psphotFindDetections (config, view, true)) { // pass 1 // this only happens if we had an error in psphotFindDetections @@ -66,25 +67,33 @@ // construct sources and measure basic stats (saved on detections->newSources) + // only run this on detections from the input images, not chisq image if (!psphotSourceStats (config, view, true)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); return psphotReadoutCleanup (config, view); } - if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, view); + + // *** generate the objects (which unify the sources from the different images) + psArray *objects = psphotMatchSources (config, view); + + // construct sources for the newly-generated sources (from other images) + if (!psphotSourceStats (config, view, false)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); } // find blended neighbors of very saturated stars (detections->newSources) - if (!psphotDeblendSatstars (config, view)) { - psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - return psphotReadoutCleanup (config, view); - } + // if (!psphotDeblendSatstars (config, view)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + // return psphotReadoutCleanup (config, view); + // } // mark blended peaks PS_SOURCE_BLEND (detections->newSources) - if (!psphotBasicDeblend (config, view)) { - psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); - return psphotReadoutCleanup (config, view); - } + // if (!psphotBasicDeblend (config, view)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); + // return psphotReadoutCleanup (config, view); + // } // classify sources based on moments, brightness + // only run this on detections from the input images, not chisq image if (!psphotRoughClass (config, view)) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); @@ -92,4 +101,5 @@ } // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) + // only run this on detections from the input images, not chisq image if (!psphotImageQuality (config, view)) { // pass 1 psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); @@ -121,9 +131,7 @@ psphotMergeSources (config, view); - // *** generate the objects (which unify the sources from the different images) - // pmArray *objects = psphotMatchSources (config, view); - // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - // psphotFitSourcesLinearStack (config, objects, FALSE); + psphotFitSourcesLinearStack (config, objects, FALSE); + psFree (objects); // identify CRs and extended sources @@ -153,4 +161,7 @@ psphotSourceFreePixels (config, view); + // remove chisq image from config->file:PSPHOT.INPUT (why?) + psphotStackRemoveChisqFromInputs(config); + // create the exported-metadata and free local data return psphotReadoutCleanup (config, view); Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSubtractBackground.c =================================================================== --- branches/eam_branches/stackphot.20100406/psphot/src/psphotSubtractBackground.c (revision 27628) +++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSubtractBackground.c (revision 27649) @@ -23,5 +23,11 @@ pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest assert (modelFile); - pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa); + + pmReadout *model = NULL; + if (modelFile->mode == PM_FPA_MODE_INTERNAL) { + model = modelFile->readout; + } else { + model = pmFPAviewThisReadout (view, modelFile->fpa); + } assert (model);