Changeset 27933 for branches/pap/pswarp/src/pswarpLoop.c

Timestamp:

May 12, 2010, 12:56:07 PM (16 years ago)

Author:

Paul Price

Message:

Implementing better Jacobian correction. Previously, was calculating Jacobian from the projections (tangent plane --> sky), which does not account for scales in the transformations from the pixels to the tangent plane. Now using the local linear mapping to provide the Jacobian, so that it's applied on small areas. Also using the (square root of the) Jacobian to fix the plate scale of the covariance (instead of a single value for the entire focal plane) so the noise model should be more accurate as well.

File:

• branches/pap/pswarp/src/pswarpLoop.c (modified) (1 diff)

branches/pap/pswarp/src/pswarpLoop.c

@@ -273 +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);
-    }
-
-    // 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));
-        }
-        psKernel *covar = psImageCovarianceScale(output->covariance, outScale / inScale); // Scaled covariance
-        psFree(output->covariance);
-        output->covariance = covar;
     }