IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19395


Ignore:
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.

Location:
trunk/pswarp/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/pswarp/src/pswarp.h

    r18839 r19395  
    3535    int nXpts, nYpts;                   // number of x,y samples in the grid
    3636    int nXpix, nYpix;                   // x,y spacing in src image pixels of grid samples
    37     double xMin,  yMin;                 // coordinate of first grid sample
     37    double xMin,  yMin;                 // coordinate of first grid sample
    3838} pswarpMapGrid;
    3939
    … …  
    5151
    5252    // output values for this tile
    53     long goodPixels;
     53    long goodPixels;                    // Number of good pixels
     54    int xMin, xMax, yMin, yMax;         // Bounds of tile
    5455} pswarpTransformTileArgs;
    5556

  • trunk/pswarp/src/pswarpTransformReadout.c

    r19391 r19395  
    1 # include "pswarp.h"
     1#include "pswarp.h"
    22
    33// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
    … …  
    55{
    66    // XXX this implementation currently ignores the use of the region
    7     psImage *region = NULL;
     7    psImage *region = NULL;             // Region to transform
    88
    9     psTimerStart ("warp");
     9    psTimerStart("warp");
    1010
    1111    // Get warp parameters
    12     bool mdok;
    13     int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX");
    14     int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY");
    15     psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE");
     12    bool mdok;                          // Status of MD lookup
     13    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); // Number of grid points in x
     14    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); // Number of grid points in y
     15    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
     16                                                                   "INTERPOLATION.MODE"); // Mode for interp
    1617
    1718    // load the recipe
    18     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     19    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE);
    1920    psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
    2021
    … …  
    2324    psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config);
    2425    psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config);
    25     psAssert (mdok, "MASK.INPUT was not defined");
     26    psAssert(mdok, "MASK.INPUT was not defined");
    2627
    27     int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
    28     if (!mdok) nThreads = 0;
    29 
    30     // Flux fraction for "poor"
    31     float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC");
     28    int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
     29    if (!mdok) {
     30        nThreads = 0;
     31    }
     32    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor"
    3233
    3334    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
    3435    // output coordinates to input coordinates
    35     pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
     36    pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY);
    3637
    3738    // XXX optionally modify the grid based on this result and force the maxError < XXX
    38     double maxError = pswarpMapGridMaxError (grid);
    39     psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);
     39    double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid
     40    psLogMsg("pswarp", 3, "maximum error using this grid sampling: %lf\n", maxError);
     41
     42    // Get range of interest
     43    int xOutMin, xOutMax, yOutMin, yOutMax; // Output pixel range
     44    pswarpMatchRange(&xOutMin, &yOutMin, &xOutMax, &yOutMax, input, output);
     45
     46    // Check the range of the grid coordinates
     47#define CHECK_GRID_RANGE() { \
     48        xGridMin = PS_MIN(xGridMin, xGrid); \
     49        xGridMax = PS_MAX(xGridMax, xGrid); \
     50        yGridMin = PS_MIN(yGridMin, yGrid); \
     51        yGridMax = PS_MAX(yGridMax, yGrid); \
     52    }
     53
     54    int xGridMin = nGridX, xGridMax = 0, yGridMin = nGridY, yGridMax = 0; // Grid range
     55    int xGrid, yGrid;                   // Grid coordinates
     56    pswarpMapGridSetGrid(grid, xOutMin, yOutMin, &xGrid, &yGrid);
     57    CHECK_GRID_RANGE();
     58    pswarpMapGridSetGrid(grid, xOutMin, yOutMax, &xGrid, &yGrid);
     59    CHECK_GRID_RANGE();
     60    pswarpMapGridSetGrid(grid, xOutMax, yOutMin, &xGrid, &yGrid);
     61    CHECK_GRID_RANGE();
     62    pswarpMapGridSetGrid(grid, xOutMax, yOutMax, &xGrid, &yGrid);
     63    CHECK_GRID_RANGE();
    4064
    4165    // Interpolation options : move these from the arguments to explicit assignments
    42     psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image, input->weight, input->mask, maskIn, NAN, NAN, maskBad, maskPoor, poorFrac);
     66    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image,
     67                                                                       input->weight, input->mask, maskIn,
     68                                                                       NAN, NAN, maskBad, maskPoor, poorFrac);
    4369
    4470    if (input->weight && !output->weight) {
    … …  
    5177    }
    5278
    53     static int maskNum = 0;
    54     {
    55         psString name = NULL;
    56         psStringAppend(&name, "mask_%d.fits", maskNum++);
    57         psFits *fits = psFitsOpen(name, "w");
    58         psFree(name);
    59         psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
    60         psFitsClose(fits);
    61     }
    62 
    63     // total number of good pixels across all tiles (summed below)
    64     int goodPixels = 0;
    65 
    6679    // create jobs and supply them to the threads
    67     for (int gridY = 0; gridY < grid->nYpts; gridY++) {
    68         for (int gridX = 0; gridX < grid->nXpts; gridX++) {
    69 
     80    for (int gridY = yGridMin; gridY < yGridMax; gridY++) {
     81        for (int gridX = xGridMin; gridX < xGridMax; gridX++) {
    7082            pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
    71 
    72             // these items are just views to the data; they are not freed with args
    73             args->input = input;
    74             args->output = output;
    75             args->grid = grid;
    76             args->interp = interp;
    77             args->region = region;
     83            args->input = psMemIncrRefCounter(input);
     84            args->output = psMemIncrRefCounter(output);
     85            args->grid = psMemIncrRefCounter(grid);
     86            args->interp = psMemIncrRefCounter(interp);
     87            args->region = psMemIncrRefCounter(region);
    7888
    7989            args->gridX = gridX;
    … …  
    8393            // allocate a job
    8494            psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
    85 
    86             // construct the arguments for this job
    87             // job is pswarpTransformTile (gridX, gridY);
    88             psArrayAdd (job->args, 1, args);
    89             // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);
    90 
    91             // call: pswarpTransformTile (args);
    92             if (!psThreadJobAddPending (job)) {
     95            psArrayAdd(job->args, 1, args);
     96            if (!psThreadJobAddPending(job)) {
    9397                psError(PS_ERR_UNKNOWN, false, "Unable to warp image.");
    9498                return false;
    … …  
    109113    // we have only supplied one type of job, so we can assume the types here
    110114    psThreadJob *job = NULL;
     115    int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
     116    int goodPixels = 0;                 // total number of good pixels across all tiles
    111117    while ((job = psThreadJobGetDone()) != NULL) {
    112118        if (job->args->n < 1) {
    … …  
    116122            // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
    117123            goodPixels += args->goodPixels;
     124            xMin = PS_MIN(args->xMin, xMin);
     125            xMax = PS_MAX(args->xMax, xMax);
     126            yMin = PS_MIN(args->yMin, yMin);
     127            yMax = PS_MAX(args->yMax, yMax);
    118128        }
    119         psFree (job);
     129        psFree(job);
    120130    }
    121131    psFree(interp);
    122132    psFree(grid);
    123133
     134    if (xMin < xMax && yMin < yMax) {
     135        psTrace("pswarp.transform", 1, "Bounds [%d:%d,%d:%d]\n", xMin, xMax, yMin, yMax);
     136    } else {
     137        psTrace("pswarp.transform", 1, "No overlap\n");
     138    }
    124139
    125140    // Store the variance factor and number of good pixels
    126141    if (goodPixels > 0) {
    127142        // Variance factor: large factor --> small scale
    128         float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp);
     143        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0,
     144                                                           input->image->numRows / 2.0 + input->image->row0,
     145                                                           interp);
    129146        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
    130147        if (vfItem) {
    … …  
    145162
    146163    if (goodPixels > 0) {
    147         if (!pswarpTransformSources (output, input, config)) {
     164        if (!pswarpTransformSources(output, input, config)) {
    148165            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    149166            return false;
    150167        }
     168
     169        // Data is only written out if there are good pixels
     170        output->data_exists = true;
     171        output->parent->data_exists = true;
     172        output->parent->parent->data_exists = true;
    151173    }
    152174
    153     // XXX should we not write anything out if there are no good pixels?
    154     output->data_exists = true;
    155     output->parent->data_exists = true;
    156     output->parent->parent->data_exists = true;
    157 
    158     {
    159         psString name = NULL;
    160         psStringAppend(&name, "mask_%d.fits", maskNum++);
    161         psFits *fits = psFitsOpen(name, "w");
    162         psFree(name);
    163         psFitsWriteImage(fits, NULL, output->mask, 0, NULL);
    164         psFitsClose(fits);
    165     }
    166 
    167     psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
     175    psLogMsg("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
    168176
    169177    return true;

  • trunk/pswarp/src/pswarpTransformTile.c

    r19149 r19395  
    1 # include "pswarp.h"
     1#include "pswarp.h"
    22
    3 void pswarpTransformTileArgsFree (pswarpTransformTileArgs *args) {
    4     if (!args) return;
     3static void transformTileArgsFree(pswarpTransformTileArgs *args)
     4{
     5    psFree(args->input);
     6    psFree(args->output);
     7    psFree(args->grid);
     8    psFree(args->interp);
     9    psFree(args->region);
    510    return;
    611}
    712
    8 pswarpTransformTileArgs *pswarpTransformTileArgsAlloc() {
    9 
    10     pswarpTransformTileArgs *args = (pswarpTransformTileArgs *)psAlloc(sizeof(pswarpTransformTileArgs));
    11     psMemSetDeallocator(args, (psFreeFunc)pswarpTransformTileArgsFree);
     13pswarpTransformTileArgs *pswarpTransformTileArgsAlloc()
     14{
     15    pswarpTransformTileArgs *args = psAlloc(sizeof(pswarpTransformTileArgs));
     16    psMemSetDeallocator(args, (psFreeFunc)transformTileArgsFree);
    1217
    1318    args->input = NULL;
    … …  
    2126
    2227    args->goodPixels = 0;
     28    args->xMin = PS_MAX_S32;
     29    args->xMax = PS_MIN_S32;
     30    args->yMin = PS_MAX_S32;
     31    args->yMax = PS_MIN_S32;
    2332
    2433    return args;
    2534}
    2635
    27 bool pswarpTransformTile (pswarpTransformTileArgs *args) {
     36bool pswarpTransformTile(pswarpTransformTileArgs *args)
     37{
     38    psImage *inImage = args->input->image; // Input image
     39    psImage *outImage = args->output->image; // Output image
    2840
    29     // int inCol0 = args->input->image->col0;
    30     // int inRow0 = args->input->image->row0;
    31     int inNcol = args->input->image->numCols;
    32     int inNrow = args->input->image->numRows;
     41    int inNumCols = inImage->numCols, inNumRows = inImage->numRows; // Size of input image
     42    int outNumCols = outImage->numCols, outNumRows = outImage->numRows; // Size of output image
     43    int outCol0 = outImage->col0, outRow0 = outImage->row0; // Offset of output image
    3344
    34     int outCol0 = args->output->image->col0;
    35     int outRow0 = args->output->image->row0;
    36     int outNcol = args->output->image->numCols;
    37     int outNrow = args->output->image->numRows;
     45    psPlane minPt, maxPt;               // Minimum and maximum points for this tile
     46    pswarpMapGridCoordRange(args->grid, args->gridX, args->gridY, &minPt, &maxPt);
    3847
    39     // get the coordinate range for this grid tile
    40     psPlane minPt, maxPt;
    41     pswarpMapGridCoordRange (args->grid, args->gridX, args->gridY, &minPt, &maxPt);
    42 
    43     psF32 **outImageData     = (args->output->image)  ? args->output->image->data.F32 : NULL;
     48    // Dereference images for convenience
     49    psF32 **outImageData     = args->output->image->data.F32;
    4450    psF32 **outVarData       = (args->output->weight) ? args->output->weight->data.F32 : NULL;
    4551    psMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_MASK_DATA : NULL;
    4652    psMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_MASK_DATA : NULL;
    4753
    48     pswarpMap *map = args->grid->maps[args->gridX][args->gridY];
    49 
    50     psImage *region = args->region;
    51 
    52     double xInRaw, yInRaw;
    53 
    54     // output values for this pixel
    55     double imageValue;
    56     double varValue;
    57     psMaskType maskValue;
     54    pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; // Map for this tile
     55    psImage *region = args->region;     // Region to transform
    5856
    5957    // Bounds for iteration
    6058    int xMin = PS_MAX(minPt.x, 0);
    61     int xMax = PS_MIN(maxPt.x, outNcol);
     59    int xMax = PS_MIN(maxPt.x, outNumCols);
    6260    int yMin = PS_MAX(minPt.y, 0);
    63     int yMax = PS_MIN(maxPt.y, outNrow);
    64 
     61    int yMax = PS_MIN(maxPt.y, outNumRows);
    6562
    6663    // Iterate over the output image pixels (parent frame)
    … …  
    7067
    7168            // Only transform those pixels requested
    72             if (region && region->data.U8[y][x]) continue;
     69            if (region && region->data.U8[y][x]) {
     70                continue;
     71            }
    7372
    7473            // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
    7574            // both are in the parent frames of the input and output images.
    76             pswarpMapApply (&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
    77 
    78             double xIn = xInRaw - outCol0;      // Position on input image
    79             double yIn = yInRaw - outRow0;      // Position on input image
    80 
    81             if (xIn < 0) continue;
    82             if (yIn >= inNcol) continue;
    83             if (yIn < 0) continue;
    84             if (yIn >= inNrow) continue;
    85 
    86             goodPixels++;
     75            double xInRaw, yInRaw;      // Input raw pixel coordinates
     76            pswarpMapApply(&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
     77            double xIn = xInRaw - outCol0, yIn = yInRaw - outRow0; // Position on input image
     78            if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
     79                continue;
     80            }
    8781
    8882            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
    89             maskValue = inMaskData ? inMaskData[(int)yIn][(int)xIn] : 0;
     83            double imageValue, varValue; // Value of image and variance map
     84            psMaskType maskValue = inMaskData ? inMaskData[(int)yIn][(int)xIn] : 0; // Value of mask
    9085            if (!psImageInterpolate(&imageValue, &varValue, &maskValue, xIn, yIn, args->interp)) {
    9186                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    … …  
    9388            }
    9489
    95             int xOut = x - outCol0;     // Position on output image
    96             int yOut = y - outRow0;     // Position on output image
     90            int xOut = x - outCol0, yOut = y - outRow0; // Position on output image
    9791
    98             // not all images need be transformed
    9992            if (outImageData) {
    10093                outImageData[yOut][xOut] = imageValue;
    … …  
    10699                outMaskData[yOut][xOut] = maskValue;
    107100            }
     101
     102            goodPixels++;
    108103        }
    109104    }
     105
    110106    args->goodPixels = goodPixels;
     107    args->xMin = xMin;
     108    args->xMax = xMax;
     109    args->yMin = yMin;
     110    args->yMax = yMax;
     111
    111112    return true;
    112113}
Note: See TracChangeset for help on using the changeset viewer.