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

Timestamp:

Jan 15, 2006, 8:30:16 AM (21 years ago)

Author:

eugene

Message:

API cleanup, removed old test files, some re-organization, changed FLT/GAL to EXT

File:

• trunk/psphot/src/psphotEnsemblePSF.c (modified) (7 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) { 
@@ -107 +61 @@
         }
 
-        // XXX EAM : add option to use FLT or PSF form
+        // XXX EAM : add option to use EXT or PSF form
         // use the source moments, etc to guess basic model parameters
-        pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); 
+        pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 
         if (inSource->mode &  PM_SOURCE_SATSTAR) {
-            modelFLT->params->data.F32[2] = inSource->moments->x;
-            modelFLT->params->data.F32[3] = inSource->moments->y;
+            modelEXT->params->data.F32[2] = inSource->moments->x;
+            modelEXT->params->data.F32[3] = inSource->moments->y;
         } else {
             // peak-up on peak (for non-sat objects)
@@ -123 +77 @@
 
             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);
+                     modelEXT->params->data.F32[2], modelEXT->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;
+                modelEXT->params->data.F32[2] = min.x + ix;
+                modelEXT->params->data.F32[3] = min.y + iy;
             }
             psFree (bicube);
@@ -134 +88 @@
 
         // set PSF parameters for this model
-        pmModel *model = pmModelFromPSF (modelFLT, psf);
-        psFree (modelFLT);
+        pmModel *model = pmModelFromPSF (modelEXT, psf);
+        psFree (modelEXT);
 
         // save the original coords
@@ -185 +139 @@
 
         // diagonal element (auto-cross-product)
-        f = psphotCrossProduct (Mi, Mi);
+        f = pmSourceCrossProduct (Mi, Mi);
         psSparseMatrixElement (sparse, i, i, f);
 
         // find the image x model value
-        f = psphotCrossProduct (Fi, Mi);
+        f = pmSourceCrossProduct (Fi, Mi);
         psSparseVectorElement (sparse, i, f);
 
@@ -203 +157 @@
             
             // got an overlap; calculate cross-product and add to output array
-            f = psphotCrossProduct (Mi, Mj);
+            f = pmSourceCrossProduct (Mi, Mj);
             psSparseMatrixElement (sparse, j, i, f); 
         }
@@ -235 +189 @@
     return true;
 }
-
-
-float psphotCrossProduct (pmSource *Mi, pmSource *Mj) {
-
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-    float flux;
-
-    psImage *Pi = Mi->pixels;
-    psImage *Pj = Mj->pixels;
-
-    psImage *Ti = Mi->mask;
-    psImage *Tj = Mj->mask;
-
-    Xs = PS_MAX (Pi->col0, Pj->col0);
-    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
-    
-    Ys = PS_MAX (Pi->row0, Pj->row0);
-    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
-    
-    xIs = Xs - Pi->col0;
-    xJs = Xs - Pj->col0;
-    yIs = Ys - Pi->row0;
-    yJs = Ys - Pj->row0;
-
-    xIe = Xe - Pi->col0;
-    yIe = Ye - Pi->row0;
-
-    flux = 0;
-    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
-        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.U8[yi][xi]) continue;
-            if (Tj->data.U8[yj][xj]) continue;
-            flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
-        }
-    }
-    return (flux);
-}