Index: branches/eam_branches/ipp-20140423/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20140423/psphot/src/psphot.h (revision 36827) +++ branches/eam_branches/ipp-20140423/psphot/src/psphot.h (revision 36828) @@ -432,10 +432,8 @@ typedef struct { - float fRmajorMin; - float fRmajorMax; - float fRmajorDel; - float fRminorMin; - float fRminorMax; - float fRminorDel; + float Q; + float NSigma; + float clampSN; + int extModelType; } psphotGalaxyShapeOptions; @@ -537,6 +535,7 @@ bool psphotGalaxyShape_Threaded (psThreadJob *job); bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize); -bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults); +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults); psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc(); +bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *defaultOptions); bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule); Index: branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummary.c =================================================================== --- branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummary.c (revision 36827) +++ branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummary.c (revision 36828) @@ -186,13 +186,27 @@ psFitsClose(fits); - // Convert to an array indexed by ID - psArray *sortedTable = psArrayAlloc(4*inTable->n); + // Convert to a set of arrays indexed by model type + 1 + // which contain pointers to arrays indexed by ID +#define MAX_MODEL_TYPE 10 + psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1); for (int i=0; in; i++) { psMetadata *row = inTable->data[i]; psS32 ID = psMetadataLookupS32(&status, row, "ID"); + psS32 modelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); + // XXX: need to use the lookup table functions to be ready for changes in the model type numbers + if (modelType+1 >= MAX_MODEL_TYPE) { + psError(PSPHOT_ERR_CONFIG, false, "found modelType %d max allowed is %d", modelType, MAX_MODEL_TYPE); + return false; + } + psArray *sortedTable = sortedTables->data[modelType+1]; + if (!sortedTable) { + sortedTable = psArrayAlloc(4*inTable->n); + sortedTables->data[modelType+1] = sortedTable; + // dont' free sortedTable the array of tables gets our reference + } if (ID >= sortedTable->n) { sortedTable = psArrayRealloc(sortedTable, 2*ID); // Why doesn't psArrayRealloc do this????? - sortedTable->n = sortedTable->nalloc; + sortedTable->n = sortedTable->nalloc; } if (sortedTable->data[ID]) { @@ -202,7 +216,7 @@ sortedTable->data[ID] = psMemIncrRefCounter(row); } - psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLES", PS_META_REPLACE, "cff tables", sortedTables); psFree(inTable); - psFree(sortedTable); + psFree(sortedTables); return true; Index: branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummaryReadout.c (revision 36827) +++ branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummaryReadout.c (revision 36828) @@ -3,23 +3,31 @@ typedef struct { + psF32 Q; + psF32 NSigma; + psF32 clampSN; +#ifdef notdef + // these are per object now int numTrials; - psF32 fRmajorMin; - psF32 fRmajorMax; - psF32 fRmajorDel; - psF32 fRminorMin; - psF32 fRminorMax; - psF32 fRminorDel; - psVector *rMajor; - psVector *rMinor; - psArray *zeroPt; - psArray *exptime; - psArray *cffTable; -} galaxyShapeOptions; - -static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); + psF64 fRmajorMin; + psF64 fRmajorMax; + psF64 fRmajorDel; + psF64 fRminorMin; + psF64 fRminorMax; + psF64 fRminorDel; + psVector *fRmajor; + psVector *fRminor; +#endif + psArray *zeroPt; // zero points for each input used with exptime to scale flux + psArray *exptime; // exposure times for each input + psArray *cffTables; // one for each model type index is (model_type + 1) entry 0 is no extended model (star) +} galaxyShapeSummaryOptions; + + + +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options); static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source); -static bool setOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors); -static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe); +static bool setOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors); +static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe); @@ -47,13 +55,10 @@ psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", outputExptime); - // don't think this one matters -// psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", - // outputExptime); // Create objects from the various input's sources // loop over the available readouts -// psArray *readouts = psArrayAlloc(num); + psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32); - galaxyShapeOptions options; + galaxyShapeSummaryOptions options; psArray *objects = NULL; for (int index = 0; index < num; index++) { @@ -72,6 +77,6 @@ return false; } - options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE"); - if (!options.cffTable) { + options.cffTables = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLES"); + if (!options.cffTables) { psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments."); return false; @@ -91,6 +96,4 @@ psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime; fluxScaleFactor->data.F32[index] = scaleFactor; - -// readouts->data[index] = readout; // find detections @@ -177,15 +180,27 @@ } -static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options) { +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options) { pmSource *outSrc = NULL; - psVector *sumWeightedFlux = NULL; - psVector *sumInvSig2 = NULL; - psVector *numerator = NULL; -// psVector *tmp = NULL; - psF32 totalNPix = 0; - long vectorLength = 0; - + + pmSource *source0 = obj->sources->data[0]; + int ID = source0->seq; + + // Find the parameters from the cff tables. nModels is the number of extended model types that + // this source has an entry for + int nModels = 0; + // skip table 0 which is modelType -1 not extended + for (int i=1; icffTables->n; i++) { + psArray *table = options->cffTables->data[i]; + if (!table) continue; + if (ID >= table->n) continue; + if (table->data[ID]) { + nModels++; + } + } + + // Loop over sources for this object. Start output source for first input that has + // good pixel frac for (int i=0; i < obj->sources->n; i++) { pmSource *source = obj->sources->data[i]; @@ -212,87 +227,193 @@ outSrc->extpars = psMemIncrRefCounter(source->extpars); } - if (source->modelEXT) { - outSrc->modelEXT = psMemIncrRefCounter(source->modelEXT); - } - if (source->modelFits) { - outSrc->modelFits = psMemIncrRefCounter(source->modelFits); - } - } - - if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) { - if (numerator == NULL) { - vectorLength = source->galaxyFits->chisq->n; - // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32); - sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(sumWeightedFlux, 0.0); - sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(sumInvSig2, 0.0); - numerator = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(numerator, 0.0); - } - - psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d", - vectorLength, options->numTrials); - psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld", - source->galaxyFits->chisq->n, vectorLength); - - - psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + } + if (nModels && source->modelFits && source->modelFits->n) { + if (!outSrc->modelFits) { + outSrc->modelFits = psArrayAllocEmpty(nModels); + } + for (int iModel = 0; iModel < source->modelFits->n; iModel++) { + pmModel *newModel = source->modelFits->data[iModel]; + pmModel *outModel = NULL; + for (int jModel = 0; jModel < outSrc->modelFits->n; jModel++) { + outModel = outSrc->modelFits->data[jModel]; + if (newModel->type == outModel->type) { + // already have an entry for this type + break; + } + outModel = NULL; + } + if (!outModel) { + // no model of this type in output source yet. Add this one. + // Note: the parameters that we need (position, angle, index, if applicable) + // are the same for all sources so copying any of them is fine + psArrayAdd(outSrc->modelFits, 1, newModel); + } + } + } + } + + if (!outSrc) { + // no good measurements for this source (how?) + return NULL; + } + + // now loop over the model fits and summarize the galaxy shape measurements for each + for (int iModel=0; outSrc->modelFits && iModel < outSrc->modelFits->n; iModel++) { + pmModel *model = outSrc->modelFits->data[iModel]; + int modelType = model->type; + psAssert(modelType >= 0 && modelType+1 < options->cffTables->n, "model type out of range"); + + psArray *cffTable = options->cffTables->data[modelType+1]; + psAssert(cffTable, "missing cff table?"); + + psMetadata *cffRow = cffTable->data[outSrc->seq]; + if (!cffRow) { + psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d model %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + bool mdok; + psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR"); + if (!mdok) { + psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d type %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR"); + if (!mdok) { + psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d type %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + + psVector *sumWeightedFlux = NULL; + psVector *sumInvSig2 = NULL; + psVector *numerator = NULL; + psVector *fRmajorVec = NULL; + psVector *fRminorVec = NULL; + psF32 totalNPix = 0; + long vectorLength = 0; + psF32 fRmajorMin = 0; + psF32 fRmajorMax = 0; + psF32 fRmajorDel = 0; + psF32 fRminorMin = 0; + psF32 fRminorMax = 0; + psF32 fRminorDel = 0; + for (int i=0; i < obj->sources->n; i++) { + pmSource *source = obj->sources->data[i]; + + if (source->pixWeightNotPoor < .9) continue; + + // The only parameters that we are summarizing currently are galaxy shapes. + if (!source->modelFits) continue; + if (!source->galaxyFits) continue; + + // XXX: put this into a function + pmSourceGalaxyFits *galaxyFits = NULL; + for (int j=0; jgalaxyFits->n; j++) { + galaxyFits = source->galaxyFits->data[j]; + if (galaxyFits && galaxyFits->modelType == modelType) break; + galaxyFits = NULL; + } + + if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) { + if (numerator == NULL) { + // first source with galaxyFits allocate accumulators + vectorLength = galaxyFits->chisq->n; + sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumWeightedFlux, 0.0); + sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumInvSig2, 0.0); + numerator = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(numerator, 0.0); + + // Initialize the fractional radius vectors + // save these to initialize the output galaxy fits + fRmajorMin = galaxyFits->fRmajorMin; + fRmajorMax = galaxyFits->fRmajorMax; + fRmajorDel = galaxyFits->fRmajorDel; + fRminorMin = galaxyFits->fRminorMin; + fRminorMax = galaxyFits->fRminorMax; + fRminorDel = galaxyFits->fRminorDel; + fRmajorVec = psVectorAlloc(vectorLength, PS_TYPE_F32); + fRminorVec = psVectorAlloc(vectorLength, PS_TYPE_F32); + int k = 0; + for (float fRmajor = fRmajorMin; fRmajor < fRmajorMax + 0.5*fRmajorDel; fRmajor += fRmajorDel) { + for (float fRminor = fRminorMin; fRminor < fRminorMax + 0.5*fRminorDel; fRminor += fRminorDel) { + fRminorVec->data.F32[k] = fRminor; + fRmajorVec->data.F32[k] = fRmajor; + k++; + } + } + psAssert(k == vectorLength, "mismatched vectors"); + } + + // Die if the lengths of the vectors is not the same in all sources +#ifdef notdef + psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d", + vectorLength, options->numTrials); +#endif + psAssert(galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld", + galaxyFits->chisq->n, vectorLength); + + psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + + totalNPix += galaxyFits->nPix; + + for (int k = 0; k < vectorLength; k++) { + psF32 chisq = galaxyFits->chisq->data.F32[k]; + psF32 flux = galaxyFits->Flux->data.F32[k] * scaleFactor; + psF32 dFlux = galaxyFits->dFlux->data.F32[k] * scaleFactor; + + numerator->data.F32[k] += chisq * galaxyFits->nPix; + + psF32 invSig2 = 1.0 / (dFlux * dFlux); + sumInvSig2->data.F32[k] += invSig2; + + sumWeightedFlux->data.F32[k] += flux * invSig2; + } + } + } + + if (vectorLength) { + // allocate galaxyFits for the output source + if (!outSrc->galaxyFits) { + outSrc->galaxyFits = psArrayAllocEmpty(1); + } + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + psArrayAdd(outSrc->galaxyFits, 1, galaxyFits); + psFree(galaxyFits); + + galaxyFits->nPix = totalNPix; + galaxyFits->modelType = modelType; + psVector *fluxVec = + galaxyFits->Flux = psVectorRecycle(galaxyFits->Flux, vectorLength, PS_TYPE_F32); + psVector *dFluxVec = + galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32); + psVector *chisqVec = + galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32); + + galaxyFits->fRmajorMin = fRmajorMin; + galaxyFits->fRmajorMax = fRmajorMax; + galaxyFits->fRmajorDel = fRmajorDel; + galaxyFits->fRminorMin = fRminorMin; + galaxyFits->fRminorMax = fRminorMax; + galaxyFits->fRminorDel = fRminorDel; + + // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq + + int min_k = -1; + psF32 minChisq = NAN; for (int k = 0; k < vectorLength; k++) { - - psF32 chisq = source->galaxyFits->chisq->data.F32[k]; - psF32 flux = source->galaxyFits->Flux->data.F32[k] * scaleFactor; - psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor; - // for chisq - // Numerator = sum( nPix[i] * chisq[i] ) - // denominator will be sum( nPix ) - - -// tmp = (psVector *) psBinaryOp(tmp, - // source->galaxyFits->chisq, "*", psScalarAlloc(source->galaxyFits->nPix, PS_TYPE_F32)); -// numerator = (psVector *) psBinaryOp(numerator, numerator, "+", tmp); - - numerator->data.F32[k] += chisq * source->galaxyFits->nPix; - totalNPix += source->galaxyFits->nPix; - - // sumInvSig2 = sum( 1 / dFlux**2 ) -// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->dFlux, "*", source->galaxyFits->dFlux); -// tmp = (psVector *) psBinaryOp(tmp, psScalarAlloc(1.0, PS_TYPE_F32), "/", tmp); -// sumInvSig2 = (psVector *) psBinaryOp(sumInvSig2, sumInvSig2, "+", tmp); -// - psF32 invSig2 = 1.0 / (dFlux * dFlux); - sumInvSig2->data.F32[k] += invSig2; - - // sumWeightedFlux = sum ( Flux / dFlux**2 ) -// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp); -// sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp); - sumWeightedFlux->data.F32[k] += flux * invSig2; - } - } - } - - if (outSrc) { - if (vectorLength) { - outSrc->galaxyFits = pmSourceGalaxyFitsAlloc(); - outSrc->galaxyFits->nPix = totalNPix; - outSrc->galaxyFits->Flux = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32); - outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32); - outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32); - for (int k = 0; k < vectorLength; k++) { - outSrc->galaxyFits->Flux->data.F32[k] = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k]; - outSrc->galaxyFits->dFlux->data.F32[k] = 1.0 / sumInvSig2->data.F32[k]; - outSrc->galaxyFits->chisq->data.F32[k] = numerator->data.F32[k] / totalNPix; - -#ifdef nodef - outSrc->galaxyFits->Flux = (psVector *) psBinaryOp(outSrc->galaxyFits->Flux, - sumWeightedFlux, "/", sumInvSig2); - - outSrc->galaxyFits->dFlux = (psVector *) psBinaryOp(outSrc->galaxyFits->dFlux, - psScalarAlloc(1.0, PS_TYPE_F32), "/", sumInvSig2); - - outSrc->galaxyFits->chisq = (psVector *) psBinaryOp(outSrc->galaxyFits->chisq, - numerator, "/", psScalarAlloc(totalNPix, PS_TYPE_F32)); -#endif + fluxVec->data.F32[k] = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k]; + dFluxVec->data.F32[k] = 1.0 / sumInvSig2->data.F32[k]; + + psF32 thischisq = chisqVec->data.F32[k] = numerator->data.F32[k] / totalNPix; + + if (isfinite(thischisq) && (!isfinite(minChisq) || thischisq < minChisq)) { + min_k = k; + minChisq = thischisq; + } } @@ -300,52 +421,230 @@ psFree(sumInvSig2); psFree(sumWeightedFlux); -// psFree(tmp); - - int min_j = -1; - psF32 minChisq = NAN; - psVector *chisq = outSrc->galaxyFits->chisq; - for (int j=0; j < chisq->n; j++) { - psF32 thischisq = chisq->data.F32[j]; - if (isfinite(thischisq) && (!isfinite(minChisq) || thischisq < minChisq)) { - min_j = j; - minChisq = thischisq; + + if (min_k >= 0 && isfinite(minChisq)) { + // copy the best fit params to the model + // fractional radii with the lowest chisq + psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->type); + + // examine the params for the trial with minimum chisq. + bool fitMajor = true; + bool fitMinor = true; + psF64 fRmajorBest = fRmajorVec->data.F32[min_k]; + if ((fabs(fRmajorBest - galaxyFits->fRmajorMin) < galaxyFits->fRmajorDel) || + (fabs(fRmajorBest - galaxyFits->fRmajorMax) < galaxyFits->fRmajorDel)) { + fitMajor = false; } - } - if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) { - // copy the best fit params to the model - psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type); - if (outSrc->seq >= options->cffTable->n) { - psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld", - outSrc->seq, options->cffTable->n); - return NULL; + psF64 fRminorBest = fRminorVec->data.F32[min_k]; + if ((fabs(fRminorBest - galaxyFits->fRminorMin) < galaxyFits->fRminorDel) || + (fabs(fRminorBest - galaxyFits->fRminorMax) < galaxyFits->fRminorDel)) { + fitMinor = false; } - psMetadata *row = options->cffTable->data[outSrc->seq]; - if (!row) { - psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq); - psFree(outSrc); - return NULL; + // If either major or minor is at one of the limits do not fit report the minimum value + bool useFit = fitMajor && fitMinor; + + // save the flux and dFlux values from entry with lowest chisq + psF64 fluxBest = fluxVec->data.F32[min_k]; + psF64 dFluxBest = dFluxVec->data.F32[min_k]; + psF64 dFlux0 = NAN; + psF64 flux0 = NAN; + + if (useFit) { + #define NUM_TRIALS_INIT 9 + psVector *major = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *minor = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *flux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + + // XXX: use a recipe parameter instead of 2.2 + psF64 maxDeltaMaj = 2.2 * galaxyFits->fRmajorDel; + psF64 maxDeltaMin = 2.2 * galaxyFits->fRminorDel; + + psF64 majorMin = NAN; + psF64 majorMax = NAN; + psF64 minorMin = NAN; + psF64 minorMax = NAN; + psF64 chisqMin = NAN; + psF64 chisqMax = NAN; + for (int k = 0; k < vectorLength; k++) { + + if (fabs(fRmajorVec->data.F32[k] - fRmajorBest) < maxDeltaMaj && + fabs(fRminorVec->data.F32[k] - fRminorBest) < maxDeltaMin) { + + if (isfinite(chisqVec->data.F32[k]) && + isfinite(fluxVec->data.F32[k]) && + isfinite(dFluxVec->data.F32[k])) { + + // compute major and minor radius vectors from nominal and trial fractions + // also find the ranges in the vectors + psF64 thisMajor = R_MAJOR * fRmajorVec->data.F32[k]; + if (!isfinite(majorMin) || thisMajor < majorMin) { + majorMin = thisMajor; + } + if (!isfinite(majorMax) || thisMajor > majorMax) { + majorMax = thisMajor; + } + psVectorAppend(major, thisMajor); + + psF64 thisMinor = R_MINOR * fRminorVec->data.F32[k]; + if (!isfinite(minorMin) || thisMinor < minorMin) { + minorMin = thisMinor; + } + if (!isfinite(minorMax) || thisMinor > minorMax) { + minorMax = thisMinor; + } + psVectorAppend(minor, thisMinor); + + psF64 thisChisq = chisqVec->data.F32[k]; + if (!isfinite(chisqMin) || thisChisq < chisqMin) { + chisqMin = thisChisq; + } + if (!isfinite(chisqMax) || thisChisq > chisqMax) { + chisqMax = thisChisq; + } + psVectorAppend(chisq, thisChisq); + + psVectorAppend(flux, fluxVec->data.F32[k]); + psVectorAppend(dFlux, dFluxVec->data.F32[k]); + } + } + } + + + // see if we ever get too few good values (haven't seen this happen) + if (major->n < NUM_TRIALS_INIT) { + fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n", + major->n, NUM_TRIALS_INIT, outSrc->seq); + } + + // rescale data - this helps avoid precision errors. + psF64 majorRange = majorMax - majorMin; + if (majorRange == 0) { + majorRange = 1; + } + psF64 minorRange = minorMax - minorMin; + if (minorRange == 0) { + minorRange = 1; + } + psF64 chisqRange = chisqMax - chisqMin; + if (chisqRange == 0) { + chisqRange = 1; + } + for (int k = 0; k < major->n; k++) { + major->data.F64[k] = (major->data.F64[k] - majorMin) / majorRange; + minor->data.F64[k] = (minor->data.F64[k] - minorMin) / minorRange; + chisq->data.F64[k] = (chisq->data.F64[k] - chisqMin) / chisqRange; + } + + // Fit chisq versus rMajor and rMinor + psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + psPolynomial2D *fluxErrorFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET; + chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET; + chisqFit->coeffMask[1][2] = PS_POLY_MASK_SET; + fluxFit->coeffMask[2][2] = PS_POLY_MASK_SET; + fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET; + fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[2][2] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[2][1] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[1][2] = PS_POLY_MASK_SET; + + bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor); + + if (goodFit) { + // find the axes for the minimum chisq from the formula + psF64 **coeff = chisqFit->coeff; + psF64 D = 4*coeff[2][0]*coeff[0][2] - coeff[1][1]*coeff[1][1]; + psF64 major0 = ((coeff[1][1]*coeff[0][1] - 2*coeff[0][2]*coeff[1][0]) / D) * majorRange + majorMin; + psF64 minor0 = ((coeff[1][1]*coeff[1][0] - 2*coeff[2][0]*coeff[0][1]) / D) * minorRange + minorMin; + + // estimated chisq at the minimum + psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin; + + // now fit the flux ... + bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor); + + // .. and compute flux at the minimum chisq position + if (goodFluxFit) { + flux0 = psPolynomial2DEval(fluxFit, major0, minor0); + } else { + flux0 = fluxBest; + } + + // .. and compute dFlux at the minimum chisq position + bool goodFluxErrorFit = psVectorFitPolynomial2D(fluxErrorFit, NULL, 0xFF, dFlux, NULL, major, minor); + if (goodFluxErrorFit) { + dFlux0 = psPolynomial2DEval(fluxErrorFit, major0, minor0); + } else { + dFlux0 = dFluxBest; + } + +#ifdef PRINTVALS + fprintf (stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %4.2f | %6.3f %6.3f %6.3f %4.2f | %7.4f %7.4f %7.1f %7.1f\n", + outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest*R_MAJOR, major0, major0/R_MAJOR, R_MINOR, fRminorBest*R_MINOR, minor0, minor0/R_MINOR, minChisq, chisq0, fluxBest, flux0); +#endif + + + axes.major = major0; + axes.minor = minor0; + model->chisq = chisq0; +#ifdef DUMPVECTORS + char fn[80]; + + sprintf(fn, "vectors/%s.%05d.txt", goodFit ? "g" : "b", outSrc->seq); + + FILE *f = fopen(fn, "w"); + fprintf(f, "#major minor chisq i\n"); + for (int k = 0; k < chisq->n; k++) { + fprintf(f, "%10.6f %10.6f %10.6f %5d\n", + major->data.F64[k]/R_MAJOR, minor->data.F64[k]/R_MINOR, chisq->data.F64[k], k); + } + fclose(f); +#endif + } else { + model->flags |= PM_MODEL_STATUS_NONCONVERGE; + +#ifdef PRINTVALS + fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f bad fit\n", + outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); +#endif + psErrorClear(); + // psFree(outSrc->galaxyFits); + + } + psFree(chisqFit); + psFree(fluxFit); + psFree(fluxErrorFit); + psFree(major); + psFree(minor); + psFree(chisq); + psFree(flux); + psFree(dFlux); + } else { + // No fit + // Set the values based on the trial with the best chisq + axes.major = R_MAJOR * fRmajorBest; + axes.minor = R_MINOR * fRminorBest; + model->chisq = minChisq; + flux0 = fluxBest; + dFlux0 = dFluxBest; +#ifdef PRINTVALS + fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f skip fit\n", + outSrc->seq, min_k, 0L, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); +#endif } - bool mdok; - psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR"); - if (!mdok) { - psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq); - psFree(outSrc); - return NULL; - } - psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR"); - if (!mdok) { - psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq); - psFree(outSrc); - return NULL; - } - axes.major = rMajor * options->rMajor->data.F32[min_j]; - axes.minor = rMinor * options->rMinor->data.F32[min_j]; - pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type); - outSrc->modelEXT->chisq = minChisq; - outSrc->modelEXT->mag = -2.5 * log10(outSrc->galaxyFits->Flux->data.F32[min_j]); - outSrc->modelEXT->magErr = 1.0 / - (outSrc->galaxyFits->Flux->data.F32[min_j]/outSrc->galaxyFits->dFlux->data.F32[min_j]); // 1 / SN - } - } + // now save the model parameters in the model structure + + pmPSF_AxesToModel (model->params->data.F32, axes, modelType); + + model->mag = -2.5 * log10(flux0); + model->magErr = dFlux0 / flux0; // 1 / SN + // XXX: should there be a different flag if we didn't do a fit of the chisq + model->flags |= PM_MODEL_STATUS_FITTED; + } + } + psFree(fRminorVec); + psFree(fRmajorVec); } @@ -360,11 +659,20 @@ } -static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) { +static bool setOptions(galaxyShapeSummaryOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) { bool status; - bool useAnalysis; - + bool useAnalysis; // fall back to recipe if we dont' find values in analysis. Probably should no longer do this + +#ifdef notdef psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN"); +#endif + psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_Q"); psMetadata *md = useAnalysis ? readout->analysis : recipe; + GETVAL(Q, "GALAXY_SHAPES_Q"); + GETVAL(NSigma, "GALAXY_SHAPES_NSIGMA"); + GETVAL(clampSN, "GALAXY_SHAPES_CLAMP_SN"); + +#ifdef notdef + // these are per object now GETVAL(fRmajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN"); GETVAL(fRmajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX"); @@ -378,6 +686,6 @@ if (makeVectors) { - opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); - opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->fRminor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->fRmajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); int i = 0; for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; @@ -385,6 +693,6 @@ for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { - opt->rMinor->data.F32[i] = fRminor; - opt->rMajor->data.F32[i] = fRmajor; + opt->fRminor->data.F32[i] = fRminor; + opt->fRmajor->data.F32[i] = fRmajor; i++; } @@ -392,7 +700,8 @@ psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials); } else { - opt->rMinor = NULL; - opt->rMajor = NULL; - } + opt->fRminor = NULL; + opt->fRmajor = NULL; + } +#endif return true; @@ -405,6 +714,6 @@ } -static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe) { - galaxyShapeOptions thisReadoutsOptions; +static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe) { + galaxyShapeSummaryOptions thisReadoutsOptions; if (!setOptions(&thisReadoutsOptions, readout, recipe, false)) { @@ -412,4 +721,10 @@ return false; } + CHECKVAL(options, thisReadoutsOptions, Q, "mismatched Q"); + CHECKVAL(options, thisReadoutsOptions, NSigma, "mismatched NSIGMA"); + CHECKVAL(options, thisReadoutsOptions, clampSN, "mismatched cleampSN"); + +#ifdef notdef + // these are per object now CHECKVAL(options, thisReadoutsOptions, numTrials, "mismatched number of trials") CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin") @@ -419,4 +734,5 @@ CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax") CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel") +#endif return true; Index: branches/eam_branches/ipp-20140423/psphot/src/psphotGalaxyShape.c =================================================================== --- branches/eam_branches/ipp-20140423/psphot/src/psphotGalaxyShape.c (revision 36827) +++ branches/eam_branches/ipp-20140423/psphot/src/psphotGalaxyShape.c (revision 36828) @@ -80,9 +80,21 @@ assert (markVal); - // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; - // bit-mask to mark pixels not used in analysis psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc(); + opt->Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); + psAssert (status, "missing GALAXY_SHAPES_Q"); + opt->NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); + psAssert (status, "missing GALAXY_SHAPES_NSIGMA"); + opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); + psAssert (status, "missing GALAXY_SHAPES_NSIGMA"); + psString modelTypeToSave = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + if (modelTypeToSave && strcmp(modelTypeToSave, "ALL") && strcmp(modelTypeToSave, "BEST")) { + opt->extModelType = pmModelClassGetType(modelTypeToSave); + } else { + opt->extModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + } + +#ifdef notdef opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5; opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0; @@ -91,4 +103,5 @@ opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0; opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1; +#endif // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) @@ -238,89 +251,97 @@ bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) { - // we use the 0th model to define the initial guess shape - pmModel *model = source->modelFits->data[0]; - if (!model) return false; - - pmModelType modelType = model->type; - - // we are using fitOptions->mode : be sure this makes sense - pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); - if (!pcm) return false; - - // we are fitting only PM_PAR_I0; the shape elements are generated from a grid - psF32 *PAR = pcm->modelConv->params->data.F32; - - // double check that the guess is carried along... - - // I have some source guess (e0, e1, e2) - psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); - - if (!source->galaxyFits) { - source->galaxyFits = pmSourceGalaxyFitsAlloc(); - } - - float fRmajorBest = NAN; - float fRminorBest = NAN; - float chisqBest = NAN; - for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) { - for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { - - psEllipseAxes testAxes = guessAxes; - testAxes.major = guessAxes.major * fRmajor; - testAxes.minor = guessAxes.minor * fRminor; - - pmPSF_AxesToModel (PAR, testAxes, modelType); - - psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true); - - int i = source->galaxyFits->chisq->n - 1; - float flux = source->galaxyFits->Flux->data.F32[i]; - if (isfinite(flux)) { - float thisChisq = source->galaxyFits->chisq->data.F32[i]; - if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) { - chisqBest = thisChisq; - fRmajorBest = fRmajor; - fRminorBest = fRminor; + + for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) { + pmModel *model = source->modelFits->data[iModel]; + if (!model) return false; + + pmModelType modelType = model->type; + model->flags = PM_MODEL_STATUS_NONE; + + // we are using fitOptions->mode : be sure this makes sense + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (!pcm) return false; + + // we are fitting only PM_PAR_I0; the shape elements are generated from a grid + psF32 *PAR = pcm->modelConv->params->data.F32; + + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + + // set the options for this source + if (!psphotGalaxyShapeOptionsSet(source, galaxyFits, opt)) { + // this source won't work + psFree(pcm); + psFree(galaxyFits); + continue; + } + + if (!source->galaxyFits) { + source->galaxyFits = psArrayAllocEmpty(1); + } + psArrayAdd(source->galaxyFits, 1, galaxyFits); + psFree(galaxyFits); + + // I have some source guess (e0, e1, e2) + psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); + + float fRmajorBest = NAN; + float fRminorBest = NAN; + float chisqBest = NAN; + for (float fRmajor = galaxyFits->fRmajorMin; fRmajor < galaxyFits->fRmajorMax + 0.5*galaxyFits->fRmajorDel; fRmajor += galaxyFits->fRmajorDel) { + for (float fRminor = galaxyFits->fRminorMin; fRminor < galaxyFits->fRminorMax + 0.5*galaxyFits->fRminorDel; fRminor += galaxyFits->fRminorDel) { + + psEllipseAxes testAxes = guessAxes; + testAxes.major = guessAxes.major * fRmajor; + testAxes.minor = guessAxes.minor * fRminor; + + pmPSF_AxesToModel (PAR, testAxes, modelType); + + psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true); + + int i = galaxyFits->chisq->n - 1; + float flux = galaxyFits->Flux->data.F32[i]; + if (isfinite(flux)) { + float thisChisq = galaxyFits->chisq->data.F32[i]; + if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) { + chisqBest = thisChisq; + fRmajorBest = fRmajor; + fRminorBest = fRminor; + } } + // reset I0 to avoid potential problems on the next iteration + PAR[PM_PAR_I0] = 1.0; } - // reset I0 to avoid potential problems on the next iteration - PAR[PM_PAR_I0] = 1.0; - } - } - -#define SAVE_BEST_MODEL -#ifdef SAVE_BEST_MODEL - // Save model with smallest chisq - if (isfinite(chisqBest)) { -#else - // Save model with nominal parameters - { - fRmajorBest = 1; - fRminorBest = 1; + } + + if (isfinite(chisqBest)) { + // now save the best fitting model as the source's extended model ... + psEllipseAxes testAxes = guessAxes; + + // ... unless this macro is defined +#ifndef SAVE_NOMINAL_MODEL + testAxes.major = guessAxes.major * fRmajorBest; + testAxes.minor = guessAxes.minor * fRminorBest; #endif - // now save the best fitting model as the source's extended model - psEllipseAxes testAxes = guessAxes; - testAxes.major = guessAxes.major * fRmajorBest; - testAxes.minor = guessAxes.minor * fRminorBest; - - pmPSF_AxesToModel (PAR, testAxes, modelType); - - psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false); - - // Replace modelEXT with this model, if the model is good - if (isfinite(PAR[PM_PAR_I0])) { - psFree (source->modelEXT); - - source->modelEXT = psMemIncrRefCounter (pcm->modelConv); - source->type = PM_SOURCE_TYPE_EXTENDED; - source->mode |= PM_SOURCE_MODE_EXTMODEL; - source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; - - // cache the model flux - pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); - } - } - - psFree (pcm); + + pmPSF_AxesToModel (PAR, testAxes, modelType); + + psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false); + + // Replace modelEXT with the best model from the first of the model fits if one of them is good + if (isfinite(PAR[PM_PAR_I0]) && modelType == opt->extModelType) { + psFree (source->modelEXT); + + source->modelEXT = psMemIncrRefCounter (pcm->modelConv); + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + + // cache the model flux + pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); + } + } + + psFree (pcm); + } return true; } @@ -328,5 +349,5 @@ // fit the given model to the source and find chisq & normalization // XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta? -bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults) { +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -384,8 +405,9 @@ if (saveResults) { - psVectorAppend (source->galaxyFits->Flux, flux); - psVectorAppend (source->galaxyFits->dFlux, dflux); - psVectorAppend (source->galaxyFits->chisq, Chisq); - source->galaxyFits->nPix = nPix; + psVectorAppend (galaxyFits->Flux, flux); + psVectorAppend (galaxyFits->dFlux, dflux); + psVectorAppend (galaxyFits->chisq, Chisq); + galaxyFits->nPix = nPix; + galaxyFits->modelType = pcm->modelConv->type; } @@ -404,12 +426,28 @@ psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree); - opt->fRmajorMin = 0.5; - opt->fRmajorMax = 2.0; - opt->fRmajorDel = 0.1; - opt->fRminorMin = 0.5; - opt->fRminorMax = 2.0; - opt->fRminorDel = 0.1; - return opt; } +bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *options) +{ + // XXX: put these in recipe + // doesn't make sense to use too small of a spacing + float clampedSN = source->extSN < options->clampSN ? source->extSN : options->clampSN; + + float f_del = options->Q / clampedSN; + float f_min = 1 - options->NSigma * f_del; + float f_max = 1 + options->NSigma * f_del; + + // if f_min goes negative skip this object + if (f_min < 0) return false; + + galaxyFits->fRmajorMin = f_min; + galaxyFits->fRmajorMax = f_max; + galaxyFits->fRmajorDel = f_del; + galaxyFits->fRminorMin = f_min; + galaxyFits->fRminorMax = f_max; + galaxyFits->fRminorDel = f_del; + + return true; +} + Index: branches/eam_branches/ipp-20140423/psphot/src/psphotKronIterate.c =================================================================== --- branches/eam_branches/ipp-20140423/psphot/src/psphotKronIterate.c (revision 36827) +++ branches/eam_branches/ipp-20140423/psphot/src/psphotKronIterate.c (revision 36828) @@ -532,4 +532,6 @@ source->moments->KronFlux = NAN; source->moments->KronFluxErr = NAN; + source->moments->KronFinner = NAN; + source->moments->KronFouter = NAN; return false; }