Index: trunk/psLib/src/imageops/psImageMapFit.c =================================================================== --- trunk/psLib/src/imageops/psImageMapFit.c (revision 37808) +++ trunk/psLib/src/imageops/psImageMapFit.c (revision 37811) @@ -286,9 +286,16 @@ // fprintf (stderr, "Total: %f\n", Total); + double MaxPivot = 0.0; + for (int i = 0; i < Nx*Ny; i++) { + MaxPivot = PS_MAX(MaxPivot, fabs(A->data.F32[i][i])); + // fprintf (stderr, "piv, max: %f : %f\n", A->data.F32[i][i], MaxPivot); + } + // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0 psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8); psVectorInit (Empty, 0); + double MinPivot = 0.025*MaxPivot; for (int i = 0; i < Nx*Ny; i++) { - if (A->data.F32[i][i] == 0.0) { + if (fabs(A->data.F32[i][i]) < MinPivot) { Empty->data.S8[i] = 1; for (int j = 0; j < Nx*Ny; j++) { Index: trunk/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 37808) +++ trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 37811) @@ -45,4 +45,6 @@ #include "pmSourceVisual.h" +bool pmPSF_DataDump (char *filename, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *mask); + /***************************************************************************** pmPSFFromPSFtry (psfTry): build a PSF model from a collection of source->modelEXT entries @@ -356,2 +358,21 @@ } +bool pmPSF_DataDump (char *filename, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *mask) { + + + FILE *f = fopen (filename, "w"); + if (!f) return false; + + for (int i = 0; i < x->n; i++) { + + fprintf (f, "%6.1f %6.1f : %5.2f %5.2f %5.2f : %2d\n", + x->data.F32[i], + y->data.F32[i], + e0->data.F32[i], + e1->data.F32[i], + e2->data.F32[i], + mask->data.U8[i]); + } + fclose (f); + return true; +} Index: trunk/psModules/src/objects/pmTrend2D.c =================================================================== --- trunk/psModules/src/objects/pmTrend2D.c (revision 37808) +++ trunk/psModules/src/objects/pmTrend2D.c (revision 37811) @@ -211,4 +211,17 @@ // XXX need to add the API which adjusts the scale status = psImageMapClipFit(pGoodFit, trend->map, trend->stats, mask, maskVal, x, y, f, df); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "failed to build PSF model map"); + return false; + } + + // the psf model map can have nan pixels: repair this by extrapolation / interpolation + // XXX TEST: p_psImagePrint(0, trend->map->map, "before"); + status = psImageMapRepair (trend->map->map); + // XXX TEST: p_psImagePrint(0, trend->map->map, "after"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "failed to repair PSF model map"); + return false; + } break; @@ -217,4 +230,57 @@ } return status; +} + +bool psImageMapRepair (psImage *map) { + + + // XXX why is my repair not working?? + + // patch over bad regions (use average of 8 possible neighbor pixels) + // XXX consider testing all pixels against the 8 neighbors and replacing outliers... + double Count = 0; // number of good pixels + double Value = 0; // sum of good pixel's value + for (int iy = 0; iy < map->numRows; iy++) { + for (int ix = 0; ix < map->numCols; ix++) { + if (isfinite(map->data.F32[iy][ix])) { + Value += map->data.F32[iy][ix]; + Count++; + continue; + } + + double value = 0; + double count = 0; + for (int jy = iy - 1; jy <= iy + 1; jy++) { + if (jy < 0) continue; + if (jy >= map->numRows) continue; + for (int jx = ix - 1; jx <= ix + 1; jx++) { + if (!jx && !jy) continue; + if (jx < 0) continue; + if (jx >= map->numCols) continue; + if (!isfinite(map->data.F32[jy][jx])) continue; + value += map->data.F32[jy][jx]; + count += 1.0; + } + } + if (count > 0) { + // psLogMsg ("psphot", PS_LOG_DETAIL, "patching image map %d, %d: %f (%d pts)\n", ix, iy, (value / count), (int) count); + map->data.F32[iy][ix] = value / count; + } + } + } + if (Count == 0) { + psError(PS_ERR_UNKNOWN, true, "failed to repair PSF model map"); + return false; + } + Value /= Count; + + // patch over remaining bad regions (use global average) + for (int iy = 0; iy < map->numRows; iy++) { + for (int ix = 0; ix < map->numCols; ix++) { + if (!isnan(map->data.F32[iy][ix])) continue; + map->data.F32[iy][ix] = Value; + } + } + return true; } Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 37808) +++ trunk/psphot/src/psphotChoosePSF.c (revision 37811) @@ -476,6 +476,6 @@ psFree (modelPSF); - float fwhmtest = pmPSFtoFWHM(psf, xc, yc); - fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + // float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); } } @@ -630,6 +630,6 @@ psFree (modelPSF); - float fwhmtest = pmPSFtoFWHM(psf, xc, yc); - fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + // float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); } Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 37808) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 37811) @@ -38,11 +38,11 @@ /*** for the moment, this test code : it is not thread safe ***/ -int Nall = 0; -int Nskip1 = 0; -int Nskip2 = 0; -int Nskip3 = 0; -int Nskip4 = 0; -int Nskip5 = 0; -int Nskip6 = 0; +static int Nall = 0; +static int Nskip1 = 0; +static int Nskip2 = 0; +static int Nskip3 = 0; +static int Nskip4 = 0; +static int Nskip5 = 0; +static int Nskip6 = 0; # define SKIP(VALUE) { VALUE++; continue; } Index: trunk/psphot/src/psphotGalaxyParams.c =================================================================== --- trunk/psphot/src/psphotGalaxyParams.c (revision 37808) +++ trunk/psphot/src/psphotGalaxyParams.c (revision 37811) @@ -1,5 +1,6 @@ # include "psphotInternal.h" -#define USE_FOOTPRINTS 1 +# define PRINT_STUFF 0 +# define USE_FOOTPRINTS 1 // measure some properties on the exended objects @@ -40,19 +41,19 @@ } -#ifdef notdef +// # define COUNT_SKIPS +# ifdef COUNT_SKIPS /*** for the moment, this test code : it is not thread safe ***/ -int Nall = 0; -int Nskip1 = 0; -int Nskip2 = 0; -int Nskip3 = 0; -int Nskip4 = 0; -int Nskip5 = 0; -int Nskip6 = 0; - +static int Nskip1 = 0; +static int Nskip2 = 0; +static int Nskip3 = 0; +static int Nskip4 = 0; +static int Nskip5 = 0; +static int Nskip6 = 0; +static int Nskip7 = 0; +static int Nskip8 = 0; # define SKIP(VALUE) { VALUE++; continue; } - -#endif //notdef - -#define SKIP(VALUE) {continue;} +# else +# define SKIP(VALUE) {continue;} +# endif //notdef // aperture-like measurements for extended sources @@ -104,8 +105,8 @@ psImage *footprintImage = NULL; -#ifdef USE_FOOTPRINTS +#if USE_FOOTPRINTS footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32); psImageInit(footprintImage, 0); - pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, true); + pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, false); #endif @@ -207,5 +208,5 @@ // # if (PS_TRACE_ON) -# if (0) +# ifdef COUNT_SKIPS fprintf (stderr, "ext analysis skipped @ 1 : %d\n", Nskip1); fprintf (stderr, "ext analysis skipped @ 2 : %d\n", Nskip2); @@ -214,4 +215,6 @@ fprintf (stderr, "ext analysis skipped @ 5 : %d\n", Nskip5); fprintf (stderr, "ext analysis skipped @ 6 : %d\n", Nskip6); + fprintf (stderr, "ext analysis skipped @ 7 : %d\n", Nskip7); + fprintf (stderr, "ext analysis skipped @ 8 : %d\n", Nskip8); #endif @@ -221,4 +224,7 @@ // psphotVisualShowPetrosians (sources); + + // XXX TEST: SAVE IMAGE: + // psphotSaveImage (NULL, readout->image, "resid.fits"); return true; @@ -240,5 +246,5 @@ psMetadata *recipe = job->args->data[3]; -#ifndef USE_FOOTPRINTS +#if (!USE_FOOTPRINTS) float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); #endif @@ -256,5 +262,5 @@ // don't try and use bad models pmModelStatus badModel = PM_MODEL_STATUS_NONE; - badModel |= PM_MODEL_STATUS_NONCONVERGE; +// badModel |= PM_MODEL_STATUS_NONCONVERGE; this causes most objects to get unmeasured badModel |= PM_MODEL_STATUS_BADARGS; badModel |= PM_MODEL_STATUS_OFFIMAGE; @@ -291,4 +297,13 @@ // Nall ++; + + bool testObject = FALSE; + // testObject = testObject || ((fabs(source->peak->xf - 5788) < 10) && (fabs(source->peak->yf - 5711) < 10)); + // testObject = testObject || ((fabs(source->peak->xf - 96) < 10) && (fabs(source->peak->yf - 5550) < 10)); + // testObject = testObject || ((fabs(source->peak->xf - 300) < 10) && (fabs(source->peak->yf - 5988) < 10)); + if (testObject) { + fprintf (stderr, "test object: %f, %f\n", source->peak->xf, source->peak->yf); + testObject = TRUE; + } // rules for measuring petrosian parameters for specific objects are set in @@ -403,5 +418,5 @@ #endif -#ifdef USE_FOOTPRINTS +#if (USE_FOOTPRINTS) psS32 **vFootprint = footprintImage->data.S32; #else @@ -436,12 +451,18 @@ } -#ifdef USE_FOOTPRINTS - if (vFootprint[yImage][xImage]) { -#else - if (fabs(pixel) >= threshold) { -#endif +#if (USE_FOOTPRINTS) + // this is a pixel in the object + if (vFootprint[yImage][xImage] == source->peak->footprint->id) { SET_MARK(MARK_SOURCE); continue; + } + // this is a pixel in another object + if (vFootprint[yImage][xImage] > 0.0) continue; +#else + // this is a pixel in an object (this or the other) + if (fabs(pixel) >= threshold) { + continue } +#endif // we have a background pixel. // Check whether it's mirror is a background pixel as well. @@ -456,5 +477,5 @@ if (isnan(mirrorPixel)) continue; -#ifdef USE_FOOTPRINTS +#if (USE_FOOTPRINTS) if (vFootprint[yFlip + row0][xFlip + col0] == 0) { #else @@ -581,10 +602,14 @@ source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA; -#ifdef PRINT_STUFF +#if (PRINT_STUFF) if (Ngood % 40 == 1) { - printf(" id sumI rHL rad sumI gRT gRA sumRt sumBt numPixels numBGPixels\n"); - } - printf("%5d %5.2e %3.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n", - source->id, sumI, source->extpars->ghalfLightRadius, radius, source->extpars->gRT, source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n); + fprintf (stderr, " id sumI rHL rad sumI gRT gRA sumRt sumBt numPixels numBGPixels\n"); + } + fprintf (stderr, "%5d %5.1f %5.1f : %5.2e %4.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n", + source->id, source->peak->xf, source->peak->yf, + sumI, source->extpars->ghalfLightRadius, + radius, + source->extpars->gRT, + source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n); #endif @@ -592,5 +617,11 @@ restoreModel: - if ((extModel != sersicModel) || !isExtended) { + if (testObject) { + fprintf (stderr, "done with %d\n", source->id); + } + +# define USE_BEST_FIT TRUE + // XXX TEST: Keep the sersic model regardless + if (USE_BEST_FIT && ((extModel != sersicModel) || !isExtended)) { // add the sersic model back in if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {