Changeset 18753 for branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout.c

Timestamp:

Jul 27, 2008, 12:41:58 PM (18 years ago)

Author:

eugene

Message:

turn on threading; add -threads argument, split out pswarpTransformSources

File:

• branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout.c (modified) (1 diff)

branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout.c

@@ -1 +1 @@
 # include "pswarp.h"
 
-bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) {
-
+// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
+bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config)
+{
     // XXX this implementation currently ignores the use of the region
     psImage *region = NULL;
-    pmCell *cell = NULL;
 
-    // select the current recipe
-    // psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    psTimerStart ("warp");
 
-    int outNx = output->image->numCols;
-    int outNy = output->image->numRows;
+    // 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");
 
-    psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
-    psPlane *outPix = psPlaneAlloc();   // Coordinates on the output detector
+    // load the recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
 
-    psPlane *FP = psPlaneAlloc();       // Coordinates on the focal plane
-    psPlane *TP = psPlaneAlloc();       // Coordinates on the tangent plane
-    psSphere *sky = psSphereAlloc();    // Coordinates on the sky
+    // output mask bits
+    psMaskType maskIn   = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); 
+    psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); 
+    psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config); 
+    psAssert (mdok, "MASK.INPUT was not defined");
 
-    cell = input->parent;
-    pmChip *chipInput = cell->parent;
-    pmFPA *fpaInput = chipInput->parent;
+    int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
+    if (!mdok) nThreads = 0;
 
-    cell = output->parent;
-    pmChip *chipOutput = cell->parent;
-    pmFPA *fpaOutput = chipOutput->parent;
+    // Flux fraction for "poor"
+    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); 
 
-    psF32 **outData = output->image->data.F32;
-    psImage *inImage = input->image;
+    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
+    // output coordinates to input coordinates
+    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
 
-    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage,
-                                                                       NULL, NULL, 0, NAN, NAN, 0, 0, 0.0);
+    // 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);
 
-    // Iterate over the output image pixels
-    for (int y = 0; y < outNy; y++) {
-        for (int x = 0; x < outNx; x++) {
-            // Only transform those pixels requested
-            if (region && region->data.U8[y][x]) continue;
+    // 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);
 
-            // XXX double check this 1/2 pixel offset
-            outPix->x = (double)x + 0.5;
-            outPix->y = (double)y + 0.5;
+    if (input->weight && !output->weight) {
+        output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
+        psImageInit(output->weight, NAN);
+    }
+    if ((input->mask || maskPoor || maskBad) && !output->mask) {
+        output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
+        psImageInit(output->mask, maskBad);
+    }
 
-            psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
-            psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
-            psDeproject (sky, TP, fpaOutput->toSky);
+    // total number of good pixels across all tiles (summed below)
+    int goodPixels = 0;
 
-            psProject (TP, sky, fpaInput->toSky);
-            psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
-            psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
+    // create jobs and supply them to the threads
+    for (int gridY = 0; gridY < grid->nYpts; gridY++) {
+        for (int gridX = 0; gridX < grid->nXpts; gridX++) {
 
-            // XXX get interpolation method from the recipe
-            double value;
-            if (!psImageInterpolate(&value, NULL, NULL, inPix->x, inPix->y, interp)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
-                psFree(interp);
-                psFree(inPix);
-                psFree(outPix);
-                psFree(FP);
-                psFree(TP);
-                psFree(sky);
-                return false;
-            }
+            pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
 
-            outData[y][x] = value;
-            // modify zero and scale?
+            // 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->gridX = gridX;
+            args->gridY = gridY;
+            args->goodPixels = 0;
+
+            if (nThreads) {
+                // allocate a job
+                psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE", 0);
+
+                // 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);
+                psThreadJobAddPending (job);
+            } else {
+                pswarpTransformTile (args);
+                goodPixels += args->goodPixels;
+            }
+            psFree (args);
+        }
+    }
+
+    // wait for the threads to finish and manage results
+    if (nThreads) {
+        // wait here for the threaded jobs to finish
+        if (!psThreadPoolWait ()) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+            return false;
+        }
+        fprintf (stderr, "success for threaded jobs\n");
+
+        // each job records its own goodPixel values; sum them here
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                pswarpTransformTileArgs *args = job->args->data[0];
+                // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
+                goodPixels += args->goodPixels;
+            }
+            psFree (job);
+        }
+    }
+    psFree(interp);
+    psFree(grid);
+
+
+    // 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); 
+        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
+        if (vfItem) {
+            psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX);
+            psAssert(goodpixItem, "It should be where we left it!");
+            psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64,
+                     "Should be the type we said.");
+
+            vfItem->data.F32 += varFactor * goodPixels;
+            goodpixItem->data.S64 += goodPixels;
+        } else {
+            psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0,
+                             "Variance factor weighted by the good pixels", varFactor * goodPixels);
+            psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0,
+                             "Number of good pixels", goodPixels);
         }
     }
-    psFree(interp);
-    psFree(inPix);
-    psFree(outPix);
-    psFree(FP);
-    psFree(TP);
-    psFree(sky);
+
+    if (goodPixels > 0) {
+        if (!pswarpTransformSources (output, input, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+            return false;
+        }
+    }
+
+    // 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;
+
+    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
 
     return true;