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

Timestamp:

Feb 5, 2009, 5:03:33 PM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:

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

trunk/psphot/src/psphotMakeResiduals.c

@@ -58 +58 @@
     // for each input source:
     // - construct a residual image, renormalized
-    // - construct a renormalized weight image
+    // - construct a renormalized variance image
     // - construct a new mask image
 
     // construct the output residual table (Nx*DX,Ny*DY)
     // for each output pixel:
-    // - construct a histogram of the values & weights (interpolate to the common pixel coordinate)
+    // - construct a histogram of the values & variances (interpolate to the common pixel coordinate)
     // - measure the robust median & sigma
     // - reject (mask) input pixels which are outliers
     // - re-measure the robust median & sigma
-    // - set output pixel, weight, and mask
+    // - set output pixel, variance, and mask
 
     // these mask values do not correspond to the recipe values: they
@@ -76 +76 @@
     // psVectorStats
 
-    const psImageMaskType badMask     = 0x01;   // mask bits
-    const psImageMaskType poorMask    = 0x02;   // from psImageInterpolate
-    const psImageMaskType clippedMask = 0x04;   // mask bit set for clipped values
+    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;
 
@@ -101 +101 @@
 
         psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
-        psImage *weight = psImageCopy (NULL, source->weight,   PS_TYPE_F32);
+        psImage *variance = psImageCopy (NULL, source->variance,   PS_TYPE_F32);
         psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK);
         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
 
-        // re-normalize image and weight
+        // re-normalize image and variance
         float Io = model->params->data.F32[PM_PAR_I0];
         psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
-        psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
-
-        // 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);
+        psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
+
+        // we interpolate the image and variance - include the mask or not?
+        // XXX why not the mask?
+        // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0);
+        psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0);
         psArrayAdd (input, 100, interp);
 
@@ -128 +128 @@
         psFree (mask);
         psFree (image);
-        psFree (weight);
+        psFree (variance);
         psFree (interp);
     }
@@ -176 +176 @@
                     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;
+                  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.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
@@ -204 +204 @@
 
             // mark input pixels which are more than N sigma from the median
-            int nKeep = 0;
+            int nKeep = 0;
             for (int i = 0; i < fluxes->n; i++) {
                 float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
@@ -210 +210 @@
                 float swing = fabs(delta) / sigma;
 
-                // mask pixels which are out of range
+                // mask pixels which are out of range
                 if (swing > nSigma) {
                     fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = clippedMask;
                 }
-                if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
+                if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
             }
 
@@ -225 +225 @@
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
-                //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
+                //resid->variance->data.F32[oy][ox] = fluxStats->sampleStdev;
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {
@@ -231 +231 @@
                 }
 
-                // 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]);
+                // 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 {
@@ -268 +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.PM_TYPE_RESID_MASK_DATA[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.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
-                //resid->weight->data.F32[oy][ox] = XXX;
+                //resid->variance->data.F32[oy][ox] = XXX;
             }
         }
@@ -301 +301 @@
       psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
       psphotSaveImage (NULL, resid->Ry,     "resid.ry.fits");
-      psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
+      psphotSaveImage (NULL, resid->variance, "resid.wt.fits");
       psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
     }