Changeset 27933 for branches/pap/pswarp/src/pswarpTransformReadout.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/pswarpTransformReadout.c (modified) (3 diffs)

branches/pap/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) {