Changeset 21183 for trunk/psphot/src/psphotMakeResiduals.c

Timestamp:

Jan 26, 2009, 8:40:07 PM (17 years ago)

Author:

eugene

Message:

incorporating changes from 16bit mask upgrades (eam_branch_20081230)

File:

• trunk/psphot/src/psphotMakeResiduals.c (modified) (13 diffs)

trunk/psphot/src/psphotMakeResiduals.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal) {
+bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal) {
 
     bool status, isPSF;
     double flux, dflux;
-    psU8 mflux;
+    psImageMaskType mflux;
 
     psTimerStart ("psphot.residuals");
@@ -69 +69 @@
     // - set output pixel, weight, and mask
 
-    // XXX need to set these correctly!!
-    const int badMask = 1;              // mask bits
-    const int poorMask = 2;             //       from psImageInterpolate
-    const int clippedMask = 4;          // mask bit set for clipped values
+    // these mask values do not correspond to the recipe values: they
+    // are not propagated to images: they just need to fit in an 8-bit
+    // value.  they are supplied to psImageInterpolate, which takes a
+    // psImageMaskType; the mask portion of the result from
+    // psImageInterpolate is supplied to fmasks, which is then used by
+    // psVectorStats
+
+    const psImageMaskType badMask     = 0x01;   // mask bits
+    const psImageMaskType poorMask    = 0x02;   // from psImageInterpolate
+    const psImageMaskType clippedMask = 0x04;   // mask bit set for clipped values
+    const psVectorMaskType fmaskVal = badMask | poorMask | clippedMask;
 
     // determine the maximum image size from the input sources
@@ -95 +102 @@
         psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
         psImage *weight = psImageCopy (NULL, source->weight,   PS_TYPE_F32);
-        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_U8);
+        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK);
         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
 
@@ -103 +110 @@
         psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
 
-        // we will interpolate the image and weight - include the mask or not?
-        // XXX consider better values for the mask bits
-        psImageInterpolation *interp =
-            psImageInterpolationAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0);
-        psArrayAdd (input,  100, interp);
+        // we interpolate the image and weight - include the mask or not?
+        // XXX why not the mask?
+        // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0);
+        psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0);
+        psArrayAdd (input, 100, interp);
 
         // save the X,Y position for future reference
@@ -131 +138 @@
     psVector *fluxes  = psVectorAlloc (input->n, PS_TYPE_F32);
     psVector *dfluxes = psVectorAlloc (input->n, PS_TYPE_F32);
-    psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_U8);
+    psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK);
 
     // statistic to use to determine baseline for clipping
@@ -164 +171 @@
                 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;
-                    fmasks->data.U8[i] = 1;
-                    // fprintf (stderr, "off image: %f %f : %f %f\n", ix, iy, flux, dflux);
-                }
-                fluxes->data.F32[i] = flux;
-                dfluxes->data.F32[i] = dflux;
-                fmasks->data.U8[i] = mflux;
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
+                } else {
+                  fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = mflux; // XXX is mflux IMAGE or VECTOR type?
+                }
+                fluxes->data.F32[i] = flux;
+                dfluxes->data.F32[i] = dflux;
                 if (isnan(flux)) {
-                    fmasks->data.U8[i] = 1;
-                }
-                if (fmasks->data.U8[i] == 0) {
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
+                }
+                if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) {
                     nGoodPixel ++;
                 }
             }
 
-            // skip pixels which are off the image...
+            // skip pixels with insufficient data
             bool validPixel = (SPATIAL_ORDER == 0) ? (nGoodPixel > 1) : (nGoodPixel > 3);
             if (!validPixel) {
@@ -186 +194 @@
                 resid->Rx->data.F32[oy][ox] = 0.0;
                 resid->Ry->data.F32[oy][ox] = 0.0;
-                resid->mask->data.U8[oy][ox] = 1;
+                resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
                 continue;
             }
@@ -192 +200 @@
             // measure the robust median to determine a baseline reference value
             *fluxClip = *fluxClipDef;
-            psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
+            psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal);
             psErrorClear();             // clear (ignore) any outstanding errors
 
@@ -202 +210 @@
                 float swing = fabs(delta) / sigma;
 
-                // make this a user option
+                // mask pixels which are out of range
                 if (swing > nSigma) {
-                    fmasks->data.U8[i] = clippedMask;
-                }
-                if (!fmasks->data.U8[i]) nKeep++;
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = clippedMask;
+                }
+                if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
             }
 
@@ -212 +220 @@
                 // measure the desired statistic on the unclipped pixels
                 *fluxStats = *fluxStatsDef;
-                psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
+                psVectorStats (fluxStats, fluxes, NULL, fmasks, fmaskVal);
                 psErrorClear();         // clear (ignore) any outstanding errors
 
@@ -220 +228 @@
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {
-                  resid->mask->data.U8[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.U8[oy][ox]);
+                  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 {
@@ -230 +238 @@
                 psVectorInit(B, 0.0);
                 for (int i = 0; i < fluxes->n; i++) {
-                    if (fmasks->data.U8[i]) continue;
+                    if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
                     B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
                     B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
@@ -260 +268 @@
 
                 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.U8[oy][ox]);
+                // 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.U8[oy][ox] = 1;
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
                 //resid->weight->data.F32[oy][ox] = XXX;