Changeset 27840 for branches/simtest_nebulous_branches/psphot/src/psphotMakeResiduals.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotMakeResiduals.c (modified) (8 diffs)

branches/simtest_nebulous_branches/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/simtest_nebulous_branches/psphot/src/psphotMakeResiduals.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+# define RESIDUAL_SOFTENING 0.005 
 
 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal) {
@@ -31 +33 @@
 
     float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
+    PS_ASSERT (status, false);
+
+    float radiusMax = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.RADIUS");
     PS_ASSERT (status, false);
 
@@ -171 +176 @@
                 bool offImage = false;
                 if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
-                    // fprintf (stderr, "off image: %f %f : %f %f\n", ix, iy, flux, dflux);
                     // This pixel is off the image
                     offImage = true;
@@ -179 +183 @@
                 }
                 fluxes->data.F32[i] = flux;
-                dfluxes->data.F32[i] = dflux;
+                dfluxes->data.F32[i] = hypot(dflux, RESIDUAL_SOFTENING);
                 if (isnan(flux)) {
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
+                }
+                if (isnan(dflux)) {
                     fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
                 }
@@ -234 +241 @@
                 }
 
+                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+                if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                  continue;
+                }
+
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
@@ -248 +261 @@
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
-
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
-
             } else {
                 assert (SPATIAL_ORDER == 1);
+
+                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+                if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                  continue;
+                }
+
                 psImageInit(A, 0.0);
                 psVectorInit(B, 0.0);
@@ -275 +292 @@
 
                 if (!psMatrixGJSolve(A, B)) {
-                    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals",
-                            oy, ox);
-                    psFree(resid); resid = NULL;
-                    break;
+                    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                    psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox);
+                    continue;
                 }
 
@@ -286 +302 @@
 
                 float dRo = sqrt(A->data.F32[0][0]);
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n",
-                // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
-                }
-                //resid->variance->data.F32[oy][ox] = XXX;
+                  resid->Ro->data.F32[oy][ox] = 0.0;
+                  resid->Rx->data.F32[oy][ox] = 0.0;
+                  resid->Ry->data.F32[oy][ox] = 0.0;
+                }
             }
         }

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack