Changeset 5828 for trunk/psphot/src/psphotEnsemblePSF.c

Timestamp:

Dec 22, 2005, 4:20:27 PM (21 years ago)

Author:

eugene

Message:

various cleanups, fixes, additions...

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (5 diffs)

trunk/psphot/src/psphotEnsemblePSF.c

@@ -1 +1 @@
 # include "psphot.h"
+
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
+
+    int ix = x - image->col0;
+    int iy = y - image->row0;
+
+    psF32 *Fm = &image->data.F32[iy - 1][ix];
+    psF32 *Fo = &image->data.F32[iy + 0][ix];
+    psF32 *Fp = &image->data.F32[iy + 1][ix];
+
+    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
+    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
+    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
+    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
+    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
+    poly->mask[2][2] = 1;
+    poly->mask[1][2] = 1;
+    poly->mask[2][1] = 1;
+
+    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
+
+    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
+    poly->coeff[0][1] = (Fyp - Fym)/6.0;
+    
+    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
+    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
+    
+    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
+    
+    return (poly);
+}
+
+psPlane psImageBicubeMin (psPolynomial2D *poly) {
+
+    psPlane min;
+
+    min.xErr = min.yErr = 0;
+
+    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
+
+    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
+    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
+    return (min);
+}
 
 bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
@@ -27 +73 @@
     index->n = 0;
 
+    bool UseAnalysisRegion = false;
+    psRegion AnalysisRegion;
+    char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION");
+    if (status) {
+        UseAnalysisRegion = true;
+        AnalysisRegion = psRegionFromString (region);
+        psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
+                  AnalysisRegion.x0, AnalysisRegion.y0, 
+                  AnalysisRegion.x1, AnalysisRegion.y1);
+    }
+
     for (int i = 0; i < sources->n; i++) {
         pmSource *inSource = sources->data[i];
@@ -35 +92 @@
         if (inSource->type == PM_SOURCE_DEFECT) continue; 
         if (inSource->type == PM_SOURCE_SATURATED) continue;
+
+        if (UseAnalysisRegion) {
+            if (inSource->moments->x < AnalysisRegion.x0) continue;
+            if (inSource->moments->y < AnalysisRegion.y0) continue;
+            if (inSource->moments->x > AnalysisRegion.x1) continue;
+            if (inSource->moments->y > AnalysisRegion.y1) continue;
+        }
 
         pmSource *otSource = pmSourceAlloc ();
@@ -51 +115 @@
             modelFLT->params->data.F32[2] = inSource->moments->x;
             modelFLT->params->data.F32[3] = inSource->moments->y;
-        }
-        // XXX EAM : add option to peak-up on peak (for non-sat objects)
+        } else {
+            // peak-up on peak (for non-sat objects)
+
+            int ix = inSource->peak->x;
+            int iy = inSource->peak->y;
+
+            psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
+            psPlane min = psImageBicubeMin (bicube);
+
+            psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
+                     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
+            
+            // if min point is too deviant, keep the old value
+            if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+                modelFLT->params->data.F32[2] = min.x + ix;
+                modelFLT->params->data.F32[3] = min.y + iy;
+            }
+            psFree (bicube);
+        }
 
         // set PSF parameters for this model
@@ -149 +230 @@
         // subtract object
         pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
+        Fi->mode |= PM_SOURCE_SUBTRACTED;
+        Fi->mode |= PM_SOURCE_TEMPSUB;
     }
 
     // XXX EAM : need to free up many things here
 
-    psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f\n", psTimerMark ("psphot"));
     return true;
 }