Index: branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c (revision 36153) +++ branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c (revision 36154) @@ -72,4 +72,7 @@ assert (modelType > -1); + // how to deal with this? + int galaxyModelType = pmModelClassGetType ("PS_MODEL_EXP"); + // 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 @@ -87,9 +90,9 @@ } - // here are the things we read from the FITS table - - unsigned int ID = psMetadataLookupU32 (&status, row, "IPP_IDET"); - float X = psMetadataLookupF32 (&status, row, "X_PSF"); - float Y = psMetadataLookupF32 (&status, row, "Y_PSF"); + // here are the things we read from the FITS table (XXX modify names if needed) + + unsigned int ID = psMetadataLookupU32 (&status, row, "IPP_IDET"); + float X = psMetadataLookupF32 (&status, row, "X_PSF"); + float Y = psMetadataLookupF32 (&status, row, "Y_PSF"); float apRadius = psMetadataLookupS32 (&status, row, "AP_RADIUS"); @@ -97,4 +100,5 @@ float petRadius = psMetadataLookupU32 (&status, row, "PETRO_RADIUS"); bool fitGalaxy = psMetadataLookupU32 (&status, row, "FIT_GALAXY"); + bool psfStar = psMetadataLookupU32 (&status, row, "PSF_STAR"); float Rmajor = psMetadataLookupU32 (&status, row, "R_MAJOR"); @@ -108,11 +112,16 @@ source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS + 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; + source->seq = ID; + PAR[PM_PAR_XPOS] = X; + PAR[PM_PAR_YPOS] = Y; dPAR[PM_PAR_XPOS] = 0.0; @@ -129,7 +138,11 @@ source->psfMag = 0.0; - source->psfMagErr = 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; @@ -138,14 +151,46 @@ float peakFlux = 1.0; - - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + source->peak = pmPeakAlloc(X, Y, peakFlux, PM_PEAK_LONE); + source->peak->xf = X; // more accurate position + source->peak->yf = Y; // more accurate position + source->peak->dx = 0.0; + source->peak->dy = 0.0; source->peak->rawFlux = peakFlux; source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - assert (status); + + source->moments = pmMomentsAlloc (); + source->moments->Mx = X; + source->moments->My = Y; + source->moments->Mrf = kronRadius; + source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + // XXX does the above cause a problem with pmSourceMoments (do we run it?) + + if (isfinite(petRadius)) { + source->extpars = pmSourceExtendedParsAlloc (); + source->extpars->petrosianRadius = petRadius; + } + + if (fitGalaxy) { + source->modelFits = psArrayAllocEmpty (models->list->n); + pmModel *model = pmModelAlloc(modelType); + 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; + + psEllipseAxes galAxes; + galAxes.major = Rmajor; + galAxes.minor = Rminor; + galAxes.theta = theta; // XXX degrees or radians? + + pmPSF_AxesToModel (PAR, guessAxes, galModelType); + xPAR[PM_PAR_7] = Sindex; + + psArrayAdd (source->modelFits, 1, model); + psFree (model); + + } sources->data[i] = source; Index: branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36153) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c (revision 36154) @@ -210,8 +210,14 @@ if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; - // XXX where are we storing the supplied kron radius? + // psphotSetRadiusMomentsExact sets the radius based on Mrf if (!isfinite(source->moments->Mrf)) continue; + // modelFits is allocated if a galaxy fit is requested + if (!source->modelFits) continue; + + psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs + // skip saturated stars modeled with a radial profile + // XXX worry about this at some point.. if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; @@ -221,11 +227,4 @@ } - // XXX this is not right (use same concept as general source fit) - float windowRadius = 6.0*source->moments->Mrf ; - - // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); - psAssert (source->pixels, "WTF?"); - // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr // do the following for a set of shapes (Ex,Ey) @@ -245,15 +244,18 @@ pmModelType modelType = pmModelClassGetType ("PS_MODEL_EXP"); - pmModel *model = pmModelAlloc(modelType); - if (!model) { - return NULL; - } - + pmModel *model = source->modelFits->data[0]; + if (!model) return false; + + // we are using fitOptions->mode : be sure this makes sense pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + // 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) - psEllipsePol guessPol; + psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); + psEllipsePol guessPol = psEllipseAxesToPol (guessAxes); // I am going to retain the value of e0, and grid search around e1,e2 Index: branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c (revision 36153) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c (revision 36154) @@ -105,5 +105,5 @@ { pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF"); - if (!readoutCMF) goto loadTXT; + if (!readoutCMF) goto loadCFF; extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS"); @@ -111,4 +111,31 @@ for (int i = 0; i < extCMF->allSources->n; i++) { pmSource *source = extCMF->allSources->data[i]; + source->mode |= PM_SOURCE_MODE_EXTERNAL; + + // the supplied peak flux needs to be re-normalized + source->peak->rawFlux = 1.0; + source->peak->smoothFlux = 1.0; + source->peak->detValue = 1.0; + + // drop the loaded source modelPSF + psFree (source->modelPSF); + source->modelPSF = NULL; + source->imageID = index; + + psArrayAdd (detections->newSources, 100, source); + } + } + } + +loadCFF: + // load data from input CFF file: + { + pmReadout *readoutCFF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CFF"); + if (!readoutCFF) goto loadTXT; + + extCFF = psMetadataLookupPtr (NULL, readoutCFF->analysis, "PSPHOT.DETECTIONS"); + if (extCFF) { + for (int i = 0; i < extCFF->allSources->n; i++) { + pmSource *source = extCFF->allSources->data[i]; source->mode |= PM_SOURCE_MODE_EXTERNAL; Index: branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c =================================================================== --- branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c (revision 36153) +++ branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c (revision 36154) @@ -187,4 +187,7 @@ PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->variance, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->extpars, false); // the peak position is less accurate but less subject to extreme deviations @@ -202,5 +205,5 @@ // Calculate the Petro magnitude (make this block optional?) // XXX set the aperture here - float radPetro = 2.0*source->moments->Mrf; + float radPetro = PETROSIAN_RADII*source->extpars->petrosianRadius; float radPetro2 = radPetro*radPetro; @@ -242,6 +245,7 @@ // return the flux and error to parameters? - // source->moments->PetroFlux = Sum; - // source->moments->PetroFluxErr = sqrt(Var); + source->extpars->petrosianFlux = Sum; + source->extpars->petrosianFluxErr = sqrt(Var); + source->extpars->petrosianFill = nPetroPix / (M_PI * radPetro2); fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var));