Changeset 28006 for trunk/pswarp/src

Timestamp:

May 18, 2010, 2:34:09 PM (16 years ago)

Author:

eugene

Message:

merging changes from PAP branch with modifications to the covariance matrix calculation

Location:

trunk

Files:

• . (modified) (1 prop)
• pswarp/src/pswarp.h (modified) (2 diffs)
• pswarp/src/pswarpLoop.c (modified) (2 diffs)
• pswarp/src/pswarpTransformReadout.c (modified) (3 diffs)
• pswarp/src/pswarpTransformTile.c (modified) (7 diffs)

trunk/pswarp/src/pswarp.h

@@ -29 +29 @@
 #define PSASTRO_RECIPE "PSASTRO" ///< Name of the recipe to use
 
-#define PSWARP_ANALYSIS_VARFACTOR "PSWARP.VARFACTOR" ///< Name for variance factor in analysis metadata
 #define PSWARP_ANALYSIS_GOODPIX   "PSWARP.GOODPIX" ///< Name for number of good pixels in analysis metadata
 #define PSWARP_ANALYSIS_COVARIANCES "PSWARP.COVARIANCES" ///< Name for covariance matrices on analysis MD
+#define PSWARP_ANALYSIS_JACOBIAN "PSWARP.JACOBIAN" ///< Name for Jacobian on analysis MD
 
 /**
@@ -70 +70 @@
     int xMin, xMax, yMin, yMax;         ///< Bounds of tile
     psKernel *covariance;               ///< Covariance matrix
+    double jacobian;                    ///< (Square root of) local Jacobian
 } pswarpTransformTileArgs;
 

trunk/pswarp/src/pswarpLoop.c

@@ -267 +267 @@
     }
 
-    // Correct image for change in the plate scale
-    {
-        psAssert(input && input->fpa && input->fpa->toSky, "Require astrometry for input");
-        psAssert(outFPA && outFPA && outFPA->toSky, "Require astrometry for output");
-
-        double inScale = input->fpa->toSky->Xs + input->fpa->toSky->Ys; // Plate scale for input
-        double outScale = outFPA->toSky->Xs + outFPA->toSky->Ys; // Plate scale for output
-        float correction = PS_SQR(outScale / inScale); // Correction factor to apply to image
-        psLogMsg("pswarp", PS_LOG_INFO, "Correcting flux by %f to account for pixel scales", correction);
-
-        psBinaryOp(output->image, output->image, "*", psScalarAlloc(correction, PS_TYPE_F32));
-        if (output->variance) {
-            psBinaryOp(output->variance, output->variance, "*",
-                       psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
-        }
-    }
-
-
     // Set covariance matrix for output
     {
@@ -291 +273 @@
         psAssert(covariances, "Should be there");
         psArray *covars = psListToArray(covariances); // Array of covariance matrices
-        output->covariance = psImageCovarianceAverage(covars);
+        psKernel *covar = psImageCovarianceAverage(covars);
         psFree(covars);
         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
+
+        // Correct covariance matrix scale for the mean (square root of the) Jacobian
+        double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+        int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
+        jacobian /= goodPixels;
+        output->covariance = psImageCovarianceScale(covar, jacobian);
+        psFree(covar);
+
         psImageCovarianceTransfer(output->variance, output->covariance);
     }

trunk/pswarp/src/pswarpTransformReadout.c

@@ -148 +148 @@
     psThreadJob *job = NULL;
     int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
-    int goodPixels = 0;                 // total number of good pixels across all tiles
+    int goodPixels = psMetadataLookupS32(&mdok, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Number of pixels
     psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
                                               PSWARP_ANALYSIS_COVARIANCES); // Collection of covar. matrices
@@ -157 +157 @@
         psFree(covariances);            // Drop reference; still have the copy on the analysis metadata
     }
+    double jacobian = psMetadataLookupF64(&mdok, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+    if (!isfinite(jacobian)) {
+        jacobian = 0.0;
+    }
+
     while ((job = psThreadJobGetDone()) != NULL) {
         if (job->args->n < 1) {
@@ -171 +176 @@
                 psListAdd(covariances, PS_LIST_TAIL, args->covariance);
             }
+            if (args->goodPixels > 0 && isfinite(args->jacobian)) {
+                jacobian += args->jacobian * args->goodPixels;
+            }
         }
         psFree(job);
     }
     psFree(grid);
+
+    psMetadataAddS32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, PS_META_REPLACE,
+                     "Number of good pixels", goodPixels);
+    psMetadataAddF64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_JACOBIAN, PS_META_REPLACE,
+                     "Jacobian of transformation", jacobian);
 
     if (xMin < xMax && yMin < yMax) {

trunk/pswarp/src/pswarpTransformTile.c

@@ -44 +44 @@
     args->yMax = PS_MIN_S32;
     args->covariance = NULL;
+    args->jacobian = NAN;
 
     return args;
@@ -50 +51 @@
 bool pswarpTransformTile(pswarpTransformTileArgs *args)
 {
-    psImage *inImage = args->input->image; ///< Input image
-    psImage *outImage = args->output->image; ///< Output image
+    pmReadout *input = args->input;        // Input readout
+    psImage *inImage = input->image, *inMask = input->mask; // Input images
+    pmReadout *output = args->output;      // Output readout
+    psImage *outImage = output->image, *outVariance = output->variance, *outMask = output->mask; // Outputs
 
     int inNumCols = inImage->numCols, inNumRows = inImage->numRows; ///< Size of input image
@@ -61 +64 @@
 
     // Dereference images for convenience
-    psF32 **outImageData     = args->output->image->data.F32;
-    psF32 **outVarData       = (args->output->variance) ? args->output->variance->data.F32 : NULL;
-    psImageMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
-    psImageMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psF32 **outImageData     = outImage->data.F32;
+    psF32 **outVarData       = outVariance ? outVariance->data.F32 : NULL;
+    psImageMaskType **outMaskData = outMask ? outMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+    psImageMaskType **inMaskData  = inMask ? inMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
 
     pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; ///< Map for this tile
@@ -74 +77 @@
     int yMin = PS_MAX(minPt.y, 0);
     int yMax = PS_MIN(maxPt.y, outNumRows);
+
+    double jacobian = map->Xx * map->Yy - map->Yx * map->Xy; // Jacobian of transformation
+    double jacobian2 = PS_SQR(jacobian);                     // Square Jacobian
 
     // Iterate over the output image pixels (parent frame)
@@ -87 +93 @@
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
-            double xInRaw, yInRaw;      // Input raw pixel coordinates
-            pswarpMapApply(&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
-            double xIn = xInRaw - outCol0, yIn = yInRaw - outRow0; // Position on input image
+            double xIn, yIn;            // Input pixel coordinates
+            pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
             if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
                 continue;
@@ -105 +110 @@
 
             if (outImageData) {
-                outImageData[yOut][xOut] = imageValue;
+                outImageData[yOut][xOut] = imageValue * jacobian;
             }
             if (outVarData) {
-                outVarData[yOut][xOut] = varValue;
+                outVarData[yOut][xOut] = varValue * jacobian2;
             }
             if (outMaskData) {
@@ -122 +127 @@
         double xIn, yIn;                // Position of interest on input
         pswarpMapApply(&xIn, &yIn, map, xOut + 0.5, yOut + 0.5);
-        // XXX Why are we subtracting the *output* col0,row0 from the *input* coordinates?
-        // I expect these are zero, so that it makes no difference, but it's distracting.
-        xIn -= outCol0;
-        yIn -= outRow0;
         psKernel *kernel = psImageInterpolationKernel(xIn, yIn, args->interp->mode); // Interpolation kernel
         args->covariance = psImageCovarianceCalculate(kernel, args->input->covariance);
+        args->jacobian = sqrt(jacobian);
         psFree(kernel);
     }