Changeset 19881

Timestamp:

Oct 3, 2008, 11:00:18 AM (18 years ago)

Author:

eugene

Message:

converted pmMoments elements to Mxx, Myy, Mxxx, etc; using pmPeak position in most cases instead of moments

Location:

trunk/psphot/src

Files:

• models/pmModel_STRAIL.c (modified) (3 diffs)
• models/pmModel_TEST1.c (modified) (1 diff)
• psphotBlendFit.c (modified) (4 diffs)
• psphotChoosePSF.c (modified) (1 diff)
• psphotExtendedSourceAnalysis.c (modified) (5 diffs)
• psphotExtendedSourceFits.c (modified) (4 diffs)
• psphotFakeSources.c (modified) (1 diff)
• psphotGuessModels.c (modified) (1 diff)
• psphotModelTest.c (modified) (1 diff)
• psphotMosaicSubimage.c (modified) (1 diff)
• psphotSourceFits.c (modified) (3 diffs)
• psphotSummaryPlots.c (modified) (3 diffs)

trunk/psphot/src/models/pmModel_STRAIL.c

@@ -478 +478 @@
 {
     pmMoments *Smoments = source->moments;
+    pmPeak    *peak    = source->peak;
     psF32     *params  = model->params->data.F32;
 
@@ -484 +485 @@
     psEllipseMoments moments;
 
-    moments.x2 = PS_SQR(Smoments->Sx);
-    moments.y2 = PS_SQR(Smoments->Sy);
-    moments.xy = Smoments->Sxy;
+    moments.x2 = Smoments->Mxx;
+    moments.y2 = Smoments->Myy;
+    moments.xy = Smoments->Mxy;
     //sometimes these moment inputs are zero...why?
 
@@ -494 +495 @@
     shape = psEllipseAxesToShape(axes);
 
-    params[0] = Smoments->Sky;
-    params[1] = Smoments->Peak - Smoments->Sky;
-    params[2] = Smoments->x;
-    params[3] = Smoments->y;
+    params[PM_PAR_SKY]  = Smoments->Sky;
+    params[PM_PAR_I0]   = peak->flux;
+    params[PM_PAR_XPOS] = peak->xf;
+    params[PM_PAR_YPOS] = peak->yf;
+
     params[7] = 2 * axes.major;
     params[8] = axes.theta;

trunk/psphot/src/models/pmModel_TEST1.c

@@ -126 +126 @@
 {
     pmMoments *moments = source->moments;
+    pmPeak    *peak    = source->peak;
     psF32     *PAR  = model->params->data.F32;
 
+    psEllipseMoments emoments;
+    emoments.x2 = moments->Mxx;
+    emoments.y2 = moments->Myy;
+    emoments.xy = moments->Mxy;
+
+    // force the axis ratio to be < 20.0
+    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
+    psEllipseShape shape = psEllipseAxesToShape (axes);
+
     PAR[PM_PAR_SKY] = moments->Sky;
-    PAR[PM_PAR_I0] = moments->Peak - moments->Sky;
-    PAR[PM_PAR_XPOS] = moments->x;
-    PAR[PM_PAR_YPOS] = moments->y;
-    PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx);
-    PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy);
+    PAR[PM_PAR_I0]   = peak->flux;
+    PAR[PM_PAR_XPOS] = peak->xf;
+    PAR[PM_PAR_YPOS] = peak->yf;
+    PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
+    PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy);
+    PAR[PM_PAR_SXY] = shape.sxy;
     PAR[PM_PAR_SXY] = 0.0;  // XXX we can get this right if we do the integral
 

trunk/psphot/src/psphotBlendFit.c

@@ -55 +55 @@
 
         // limit selection to some SN limit
-        // XXX this should use peak?
-        if (source->moments == NULL) continue;
-        if (source->moments->SN < FIT_SN_LIM) continue;
+        if (source->peak->SN < FIT_SN_LIM) continue;
 
         // XXX this should use peak?
-        if (source->moments->x < AnalysisRegion.x0) continue;
-        if (source->moments->y < AnalysisRegion.y0) continue;
-        if (source->moments->x > AnalysisRegion.x1) continue;
-        if (source->moments->y > AnalysisRegion.y1) continue;
+        if (source->peak->xf < AnalysisRegion.x0) continue;
+        if (source->peak->yf < AnalysisRegion.y0) continue;
+        if (source->peak->xf > AnalysisRegion.x1) continue;
+        if (source->peak->yf > AnalysisRegion.y1) continue;
 
         // if model is NULL, we don't have a starting guess
@@ -95 +93 @@
             if (psphotFitBlob (readout, source, sources, psf, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
-                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
                 continue;
@@ -102 +100 @@
             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_STAR;
-                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
                 Npsf ++;
                 continue;
@@ -108 +106 @@
         }
 
-        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y);
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
         Nfail ++;
 

trunk/psphot/src/psphotChoosePSF.c

@@ -474 +474 @@
         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
 
-        // moments->Sx,Sy are roughly sigma_x,y
-        moments.x2 = PS_SQR(source->moments->Sx);
-        moments.y2 = PS_SQR(source->moments->Sy);
-        moments.xy = source->moments->Sxy;
+        moments.x2 = source->moments->Mxx;
+        moments.y2 = source->moments->Myy;
+        moments.xy = source->moments->Mxy;
 
         // limit axis ratio < 20.0

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -71 +71 @@
                 // all measurements below require the radial profile; skip them all
                 // re-subtract the object, leave local sky
-                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
                 source->mode |= PM_SOURCE_MODE_SUBTRACTED;
@@ -81 +81 @@
         if (doIsophotal) {
             if (!psphotIsophotal (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             } else {
-                psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 Nisophot ++;
             }
@@ -91 +91 @@
         if (doPetrosian) {
             if (!psphotPetrosian (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             } else {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 Npetro ++;
             }
@@ -101 +101 @@
         if (doKron) {
             if (!psphotKron (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             } else {
-                psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+                psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 Nkron ++;
             }
@@ -114 +114 @@
                 return false;
             } 
-            psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             Nannuli ++;
         }

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -163 +163 @@
               modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
               if (!modelFit) {
-                  psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->x, source->moments->y);
+                  psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
                   continue;
               }
-              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
               Nconvolve ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -176 +176 @@
               modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal);
               if (!modelFit) {
-                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->x, source->moments->y);
+                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
                   continue;
               }
               pmSourceCacheModel (source, maskVal);
-              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq);
+              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq);
               Nplain ++;
               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
@@ -217 +217 @@
         if (minModel == -1) {
           // re-subtract the object, leave local sky
-          psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y);
+          psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->Mx, source->moments->My);
 
           // replace original model, subtract it
@@ -257 +257 @@
         psFree (modelFluxes);
 
-        psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
-        psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+        psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
+        psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 
         if (savePics) {

trunk/psphot/src/psphotFakeSources.c

@@ -10 +10 @@
         pmSource *source = pmSourceAlloc ();
         source->moments = pmMomentsAlloc ();
-        source->moments->x = 10;
-        source->moments->y = 10;
-        source->moments->Sx = 1;
-        source->moments->Sy = 1;
-        source->moments->Sxy = 0;
+        source->moments->Mx = 10;
+        source->moments->My = 10;
+        source->moments->Mxx = 1;
+        source->moments->Myy = 1;
+        source->moments->Mxy = 0;
         source->moments->Sum = 1000;
         source->moments->Peak = 100;

trunk/psphot/src/psphotGuessModels.c

@@ -59 +59 @@
             // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
             if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
-                modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->x;
-                modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y;
+                modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
+                modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
             } else {
                 modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;

trunk/psphot/src/psphotModelTest.c

@@ -143 +143 @@
     source->peak->value = source->moments->Peak;
 
-    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y);
-    fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Sx, source->moments->Sy, source->moments->Sxy);
+    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
+    fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
 
     psEllipseMoments moments;
-    moments.x2 = source->moments->Sx;
-    moments.y2 = source->moments->Sy;
-    moments.xy = source->moments->Sxy;
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
     psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
 

trunk/psphot/src/psphotMosaicSubimage.c

@@ -38 +38 @@
                 newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak;
             } else {
-                newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak;
+                newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix];
             }
         }

trunk/psphot/src/psphotSourceFits.c

@@ -218 +218 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
-    if (source->moments->SN < EXT_MIN_SN) return false;
+    if (source->peak->SN < EXT_MIN_SN) return false;
 
     // recalculate the source moments using the larger extended-source moments radius
@@ -331 +331 @@
     maskVal |= markVal;
 
+    // XXX this is really poor: if we don't have moments for the source, we have no guess.
+    // force the measurement?
+    psAssert (source->moments, "moments are re-calculated for any extended source");
+        
     // make a guess at the position of the two sources
-    moments.x2 = source->moments->Sx;
-    moments.y2 = source->moments->Sy;
-    moments.xy = source->moments->Sxy;
+    moments.x2 = source->moments->Mxx;
+    moments.xy = source->moments->Mxy;
+    moments.y2 = source->moments->Myy;
     axes = psEllipseMomentsToAxes (moments, 20.0);
 
@@ -380 +384 @@
     psphotCheckRadiusEXT (readout, source, EXT, markVal);
 
-    if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) {
-        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy);
+    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
+        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
     }
 

trunk/psphot/src/psphotSummaryPlots.c

@@ -61 +61 @@
         if (source->moments == NULL) continue;
     
-        xFaint->data.F32[nF] = source->moments->Sx;
-        yFaint->data.F32[nF] = source->moments->Sy;
+        xFaint->data.F32[nF] = source->moments->Mxx;
+        yFaint->data.F32[nF] = source->moments->Myy;
         nF++;
     
@@ -68 +68 @@
         if (source->moments->SN < 25) continue;
 
-        xBright->data.F32[nB] = source->moments->Sx;
-        yBright->data.F32[nB] = source->moments->Sy;
+        xBright->data.F32[nB] = source->moments->Mxx;
+        yBright->data.F32[nB] = source->moments->Myy;
         nB++;
     }
@@ -160 +160 @@
 
         psEllipseMoments moments;
-        moments.x2 = source->moments->Sx;
-        moments.y2 = source->moments->Sy;
-        moments.xy = source->moments->Sxy;
+        moments.x2 = source->moments->Mxx;
+        moments.xy = source->moments->Mxy;
+        moments.y2 = source->moments->Myy;
 
         psEllipseShape shape;