Index: /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.c (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.c (revision 36829) @@ -835,4 +835,5 @@ psImageMaskType goodMask, // Value for good pixels bool safe, // Safe combination? + int nminpix, // Minimum number of input per pixel float invTotalWeight // Inverse of total weight for all inputs ) @@ -854,5 +855,11 @@ CHECKPIX(x, y, "bad vs good : %x %x %x\n", maskValue, badMask, blankMask); - switch (num) { + //MEH -- hackish adding of lower limit for input per pixel + int numN = num; + if (num < nminpix) { + CHECKPIX(x, y, "Nmin (%d) inputs (%d) to combine, pixel %d,%d is manually set bad\n", nminpix, numN, x, y); + numN = 0; + } + switch (numN) { case 0: { // Nothing to combine: it's bad @@ -1518,4 +1525,5 @@ bool useVariance, bool safe, + int nminpix, bool rejection) { @@ -1692,5 +1700,5 @@ psImageMaskType goodMask = 0; // OR of mask bits in all good input pixels combineExtract(&num, &suspect, &badMask, &goodMask, buffer, combinedImage, combinedMask, combinedVariance, input, weights, exps, addVariance, reject, x, y, badMaskBits, suspectMaskBits); - combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, totalExpWeight); + combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight, num, buffer, x, y, blankMaskBits, badMask, goodMask, safe, nminpix, totalExpWeight); if (iter > 0) { Index: /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.h (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/imcombine/pmStack.h (revision 36829) @@ -61,4 +61,5 @@ bool useVariance, ///< Use variance values for rejection? bool safe, ///< Play safe with small numbers of input pixels (mask if N <= 2)? + int nminpix, ///< Minimum number input per pixel to combine bool rejectInspect ///< Reject pixels instead of marking them for inspection? ); Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.c (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.c (revision 36829) @@ -162,4 +162,5 @@ source->sky = NAN; source->skyErr = NAN; + source->extSN = NAN; source->region = psRegionSet(NAN, NAN, NAN, NAN); Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.h (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.h (revision 36829) @@ -113,4 +113,5 @@ float sky; ///< The sky at the center of the object float skyErr; ///< The sky error at the center of the object + float extSN; ///< for externally supplied source the kron signal to noise (used by full force) psRegion region; ///< area on image covered by selected pixels @@ -119,5 +120,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections - pmSourceGalaxyFits *galaxyFits; ///< fits to galaxy models (psphotFullForce only) + psArray *galaxyFits; ///< fits to galaxy models (psphotFullForce only) pmSourceLensing *lensingOBJ; ///< lensing moments parameters (per object) pmSourceLensing *lensingPSF; ///< lensing moments parameters (psf, interpolated) Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceExtendedPars.h (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceExtendedPars.h (revision 36829) @@ -83,8 +83,16 @@ typedef struct { + int modelType; psVector *Flux; psVector *dFlux; psVector *chisq; - int nPix; + int nPix; + bool reducedTrials; + float fRmajorMin; + float fRmajorMax; + float fRmajorDel; + float fRminorMin; + float fRminorMax; + float fRminorDel; } pmSourceGalaxyFits; Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.c (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.c (revision 36829) @@ -1202,5 +1202,11 @@ break; - case PM_FPA_FILE_CFF: + case PM_FPA_FILE_CFF: { + // determine the output table format + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } // read in header, if not yet loaded hdu = pmFPAviewThisHDU (view, file->fpa); @@ -1244,5 +1250,5 @@ } - sources = pmSourcesRead_CFF(file->fits, hdu->header); + sources = pmSourcesRead_CFF(file->fits, hdu->header, recipe); psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname); @@ -1250,4 +1256,5 @@ psFree (dataname); psFree (tableHeader); + } break; Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.h (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.h (revision 36829) @@ -55,5 +55,5 @@ psArray *pmSourcesReadCMP (char *filename, psMetadata *header); -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe); bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe); Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CFF.c (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CFF.c (revision 36829) @@ -53,6 +53,6 @@ #include "pmSourceOutputs.h" -// read in a readout from the fits file -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header) +// read in sources readout from a cff fits file +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -75,12 +75,42 @@ PS_ASSERT_INT_NONNEGATIVE(modelType, NULL); + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + + psString modelForce = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + psF32 forceDevSersicMin = NAN; + bool forceAll = false; + int modelTypeForce = 0; + if (!strcmp(modelForce, "ALL")) { + forceAll = true; + } else if (!strcmp(modelForce, "PS_MODEL_SERSIC")) { + modelTypeForce = sersicModelType; + forceDevSersicMin = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + } else { + modelTypeForce = pmModelClassGetType(modelForce); + PS_ASSERT_INT_NONNEGATIVE(modelTypeForce, NULL); + } + + // skip extended model types for likely stars + // max value of KronMag - psfMag to keep ... + psF32 starCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_MAGDIFF_MAX"); + if (!status) { + starCut = 0; + } + // ... unless SN is less than this value + psF32 SNMinForCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_CUT_SN_MIN"); + if (!status) { + SNMinForCut = 10; + } + // We get the size of the table, and allocate the array of sources first because the table // is large and ephemeral --- when the table gets blown away, whatever is allocated after // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { + long numRows = psFitsTableSize(fits); // Number of rows in table + psArray *sources = psArrayAllocEmpty(numRows); // Array of sources, to return + + // convert the table to the pmSource entries + pmSource *source = NULL; + for (int i = 0; i < numRows; i++) { psMetadata *row = psFitsReadTableRow(fits, i); // Table row if (!row) { @@ -108,4 +138,5 @@ float kronRadius = psMetadataLookupF32 (&status, row, "KRON_RADIUS"); float petRadius = psMetadataLookupF32 (&status, row, "PETRO_RADIUS"); + float SN = psMetadataLookupF32 (&status, row, "SN"); bool fitGalaxy = psMetadataLookupU8 (&status, row, "FIT_GALAXY"); bool psfStar = psMetadataLookupU8 (&status, row, "PSF_STAR"); @@ -114,7 +145,11 @@ float Rminor = psMetadataLookupF32 (&status, row, "R_MINOR"); float theta = psMetadataLookupF32 (&status, row, "THETA"); + float chisq = psMetadataLookupF32 (&status, row, "CHISQ"); + float nDOF = psMetadataLookupF32 (&status, row, "NDOF"); + float magDiff = psMetadataLookupF32 (&status, row, "MAG_DIFF"); + psS16 modelFlags = psMetadataLookupS32 (&status, row, "MODEL_FLAGS"); int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); - if (status) { + if (status && galaxyModelType >= 0) { galaxyModelType = pmModelClassGetLocalType(galaxyModelType); } else { @@ -123,79 +158,112 @@ float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; -// RoughClass wants source type to be unknown -// source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - source->type = PM_SOURCE_TYPE_UNKNOWN; - - // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and - // stars that are masked - if (psfStar) { - source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; - } - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - source->seq = ID; - PAR[PM_PAR_XPOS] = X; - PAR[PM_PAR_YPOS] = Y; - - dPAR[PM_PAR_XPOS] = 0.0; - dPAR[PM_PAR_YPOS] = 0.0; - - PAR[PM_PAR_SKY] = 0.0; - dPAR[PM_PAR_SKY] = 0.0; - - PAR[PM_PAR_I0] = 1.0; - dPAR[PM_PAR_I0] = 0.0; - - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - source->psfMag = 0.0; - source->psfMagErr = 0.0; - source->apMag = 0.0; - source->apRadius = apRadius; - - // we generate a somewhat fake PSF model here -- - // in most (all?) contexts, we will replace this with a measured psf model - // elsewhere - axes.major = 1.0; - axes.minor = 1.0; - axes.theta = 0.0; - pmPSF_AxesToModel (PAR, axes, modelType); - - // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second - source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); - source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here - source->peak->dx = 0.0; - source->peak->dy = 0.0; - - source->moments = pmMomentsAlloc (); - source->moments->Mx = X; - source->moments->My = Y; - source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment - - // Don't mark the moments as measured because that causes many fields to be left blank. - // The moments code knows not to change the position or the Mrf for external sources - // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; - - if (isfinite(petRadius)) { - source->extpars = pmSourceExtendedParsAlloc (); - source->extpars->petrosianRadius = petRadius; - } - + if (!source || ID != source->seq) { + if (source) { + psArrayAdd (sources, 1, source); + psFree(source); + } + source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + // RoughClass wants source type to be unknown + // source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + source->type = PM_SOURCE_TYPE_UNKNOWN; + + // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and + // stars that are masked + if (psfStar) { + source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; + } + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + source->seq = ID; + + PAR[PM_PAR_XPOS] = X; + PAR[PM_PAR_YPOS] = Y; + + dPAR[PM_PAR_XPOS] = 0.0; + dPAR[PM_PAR_YPOS] = 0.0; + + PAR[PM_PAR_SKY] = 0.0; + dPAR[PM_PAR_SKY] = 0.0; + + PAR[PM_PAR_I0] = 1.0; + dPAR[PM_PAR_I0] = 0.0; + + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + source->psfMag = 0.0; + source->psfMagErr = 0.0; + source->apMag = 0.0; + source->apRadius = apRadius; + + // we generate a somewhat fake PSF model here -- + // in most (all?) contexts, we will replace this with a measured psf model + // elsewhere + axes.major = 1.0; + axes.minor = 1.0; + axes.theta = 0.0; + pmPSF_AxesToModel (PAR, axes, modelType); + + // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second + source->peak = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE); + source->peak->xf = X; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->yf = Y; // pmPeakAlloc converts X,Y to int, so reset here + source->peak->dx = 0.0; + source->peak->dy = 0.0; + + source->extSN = SN; + + source->moments = pmMomentsAlloc (); + source->moments->Mx = X; + source->moments->My = Y; + source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment + + // Don't mark the moments as measured because that causes many fields to be left blank. + // The moments code knows not to change the position or the Mrf for external sources + // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + + if (isfinite(petRadius)) { + source->extpars = pmSourceExtendedParsAlloc (); + source->extpars->petrosianRadius = petRadius; + } + + } + bool saveExtModelParams = false; if (fitGalaxy && galaxyModelType >= 0) { - source->modelFits = psArrayAllocEmpty (1); - pmModel *model = pmModelAlloc(galaxyModelType); - psF32 *xPAR = model->params->data.F32; - - xPAR[PM_PAR_SKY] = 0.0; - xPAR[PM_PAR_I0] = 1.0; - xPAR[PM_PAR_XPOS] = X; + // skip likely stars + if (magDiff < starCut || SN < SNMinForCut) { + if (forceAll) { + saveExtModelParams = true; + } else if (galaxyModelType == modelTypeForce) { + // This is the model type that we are looking for + // proceed + saveExtModelParams = true; + } else if (modelTypeForce == sersicModelType && galaxyModelType == devModelType) { + // We're doing sersic models, if sersic index is greater than the recipe's minimum index + // do dev model as well + if (isfinite(forceDevSersicMin) || Sindex >= forceDevSersicMin) { + saveExtModelParams = true; + } + } else { + // not interested in this model + } + } + } + + if (saveExtModelParams) { + if (!source->modelFits) { + source->modelFits = psArrayAllocEmpty (1); + } + pmModel *model = pmModelAlloc(galaxyModelType); + psF32 *xPAR = model->params->data.F32; + + xPAR[PM_PAR_SKY] = 0.0; + xPAR[PM_PAR_I0] = 1.0; + xPAR[PM_PAR_XPOS] = X; xPAR[PM_PAR_YPOS] = Y; @@ -210,24 +278,19 @@ } + model->chisq = chisq; + model->nDOF = nDOF; + model->flags = modelFlags; + psArrayAdd (source->modelFits, 1, model); -#ifdef notyet - // XXX: set source->modelEXT to this model and flag as extended so that we have a better - // shot at subtracting extended sources? - - // This doesn't work right. Need to do some more work to flesh out the model before it can be - // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to - // First find the best extended models and then do everything else. - source->modelEXT = psMemIncrRefCounter(model); - // is this safe? (no see above) - source->type = PM_SOURCE_TYPE_EXTENDED; - source->mode |= PM_SOURCE_MODE_EXTMODEL; -#endif - psFree (model); } - sources->data[i] = source; psFree(row); + } + if (source) { + // close out last source + psArrayAdd (sources, 1, source); + psFree(source); } @@ -250,12 +313,30 @@ pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + pmModelType selectedModelType = -1; + bool chooseBest = false; + bool chooseAll = false; psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF"); - pmModelType selectedModelType = -1; + if (mdok && modelToChoose != NULL) { - if (strcmp(modelToChoose, "BEST")) { + if (!strcmp(modelToChoose, "BEST")) { + chooseBest = true; + } else if (!strcmp(modelToChoose, "ALL")) { + chooseAll = true; + } else if (strcmp(modelToChoose, "PS_MODEL_SERSIC")) { + // We have selected a model type other than Sersic. + // Save it's type for use below. Sersic is handled specially selectedModelType = pmModelClassGetType(modelToChoose); } } + + // minimum sersic index to force devModel + psF32 sersicMinDev = psMetadataLookupF32(&mdok, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + if (!mdok) { + sersicMinDev = NAN; + } + + sources = psArraySort (sources, pmSourceSortBySeq); for (int i = 0; i < sources->n; i++) { @@ -263,88 +344,197 @@ pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - psF32 xPos, yPos, flux, rMajor, rMinor, theta; - psS32 modelType = 0; + #define MAX_ROWS_PER_SRC 10 + psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC]; + psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC]; + psF32 chisq[MAX_ROWS_PER_SRC], nDOF[MAX_ROWS_PER_SRC]; + psS32 modelFlags[MAX_ROWS_PER_SRC]; + psS32 modelType[MAX_ROWS_PER_SRC]; + psF32 sersicIndex = NAN; bool fitGalaxy = false; bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; - psF32 sersicIndex = 0; - // For now only perform galaxy fits on extended objects - if (source->modelEXT == NULL) { - pmModel *model = source->modelPSF; - if (model == NULL) continue; - psF32 *PAR = model->params->data.F32; - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(source->psfFlux)) { - continue; - } - - xPos = model->params->data.F32[PM_PAR_XPOS]; - yPos = model->params->data.F32[PM_PAR_YPOS]; - flux = source->psfFlux; - rMajor = 0; - rMinor = 0; - theta = 0; - } else { - // Find the model with the selected type. If selected type is -1 choose the one selected ad - // modelEXT which was the best - int iModel = -1; - if (source->modelEXT) { - pmModelType ext_model_type = selectedModelType != -1 ? selectedModelType : source->modelEXT->type; + int n_rows = 0; + + psF32 kronFlux = source->moments->KronFlux; + psF32 SN = NAN; + psF32 magDiff = NAN; + if (isfinite(kronFlux) && isfinite(source->moments->KronFluxErr) && isfinite(source->psfMag)) { + SN = kronFlux/source->moments->KronFluxErr; + // kronMag - psfMag for use as star/glaxy separator + magDiff = -2.5 * log10(kronFlux) - source->psfMag ; + } + + // start with psf model + pmModel *model = source->modelPSF; + if (model == NULL) continue; + psF32 *PAR = model->params->data.F32; + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || + !isfinite(source->psfFlux)) { + continue; + } + + xPos[0] = model->params->data.F32[PM_PAR_XPOS]; + yPos[0] = model->params->data.F32[PM_PAR_YPOS]; + flux[0] = source->psfFlux; + rMajor[0] = 0; + rMinor[0] = 0; + theta[0] = 0; + modelType[0] = -1; + sersicIndex = NAN; + chisq[0] = NAN; + nDOF[0] = NAN; + modelFlags[0] = 0; + + if (source->modelFits != NULL) { + // figure out which models to use based on recipe paramters + if (chooseAll) { + // Save parameters for all valid extended models + + // but make sure we aren't going to overflow our arrays + assert (source->modelFits->n < MAX_ROWS_PER_SRC); + for (int j=0; jmodelFits->n; j++) { - pmModel *aModel = source->modelFits->data[j]; - if (aModel->type == ext_model_type) { - iModel = j; - break; + pmModel *model = source->modelFits->data[j]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[n_rows] = PAR[PM_PAR_XPOS]; + yPos[n_rows] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[n_rows] = axes.major; + rMinor[n_rows] = axes.minor; + theta[n_rows] = axes.theta*PS_DEG_RAD; + flux[n_rows] = pow(10.0, -0.4*model->mag); + modelType[n_rows] = model->type; + chisq[n_rows] = model->chisq; + nDOF[n_rows] = model->nDOF; + modelFlags[n_rows] = model->flags; + fitGalaxy = true; + if (model->params->n == 8) { + // this will save the sersic index in all models but that might + // be useful + sersicIndex = 0.5 / PAR[PM_PAR_7]; + } else { + sersicIndex = NAN; + } + + n_rows++; + } + } + } else { + int jModelSersic = -1; + int jModelDev = -1; + int jModelSelected = -1; + psF32 minChisq = NAN; + for (int j=0; jmodelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + if (chooseBest) { + // choose the model with lowest chisq + if (isfinite(model->chisq) && (!isfinite(minChisq) || model->chisq < minChisq)) { + jModelSelected = j; + minChisq = model->chisq; + } + } else { + // find the index of models of interest + if (model->type == selectedModelType) { + jModelSelected = j; + } else if (model->type == sersicModelType) { + jModelSersic = j; + } else if (model->type == devModelType) { + jModelDev = j; + } + } + } + if (jModelSelected >= 0 || jModelSersic >= 0) { + // If a specific non-sersic model we take paramers from that one. + // Otherwise we do the sersic model. + pmModel *model = jModelSelected >= 0 ? source->modelFits->data[jModelSelected] : + source->modelFits->data[jModelSersic]; + psF32 *PAR = model->params->data.F32; + + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[0] = PAR[PM_PAR_XPOS]; + yPos[0] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[0] = axes.major; + rMinor[0] = axes.minor; + theta[0] = axes.theta*PS_DEG_RAD; + flux[0] = pow(10.0, -0.4*model->mag); + modelType[0] = model->type; + chisq[0] = model->chisq; + nDOF[0] = model->nDOF; + modelFlags[0] = model->flags; + fitGalaxy = true; + if (model->type == sersicModelType) { + PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); + sersicIndex = 0.5 / PAR[PM_PAR_7]; + } + + n_rows = 1; + + // Unless a specific non-sersic model type was selected do dev model for sources with + // sersic index above the recipe limit. + if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) && + sersicIndex > sersicMinDev) { + + model = source->modelFits->data[jModelDev]; + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY]) && isfinite(PAR[PM_PAR_SXY]) && + isfinite(model->mag)) { + + xPos[1] = PAR[PM_PAR_XPOS]; + yPos[1] = PAR[PM_PAR_YPOS]; + + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + rMajor[1] = axes.major; + rMinor[1] = axes.minor; + theta[1] = axes.theta*PS_DEG_RAD; + flux[1] = pow(10.0, -0.4*model->mag); + modelType[1] = model->type; + sersicIndex = NAN; + chisq[1] = model->chisq; + nDOF[1] = model->nDOF; + modelFlags[1] = model->flags; + n_rows = 2; + } + } } } } - if (iModel == -1) { - // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits? - // XXX: Should this be an assert? - continue; - } - pmModel *model = source->modelFits->data[iModel]; - psF32 *PAR = model->params->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { - // XXX: do we need to handle this type - continue; - } else { - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY]) || - !isfinite(model->mag)) { - // bad model - continue; - } - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); - rMajor = axes.major; - rMinor = axes.minor; - theta = axes.theta*PS_DEG_RAD; - flux = pow(10.0, -0.4*model->mag); - fitGalaxy = true; - } - modelType = model->type; - if (modelType == sersicModelType) { - PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); - sersicIndex = 0.5 / PAR[PM_PAR_7]; - } - } - psMetadata *row = psMetadataAlloc(); - psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux/exptime); - psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); - psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); - psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); - psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); - psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor); - psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor); - psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType); - psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); - - psArrayAdd(table, 100, row); - psFree(row); + } + if (n_rows == 0) { + // didn't find a suitable extended model we have one row with the psf parameters + n_rows = 1; + } + + for (int j = 0; j < n_rows; j++) { + psMetadata *row = psMetadataAlloc(); + psMetadataAddU32 (row, PS_LIST_TAIL, "ID", 0, "IPP detection identifier", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X", 0, "x coordinate", xPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y", 0, "y coordinate", yPos[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux[j]/exptime); + psMetadataAddF32 (row, PS_LIST_TAIL, "SN", 0, "kron flux signal to noise", SN); + psMetadataAddF32 (row, PS_LIST_TAIL, "MAG_DIFF", 0, "psf mag - kron mag", magDiff); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); + psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS", 0, "Petrosian Radius", source->extpars ? source->extpars->petrosianRadius : NAN); + psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY", 0, "source has xfit", fitGalaxy); + psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR", 0, "source was psf star", psfStar); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR", 0, "radius of major axis", rMajor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR", 0, "radius of minor axis", rMinor[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); + psMetadataAddF32 (row, PS_LIST_TAIL, "CHISQ", 0, "chisq", chisq[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NDOF", 0, "n degrees of freedom", nDOF[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_FLAGS", 0, "model flags", modelFlags[j]); + + psArrayAdd(table, 100, row); + psFree(row); + } } Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36828) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36829) @@ -482,6 +482,8 @@ } +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif table = psArrayAllocEmpty (sources->n); @@ -821,6 +823,8 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif float magOffset; @@ -942,5 +946,5 @@ // EAM : adding for PV2 outputs: - @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0); + @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); @@ -1109,4 +1113,5 @@ model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ"); model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF"); + model->flags = psMetadataLookupS16(&status, row, "EXT_FLAGS"); // EXT_MODEL_TYPE gives the model chosen by psphot as the best. @@ -1225,6 +1230,8 @@ psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); +#ifdef SORT_EXTENSIONS_BY_FLUX // let's write these out in S/N order sources = psArraySort (sources, pmSourceSortByFlux); +#endif table = psArrayAllocEmpty (sources->n); @@ -1439,8 +1446,15 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); - psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE"); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); + psF32 Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q); + + psF32 NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma); + + psF32 clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); + psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN); + + // They are probably already in this order but ... + sources = psArraySort (sources, pmSourceSortBySeq); psArray *table = psArrayAllocEmpty (sources->n); @@ -1459,28 +1473,48 @@ if (source->galaxyFits == NULL) continue; - pmModel *model = source->modelFits->data[0]; - if (!model) return false; - - // X,Y coordinates are stored with the model parameters - psF32 *PAR = model->params->data.F32; - - psMetadata *row = psMetadataAlloc (); - - // we write out the x,y positions so people can link to the psf either way (position or ID) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT", 0, "model x coordinate", PAR[PM_PAR_XPOS]); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT", 0, "model y coordinate", PAR[PM_PAR_YPOS]); - psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX", 0, "number of pixels for fits", source->galaxyFits->nPix); - - psVector *Flux = source->galaxyFits->Flux; - psVector *dFlux = source->galaxyFits->dFlux; - psVector *chisq = source->galaxyFits->chisq; - - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX", PS_META_REPLACE, "normalization for galaxy flux", Flux); - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux); - psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ", PS_META_REPLACE, "galaxy fit chisq", chisq); - - psArrayAdd (table, 100, row); - psFree (row); + for (int iModel = 0; iModel < source->modelFits->n; iModel++) { + pmModel *model = source->modelFits->data[iModel]; + if (!model) continue; + + pmSourceGalaxyFits *galaxyFits = NULL; + for (int iFit = 0; iFit < source->galaxyFits->n; iFit++) { + galaxyFits = source->galaxyFits->data[iFit]; + if (model->type == galaxyFits->modelType) break; + galaxyFits = NULL; + } + + if (!galaxyFits) continue; + + // X,Y coordinates are stored with the model parameters + psF32 *PAR = model->params->data.F32; + + psMetadata *row = psMetadataAlloc (); + + // we write out the x,y positions so people can link to the psf either way (position or ID) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", galaxyFits->modelType); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT", 0, "model x coordinate", PAR[PM_PAR_XPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT", 0, "model y coordinate", PAR[PM_PAR_YPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX", 0, "number of pixels for fits", galaxyFits->nPix); + // psMetadataAddS32 (row, PS_LIST_TAIL, "REDUCED_TRIALS", 0, "source has reduced number of trials", galaxyFits->reducedTrials); + + psVector *Flux = galaxyFits->Flux; + psVector *dFlux = galaxyFits->dFlux; + psVector *chisq = galaxyFits->chisq; + + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX", PS_META_REPLACE, "normalization for galaxy flux", Flux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ", PS_META_REPLACE, "galaxy fit chisq", chisq); + + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MIN", 0, "fractional major axis min", galaxyFits->fRmajorMin); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MAX", 0, "fractional major axis max", galaxyFits->fRmajorMax); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_DEL", 0, "fractional major axis max", galaxyFits->fRmajorDel); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MIN", 0, "fractional minor axis min", galaxyFits->fRminorMin); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MAX", 0, "fractional minor axis max", galaxyFits->fRminorMax); + psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_DEL", 0, "fractional minor axis max", galaxyFits->fRminorDel); + + psArrayAdd (table, 100, row); + psFree (row); + } } @@ -1520,4 +1554,20 @@ return false; } + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + psF32 Q = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_Q"); + // XXX: turn this into an assert once we're done + if (status) { + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q); + + psF32 NSigma = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_NSIGMA"); + psAssert(status, "missing GALAXY_SHAPES_NSIGMA"); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma); + + psF32 clampSN = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_CLAMP_SN"); + psAssert(status, "missing GALAXY_SHAPES_CLAMP_SN"); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN); + } for (long i = 0; i < numSources; i++) { @@ -1541,4 +1591,5 @@ } + int modelType = psMetadataLookupS32 (&status, row, "MODEL_TYPE"); psVector *Flux = psMetadataLookupVector(&status, row, "GAL_FLUX"); psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR"); @@ -1546,14 +1597,29 @@ if (Flux && Flux->n > 0) { - psFree(source->galaxyFits); - source->galaxyFits = pmSourceGalaxyFitsAlloc(); - source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX"); - - psFree(source->galaxyFits->Flux); - source->galaxyFits->Flux = psMemIncrRefCounter(Flux); - psFree(source->galaxyFits->dFlux); - source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux); - psFree(source->galaxyFits->chisq); - source->galaxyFits->chisq = psMemIncrRefCounter(chisq); + if (!source->galaxyFits) { + source->galaxyFits = psArrayAllocEmpty(1); + } + + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + + psArrayAdd(source->galaxyFits, 1, galaxyFits); + + psFree(galaxyFits); + galaxyFits->modelType = modelType; + galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX"); + + galaxyFits->fRmajorMin = psMetadataLookupF32(&status, row, "FR_MAJOR_MIN"); + galaxyFits->fRmajorMax = psMetadataLookupF32(&status, row, "FR_MAJOR_MAX"); + galaxyFits->fRmajorDel = psMetadataLookupF32(&status, row, "FR_MAJOR_DEL"); + galaxyFits->fRminorMin = psMetadataLookupF32(&status, row, "FR_MINOR_MIN"); + galaxyFits->fRminorMax = psMetadataLookupF32(&status, row, "FR_MINOR_MAX"); + galaxyFits->fRminorDel = psMetadataLookupF32(&status, row, "FR_MINOR_DEL"); + + psFree(galaxyFits->Flux); + galaxyFits->Flux = psMemIncrRefCounter(Flux); + psFree(galaxyFits->dFlux); + galaxyFits->dFlux = psMemIncrRefCounter(dFlux); + psFree(galaxyFits->chisq); + galaxyFits->chisq = psMemIncrRefCounter(chisq); }