Changeset 28484 for branches/pap/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Jun 24, 2010, 2:59:09 PM (16 years ago)

Author:

Paul Price

Message:

Merging trunk in advance of reintegrating into trunk.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSourcePhotometry.c (modified) (14 diffs)

branches/pap/psModules/src/objects/pmSourcePhotometry.c

@@ -107 +107 @@
         // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
         double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
-        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
-        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
-        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
-        source->psfMag = -2.5*log10(source->psfFlux);
+        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
+        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
+        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
+        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
+        source->psfMag = -2.5*log10(source->psfFlux);
     } else {
         status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
-        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
+        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     }
 
@@ -120 +120 @@
     if (source->modelFits) {
         bool foundEXT = false;
-        for (int i = 0; i < source->modelFits->n; i++) {
-            pmModel *model = source->modelFits->data[i];
-            status = pmSourcePhotometryModel (&model->mag, NULL, model);
-            if (model == source->modelEXT) foundEXT = true;
-        }
-        if (foundEXT) {
-            source->extMag = source->modelEXT->mag;
-        } else {
-            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
-        }
+        for (int i = 0; i < source->modelFits->n; i++) {
+            pmModel *model = source->modelFits->data[i];
+            status = pmSourcePhotometryModel (&model->mag, NULL, model);
+            if (model == source->modelEXT) foundEXT = true;
+        }
+        if (foundEXT) {
+            source->extMag = source->modelEXT->mag;
+        } else {
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
+        }
     } else {
-        if (source->modelEXT) {
-            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
-        }
+        if (source->modelEXT) {
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
+        }
     }
 
@@ -204 +204 @@
         }
         if (mode & PM_SOURCE_PHOT_APCORR) {
-            // XXX this should be removed -- we no longer fit for the 'sky bias'
+            // XXX this should be removed -- we no longer fit for the 'sky bias'
             rflux   = pow (10.0, 0.4*source->psfMag);
             source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
@@ -236 +236 @@
     flux = model->modelFlux (model->params);
     if (flux > 0) {
-        mag = -2.5*log10(flux);
+        mag = -2.5*log10(flux);
     }
     if (fitMag) {
-        *fitMag = mag;
+        *fitMag = mag;
     }
     if (fitFlux) {
-        *fitFlux = flux;
+        *fitFlux = flux;
     }
 
@@ -380 +380 @@
 
     if (source->diffStats == NULL) {
-        source->diffStats = pmSourceDiffStatsAlloc();
+        source->diffStats = pmSourceDiffStatsAlloc();
     }
 
@@ -388 +388 @@
     int   nMask = 0;
     int   nBad  = 0;
-    
+
     psImage *flux     = source->pixels;
     psImage *variance = source->variance;
@@ -394 +394 @@
 
     for (int iy = 0; iy < flux->numRows; iy++) {
-        for (int ix = 0; ix < flux->numCols; ix++) {
+        for (int ix = 0; ix < flux->numCols; ix++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
-                nMask ++;
-                continue;
-            }
-
-            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
-
-            if (SN > +FLUX_LIMIT) { 
-                nGood ++;
-                fGood += fabs(flux->data.F32[iy][ix]);
-            }
-
-            if (SN < -FLUX_LIMIT) { 
-                nBad ++;
-                fBad += fabs(flux->data.F32[iy][ix]);
-            }
-        }
+                nMask ++;
+                continue;
+            }
+
+            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+
+            if (SN > +FLUX_LIMIT) {
+                nGood ++;
+                fGood += fabs(flux->data.F32[iy][ix]);
+            }
+
+            if (SN < -FLUX_LIMIT) {
+                nBad ++;
+                fBad += fabs(flux->data.F32[iy][ix]);
+            }
+        }
     }
 
     source->diffStats->nGood      = nGood;
-    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);          
-    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);          
-    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
+    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);
+    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);
+    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);
     source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
 
@@ -628 +628 @@
         }
     }
-    model->nPix = Npix; 
+    model->nPix = Npix;
     model->nDOF = Npix - 1;
     model->chisq = dC;
@@ -636 +636 @@
 
 
-double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
+double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -652 +652 @@
     for (int yi = 0; yi < Pi->numRows; yi++) {
         for (int xi = 0; xi < Pi->numCols; xi++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
                 continue;
             if (!unweighted_sum) {
@@ -684 +684 @@
 }
 
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -727 +727 @@
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
-                continue;
-            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+                continue;
+            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
                 continue;
 
@@ -746 +746 @@
 }
 
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -789 +789 @@
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi])
-                continue;
-            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj])
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+                continue;
+            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
                 continue;