Changeset 5654 for trunk/psphot/src/psphotSparseMatrix.c

Timestamp:

Dec 1, 2005, 9:55:16 AM (21 years ago)

Author:

eugene

Message:

adding ensemble photometry

File:

• trunk/psphot/src/psphotSparseMatrix.c (modified) (4 diffs)

trunk/psphot/src/psphotSparseMatrix.c

@@ -1 +1 @@
 # include "psphot.h"
 
-psSparse *psphotStarOverlaps (psArray *sources) {
+bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    psVector *params = model->params;
+    bool  status;
+    float x;
+    float y;
+    float Sky;
+    int   Nfit = 0;
 
-    // source analysis is done in S/N order (brightest first)
+    psTimerStart ("psphot");
+
+    // source analysis is done in spatial order
     sources = psArraySort (sources, psphotSortByY);
 
-    // fill out the models for each object
-    // XXX need to assign the model parameters using the PSF function, with unit peak
+    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
+    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
+
+    // set the object surface-brightness limit for fitted pixels
+    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
+    psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
+
+    // vector for model input coordinates
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+
+    // this function specifies the radius at this the model hits the given flux
+    pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
+
+    // this is the function to get the model function
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    // pre-calculate all model pixels
     psArray  *pixlist = psArrayAlloc (sources->n);
     for (int k = 0; k < sources->n; k++) {
-        src = sources->data[k];
+        source = sources->data[k];
+
+        // use the source moments, etc to guess basic model parameters
+        pmModel *modelFLT = pmSourceModelGuess (source, psf->type); 
+
+        // set PSF parameters for this model
+        pmModel *model = pmModelFromPSF (modelFLT, psf);
+        // model->params->data.F32[2] = source->peak->x;
+        // model->params->data.F32[3] = source->peak->y;
+        model->params->data.F32[1] = 1.0;
+        x = model->params->data.F32[2];
+        y = model->params->data.F32[3];
+        psFree (modelFLT);
+
+        // set the fit radius based on the object flux limit and the model
+        // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
+        model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
+        if (isnan(model->radius)) {
+          psAbort ("apply_psf_model", "error in radius");
+        }
+
+        // build the model image 
         nx = src->pixels->numCols;
         ny = src->pixels->numRows;
         im = psImageAlloc (nx, ny, PS_DATA_F32);
+        mk = psImageAlloc (nx, ny, PS_DATA_U8);
+
+        // fill in the model pixel values
+        psImageKeepCircle (mk, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+        psVector *params = model->params;
         for (int i = 0; i < nx; i++) {
             for (int j = 0; j < ny; j++) {
+                if (mk->data.U8[j][i]) continue;
                 x->data.F32[0] = (float) i + pix->col0;
                 x->data.F32[1] = (float) j + pix->row0;
@@ -28 +75 @@
     }
     
-    psVector *Aij = psVectorAlloc (100, PS_DATA_F32);
-    psVector *Bfj = psVectorAlloc (100, PS_DATA_F32);
-    psVector *Si  = psVectorAlloc (100, PS_DATA_S32);
-    psVector *Sj  = psVectorAlloc (100, PS_DATA_S32);
+    // fill out the sparse matrix
+    psSparse *sparse = psSparseAlloc (sources->n, 100);
+    psVector *Aij = sparse->Aij;
+    psVector *Bfj = sparse->Bfj;
+    psVector *Si  = sparse->Si;
+    psVector *Sj  = sparse->Sj;
 
     int Ri, Rj;
@@ -42 +91 @@
         // diagonal element (auto-cross-product)
         f = psphotCrossProduct (pi, pi);
-        Aij->data.F32[Nelem] = f;
-        Si->data.S32[Nelem] = i;
-        Sj->data.S32[Nelem] = j;
-        Nelem ++;
-
-        Qii->data.F32[i] = f;
+        psSparseMatrixElement (sparse, i, i, f);
 
         // find the image x model value
         f = psphotCrossProduct (pi, pf);
-        Bfj->data.F32[i] = f;
+        psSparseVectorElement (sparse, i, f);
 
         // loop over all other stars following this one
@@ -64 +108 @@
             // got an overlap; calculate cross-product and add to output array
             f = psphotCrossProduct (pi, pj);
-            
-            Aij->data.F32[Nelem] = f;
-            Si->data.S32[Nelem] = i;
-            Sj->data.S32[Nelem] = j;
-            Nelem ++;
-
-            Aij->data.F32[Nelem] = f;
-            Si->data.S32[Nelem] = j;
-            Sj->data.S32[Nelem] = i;
-            Nelem ++;
+            psSparseMatrixElement (sparse, j, i, f); 
         }
     }
-
-    index = psVectorSortIndex (NULL, Sj);
-
-    psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32);
-    psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
-    psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
-    for (i = 0; i < Nelem; i++) {
-        j = index->data.U32[i];
-        tAij->data.F32[i] = Aij->data.F32[j];
-        tSi->data.S32[i]  = Si->data.S32[j];
-        tSj->data.S32[i]  = Sj->data.S32[j];
-    }
-    psFree (Aij);
-    psFree (Si);
-    psFree (Sj);
-
-    sparse->Aij = Aij;
-    sparse->Bfj = Bfj;
-    sparse->Qii = Qii;
-    sparse->Si = Si;
-    sparse->Sj = Sj;
-    sparse->Nelem = Nelem;
-    sparse->Nrows = sources->n;
-
+    
+    psSparseResort (sparse);
     return (match);
 }