Changeset 13064 for trunk/psModules/src/objects/pmPSFtry.c

Timestamp:

Apr 27, 2007, 12:14:08 PM (19 years ago)

Author:

magnier

Message:

converting the PSF model fits to the polariation terms (replacing Sx,Sy,Sxy)

File:

• trunk/psModules/src/objects/pmPSFtry.c (modified) (4 diffs)

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-26 01:20:29 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-27 22:14:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -352 +352 @@
         x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
         y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
+
+        // weight by the error on the source flux
+        if (dz) {
+            dz->data.F64[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0];
+        }
     }
 
@@ -360 +365 @@
     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
 
-    // skip the unfitted parameters (X, Y, Io, Sky)
+    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
+    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
+    // each source and fit this set of parameters.  These values are less tightly coupled, but
+    // are still inter-related.  The fitted values do a good job of constraining the major axis
+    // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
+    // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
+    // parameters, with the constraint that the minor axis must be greater than a minimum
+    // threshold.
+
+    // convert the measured source shape paramters to polarization terms
+    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *e1   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *e2   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < psfTry->sources->n; i++) {
+        pmSource *source = psfTry->sources->data[i];
+        if (source->modelEXT == NULL) continue;
+
+        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
+
+        e0->data.F64[i] = pol.e0;
+        e1->data.F64[i] = pol.e1;
+        e2->data.F64[i] = pol.e2;
+    }
+
+    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
+    // This way, the parameters masked by one of the fits will be applied to the others
+    for (int i = 0; i < stats->clipIter; i++) {
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
+    }
+
+    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
     for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (i == PM_PAR_SKY)
-            continue;
-        if (i == PM_PAR_I0)
-            continue;
-        if (i == PM_PAR_XPOS)
-            continue;
-        if (i == PM_PAR_YPOS)
-            continue;
+        switch (i) {
+          case PM_PAR_SKY:
+          case PM_PAR_I0:
+          case PM_PAR_XPOS:
+          case PM_PAR_YPOS:            
+          case PM_PAR_SXX:            
+          case PM_PAR_SYY:            
+          case PM_PAR_SXY:            
+            continue;
+          default:
+            break;
+        }
 
         // select the per-object fitted data for this parameter
@@ -376 +417 @@
             if (source->modelEXT == NULL)
                 continue;
-
             z->data.F64[j] = source->modelEXT->params->data.F32[i];
-            if (applyWeight) {
-                dz->data.F64[j] = source->modelEXT->dparams->data.F32[i];
-            }
-
-            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
-            if (i == PM_PAR_SXY) {
-                z->data.F64[j] = pmPSF_SXYfromModel (source->modelEXT->params->data.F32);
-            }
         }
 
         // fit the collection of measured parameters to the PSF 2D model
         // the mask is carried from previous steps and updated with this operation
-        // the weight is either the parameter error or NULL, depending on 'applyWeights'
+        // the weight is either the flux error or NULL, depending on 'applyWeights'
         if (!psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);