Index: trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- trunk/psModules/src/camera/pmFPAfile.c (revision 36287) +++ trunk/psModules/src/camera/pmFPAfile.c (revision 36375) @@ -491,4 +491,7 @@ if (!strcasecmp(type, "CMF")) { return PM_FPA_FILE_CMF; + } + if (!strcasecmp(type, "CFF")) { + return PM_FPA_FILE_CFF; } if (!strcasecmp(type, "WCS")) { Index: trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- trunk/psModules/src/camera/pmFPAfile.h (revision 36287) +++ trunk/psModules/src/camera/pmFPAfile.h (revision 36375) @@ -34,4 +34,5 @@ PM_FPA_FILE_CMP, PM_FPA_FILE_CMF, + PM_FPA_FILE_CFF, PM_FPA_FILE_WCS, PM_FPA_FILE_RAW, Index: trunk/psModules/src/camera/pmFPAfileDefine.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileDefine.c (revision 36287) +++ trunk/psModules/src/camera/pmFPAfileDefine.c (revision 36375) @@ -103,7 +103,13 @@ type = psMetadataLookupStr(&status, data, "FILE.TYPE"); + if (!type) { + psError(PM_ERR_CONFIG, true, "FILE.TYPE is not defined for %s\n", name); + psFree(file); + return NULL; + } + file->type = pmFPAfileTypeFromString(type); if (file->type == PM_FPA_FILE_NONE) { - psError(PM_ERR_CONFIG, true, "FILE.TYPE is not defined for %s\n", name); + psError(PM_ERR_CONFIG, true, "FILE.TYPE %s is not registered in pmFPAfile.c:pmFPAfileTypeFromString\n", type); psFree(file); return NULL; Index: trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileIO.c (revision 36287) +++ trunk/psModules/src/camera/pmFPAfileIO.c (revision 36375) @@ -222,4 +222,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_WCS: case PM_FPA_FILE_SRCTEXT: @@ -317,4 +318,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_WCS: case PM_FPA_FILE_PSF: @@ -483,4 +485,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: status = pmFPAviewWriteObjects (view, file, config); break; @@ -567,4 +570,5 @@ case PM_FPA_FILE_PATTERN: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_WCS: case PM_FPA_FILE_PSF: @@ -643,4 +647,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_WCS: case PM_FPA_FILE_PSF: @@ -804,4 +809,5 @@ case PM_FPA_FILE_PATTERN: case PM_FPA_FILE_CMF: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_WCS: case PM_FPA_FILE_PSF: @@ -1016,4 +1022,5 @@ case PM_FPA_FILE_CMP: case PM_FPA_FILE_WCS: + case PM_FPA_FILE_CFF: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: Index: trunk/psModules/src/config/pmConfig.c =================================================================== --- trunk/psModules/src/config/pmConfig.c (revision 36287) +++ trunk/psModules/src/config/pmConfig.c (revision 36375) @@ -1860,5 +1860,5 @@ } - return psMetadataLookupMetadata(NULL, filerules, realname); + return psMetadataLookupMetadata(&mdok, filerules, realname); } Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 36287) +++ trunk/psModules/src/objects/Makefile.am (revision 36375) @@ -40,4 +40,5 @@ pmSourceIO_SX.c \ pmSourceIO_CMP.c \ + pmSourceIO_CFF.c \ pmSourceIO_SMPDATA.c \ pmSourceIO_PS1_DEV_0.c \ Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 36287) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 36375) @@ -269,5 +269,6 @@ // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width - float scale = moments->Mrf / axes.major; + // the factor of 2.3 comes from Table 1 of Graham and Driver (2005) + float scale = moments->Mrf / axes.major / 2.3; axes.major *= scale; axes.minor *= scale; Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 36287) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 36375) @@ -63,4 +63,17 @@ // 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords // values need to be pixel coords +// + +// Notes on changing kappa value from 1.70056 to 1.678 +// I'm using a functional form f(x,y) = Io exp(-kappa (r / r_e)). +// The article by Graham & Driver (2005) uses a form Ie exp(-bn [(r / r_e) -1]) +// which is equal to Ie exp(-bn (r / r_e)) exp(bn). +// Thus, my Io = Ie exp(bn) and my kappa is their bn. +// My value of kappa is 1.700, their value for bn is 1.678., so I am off by a small amount there (1.5%). + + +#define KAPPA_EXP 1.678 +#define OLD_KAPP_EXP 1.70056 + // Lax parameter limits @@ -109,5 +122,5 @@ // for EXP, we can hard-wire kappa(1): // float index = 1.0; - float kappa = 1.70056; + float kappa = KAPPA_EXP; // sqrt(z) is r @@ -318,5 +331,5 @@ // static value for EXP: - float kappa = 1.70056; + float kappa = KAPPA_EXP; // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 36287) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 36375) @@ -75,5 +75,5 @@ // Lax parameter limits -static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 }; +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.0625 }; static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 }; Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 36287) +++ trunk/psModules/src/objects/pmModel.c (revision 36375) @@ -217,9 +217,11 @@ // the options allow us to modify various aspects of the model if (mode & PM_MODEL_OP_NORM) { + // if we are including the sky, renormalizing should force use to normalized down the sky flux + params->data.F32[PM_PAR_SKY] /= params->data.F32[PM_PAR_I0]; params->data.F32[PM_PAR_I0] = 1.0; } if (!(mode & PM_MODEL_OP_SKY)) { params->data.F32[PM_PAR_SKY] = 0.0; - } + } if (mode & PM_MODEL_OP_CENTER) { params->data.F32[PM_PAR_XPOS] = image->col0 + 0.5*image->numCols; Index: trunk/psModules/src/objects/pmModelFuncs.h =================================================================== --- trunk/psModules/src/objects/pmModelFuncs.h (revision 36287) +++ trunk/psModules/src/objects/pmModelFuncs.h (revision 36375) @@ -36,11 +36,16 @@ typedef enum { - PM_MODEL_STATUS_NONE = 0x00, ///< model fit not yet attempted, no other info - PM_MODEL_STATUS_FITTED = 0x01, ///< model fit completed - PM_MODEL_STATUS_NONCONVERGE = 0x02, ///< model fit did not converge - PM_MODEL_STATUS_OFFIMAGE = 0x04, ///< model fit drove out of range - PM_MODEL_STATUS_BADARGS = 0x08, ///< model fit called with invalid args - PM_MODEL_STATUS_LIMITS = 0x10, ///< model parameters hit limits - PM_MODEL_STATUS_WEAK_FIT = 0x20, ///< model fit met loose tolerance, but not tight tolerance + PM_MODEL_STATUS_NONE = 0x000, ///< model fit not yet attempted, no other info + PM_MODEL_STATUS_FITTED = 0x001, ///< model fit completed + PM_MODEL_STATUS_NONCONVERGE = 0x002, ///< model fit did not converge + PM_MODEL_STATUS_OFFIMAGE = 0x004, ///< model fit drove out of range + PM_MODEL_STATUS_BADARGS = 0x008, ///< model fit called with invalid args + PM_MODEL_STATUS_LIMITS = 0x010, ///< model parameters hit limits + PM_MODEL_STATUS_WEAK_FIT = 0x020, ///< model fit met loose tolerance, but not tight tolerance + PM_MODEL_STATUS_NAN_CHISQ = 0x040, ///< model fit failed with a NAN chisq + PM_MODEL_SERSIC_PCM_FAIL_GUESS = 0x080, ///< sersic model fit failed on the initial moments-based guess + PM_MODEL_SERSIC_PCM_FAIL_GRID = 0x100, ///< sersic model fit failed on the grid search + PM_MODEL_PCM_FAIL_GUESS = 0x200, ///< non-sersic model fit failed on the initial moments-based guess + PM_MODEL_BEST_FIT = 0x400, ///< this model was the best fit and was subtracted } pmModelStatus; Index: trunk/psModules/src/objects/pmModelUtils.c =================================================================== --- trunk/psModules/src/objects/pmModelUtils.c (revision 36287) +++ trunk/psModules/src/objects/pmModelUtils.c (revision 36375) @@ -129,4 +129,8 @@ bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff) { + // restrict axex to 0.5 here not below + if (axes.minor < 0.2) axes.minor = 0.2; + if (axes.major < 0.2) axes.major = 0.2; + psEllipseShape shape = psEllipseAxesToShape (axes); @@ -137,10 +141,14 @@ // set the shape parameters if (useReff) { - *Sxx = PS_MAX(0.5, shape.sx); - *Syy = PS_MAX(0.5, shape.sy); + // *Sxx = PS_MAX(0.5, shape.sx); + // *Syy = PS_MAX(0.5, shape.sy); + *Sxx = shape.sx; + *Syy = shape.sy; *Sxy = shape.sxy * 2.0; } else { - *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); - *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + // *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); + // *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + *Sxx = M_SQRT2*shape.sx; + *Syy = M_SQRT2*shape.sy; *Sxy = shape.sxy; } @@ -193,4 +201,5 @@ if (!isfinite(axes.minor)) return false; if (!isfinite(axes.theta)) return false; + if (axes.major == 0) return false; // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width Index: trunk/psModules/src/objects/pmModel_CentralPixel.c =================================================================== --- trunk/psModules/src/objects/pmModel_CentralPixel.c (revision 36287) +++ trunk/psModules/src/objects/pmModel_CentralPixel.c (revision 36375) @@ -695,5 +695,6 @@ int npix = 0; - float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex); + // -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex); + float kappa = pmSersicKappa (Sindex); float rindex = 0.5 / Sindex; @@ -703,14 +704,21 @@ float delta = 1.0 / (float) Nsub; - float off = -Nsub2 * delta; - for (float ix = off; ix < 0.5; ix += delta) { - for (float iy = off; iy < 0.5; iy += delta) { - - float dX = dx + ix; - float dY = dy + iy; + // float off = -Nsub2 * delta; + + int Sx = (int) floor(dx / delta); + int Sy = (int) floor(dy / delta); + + for (int ix = -Nsub2; ix <= Nsub2; ix++) { + float dX = delta * (Sx + ix); + for (int iy = -Nsub2; iy <= Nsub2; iy++) { + float dY = delta * (Sy + iy); float z = PS_SQR(dX / Rxx) + PS_SQR(dY / Ryy) + dX * dY * Rxy; float q = pow (z, rindex); float f = exp(-kappa*q); + + // if ((ix == 0) && (iy == 0)) { + // // fprintf (stderr, "this: %f %f %f -- full : %f %f\n", z, q, f, flux, (float) npix); + // } flux += f; Index: trunk/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 36287) +++ trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 36375) @@ -135,4 +135,12 @@ } + if (min->isInteractive) { + fprintf (stderr, "%d : ", min->iter); + for (int ti = 0; ti < params->n; ti++) { + fprintf (stderr, "%f ", params->data.F32[ti]); + } + fprintf (stderr, " : %f\n", min->value); + } + char key[10]; // used for interactive responses bool testValue = false; @@ -140,5 +148,5 @@ // set a new guess for Alpha, Beta, Params if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) { - if (min->isInteractive) { + if (false && min->isInteractive) { fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] "); if (!fgets(key, 8, stdin)) { @@ -167,5 +175,5 @@ } - if (min->isInteractive) { + if (false && min->isInteractive) { p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: "); fprintf (stdout, "last chisq : %f\n", min->value); @@ -473,10 +481,15 @@ # else if (pcm->use1Dgauss) { - // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): - // * the model flux is not masked - // * threading takes place above this level - pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type); - psImageSmooth_PreAlloc_F32 (pcm->modelConvFlux, pcm->smdata); - // psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma); + + if (USE_1D_CACHE) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type); + psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata); + } else { + pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type); + psImageSmooth2dCache_F32 (pcm->modelConvFlux, pcm->smdata2d); + } } else { psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT); @@ -493,10 +506,14 @@ # else if (pcm->use1Dgauss) { - // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): - // * the model flux is not masked - // * threading takes place above this level - dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); - psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata); - // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma); + if (USE_1D_CACHE) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmoothCache_F32 (dmodelConv, pcm->smdata); + } else { + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d); + } } else { psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT); @@ -514,10 +531,14 @@ if (pcm->use1Dgauss) { - // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): - // * the model flux is not masked - // * threading takes place above this level - dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); - psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata); - // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma); + if (USE_1D_CACHE) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmoothCache_F32 (dmodelConv, pcm->smdata); + } else { + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d); + } } else { psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf); @@ -541,5 +562,7 @@ static int Npass = 0; char name[128]; - snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name); + if (!pcm->use1Dgauss) { + snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name); + } snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name); snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name); @@ -579,5 +602,7 @@ float ymodel = pcm->modelConvFlux->data.F32[i][j]; - float yweight = 1.0 / source->variance->data.F32[i][j]; + + // XXXX note this point here::: + float yweight = pcm->poissonErrors ? 1.0 / source->variance->data.F32[i][j] : 1.0; float delta = ymodel - source->pixels->data.F32[i][j]; Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 36287) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 36375) @@ -43,5 +43,4 @@ # define USE_DELTA_PSF 0 -# define USE_1D_GAUSS 1 static void pmPCMdataFree (pmPCMdata *pcm) { @@ -58,5 +57,7 @@ psFree (pcm->psfFFT); psFree (pcm->constraint); + psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc + psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc return; } @@ -88,4 +89,7 @@ } + pcm->smdata = NULL; + pcm->smdata2d = NULL; + pcm->modelConv = NULL; pcm->psf = NULL; @@ -94,8 +98,10 @@ pcm->nDOF = 0; + pcm->poissonErrors = true; + // full convolution with the PSF is expensive. if we have to save time, we can do a 1D // convolution with a Gaussian approximation to the kernel pcm->use1Dgauss = false; - pcm->nsigma = 3.0; + pcm->nsigma = NAN; // this is set to something defined by the user pcm->sigma = 1.0; // this should be set to something sensible when the psf is known @@ -248,21 +254,115 @@ } +static int modelType_GAUSS = -1; +static int modelType_PS1_V1 = -1; + +// generate a Gaussian smoothing kernel for supplied sigma. sigma here does not need to match +// that used to allocate the structure, but it is recommended +bool psImageSmoothCacheKernel_PS1_V1 (psImageSmoothCacheData *smdata, float sigma, float kappa) { + // check for NULL structure elements? + + int size = smdata->Nrange; + + psFree (smdata->kernel); + smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32); + + double sum = 0.0; // Sum of Gaussian, for normalization + double factor = 1.0 / (sigma * M_SQRT2); // Multiplier for i -> z + + // PS1_V1 is a power-law with fitted linear term: + // 1 / (1 + kappa z + z^1.666) where z = (r/sigma)^2 + + // generate the kernel (not normalized) + for (int i = -size, j = 0; i <= size; i++, j++) { + float z = PS_SQR(i * factor); + sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666)); + } + + // renormalize kernel to integral of 1.0 + for (int i = 0; i < 2 * size + 1; i++) { + smdata->kernel->data.F32[i] /= sum; + } + + return true; +} + +psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) { + + psAssert (modelPSF, "psf model must be defined"); + + psEllipseAxes axes; + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); + + *sigma = NAN; + *kappa = NAN; + + // XXX need to do this more carefully + if (modelPSF->type == modelType_GAUSS) { + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; + } + if (modelPSF->type == modelType_PS1_V1) { + *sigma = 0.5 * (axes.major + axes.minor); + *kappa = PAR[PM_PAR_7]; + } + psAssert (isfinite(*sigma), "invalid model type"); + + // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector + psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (flux, *sigma, nsigma); + + if (modelPSF->type == modelType_GAUSS) { + psImageSmoothCacheKernel_Gauss (smdata, *sigma); + } + if (modelPSF->type == modelType_PS1_V1) { + psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa); + } + + return smdata; +} + +psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) { + + psAssert (modelPSF, "psf model must be defined"); + + psEllipseAxes axes; + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); + + *sigma = NAN; + *kappa = NAN; + + // XXX need to do this more carefully + if (modelPSF->type == modelType_GAUSS) { + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; + } + if (modelPSF->type == modelType_PS1_V1) { + *sigma = 0.5 * (axes.major + axes.minor); + *kappa = PAR[PM_PAR_7]; + } + psAssert (isfinite(*sigma), "invalid model type"); + + // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma); + + if (modelPSF->type == modelType_GAUSS) { + psImageSmooth2dCacheKernel_Gauss (smdata, *sigma); + } + if (modelPSF->type == modelType_PS1_V1) { + psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa); + } + + return smdata; +} + pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) { - // make sure we save a cached copy of the psf flux - pmSourceCachePSF (source, maskVal); - - // convert the cached cached psf model for this source to a psKernel - psKernel *psf = pmPCMkernelFromPSF (source, psfSize); - if (!psf) { - // NOTE: this only happens if the source is too close to an edge - model->flags |= PM_MODEL_STATUS_BADARGS; - return NULL; - } - -# if (USE_DELTA_PSF) - psImageInit (psf->image, 0.0); - psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; -# endif + modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS"); + modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1"); // count the number of unmasked pixels: @@ -298,5 +398,4 @@ if (nPix < nParams + 1) { psTrace ("psModules.objects", 4, "insufficient valid pixels\n"); - psFree (psf); psFree (constraint); model->flags |= PM_MODEL_STATUS_BADARGS; @@ -306,8 +405,9 @@ // generate PCM data storage structure pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source); - - pcm->psf = psf; pcm->modelConv = psMemIncrRefCounter(model); pcm->constraint = constraint; + + pcm->poissonErrors = fitOptions->poissonErrors; + pcm->nsigma = fitOptions->nsigma; pcm->nPix = nPix; @@ -316,21 +416,29 @@ # if (USE_1D_GAUSS) - pmModel *modelPSF = source->modelPSF; - psAssert (modelPSF, "psf model must be defined"); - - psEllipseAxes axes; - bool useReff = pmModelUseReff (modelPSF->type); - psF32 *PAR = modelPSF->params->data.F32; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); - float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); pcm->use1Dgauss = true; - pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; - pcm->nsigma = 2.0; - - pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma); + if (USE_1D_CACHE) { + pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF); + } else { + pcm->smdata2d = psImageSmooth2dCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF); + } + # else + // make sure we save a cached copy of the psf flux + pmSourceCachePSF (source, maskVal); + + // convert the cached cached psf model for this source to a psKernel + psKernel *psf = pmPCMkernelFromPSF (source, psfSize); + if (!psf) { + // NOTE: this only happens if the source is too close to an edge + model->flags |= PM_MODEL_STATUS_BADARGS; + return NULL; + } + +# if (USE_DELTA_PSF) + psImageInit (psf->image, 0.0); + psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; +# endif + pcm->psf = psf; pcm->smdata = NULL; # endif @@ -395,6 +503,29 @@ pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32); } - psFree(pcm->smdata); - pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma); + + // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain) + if (USE_1D_CACHE) { + psFree(pcm->smdata); + pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma); + + pmModel *modelPSF = source->modelPSF; + if (modelPSF->type == modelType_GAUSS) { + psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma); + } + if (modelPSF->type == modelType_PS1_V1) { + psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa); + } + } else { + psFree(pcm->smdata2d); + pcm->smdata2d = psImageSmooth2dCacheAlloc (pcm->nsigma); + + pmModel *modelPSF = source->modelPSF; + if (modelPSF->type == modelType_GAUSS) { + // psImageSmooth2dCacheKernel_Gauss (pcm->smdata2d, pcm->sigma); + } + if (modelPSF->type == modelType_PS1_V1) { + psImageSmooth2dCacheKernel_PS1_V1 (pcm->smdata2d, pcm->sigma, pcm->kappa); + } + } } @@ -403,5 +534,5 @@ // construct a realization of the source model -bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) { +bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -420,24 +551,18 @@ // convolve the model image with the PSF if (USE_1D_GAUSS) { - // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): - // * the model flux is not masked - // * threading takes place above this level - // define the Gauss parameters from the psf - pmModel *modelPSF = source->modelPSF; - psAssert (modelPSF, "psf model must be defined"); - - psEllipseAxes axes; - bool useReff = pmModelUseReff (modelPSF->type); - psF32 *PAR = modelPSF->params->data.F32; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); - float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); - - float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; - float nsigma = 2.0; - - psImageSmooth (source->modelFlux, sigma, nsigma); + float sigma = NAN; + float kappa = NAN; + + if (USE_1D_CACHE) { + psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF); + psImageSmoothCache_F32 (source->modelFlux, smdata); + psFree (smdata); + } else { + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF); + psImageSmooth2dCache_F32 (source->modelFlux, smdata); + psFree (smdata); + } + // old call: psImageSmooth (source->modelFlux, sigma, nsigma); } else { // make sure we save a cached copy of the psf flux @@ -459,5 +584,5 @@ // construct a realization of the source model -bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) { +bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -468,28 +593,23 @@ // modelFlux always has unity normalization (I0 = 1.0) - pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | PM_MODEL_OP_NORM, maskVal); // convolve the model image with the PSF if (USE_1D_GAUSS) { - // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): - // * the model flux is not masked - // * threading takes place above this level - - // define the Gauss parameters from the psf - pmModel *modelPSF = source->modelPSF; - psAssert (modelPSF, "psf model must be defined"); - - psEllipseAxes axes; - bool useReff = pmModelUseReff (modelPSF->type); - psF32 *PAR = modelPSF->params->data.F32; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); - float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); - - float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; - float nsigma = 2.0; - - psImageSmooth (source->modelFlux, sigma, nsigma); + + float sigma = NAN; + float kappa = NAN; + + if (USE_1D_CACHE) { + psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF); + psImageSmoothCache_F32 (source->modelFlux, smdata); + psFree (smdata); + } else { + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF); + psImageSmooth2dCache_F32 (source->modelFlux, smdata); + psFree (smdata); + } + // old call: psImageSmooth (source->modelFlux, sigma, nsigma); } else { // make sure we save a cached copy of the psf flux @@ -509,3 +629,2 @@ return true; } - Index: trunk/psModules/src/objects/pmPCMdata.h =================================================================== --- trunk/psModules/src/objects/pmPCMdata.h (revision 36287) +++ trunk/psModules/src/objects/pmPCMdata.h (revision 36375) @@ -14,4 +14,8 @@ /// @addtogroup Objects Object Detection / Analysis Functions /// @{ + +// XXX this is basically for testing -- when I am happy with the convolution process, I'll strip this out +# define USE_1D_CACHE 0 +# define USE_1D_GAUSS 1 /** pmPCMdata : PSF Convolved Model data storage structure @@ -36,9 +40,14 @@ int nDOF; + bool poissonErrors; + bool use1Dgauss; + float kappa; float sigma; float nsigma; - psImageSmooth_PreAlloc_Data *smdata; + // psArray *smdata; + psImageSmoothCacheData *smdata; + psImageSmooth2dCacheData *smdata2d; } pmPCMdata; @@ -96,7 +105,7 @@ bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); -bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize); +bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma); -bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize); +bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize); /// @} Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 36287) +++ trunk/psModules/src/objects/pmSource.c (revision 36375) @@ -66,4 +66,5 @@ psFree(tmp->extpars); psFree(tmp->diffStats); + psFree(tmp->galaxyFits); psFree(tmp->radialAper); psTrace("psModules.objects", 10, "---- end ----\n"); @@ -164,4 +165,5 @@ source->extpars = NULL; source->diffStats = NULL; + source->galaxyFits = NULL; source->radialAper = NULL; source->parent = NULL; @@ -690,5 +692,4 @@ // why do we recalculate moments here? // we already attempt to do this in psphotSourceStats - // pmSourceMoments (source, INNER_RADIUS); Nsatstar ++; continue; @@ -804,191 +805,4 @@ return true; } - -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -pmSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - pmSource - pmSource->peak - pmSource->pixels - pmSource->variance - pmSource->mask - -XXX: The peak calculations are done in image coords, not subImage coords. - -XXX EAM : this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons -*****************************************************************************/ -# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -/*** this been moved to pmSourceMoments.c ***/ -# if (0) -bool pmSourceMoments(pmSource *source, - psF32 radius) -{ - psTrace("psModules.objects", 10, "---- begin ----\n"); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - psF32 xPeak = source->peak->x; - psF32 yPeak = source->peak->y; - psF32 xOff = source->pixels->col0 - source->peak->x; - psF32 yOff = source->pixels->row0 - source->peak->y; - - // XXX why do I get different results for these two methods of finding Sx? - // XXX Sx, Sy would be better measured if we clip pixels close to sky - // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? - // We loop through all pixels in this subimage (source->pixels), and for each - // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk) { - vMsk++; - psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to mask: %d\n", - col, row, (int)*vMsk); - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - psF32 xDiff = col + xOff; - psF32 yDiff = row + yOff; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) { -#if 1 - psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to position: %f %f\n", - col, row, xDiff, yDiff); -#endif - continue; - } - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : check for valid S/N in pixel - // XXX EAM : should this limit be user-defined? -#if 1 - if (PS_SQR(pDiff) < wDiff) { - psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to insignificance: %f, %f\n", - col, row, pDiff, wDiff); - continue; - } -#endif - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - XY += xDiff * yWght; - X2 += xDiff * xWght; - Y2 += yDiff * yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", - numPixels, (int)(0.75*R2), Sum); - psTrace("psModules.objects", 10, "---- end (false) ----\n"); - return (false); - } - - psTrace ("psModules.objects", 4, "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - psTrace("psModules.objects", 10, "---- end(false) ----\n"); - return (false); - } - - source->moments->Mx = x + xPeak; - source->moments->My = y + yPeak; - - // XXX EAM : Sxy needs to have x*y subtracted - source->moments->Mxy = XY/Sum - x*y; - source->moments->Sum = Sum; - source->moments->SN = Sum / sqrt(Var); - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - - // XXX EAM : these values can be negative, so we need to limit the range - // XXX EAM : make the use of this consistent: should this be the second moment or sqrt? - // source->moments->Mxx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); - // source->moments->Myy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); - source->moments->Mxx = PS_MAX(X2/Sum - PS_SQR(x), 0); - source->moments->Myy = PS_MAX(Y2/Sum - PS_SQR(y), 0); - - psTrace ("psModules.objects", 4, - "sky: %f Sum: %f Mx: %f My: %f Mxx: %f Myy: %f Mxy: %f\n", - sky, Sum, source->moments->Mx, source->moments->My, - source->moments->Mxx, source->moments->Myy, source->moments->Mxy); - - psTrace("psModules.objects", 10, "---- end ----\n"); - return(true); -} -# endif // construct a realization of the source model Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 36287) +++ trunk/psModules/src/objects/pmSource.h (revision 36375) @@ -40,4 +40,6 @@ PM_SOURCE_TMPF_PETRO_KEEP = 0x0100, PM_SOURCE_TMPF_PETRO_SKIP = 0x0200, + PM_SOURCE_TMPF_EXT_FIT = 0x0400, // not just galaxies (trails as well) + PM_SOURCE_TMPF_PETRO = 0x0800, } pmSourceTmpF; @@ -117,4 +119,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections + pmSourceGalaxyFits *galaxyFits; ///< fits to galaxy models (psphotFullForce only) psArray *radialAper; ///< radial flux in circular apertures pmSource *parent; ///< reference to the master source from which this is derived Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 36375) @@ -286,2 +286,25 @@ return pars; } + +// *** pmSourceExtFitPars describes extra metadata related to an extended fit +static void pmSourceGalaxyFitsFree (pmSourceGalaxyFits *tmp) { + + psFree (tmp->Flux); + psFree (tmp->dFlux); + psFree (tmp->chisq); + + return; +} + +pmSourceGalaxyFits *pmSourceGalaxyFitsAlloc (void) { + + pmSourceGalaxyFits *tmp = (pmSourceGalaxyFits *) psAlloc(sizeof(pmSourceGalaxyFits)); + psMemSetDeallocator(tmp, (psFreeFunc) pmSourceGalaxyFitsFree); + + tmp->Flux = psVectorAllocEmpty (25, PS_TYPE_F32); + tmp->dFlux = psVectorAllocEmpty (25, PS_TYPE_F32); + tmp->chisq = psVectorAllocEmpty (25, PS_TYPE_F32); + tmp->nPix = 0; + + return tmp; +} Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 36287) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 36375) @@ -82,4 +82,11 @@ } pmSourceExtFitPars; +typedef struct { + psVector *Flux; + psVector *dFlux; + psVector *chisq; + int nPix; +} pmSourceGalaxyFits; + pmSourceRadialFlux *pmSourceRadialFluxAlloc(); bool psMemCheckSourceRadialFlux(psPtr ptr); @@ -109,4 +116,7 @@ pmSourceExtFitPars *pmSourceExtFitParsAlloc (void); +pmSourceGalaxyFits *pmSourceGalaxyFitsAlloc (void); + + /// @} # endif /* PM_SOURCE_H */ Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 36375) @@ -59,4 +59,5 @@ opt->maxTol = 1.00; opt->weight = 1.00; + opt->nsigma = 5.00; opt->maxChisqDOF = NAN; opt->poissonErrors = true; @@ -281,5 +282,11 @@ // set the model success or failure status model->flags |= PM_MODEL_STATUS_FITTED; - if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; + if (!fitStatus) { + if (isnan(myMin->value)) { + model->flags |= PM_MODEL_STATUS_NAN_CHISQ; + } else { + model->flags |= PM_MODEL_STATUS_NONCONVERGE; + } + } if (myMin->chisqConvergence) { Index: trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.h (revision 36287) +++ trunk/psModules/src/objects/pmSourceFitModel.h (revision 36375) @@ -34,4 +34,5 @@ float weight; ///< use this weight for constant-weight fits float covarFactor; ///< covariance factor for calculating the chisq + float nsigma; ///< how far out to convolve bool poissonErrors; ///< use poisson errors for fits? bool saveCovariance; Index: trunk/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitPCM.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceFitPCM.c (revision 36375) @@ -51,4 +51,6 @@ # define TIMING 0 +bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal); + bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { @@ -113,5 +115,6 @@ } else { // xxx this is wrong because it does not convolve with the psf - pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal); + pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize); + pmSourceChisqModelFlux (source, pcm->modelConv, maskVal); } if (TIMING) { t4 = psTimerMark ("pmSourceFitPCM"); } @@ -119,5 +122,12 @@ // set the model success or failure status pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED; - if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; + + if (!fitStatus) { + if (isnan(myMin->value)) { + pcm->modelConv->flags |= PM_MODEL_STATUS_NAN_CHISQ; + } else { + pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; + } + } if (myMin->chisqConvergence) { Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 36375) @@ -352,5 +352,11 @@ // set the model success or failure status model->flags |= PM_MODEL_STATUS_FITTED; - if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; + if (!fitStatus) { + if (isnan(myMin->value)) { + model->flags |= PM_MODEL_STATUS_NAN_CHISQ; + } else { + model->flags |= PM_MODEL_STATUS_NONCONVERGE; + } + } // models can go insane: reject these Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 36375) @@ -69,4 +69,5 @@ psString *xfitname, // Extension name for extended fitted measurements psString *xradname, // Extension name for radial apertures + psString *xgalname, // Extension name for galaxy shapes const pmFPAfile *file, // File of interest const pmFPAview *view // View to level of interest @@ -140,4 +141,14 @@ } *xradname = pmFPAfileNameFromRule (rule, file, view); + } + + // EXTNAME for radial apertures + if (xgalname) { + const char *rule = psMetadataLookupStr(&status, menu, "CMF.XGAL"); + if (!rule) { + psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XGAL in EXTNAME.RULES in camera.config"); + return false; + } + *xgalname = pmFPAfileNameFromRule (rule, file, view); } @@ -362,4 +373,7 @@ status &= pmSourcesWrite_##TYPE##_XRAD (file->fits, readout, sources, file->header, xradname, recipe); \ } \ + if (xgalname) { \ + status &= pmSourcesWrite_##TYPE##_XGAL (file->fits, sources, xgalname, recipe); \ + } \ } @@ -464,8 +478,11 @@ } - // if this is not TRUE, the output files only contain the psf measurements. - bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + // if none of these are TRUE, the output files only contain the psf measurements. + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool XSRC_OUTPUT = doPetrosian || doAnnuli; bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES"); + bool XGAL_OUTPUT = psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES"); // define the EXTNAME values for the different data segments: @@ -476,8 +493,10 @@ psString xfitname = NULL; psString xradname = NULL; + psString xgalname = NULL; if (!pmSourceIOextnames(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL, XFIT_OUTPUT ? &xfitname : NULL, XRAD_OUTPUT ? &xradname : NULL, + XGAL_OUTPUT ? &xgalname : NULL, file, view)) { return false; @@ -563,4 +582,7 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "XRADNAME", PS_META_REPLACE, "name of XRAD table extension", xradname); } + if (xgalname) { + psMetadataAddStr (outhead, PS_LIST_TAIL, "XGALNAME", PS_META_REPLACE, "name of XGAL table extension", xgalname); + } // these are case-sensitive since the EXTYPE is case-sensitive @@ -609,4 +631,5 @@ psFree (xfitname); psFree (xradname); + psFree (xgalname); psFree (deteffname); @@ -620,7 +643,28 @@ psFree (xfitname); psFree (xradname); + psFree (xgalname); psFree (deteffname); return false; + 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; + } + + hdu = pmFPAviewThisHDU (view, fpa); + pmConfigConformHeader(hdu->header, file->format); + psFitsWriteBlank (file->fits, hdu->header, NULL); + file->header = hdu->header; + file->wrote_phu = true; + if (!pmSourcesWrite_CFF(readout, file->fits, sources, hdu->header, recipe)) { + psError(PS_ERR_UNKNOWN, false, "failed to write CFF"); + return false; + } + break; + } + default: fprintf (stderr, "warning: type mismatch\n"); @@ -914,4 +958,16 @@ psArray *sources = NULL; pmHDU *hdu; + + // define the EXTNAME values for the different data segments: + psString headname = NULL; + psString dataname = NULL; + psString deteffname = NULL; + psString xsrcname = NULL; + psString xfitname = NULL; + psString xradname = NULL; + psString xgalname = NULL; + + psMetadata *tableHeader = NULL; + char *xtension = NULL; switch (file->type) { @@ -963,12 +1019,4 @@ hdu = pmFPAviewThisHDU (view, file->fpa); - // define the EXTNAME values for the different data segments: - psString headname = NULL; - psString dataname = NULL; - psString deteffname = NULL; - psString xsrcname = NULL; - psString xfitname = NULL; - psString xradname = NULL; - // determine the output table format. Assume if we need to output extendend source // parameters that they may exist in the input. @@ -981,8 +1029,12 @@ } - // if this is not TRUE, the output files only contain the psf measurements. - bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + // if none of these are TRUE, we only read the psf measurements + // XXX: shouldn't we look for these extensions and read the regardless of the recipe values? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool XSRC_OUTPUT = doPetrosian || doAnnuli; bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES"); + bool XGAL_OUTPUT = false; // psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES"); if (!pmSourceIOextnames(&headname, &dataname, &deteffname, @@ -990,4 +1042,5 @@ XFIT_OUTPUT ? &xfitname : NULL, XRAD_OUTPUT ? &xradname : NULL, + XGAL_OUTPUT ? &xgalname : NULL, file, view)) { return false; @@ -1033,8 +1086,8 @@ } - psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header if (!tableHeader) psAbort("cannot read table header"); - char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); if (!xtension) psAbort("cannot read table type"); if (strcmp (xtension, "BINTABLE")) { @@ -1132,4 +1185,54 @@ break; + case PM_FPA_FILE_CFF: + // read in header, if not yet loaded + hdu = pmFPAviewThisHDU (view, file->fpa); + + // look these up in the camera config? + // headrule = {CHIP.NAME}.hdr + // datarule = {CHIP.NAME}.cff + + // define the EXTNAME values for the different data segments: + headname = pmFPAfileNameFromRule("{CHIP.NAME}.hdr", file, view); + dataname = pmFPAfileNameFromRule("{CHIP.NAME}.cff", file, view); + + // advance to the IMAGE HEADER extension + if (hdu->header == NULL) { + // if the IMAGE header does not exist, we have no data for this view + if (!psFitsMoveExtNameClean (file->fits, headname)) { + readout->data_exists = false; + psFree (headname); + psFree (dataname); + return true; + } + hdu->header = psFitsReadHeader (NULL, file->fits); + } + + // advance to the table data extension + // since we have read the IMAGE header, the TABLE header should exist + if (!psFitsMoveExtName (file->fits, dataname)) { + psAbort("cannot find data extension %s in %s", dataname, file->filename); + } + + tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + // verify this is a binary table + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", dataname); + psFree(tableHeader); + return false; + } + + sources = pmSourcesRead_CFF(file->fits, hdu->header); + + psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname); + psFree (headname); + psFree (dataname); + psFree (tableHeader); + break; + default: fprintf (stderr, "warning: type mismatch\n"); Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 36375) @@ -21,4 +21,5 @@ bool pmSourcesWrite_##TYPE##_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); \ bool pmSourcesWrite_##TYPE##_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); \ + bool pmSourcesWrite_##TYPE##_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe); \ psArray *pmSourcesRead_##TYPE (psFits *fits, psMetadata *header); \ bool pmSourcesRead_##TYPE##_XSRC (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \ @@ -52,4 +53,6 @@ psArray *pmSourcesReadCMP (char *filename, psMetadata *header); +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header); +bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe); bool pmSourcesWritePSFs (psArray *sources, char *filename); Index: trunk/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36375) +++ trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36375) @@ -0,0 +1,343 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSourceSatstar.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" + +#include "pmSourceIO.h" +#include "pmSourceOutputs.h" + +// read in a readout from the fits file +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + PS_ASSERT_INT_NONNEGATIVE(modelType, NULL); + + // 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++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); + psFree(sources); + return NULL; + } + + // here are the things we read from the FITS table (XXX modify names if needed) + + unsigned int ID = psMetadataLookupU32 (&status, row, "ID"); + float X = psMetadataLookupF32 (&status, row, "X"); + assert(status); // poorly constructed cff + float Y = psMetadataLookupF32 (&status, row, "Y"); + assert(status); // poorly constructed cff + float flux = psMetadataLookupF32 (&status, row, "FLUX"); // guess at the instrumental flux + if (!isfinite(X) || !isfinite(Y) || !isfinite(flux)) { + psError(PS_ERR_IO, true, "source ID %d is invalid x: %f y: %f flux: %f", ID, X, Y, flux); + psFree(sources); + return NULL; + } + + + float apRadius = psMetadataLookupF32 (&status, row, "AP_RADIUS"); + float kronRadius = psMetadataLookupF32 (&status, row, "KRON_RADIUS"); + float petRadius = psMetadataLookupF32 (&status, row, "PETRO_RADIUS"); + bool fitGalaxy = psMetadataLookupU8 (&status, row, "FIT_GALAXY"); + bool psfStar = psMetadataLookupU8 (&status, row, "PSF_STAR"); + + float Rmajor = psMetadataLookupF32 (&status, row, "R_MAJOR"); + float Rminor = psMetadataLookupF32 (&status, row, "R_MINOR"); + float theta = psMetadataLookupF32 (&status, row, "THETA"); + + // XXX: we need to put a lookup table in the cff header to define the correspondence of the + // model type values in the cff with our models. (We want to use an interger for efficiency + // but the value for each model type is set on the organization of the the array in pmModelClass.c + // For now use the input values verbatim and trust the user that this is valid value + int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); + float Sindex = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + 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; + + 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 + source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + // XXX does the above cause a problem with pmSourceMoments (do we run it? yes we do) + + if (isfinite(petRadius)) { + source->extpars = pmSourceExtendedParsAlloc (); + source->extpars->petrosianRadius = petRadius; + } + + // XXX: should use < 0 as invalid galaxyModelType + + 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; + xPAR[PM_PAR_YPOS] = Y; + + psEllipseAxes galaxyAxes; + galaxyAxes.major = Rmajor; + galaxyAxes.minor = Rminor; + galaxyAxes.theta = theta * PS_RAD_DEG; + + pmPSF_AxesToModel (xPAR, galaxyAxes, galaxyModelType); + if (model->params->n > 7) { + xPAR[PM_PAR_7] = 0.5 / Sindex; + } + + 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); + } + + return sources; +} + +bool pmSourcesWrite_CFF(pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe) { + + char *extname = "SkyChip.cff"; + + bool mdok; + psF32 exptime = psMetadataLookupF32(&mdok, readout->parent->concepts, "CELL.EXPOSURE"); + PS_ASSERT(mdok, false); + + psArray *table = psArrayAllocEmpty(sources->n); + + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + + for (int i = 0; i < sources->n; i++) { + pmSource *thisSource = sources->data[i]; + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + psF32 xPos, yPos, flux, rMajor, rMinor, theta; + psS32 modelType = 0; + bool fitGalaxy = false; + bool psfStar = false; + psF32 sersicIndex = 0; + if (source->modelFits == 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; + psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; + rMajor = 0; + rMinor = 0; + theta = 0; + } else { + // Choose the extended model to use for this source. + // For now we use the one set as the modelEXT for this source which had the best fit. + // XXX: make this controllable by recipe + // use best (as implemented here) + // choose specific type + int iModel = -1; + if (source->modelEXT) { + pmModelType ext_model_type = source->modelEXT->type; + for (int j=0; jmodelFits->n; j++) { + pmModel *aModel = source->modelFits->data[j]; + if (aModel->type == ext_model_type) { + iModel = j; + break; + } + } + } + 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); + } + + if (!psFitsWriteTable(fits, NULL, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + // psFree(header); + + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36375) @@ -789,4 +789,6 @@ char name[64]; + pmModelType modelTypeTrail = pmModelClassGetType("PS_MODEL_TRAIL"); + // create a header to hold the output data psMetadata *outhead = psMetadataAlloc (); @@ -798,9 +800,9 @@ sources = psArraySort (sources, pmSourceSortByFlux); - @>PS1_DV2@ float magOffset; - @>PS1_DV2@ float zeroptErr; - @>PS1_DV2@ float fwhmMajor; - @>PS1_DV2@ float fwhmMinor; - @>PS1_DV2@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams @@ -823,4 +825,12 @@ @>PS1_DV2@ pmChip *chip = readout->parent->parent; + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + table = psArrayAllocEmpty (sources->n); @@ -847,5 +857,6 @@ // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; + // XXX + if (model->flags & badModel) continue; PAR = model->params->data.F32; @@ -853,6 +864,21 @@ xPos = PAR[PM_PAR_XPOS]; yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; + + // for the extended source models, we do not always fit the centroid in the non-linear fitting process + // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, + // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) + // TRAIL : X,Y are fitted + // + + // XXX this should be based on what happened, not on the model type + if (model->type == modelTypeTrail) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // this is definitely an underestimate since it does not + // account for the extent of the source + xErr = fwhmMajor * model->magErr / 2.35; + yErr = fwhmMinor * model->magErr / 2.35; + } @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; @@ -870,5 +896,4 @@ row = psMetadataAlloc (); - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" @@ -888,6 +913,10 @@ @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Magnitude using supplied calibration", model->chisq); - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Magnitude using supplied calibration", model->nDOF); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); + @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE", PS_DATA_S32, "type for chosen EXT_MODEL", source->modelEXT ? source->modelEXT->type : -1); + + // 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_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); @@ -946,5 +975,5 @@ snprintf (name, 64, "EXT_PAR_%02d", k); - + if (k < model->params->n) { psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); @@ -956,5 +985,5 @@ // optionally, write out the covariance matrix values // XXX do I need to pad this to match the biggest covar matrix? - if (model->covar) { + if (false && model->covar) { for (int iy = 0; iy < model->covar->numCols; iy++) { for (int ix = iy; ix < model->covar->numCols; ix++) { @@ -1058,4 +1087,17 @@ model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ"); + model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF"); + + // EXT_MODEL_TYPE gives the model chosen by psphot as the best. + // Putting this into the XFIT table makes 3 copies of it (one for each model) + // but since we have many fewer XFIT rows than psf rows that is cheaper than putting it + // in the psf table. + psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE"); + if (!status) { + // older cmfs don't have this column + extModelType = -1; + } + psEllipseAxes axes; axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); @@ -1072,20 +1114,40 @@ if (model->params->n > 7) { PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); - } - // read the covariance matrix - int nparams = model->params->n; - psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); - for (int y = 0; y < nparams; y++) { - for (int x = 0; x < nparams; x++) { - char name[64]; - snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); - covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + // XXX add an error: + // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); + } + + // XXX : make this depend on what is in the cmf + if (0) { + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + } + + if (modelType == extModelType) { + // The software that created this source picked this model as the best of the fits. + // Set the extModel to point to it. + // This is important for programs like psastro (skycal) so that its output cmfs + // will have valid EXT_MODEL_TYPE + psFree(source->modelEXT); + source->modelEXT = psMemIncrRefCounter(model); + source->type = PM_SOURCE_TYPE_EXTENDED; + if (0) { + // since FLAGS were read we don't need to do this + source->mode |= PM_SOURCE_MODE_EXTMODEL; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; } } - model->covar = covar; psArrayAdd(source->modelFits, 1, model); psFree(model); - psFree(row); } @@ -1209,8 +1271,8 @@ write_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); psFree (radFlux); psFree (radFluxErr); @@ -1343,2 +1405,88 @@ return true; } + +// XXX where should I record the number of columns?? +bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + bool status = false; + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { + psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); + return true; + } + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + + psArray *table = psArrayAllocEmpty (sources->n); + + for (int i = 0; i < sources->n; i++) { + + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // if we did not fit the galaxy model, modelFits will be NULL + if (source->modelFits == NULL) continue; + + // if we did not fit the galaxy model, galaxyFits will also be NULL + 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); + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} + Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 36375) @@ -713,2 +713,7 @@ return true; } + +bool pmSourcesWrite_PS1_CAL_0_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 36375) @@ -255,2 +255,7 @@ return true; } + +bool pmSourcesWrite_PS1_DEV_0_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 36375) @@ -595,2 +595,7 @@ return true; } + +bool pmSourcesWrite_PS1_DEV_1_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 36375) @@ -225,2 +225,7 @@ return true; } + +bool pmSourcesWrite_SMPDATA_XGAL(psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + return true; +} Index: trunk/psModules/src/objects/pmSourceMasks.h =================================================================== --- trunk/psModules/src/objects/pmSourceMasks.h (revision 36287) +++ trunk/psModules/src/objects/pmSourceMasks.h (revision 36375) @@ -56,4 +56,19 @@ PM_SOURCE_MODE2_DIFF_SELF_MATCH = 0x00000800, ///< positive detection match is probably this source PM_SOURCE_MODE2_SATSTAR_PROFILE = 0x00001000, ///< saturated source is modeled with a radial profile + + PM_SOURCE_MODE2_ECONTOUR_FEW_PTS = 0x00002000, ///< too few points to measure the elliptical contour + PM_SOURCE_MODE2_RADBIN_NAN_CENTER = 0x00004000, ///< radial bins failed with too many NaN center bin + PM_SOURCE_MODE2_PETRO_NAN_CENTER = 0x00008000, ///< petrosian radial bins failed with too many NaN center bin + PM_SOURCE_MODE2_PETRO_NO_PROFILE = 0x00010000, ///< petrosian not build because radial bins missing + + PM_SOURCE_MODE2_PETRO_INSIG_RATIO = 0x00020000, ///< insignificant measurement of petrosian ratio + PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN = 0x00040000, ///< petrosian ratio in the 0th bin (likely bad) + + PM_SOURCE_MODE2_EXT_FITS_RUN = 0x00080000, ///< we attempted to run extended fits on this source + PM_SOURCE_MODE2_EXT_FITS_FAIL = 0x00100000, ///< at least one of the model fits failed + PM_SOURCE_MODE2_EXT_FITS_RETRY = 0x00200000, ///< one of the model fits was re-tried with new window + PM_SOURCE_MODE2_EXT_FITS_NONE = 0x00400000, ///< ALL of the model fits failed + + } pmSourceMode2; Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 36287) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 36375) @@ -65,6 +65,9 @@ void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; } +bool pmSourceMomentsHighOrder (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal); +bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal); +bool pmSourceMomentsKronFluxes (pmSource *source, float sigma, float minSN, psImageMaskType maskVal); + // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid) - bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal) { @@ -74,18 +77,7 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // this function assumes the sky has been well-subtracted for the image - float sky = 0.0; - if (source->moments == NULL) { source->moments = pmMomentsAlloc(); } - - float Sum = 0.0; - float Var = 0.0; - float SumCore = 0.0; - float VarCore = 0.0; - float R2 = PS_SQR(radius); - float minSN2 = PS_SQR(minSN); - float rsigma2 = 0.5 / PS_SQR(sigma); // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -110,10 +102,229 @@ // of any object drops pretty quickly outside 1-2 sigmas. (The exception is bright // saturated stars, for which we need to use a very large radius here) + // NOTE: if (mode & EXTERNAL) or (mode2 & MATCHED), do not re-calculate the centroid (use peak as centroid) + // (we still call this function because it sets moments->Sum,SN,Peak,nPixels if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) { return false; } + pmSourceMomentsHighOrder (source, radius, sigma, minSN, maskVal); + + // now calculate the 1st radial moment (for kron flux) using symmetrical averaging + pmSourceMomentsRadialMoment (source, radius, minKronRadius, maskVal); + + // now calculate the kron flux values using the 1st radial moment + pmSourceMomentsKronFluxes (source, sigma, minSN, maskVal); + + psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", + source->moments->Mrf, source->moments->KronFlux, + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, + source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); + + psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f Npix: %d\n", + source->peak->xf, source->peak->yf, + source->peak->rawFlux, sqrt(source->peak->detValue), + source->moments->Mx, source->moments->My, + source->moments->Sum, + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->nPixels); + + return(true); +} + +bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { + + // First Pass: calculate the first moments (these are subtracted from the coordinates below) + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // .. etc + + float sky = 0.0; + + float peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + float Sum = 0.0; + float Var = 0.0; + float X1 = 0.0; + float Y1 = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); + + float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage + float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage + + // we are guaranteed to have a valid pixel and variance at this location (right? right?) + // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); + + // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should + // not depend on the fractional pixel location of the source. However, the aperture + // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel + // position of the expected centroid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + float yDiff = row + 0.5 - yPeak; + if (fabs(yDiff) > radius) continue; + + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row]; + + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + float xDiff = col + 0.5 - xPeak; + if (fabs(xDiff) > radius) continue; + + // radius is just a function of (xDiff, yDiff) + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + float pDiff = *vPix - sky; + float wDiff = *vWgt; + + // skip pixels below specified significance level. for a PSFs, this + // over-weights the wings of bright stars compared to those of faint stars. + // for the estimator used for extended source analysis (where the window + // function is allowed to be arbitrarily large), we need to clip to avoid + // negative second moments. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; // + if ((minSN > 0.0) && (pDiff < 0)) continue; // + + // Apply a Gaussian window function. Be careful with the window function. S/N + // weighting over weights the sky for faint sources + if (sigma > 0.0) { + float z = r2*rsigma2; + assert (z >= 0.0); + float weight = exp(-z); + + wDiff *= weight; + pDiff *= weight; + } + + Var += wDiff; + Sum += pDiff; + + float xWght = xDiff * pDiff; + float yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + } + } + + // if we have less than (1/4) of the possible pixels (in circle or box), force a retry + int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0); + + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < minPixels) || (Sum <= 0)) { + psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum); + return (false); + } + + // calculate the first moment. + float Mx = X1/Sum; + float My = Y1/Sum; + if ((fabs(Mx) > radius) || (fabs(My) > radius)) { + psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); + return (false); + } + if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) { + psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y); + } + + psTrace ("psModules.objects", 5, "id: %d, sky: %f Mx: %f My: %f Sum: %f X1: %f Y1: %f Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels); + + // add back offset of peak in primary image + // also offset from pixel index to pixel coordinate + // (the calculation above uses pixel index instead of coordinate) + // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords + + // we only update the centroid if the position is not supplied from elsewhere + bool skipCentroid = false; + skipCentroid |= (source->mode & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions + skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions + + if (skipCentroid) { + source->moments->Mx = source->peak->xf; + source->moments->My = source->peak->yf; + } else { + source->moments->Mx = Mx + xGuess; + source->moments->My = My + yGuess; + } + + source->moments->Sum = Sum; + source->moments->SN = Sum / sqrt(Var); + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + return true; +} + +float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) { + + psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (int i = 0; i < sources->n; i++) { + pmSource *src = sources->data[i]; // Source of interest + if (!src || !src->moments) { + continue; + } + + if (src->mode & PM_SOURCE_MODE_BLEND) { + continue; + } + + if (!src->moments->nPixels) continue; + + if (src->moments->SN < PSF_SN_LIM) continue; + + // XXX put in Mxx,Myy cut based on clump location + + psVectorAppend(radii, src->moments->Mrf); + } + + // find the peak in this image + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + if (!psVectorStats (stats, radii, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); + psFree(stats); + return NAN; + } + + float minRadius = stats->sampleMedian; + + psFree(radii); + psFree(stats); + return minRadius; +} + +bool pmSourceMomentsHighOrder (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal) { + + // this function assumes the sky has been well-subtracted for the image + float Sum = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); + // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates - // Xn = SUM (x - xc)^n * (z - sky) + // Xn = SUM (x - xc)^n * (z - sky) -- note that sky is 0.0 by definition here float XX = 0.0; float XY = 0.0; @@ -129,11 +340,4 @@ float YYYY = 0.0; - Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum - - // float dX = source->moments->Mx - source->peak->xf; - // float dY = source->moments->My - source->peak->yf; - // float dR = hypot(dX, dY); - // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; - // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; float Xo = source->moments->Mx; float Yo = source->moments->My; @@ -154,5 +358,4 @@ psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -173,5 +376,5 @@ if (r2 > R2) continue; - float fDiff = *vPix - sky; + float fDiff = *vPix; float pDiff = fDiff; float wDiff = *vWgt; @@ -181,8 +384,8 @@ // stars. if (PS_SQR(pDiff) < minSN2*wDiff) continue; - if ((minSN > 0.0) && (pDiff < 0)) continue; // + if ((minSN > 0.0) && (pDiff < 0)) continue; // Apply a Gaussian window function. Be careful with the window function. S/N - // weighting over weights the sky for faint sources + // weighting over-weights the sky for faint sources if (sigma > 0.0) { float z = r2 * rsigma2; @@ -230,20 +433,4 @@ XYYY += xyyy; YYYY += yyyy; - - // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) - // XXX float r = sqrt(r2); - // XXX float rf = r * fDiff; - // XXX float rh = sqrt(r) * fDiff; - // XXX float rs = fDiff; - // XXX - // XXX float rfw = r * pDiff; - // XXX float rhw = sqrt(r) * pDiff; - // XXX - // XXX RF += rf; - // XXX RH += rh; - // XXX RS += rs; - // XXX - // XXX RFW += rfw; - // XXX RHW += rhw; } } @@ -263,8 +450,11 @@ source->moments->Myyyy = YYYY/Sum; - // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging + return true; +} + +bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) { + float **vPix = source->pixels->data.F32; - float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32; psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; @@ -272,4 +462,15 @@ float RH = 0.0; float RS = 0.0; + + // centroid around which to calculate the moments + float Xo = source->moments->Mx; + float Yo = source->moments->My; + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + float R2 = PS_SQR(radius); for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -304,6 +505,6 @@ if (r2 > R2) continue; - float fDiff1 = vPix[row][col] - sky; - float fDiff2 = vPix[yFlip][xFlip] - sky; + float fDiff1 = vPix[row][col]; + float fDiff2 = vPix[yFlip][xFlip]; float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2)); @@ -329,11 +530,22 @@ kronRefRadius = MIN(radius, kronRefRadius); } - source->moments->Mrf = kronRefRadius; - - // *** now calculate the kron flux values using the 1st radial moment - - float radKinner = 1.0*kronRefRadius; - float radKron = 2.5*kronRefRadius; - float radKouter = 4.0*kronRefRadius; + + // if source is externally supplied and it already has a finite Mrf do not change it + if (! ((source->mode & PM_SOURCE_MODE_EXTERNAL) && isfinite(source->moments->Mrf))) { + source->moments->Mrf = kronRefRadius; + } + + return true; +} + +bool pmSourceMomentsKronFluxes (pmSource *source, float sigma, float minSN, psImageMaskType maskVal) { + + float **vPix = source->pixels->data.F32; + float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + float radKinner = 1.0*source->moments->Mrf; + float radKron = 2.5*source->moments->Mrf; + float radKouter = 4.0*source->moments->Mrf; int nKronPix = 0; @@ -341,7 +553,22 @@ int nInner = 0; int nOuter = 0; - Sum = Var = 0.0; + + float Sum = 0.0; + float Var = 0.0; + float SumCore = 0.0; + float VarCore = 0.0; float SumInner = 0.0; float SumOuter = 0.0; + + // centroid around which to calculate the moments + float Xo = source->moments->Mx; + float Yo = source->moments->My; + + // center of mass in subimage. Note: the calculation below uses pixel index, so we correct + // xCM, yCM from pixel coords to pixel index here. + float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + + float minSN2 = PS_SQR(minSN); // calculate the Kron flux, and related fluxes (NO symmetrical averaging) @@ -362,5 +589,5 @@ float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - float fDiff1 = vPix[row][col] - sky; + float fDiff1 = vPix[row][col]; float pDiff = fDiff1; float wDiff = vWgt[row][col]; @@ -376,5 +603,4 @@ Var += wDiff; nKronPix ++; - // if (beVerbose) fprintf (stderr, "mome: %d %d %f %f %f\n", col, row, sky, *vPix, Sum); } @@ -397,10 +623,10 @@ } // *** should I rescale these fluxes by pi R^2 / nNpix? - // XXX source->moments->KronCore = SumCore * M_PI * PS_SQR(sigma) / nCorePix; - // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix; - // XXX source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; - // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; - // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron) - PS_SQR(radKinner)) / nInner; - // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) - PS_SQR(radKron)) / nOuter; + // XXX source->moments->KronCore = SumCore * M_PI * PS_SQR(sigma) / nCorePix; + // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix; + // XXX source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; + // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; + // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron) - PS_SQR(radKinner)) / nInner; + // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) - PS_SQR(radKron)) / nOuter; source->moments->KronCore = SumCore; @@ -408,6 +634,6 @@ source->moments->KronFlux = Sum; source->moments->KronFluxErr = sqrt(Var); - source->moments->KronFinner = SumInner; - source->moments->KronFouter = SumOuter; + source->moments->KronFinner = SumInner; + source->moments->KronFouter = SumOuter; // XXX not sure I should save this here... @@ -416,198 +642,4 @@ source->moments->KronRadiusPSF = source->moments->Mrf; - psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", - source->moments->Mrf, source->moments->KronFlux, - source->moments->Mxx, source->moments->Mxy, source->moments->Myy, - source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, - source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); - - psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", - source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, source->moments->nPixels); - - return(true); -} - -bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { - - // First Pass: calculate the first moments (these are subtracted from the coordinates below) - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // .. etc - - float sky = 0.0; - - float peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - float Sum = 0.0; - float Var = 0.0; - float X1 = 0.0; - float Y1 = 0.0; - float R2 = PS_SQR(radius); - float minSN2 = PS_SQR(minSN); - float rsigma2 = 0.5 / PS_SQR(sigma); - - float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage - float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage - - // we are guaranteed to have a valid pixel and variance at this location (right? right?) - // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); - // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); - - // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should - // not depend on the fractional pixel location of the source. However, the aperture - // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel - // position of the expected centroid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - float yDiff = row + 0.5 - yPeak; - if (fabs(yDiff) > radius) continue; - - float *vPix = source->pixels->data.F32[row]; - float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row]; - - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk & maskVal) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - - float xDiff = col + 0.5 - xPeak; - if (fabs(xDiff) > radius) continue; - - // radius is just a function of (xDiff, yDiff) - float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - if (r2 > R2) continue; - - float pDiff = *vPix - sky; - float wDiff = *vWgt; - - // skip pixels below specified significance level. for a PSFs, this - // over-weights the wings of bright stars compared to those of faint stars. - // for the estimator used for extended source analysis (where the window - // function is allowed to be arbitrarily large), we need to clip to avoid - // negative second moments. - if (PS_SQR(pDiff) < minSN2*wDiff) continue; // - if ((minSN > 0.0) && (pDiff < 0)) continue; // - - // Apply a Gaussian window function. Be careful with the window function. S/N - // weighting over weights the sky for faint sources - if (sigma > 0.0) { - float z = r2*rsigma2; - assert (z >= 0.0); - float weight = exp(-z); - - wDiff *= weight; - pDiff *= weight; - } - - Var += wDiff; - Sum += pDiff; - - float xWght = xDiff * pDiff; - float yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - } - } - - // if we have less than (1/4) of the possible pixels (in circle or box), force a retry - int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0); - - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < minPixels) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum); - return (false); - } - - // calculate the first moment. - float Mx = X1/Sum; - float My = Y1/Sum; - if ((fabs(Mx) > radius) || (fabs(My) > radius)) { - psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y); - return (false); - } - if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) { - psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y); - } - - psTrace ("psModules.objects", 5, "id: %d, sky: %f Mx: %f My: %f Sum: %f X1: %f Y1: %f Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels); - - // add back offset of peak in primary image - // also offset from pixel index to pixel coordinate - // (the calculation above uses pixel index instead of coordinate) - // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords - - // we only update the centroid if the position is not supplied from elsewhere - bool skipCentroid = false; - skipCentroid |= (source->mode & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions - skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions - - if (skipCentroid) { - source->moments->Mx = source->peak->xf; - source->moments->My = source->peak->yf; - } else { - source->moments->Mx = Mx + xGuess; - source->moments->My = My + yGuess; - } - - source->moments->Sum = Sum; - source->moments->SN = Sum / sqrt(Var); - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - return true; } - -float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) { - - psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32); - - for (int i = 0; i < sources->n; i++) { - pmSource *src = sources->data[i]; // Source of interest - if (!src || !src->moments) { - continue; - } - - if (src->mode & PM_SOURCE_MODE_BLEND) { - continue; - } - - if (!src->moments->nPixels) continue; - - if (src->moments->SN < PSF_SN_LIM) continue; - - // XXX put in Mxx,Myy cut based on clump location - - psVectorAppend(radii, src->moments->Mrf); - } - - // find the peak in this image - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - - if (!psVectorStats (stats, radii, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); - psFree(stats); - return NAN; - } - - float minRadius = stats->sampleMedian; - - psFree(radii); - psFree(stats); - return minRadius; -} - Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 36287) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 36375) @@ -113,4 +113,12 @@ source->apFluxErr = NAN; + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + // XXXXXX review: // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this @@ -162,5 +170,5 @@ for (int i = 0; i < source->modelFits->n; i++) { pmModel *model = source->modelFits->data[i]; - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; + if (model->flags & badModel) continue; status = pmSourcePhotometryModel (&model->mag, NULL, model); if (model == source->modelEXT) foundEXT = true; @@ -902,4 +910,40 @@ } +bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(model, false); + + float dC = 0.0; + int Npix = 0; + + psVector *params = model->params; + psImage *image = source->pixels; + psImage *modelFlux = source->modelFlux; + psImage *mask = source->maskObj; + psImage *variance = source->variance; + + float Io = params->data.F32[PM_PAR_I0]; + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + + // skip pixels which are masked + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + + if (variance->data.F32[iy][ix] <= 0) continue; + + dC += PS_SQR (image->data.F32[iy][ix] - Io*modelFlux->data.F32[iy][ix]) / variance->data.F32[iy][ix]; + Npix ++; + } + } + model->nPix = Npix; + model->nDOF = Npix - model->nPar; + model->chisq = dC; + model->chisqNorm = dC / model->nDOF; + + return (true); +} + double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal) {