Changeset 19395 for trunk/pswarp/src/pswarpTransformReadout.c

Timestamp:

Sep 5, 2008, 12:37:59 PM (18 years ago)

Author:

Paul Price

Message:

Fixed a bug in the tile transformation that was resulting in transformed data being overwritten with untransformed data (a single character bug --- 'y' instead of 'x'!). Attempted to make transformations a bit more efficient by only transforming the tiles that overlap the image; doesn't seem to make much difference.

File:

• trunk/pswarp/src/pswarpTransformReadout.c (modified) (8 diffs)

trunk/pswarp/src/pswarpTransformReadout.c

@@ -1 @@
-# include "pswarp.h"
+#include "pswarp.h"
 
 // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
@@ -5 +5 @@
 {
     // XXX this implementation currently ignores the use of the region
-    psImage *region = NULL;
+    psImage *region = NULL;             // Region to transform
 
-    psTimerStart ("warp");
+    psTimerStart("warp");
 
     // Get warp parameters
-    bool mdok;
-    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX");
-    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY");
-    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE");
+    bool mdok;                          // Status of MD lookup
+    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); // Number of grid points in x
+    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); // Number of grid points in y
+    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
+                                                                   "INTERPOLATION.MODE"); // Mode for interp
 
     // load the recipe
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE);
     psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
 
@@ -23 +24 @@
     psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config);
     psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config);
-    psAssert (mdok, "MASK.INPUT was not defined");
+    psAssert(mdok, "MASK.INPUT was not defined");
 
-    int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
-    if (!mdok) nThreads = 0;
-
-    // Flux fraction for "poor"
-    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC");
+    int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
+    if (!mdok) {
+        nThreads = 0;
+    }
+    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor"
 
     // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
     // output coordinates to input coordinates
-    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
+    pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY);
 
     // XXX optionally modify the grid based on this result and force the maxError < XXX
-    double maxError = pswarpMapGridMaxError (grid);
-    psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);
+    double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid
+    psLogMsg("pswarp", 3, "maximum error using this grid sampling: %lf\n", maxError);
+
+    // Get range of interest
+    int xOutMin, xOutMax, yOutMin, yOutMax; // Output pixel range
+    pswarpMatchRange(&xOutMin, &yOutMin, &xOutMax, &yOutMax, input, output);
+
+    // Check the range of the grid coordinates
+#define CHECK_GRID_RANGE() { \
+        xGridMin = PS_MIN(xGridMin, xGrid); \
+        xGridMax = PS_MAX(xGridMax, xGrid); \
+        yGridMin = PS_MIN(yGridMin, yGrid); \
+        yGridMax = PS_MAX(yGridMax, yGrid); \
+    }
+
+    int xGridMin = nGridX, xGridMax = 0, yGridMin = nGridY, yGridMax = 0; // Grid range
+    int xGrid, yGrid;                   // Grid coordinates
+    pswarpMapGridSetGrid(grid, xOutMin, yOutMin, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMin, yOutMax, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMax, yOutMin, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
+    pswarpMapGridSetGrid(grid, xOutMax, yOutMax, &xGrid, &yGrid);
+    CHECK_GRID_RANGE();
 
     // Interpolation options : move these from the arguments to explicit assignments
-    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image, input->weight, input->mask, maskIn, NAN, NAN, maskBad, maskPoor, poorFrac);
+    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image,
+                                                                       input->weight, input->mask, maskIn,
+                                                                       NAN, NAN, maskBad, maskPoor, poorFrac);
 
     if (input->weight && !output->weight) {
@@ -51 +77 @@
     }
 
-    static int maskNum = 0;
-    {
-        psString name = NULL;
-        psStringAppend(&name, "mask_%d.fits", maskNum++);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
-        psFitsClose(fits);
-    }
-
-    // total number of good pixels across all tiles (summed below)
-    int goodPixels = 0;
-
     // create jobs and supply them to the threads
-    for (int gridY = 0; gridY < grid->nYpts; gridY++) {
-        for (int gridX = 0; gridX < grid->nXpts; gridX++) {
-
+    for (int gridY = yGridMin; gridY < yGridMax; gridY++) {
+        for (int gridX = xGridMin; gridX < xGridMax; gridX++) {
             pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
-
-            // these items are just views to the data; they are not freed with args
-            args->input = input;
-            args->output = output;
-            args->grid = grid;
-            args->interp = interp;
-            args->region = region;
+            args->input = psMemIncrRefCounter(input);
+            args->output = psMemIncrRefCounter(output);
+            args->grid = psMemIncrRefCounter(grid);
+            args->interp = psMemIncrRefCounter(interp);
+            args->region = psMemIncrRefCounter(region);
 
             args->gridX = gridX;
@@ -83 +93 @@
             // allocate a job
             psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
-
-            // construct the arguments for this job
-            // job is pswarpTransformTile (gridX, gridY);
-            psArrayAdd (job->args, 1, args);
-            // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);
-
-            // call: pswarpTransformTile (args);
-            if (!psThreadJobAddPending (job)) {
+            psArrayAdd(job->args, 1, args);
+            if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to warp image.");
                 return false;
@@ -109 +113 @@
     // we have only supplied one type of job, so we can assume the types here
     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
     while ((job = psThreadJobGetDone()) != NULL) {
         if (job->args->n < 1) {
@@ -116 +122 @@
             // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
             goodPixels += args->goodPixels;
+            xMin = PS_MIN(args->xMin, xMin);
+            xMax = PS_MAX(args->xMax, xMax);
+            yMin = PS_MIN(args->yMin, yMin);
+            yMax = PS_MAX(args->yMax, yMax);
         }
-        psFree (job);
+        psFree(job);
     }
     psFree(interp);
     psFree(grid);
 
+    if (xMin < xMax && yMin < yMax) {
+        psTrace("pswarp.transform", 1, "Bounds [%d:%d,%d:%d]\n", xMin, xMax, yMin, yMax);
+    } else {
+        psTrace("pswarp.transform", 1, "No overlap\n");
+    }
 
     // Store the variance factor and number of good pixels
     if (goodPixels > 0) {
         // Variance factor: large factor --> small scale
-        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp);
+        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0,
+                                                           input->image->numRows / 2.0 + input->image->row0,
+                                                           interp);
         psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
         if (vfItem) {
@@ -145 +162 @@
 
     if (goodPixels > 0) {
-        if (!pswarpTransformSources (output, input, config)) {
+        if (!pswarpTransformSources(output, input, config)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
             return false;
         }
+
+        // Data is only written out if there are good pixels
+        output->data_exists = true;
+        output->parent->data_exists = true;
+        output->parent->parent->data_exists = true;
     }
 
-    // XXX should we not write anything out if there are no good pixels?
-    output->data_exists = true;
-    output->parent->data_exists = true;
-    output->parent->parent->data_exists = true;
-
-    {
-        psString name = NULL;
-        psStringAppend(&name, "mask_%d.fits", maskNum++);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
-        psFitsClose(fits);
-    }
-
-    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
+    psLogMsg("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
 
     return true;