Index: branches/eam_branches/20090715/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25410) +++ branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25433) @@ -147,5 +147,5 @@ psphotPetrosianStats.c \ psphotPetrosianAnalysis.c \ - pmPetrosian.c + pmPetrosian.c \ psphotEfficiency.c Index: branches/eam_branches/20090715/psphot/src/pmPetrosian.h =================================================================== --- branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25410) +++ branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25433) @@ -35,12 +35,12 @@ bool psphotPetrosianFreeVectors(pmPetrosian *petrosian); -bool psphotPetrosianProfile (pmSource *source, float skynoise); -bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax); -float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax); -bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax); +bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise); +bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petro, int Nsec, float Rmax); +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); +bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax); bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian); -bool psphotEllipticalContour (pmPetrosian *petrosian); +bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian); bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise); -bool psphotPetrosianStats (pmPetrosian *petrosian); +bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian); bool psphotPetrosianSortPair (psVector *index, psVector *extra); Index: branches/eam_branches/20090715/psphot/src/psphot.h =================================================================== --- branches/eam_branches/20090715/psphot/src/psphot.h (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphot.h (revision 25433) @@ -110,4 +110,5 @@ bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type); bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); +float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal); // output functions Index: branches/eam_branches/20090715/psphot/src/psphotBlendFit.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotBlendFit.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotBlendFit.c (revision 25433) @@ -374,4 +374,8 @@ *nfail = Nfail; + // moments are modified by the fit; re-display + psphotVisualPlotMoments (recipe, sources); + psphotVisualShowResidualImage (readout); + return true; } Index: branches/eam_branches/20090715/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotEfficiency.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotEfficiency.c (revision 25433) @@ -134,7 +134,11 @@ psFree(fakeRO); psFree(magAll); + psFree(xAll); + psFree(yAll); return false; } psFree(magAll); + psFree(xAll); + psFree(yAll); psBinaryOp(ro->image, ro->image, "+", fakeRO->image); Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25433) @@ -1,5 +1,3 @@ # include "psphotInternal.h" - -// XXX consistency : theta in radians here and in calling functions // model parameters @@ -7,5 +5,5 @@ psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord); -bool psphotEllipticalContour (pmPetrosian *petrosian) { +bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian) { // use LMM to fit theta vs radius to an ellipse @@ -19,10 +17,11 @@ // arrays to hold the data to be fitted // we fit x and y vs theta in separate passes. - psArray *x = psArrayAlloc(2*radius->n); - psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32); + psArray *x = psArrayAllocEmpty(2*radius->n); + psVector *y = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32); + psVector *yErr = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32); int n = 0; for (int i = 0; i < radius->n; i++) { + if (!isfinite(radius->data.F32[i])) continue; psVector *coord = NULL; @@ -50,6 +49,14 @@ Rmax = MAX (Rmax, radius->data.F32[i]); } - assert (x->n == n); - assert (y->n == n); + x->n = n; + y->n = n; + yErr->n = n; + + if (n < 4) { + psFree (x); + psFree (y); + psFree (yErr); + return false; + } psVector *params = psVectorAlloc (3, PS_TYPE_F32); @@ -61,4 +68,5 @@ // XXX for now, no parameter masks, skip checkLimits + // XXX might help to add a limit to the angle or re-parameterize the ellipse in terms of Rxx, Ryy, Rxy // constraint->checkLimits = psastroModelBoresiteLimits; @@ -84,15 +92,19 @@ } - fprintf (stderr, "# fitted values:\n"); - fprintf (stderr, "Po: %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD); - fprintf (stderr, "Ep: %f\n", params->data.F32[PAR_EPSILON]); - fprintf (stderr, "Rm: %f\n", params->data.F32[PAR_RMIN]); - + psTrace ("psphot", 4, "# fitted values:\n"); + psTrace ("psphot", 4, "Phi: %f\n", petrosian->axes.theta*PS_DEG_RAD); + psTrace ("psphot", 4, "Rmaj: %f\n", petrosian->axes.major); + psTrace ("psphot", 4, "Rmin: %f\n", petrosian->axes.minor); + // show the results - psphotPetrosianVisualEllipticalContour (petrosian); + // psphotPetrosianVisualEllipticalContour (petrosian); psFree (x); psFree (y); psFree (yErr); + psFree (params); + psFree (covar); + psFree (myMin); + psFree (constraint); return true; Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25433) @@ -32,4 +32,5 @@ float Syy = shape.sy; + // XXX drop these two vectors? psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32); psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32); @@ -65,5 +66,8 @@ // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0); - psphotPetrosianVisualProfileByAngle (radius, flux); + // psphotPetrosianVisualProfileByAngle (radius, flux); + + psFree (radiusRaw); + psFree (fluxRaw); return true; } Index: branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c (revision 25433) @@ -239,5 +239,5 @@ psphotVisualShowResidualImage (readout); - psphotVisualShowFlags (sources); + // psphotVisualShowFlags (sources); return true; Index: branches/eam_branches/20090715/psphot/src/psphotImageQuality.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotImageQuality.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotImageQuality.c (revision 25433) @@ -279,8 +279,8 @@ #endif - psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) : %f, %f\n", + psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n", M2->n, fwhm_major, fwhm_minor); - psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n", + psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f [pixels^n]\n", vM2, dM2, vM3, dM3, vM4, dM4); Index: branches/eam_branches/20090715/psphot/src/psphotPSFConvModel.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPSFConvModel.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPSFConvModel.c (revision 25433) @@ -68,4 +68,5 @@ psVector *dparams = modelConv->dparams; + // XXX this needs to be changed to use the new psphotSetRadiusEXT function psphotCheckRadiusEXT (readout, source, modelConv, markVal); Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25433) @@ -12,4 +12,11 @@ // XXX temporary user-supplied systematic sky noise measurement (derive from background model) float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + +# if (0) + // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN + // XXX use this to set skynoise + pmReadout *backModel = psphotSelectBackground (config, view); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); +# endif // S/N limit to perform full non-linear fits @@ -33,9 +40,11 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? if (source->peak->SN < SN_LIM) continue; - if (source->extNsigma < 10.0) continue; // XXX this should not be hardwired // limit selection by analysis region @@ -50,6 +59,11 @@ } - psphotPetrosianProfile (source, skynoise); + psphotPetrosianProfile (readout, source, skynoise); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } + + psphotVisualShowResidualImage (readout); // pause and wait for user input: Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25433) @@ -7,65 +7,73 @@ // structure to something the pmRadialProfile -bool psphotPetrosianProfile (pmSource *source, float skynoise) { +bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise) { - // container to hold results from the radial profile analysis - pmPetrosian *petrosian = pmPetrosianAlloc(); + // container to hold results from the radial profile analysis + pmPetrosian *petrosian = pmPetrosianAlloc(); - int Nsec = 24; - float Rmax = 200; - float fluxMin = 0.0; - float fluxMax = source->peak->flux; + int Nsec = 24; + float Rmax = 200; + float fluxMin = 0.0; + float fluxMax = source->peak->flux; - // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux - // is measured by interpolation for small radii; for large radii, the pixels in a box - // are averaged to increase the S/N (XXX not yet done) - if (!psphotRadialProfilesByAngles (petrosian, source, Nsec, Rmax)) { - psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); - return false; - } + // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux + // is measured by interpolation for small radii; for large radii, the pixels in a box + // are averaged to increase the S/N (XXX not yet done) + if (!psphotRadialProfilesByAngles (source, petrosian, Nsec, Rmax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); + psFree (petrosian); + return false; + } - // use the radial profiles to determine the radius of a given isophote. this isophote - // is used to determine the elliptical shape of the object, so it has a relatively high - // value (nominally 50% of the peak) - if (!psphotRadiiFromProfiles (petrosian, fluxMin, fluxMax)) { - psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); - return false; - } + // use the radial profiles to determine the radius of a given isophote. this isophote + // is used to determine the elliptical shape of the object, so it has a relatively high + // value (nominally 50% of the peak) + if (!psphotRadiiFromProfiles (source, petrosian, fluxMin, fluxMax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); + psFree (petrosian); + return false; + } - // convert the isophotal radius vs angle measurements to an elliptical contour - if (!psphotEllipticalContour (petrosian)) { - psError (PS_ERR_UNKNOWN, false, "failed to measure elliptical contour"); - return false; - } + // convert the isophotal radius vs angle measurements to an elliptical contour + if (!psphotEllipticalContour (source, petrosian)) { + psLogMsg ("psphot", 3, "failed to measure elliptical contour"); + psFree (petrosian); + return false; + } - // generate a single, normalized radial profile following the elliptical contours. - // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel - if (!psphotEllipticalProfile (source, petrosian)) { - psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); - return false; - } + // generate a single, normalized radial profile following the elliptical contours. + // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel + if (!psphotEllipticalProfile (source, petrosian)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } - // integrate the radial profile for radial bins defined for the petrosian measurement: - // SB_i (r_i) where \alpha r_i < r < \beta r_i - if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) { - psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); - return false; - } + // integrate the radial profile for radial bins defined for the petrosian measurement: + // SB_i (r_i) where \alpha r_i < r < \beta r_i + if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } - // use the SB_i from above to calculate the petrosian radius and the flux within that radius - if (!psphotPetrosianStats (petrosian)) { - psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); - return false; - } + // use the SB_i from above to calculate the petrosian radius and the flux within that radius + if (!psphotPetrosianStats (source, petrosian)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } - // XXX this will only work in the psphot context, not the psphotPetrosianStudy... - // XXX add the petrosian to the pmSource structure... - psphotVisualShowPetrosian (source, petrosian); + // XXX this will only work in the psphot context, not the psphotPetrosianStudy... + // XXX add the petrosian to the pmSource structure... + // psphotVisualShowResidualImage (readout); + psphotVisualShowPetrosian (source, petrosian); - psphotPetrosianFreeVectors(petrosian); + psphotPetrosianFreeVectors(petrosian); - fprintf (stdout, "\n petrosian radius: %f\n flux: %f\n axis ratio: %f\n angle: %f\n", - petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta); + psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f", + source->peak->xf, source->peak->yf, petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta); - return true; + psFree (petrosian); + return true; } Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25433) @@ -18,10 +18,4 @@ float skyModelErrorSQ = PS_SQR(skynoise); - -# if (0) - // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN - pmReadout *backModel = psphotSelectBackground (config, view); - pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); -# endif psVector *radius = petrosian->radiusElliptical; @@ -104,9 +98,7 @@ binArea->data.F32[i] = M_PI * (rMax2 - rMin2); - if (0) { - fprintf (stderr, "%3d %5.1f %5.1f : %5.1f : %5.1f %5.1f\n", + psTrace ("psphot", 6, "%3d %5.1f %5.1f : %5.1f : %5.1f %5.1f\n", i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i], radMax->data.F32[i], radBet->data.F32[i]); - } } @@ -131,8 +123,16 @@ if (radius->data.F32[i] > Rmax) { // calculate the value for the nOut bin - psVectorStats (stats, values, NULL, NULL, 0); + float value, dvalue; + if (values->n > 0) { + psVectorStats (stats, values, NULL, NULL, 0); + value = stats->robustMedian; + dvalue = stats->robustStdev; + } else { + value = NAN; + dvalue = NAN; + } // binSB->data.F32[nOut] = stats->sampleMedian; - binSB->data.F32[nOut] = stats->robustMedian; - binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->robustStdev) / values->n + skyModelErrorSQ); + binSB->data.F32[nOut] = value; + binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ); // binSB->data.F32[nOut] = stats->fittedMean; // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ); @@ -142,8 +142,6 @@ // residual flux, but the sky from the sky model) - if (1) { - fprintf (stderr, "%3d %5.1f %5.1f : %5.1f %5.2f\n", - nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); - } + psTrace ("psphot", 5, "%3d %5.1f %5.1f : %5.1f %5.2f\n", + nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); nOut ++; @@ -170,5 +168,5 @@ petrosian->binSBstdev = binSBstdev; - psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0); + // psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0); psFree(radMin); Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25433) @@ -8,5 +8,5 @@ float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); -bool psphotPetrosianStats (pmPetrosian *petrosian) { +bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian) { float petRadius, petFlux; @@ -77,12 +77,10 @@ psVectorAppend(refRadius, binRad->data.F32[i]); - if (1) { - fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", - i, refRadius->data.F32[nOut], - binSB->data.F32[i], binSBstdev->data.F32[i], - meanSB->data.F32[nOut], meanSBerr, - petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], - fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner); - } + psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", + i, refRadius->data.F32[nOut], + binSB->data.F32[i], binSBstdev->data.F32[i], + meanSB->data.F32[nOut], meanSBerr, + petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], + fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner); // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux @@ -149,5 +147,8 @@ psFree(fluxSum); + psFree(fluxSumErr2); + psFree(refRadius); psFree(petRatio); + psFree(petRatioErr); psFree(meanSB); psFree(areaSum); Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25433) @@ -136,5 +136,5 @@ } - psphotPetrosianProfile (source, skynoise); + psphotPetrosianProfile (readout, source, skynoise); psFree (source); Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25433) @@ -19,4 +19,5 @@ static int kapa = -1; +static int kapa2 = -1; // if no valid data is supplied (NULL or n <- 0), leave limits as they were @@ -53,9 +54,9 @@ // return true; - if (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa == -1) { + // if (!pmVisualIsVisual()) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa2 == -1) { fprintf (stderr, "failure to open kapa; visual mode disabled\n"); pmVisualSetVisual(false); @@ -64,22 +65,22 @@ } - KapaClearPlots (kapa); + KapaClearPlots (kapa2); KapaInitGraph (&graphdata); - KapaSetFont (kapa, "courier", 14); + KapaSetFont (kapa2, "courier", 14); pmVisualLimitsFromVectors (&graphdata, radius, flux); - KapaSetLimits (kapa, &graphdata); - - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); - KapaSendLabel (kapa, "flux", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.style = 2; - graphdata.ptype = 0; - graphdata.size = 1.0; - KapaPrepPlot (kapa, radius->n, &graphdata); - KapaPlotVector (kapa, radius->n, radius->data.F32, "x"); - KapaPlotVector (kapa, radius->n, flux->data.F32, "y"); + KapaSetLimits (kapa2, &graphdata); + + KapaBox (kapa2, &graphdata); + KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "flux", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa2, radius->n, &graphdata); + KapaPlotVector (kapa2, radius->n, radius->data.F32, "x"); + KapaPlotVector (kapa2, radius->n, flux->data.F32, "y"); // pause and wait for user input: Index: branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25433) @@ -1,7 +1,8 @@ # include "psphotInternal.h" -// Given a source at (x,y), generate a collection of radial profiles at even angular -// separations - +// Given a source at (x,y), generate a collection of radial profiles at even angular separations + +// These functions are used to calculate the stats in a rectangle at arbitrary orientation. +// XXX Move these elsewhere (psLib?) float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta); psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta); @@ -9,5 +10,5 @@ psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); -bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) { +bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petrosian, int Nsec, float Rmax) { // we want to have an even number of sectors so we can do 180 deg symmetrizing @@ -124,4 +125,6 @@ psVector *values = psphotBoxValues (image, x, y, dL, dW, theta); + if (!values) goto escape; + if (!values->n) goto escape; psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); @@ -134,4 +137,8 @@ return value; + +escape: + psFree(values); + return NAN; } Index: branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25433) @@ -3,5 +3,5 @@ // Given the Petrosian data (radii, fluxes) determine the radius for each profile at the desisred isophote -bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax) { +bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax) { petrosian->isophotalRadii = psVectorAlloc(petrosian->theta->n, PS_TYPE_F32); @@ -10,5 +10,5 @@ psVector *radii = petrosian->radii->data[i]; psVector *fluxes = petrosian->fluxes->data[i]; - float radius = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax); + float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax); // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius); @@ -20,5 +20,5 @@ } -float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax) { +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax) { // 'flux' is a noisy sample of the galaxy radial profile at points 'radius' @@ -49,10 +49,14 @@ psVectorAppend (values, radius->data.F32[i]); } - psVectorStats (stats, values, NULL, NULL, 0); + if (values->n > 1) { + psVectorStats (stats, values, NULL, NULL, 0); - // if we have a valid range, rebin with bin size 1/2 of median radius - if (isfinite(stats->sampleMedian)) { - Rbin = MAX(1, 0.5*stats->sampleMedian); + // if we have a valid range, rebin with bin size 1/2 of median radius + if (isfinite(stats->sampleMedian)) { + Rbin = MAX(1, 0.5*stats->sampleMedian); + } } + psFree (values); + psFree (stats); } Rbin = 3; @@ -90,6 +94,12 @@ // calculate the value for the nOut bin // XXX need to fix this as well psStats (stats, values); - psVectorStats (stats, values, NULL, NULL, 0); - fluxBinned->data.F32[nOut] = stats->sampleMedian; + float value; + if (values->n > 0) { + psVectorStats (stats, values, NULL, NULL, 0); + value = stats->sampleMedian; + } else { + value = NAN; + } + fluxBinned->data.F32[nOut] = value; nOut ++; radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin; @@ -118,9 +128,10 @@ // XXX is there a macro in psLib that does this interpolation? if (i == 0) { - psWarning ("impossible condition f[0] < Fo"); - continue; - } else { - Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]); - } + psLogMsg ("psphot", 3, "bogus radial profile for ..., skipping"); + psFree (fluxBinned); + psFree (radiusBinned); + return NAN; + } + Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]); above = FALSE; } Index: branches/eam_branches/20090715/psphot/src/psphotRadiusChecks.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadiusChecks.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotRadiusChecks.c (revision 25433) @@ -86,5 +86,50 @@ // call this function whenever you (re)-define the EXT model +float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) { + + psAssert (source, "source not defined??"); + psAssert (source->peak, "peak not defined??"); + + pmPeak *peak = source->peak; + + // set the radius based on the footprint: + if (!peak->footprint) goto escape; + pmFootprint *footprint = peak->footprint; + if (!footprint->spans) goto escape; + if (footprint->spans->n < 1) goto escape; + + // find the max radius + float radius = 0.0; + for (int j = 0; j < footprint->spans->n; j++) { + pmSpan *span = footprint->spans->data[j]; + + float dY = span->y - peak->yf; + float dX0 = span->x0 - peak->xf; + float dX1 = span->x1 - peak->xf; + + radius = PS_MAX (radius, hypot(dY, dX0)); + radius = PS_MAX (radius, hypot(dY, dX1)); + } + + radius += EXT_FIT_PADDING; + if (isnan(radius)) psAbort("error in radius"); + + // redefine the pixels if needed + pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal); + return radius; + +escape: + return NAN; + // bool result = psphotCheckRadiusEXT_old (readout, source, model, markVal); + // return result; +} + +// call this function whenever you (re)-define the EXT model bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) { + + psAbort ("do not use this function"); psF32 *PAR = model->params->data.F32; Index: branches/eam_branches/20090715/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotSourceFits.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotSourceFits.c (revision 25433) @@ -223,7 +223,13 @@ if (source->peak->SN < EXT_MIN_SN) return false; + float radius = psphotSetRadiusEXT (readout, source, markVal); + + // XXX note that this changes the source moments that are published... // recalculate the source moments using the larger extended-source moments radius // at this stage, skip Gaussian windowing, and do not clip pixels by S/N - if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false; + // XXX use the deteremined value of PSF_MOMENTS_RADIUS and SIGMA and scale up? + // XXX increase the moments radius iteratively until the measured sigma / input < f? + // XXX use the footprint somehow to judge both radius and aperture? + if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false; psTrace ("psphot", 5, "trying blob...\n"); @@ -277,4 +283,5 @@ // both models failed; reject them both + // XXX -- change type flags to psf in this case and keep original moments? psFree (EXT); psFree (DBL); @@ -287,4 +294,5 @@ // save new model source->modelEXT = EXT; + source->modelEXT->radiusFit = radius; source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXTMODEL; @@ -294,5 +302,15 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + +# if (PS_TRACE_ON) psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); + if (psTraceGetLevel("psphot") >= 6) { + psLogMsg ("psphot", 1, "source 2:\n"); + for (int i = 0; i < source->modelEXT->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]); + } + } +# endif + return true; @@ -306,9 +324,10 @@ source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; source->mode |= PM_SOURCE_MODE_PAIR; + source->modelPSF->radiusFit = radius; // copy most data from the primary source (modelEXT, blends stay NULL) - // XXX use pmSourceCopy? pmSource *newSrc = pmSourceCopy (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); + newSrc->modelPSF->radiusFit = radius; // build cached models and subtract @@ -317,5 +336,19 @@ pmSourceCacheModel (newSrc, maskVal); pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal); + +# if (PS_TRACE_ON) psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]); + if (psTraceGetLevel("psphot") >= 6) { + psLogMsg ("psphot", 1, "source 1:\n"); + for (int i = 0; i < newSrc->modelPSF->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]); + } + psLogMsg ("psphot", 1, "source 2:\n"); + for (int i = 0; i < source->modelPSF->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]); + } + psphotVisualShowResidualImage (readout); + } +# endif psArrayAdd (newSources, 100, newSrc); @@ -356,5 +389,5 @@ // save the PSF model from the Ensemble fit PSF = source->modelPSF; - psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0); + // psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0); if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit"); @@ -390,6 +423,5 @@ // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit"); - - psphotCheckRadiusEXT (readout, source, EXT, markVal); + // psphotCheckRadiusEXT (readout, source, EXT, markVal); if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { @@ -401,3 +433,2 @@ return (EXT); } - Index: branches/eam_branches/20090715/psphot/src/psphotSourceSize.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25433) @@ -194,7 +194,5 @@ bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { - // select the psf stars - psArray *psfstars = psArrayAllocEmpty (100); - + // select stats from the psf stars psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -205,5 +203,4 @@ pmSource *source = sources->data[i]; if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - psArrayAdd (psfstars, 100, source); // XXX can we test if psfMag is set and calculate only if needed? @@ -229,4 +226,9 @@ options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); fprintf (stderr, "psf-sum: %f +/- %f\n", options->ApResid, options->ApSysErr); + + psFree (Ap); + psFree (ApErr); + psFree (stats); + psFree (dAp); return true; Index: branches/eam_branches/20090715/psphot/src/psphotVisual.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotVisual.c (revision 25410) +++ branches/eam_branches/20090715/psphot/src/psphotVisual.c (revision 25433) @@ -166,5 +166,5 @@ bool psphotVisualShowImage (pmReadout *readout) { - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; int kapa = psphotKapaChannel (1); @@ -2000,5 +2000,5 @@ bool psphotVisualShowResidualImage (pmReadout *readout) { - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; int myKapa = psphotKapaChannel (1); @@ -2100,7 +2100,7 @@ if (source == NULL) return true; - if (!pmVisualIsVisual()) return true; - - int kapa = psphotKapaChannel (3); + // if (!pmVisualIsVisual()) return true; + + int kapa = psphotKapaChannel (1); if (kapa == -1) return false;