IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18753


Ignore:
Timestamp:
Jul 27, 2008, 12:41:58 PM (18 years ago)
Author:
eugene
Message:

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

Location:
branches/eam_branch_20080719/pswarp/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20080719/pswarp/src/Makefile.am

    r18752 r18753  
    1717        pswarpSetMaskBits.c             \
    1818        pswarpThreadLauncher.c          \
    19         pswarpTransformReadout_Opt.c    \
    20         pswarpTransformReadout_Threaded.c   \
    21         pswarpTransformReadout_Unthreaded.c \
     19        pswarpTransformReadout.c        \
     20        pswarpTransformSources.c \
    2221        pswarpTransformTile.c           \
    2322        pswarpVersion.c           

  • branches/eam_branch_20080719/pswarp/src/pswarp.c

    r18752 r18753  
    1717    psLibInit(NULL);
    1818
    19     // create the thread pool with 4 threads, supplying our thread launcher function
    20     // XXX need to determine the number of threads from the config data
    21     psThreadPoolInit (4, &pswarpThreadLauncher);
    22    
    2319    // model inits are needed in pmSourceIO
    2420    // models defined in psphot/src/models are not available in psastro

  • branches/eam_branch_20080719/pswarp/src/pswarp.h

    r18752 r18753  
    1010#include <psmodules.h>
    1111#include <psphot.h>
     12
     13#define THREADED 1
    1214
    1315#include "pswarpErrorCodes.h"
    … …  
    3335    int nXpts, nYpts;                   // number of x,y samples in the grid
    3436    int nXpix, nYpix;                   // x,y spacing in src image pixels of grid samples
    35     int xMin,  yMin;                    // coordinate of first grid sample
     37    double xMin,  yMin;                 // coordinate of first grid sample
    3638} pswarpMapGrid;
    3739
    … …  
    5254} pswarpTransformTileArgs;
    5355
    54 bool pswarpTransformReadout_Unthreaded(pmReadout *output, pmReadout *input, pmConfig *config);
    55 bool pswarpTransformReadout_Threaded(pmReadout *output, pmReadout *input, pmConfig *config);
    56 
    5756pswarpTransformTileArgs *pswarpTransformTileArgsAlloc();
    5857bool pswarpTransformTile (pswarpTransformTileArgs *args);
    … …  
    6665void pswarpCleanup (pmConfig *config);
    6766bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);
    68 bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config);
     67bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config);
    6968
    7069bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src);
    … …  
    7574pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
    7675bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
     76bool pswarpMapGridCoordRange (pswarpMapGrid *grid, int gridX, int gridY, psPlane *min, psPlane *max);
    7777int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX);
    7878int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY);

  • branches/eam_branch_20080719/pswarp/src/pswarpArguments.c

    r18622 r18753  
    3939                         "Filename for statistics of output image", argv[N]);
    4040        psArgumentRemove(N, &argc, argv);
     41    }
     42
     43    // Number of threads
     44    if ((N = psArgumentGet(argc, argv, "-threads"))) {
     45        psArgumentRemove(N, &argc, argv);
     46        int nThreads = atoi(argv[N]);
     47        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of warp threads", nThreads);
     48        psArgumentRemove(N, &argc, argv);
     49
     50        // create the thread pool with number of desired threads, supplying our thread launcher function
     51        // XXX need to determine the number of threads from the config data
     52        psThreadPoolInit (nThreads, &pswarpThreadLauncher);
    4153    }
    4254

  • branches/eam_branch_20080719/pswarp/src/pswarpCleanup.c

    r14641 r18753  
    99    pmModelClassCleanup ();
    1010    psTimeFinalize ();
     11    psThreadPoolFinalize ();
    1112    pmConceptsDone ();
    1213    pmConfigDone ();

  • branches/eam_branch_20080719/pswarp/src/pswarpLoop.c

    r18622 r18753  
    241241                }
    242242
    243                 pswarpTransformReadout_Opt(output, readout, config);
     243                pswarpTransformReadout(output, readout, config);
    244244
    245245                pmFPAfileIOChecks(config, view, PM_FPA_AFTER);

  • branches/eam_branch_20080719/pswarp/src/pswarpMapGrid.c

    r12523 r18753  
    44// coordinates (src) to destination coordinates (dest).  we construct a grid with superpixel
    55// spacing of nXpix, nYpix.  The transformation for each grid cell is valid for the superpixel.
    6 // The grid over-fills the source image so ever source image pixel is guaranteed to have a map.
     6// The grid over-fills the source image so every source image pixel is guaranteed to have a map.
    77pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) {
    88
    … …  
    1111
    1212    // start counting from the center of the superpixels
    13     int xMin = 0.5*nXpix;
    14     int yMin = 0.5*nYpix;
     13    double xMin = 0.5*nXpix;
     14    double yMin = 0.5*nYpix;
    1515
    1616    // the map is defined for coordinates in the image parent frame.
    … …  
    4242
    4343// set the grid coordinate (gridX,gridY) for the given source image coordinate (ix,iy)
     44// XXX return true if the result is on the src image, false otherwise (???)
    4445bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) {
    4546
    4647    *gridX = 0.5 + (ix - grid->xMin) / (double) grid->nXpix;
    4748    *gridY = 0.5 + (iy - grid->yMin) / (double) grid->nYpix;
     49
     50    return true;
     51}
     52
     53// given the specified grid coordinate (gridX, gridY), return the min and max coordinates for the tile
     54bool pswarpMapGridCoordRange (pswarpMapGrid *grid, int gridX, int gridY, psPlane *min, psPlane *max) {
     55
     56    min->x = (gridX - 0.5)*grid->nXpix + grid->xMin;
     57    min->y = (gridY - 0.5)*grid->nYpix + grid->yMin;
     58
     59    max->x = min->x + grid->nXpix;
     60    max->y = min->y + grid->nYpix;
     61
    4862    return true;
    4963}

  • branches/eam_branch_20080719/pswarp/src/pswarpThreadLauncher.c

    r18752 r18753  
    11# include "pswarp.h"
     2// XXX add in a thread exit api (set a local 'exit' variable)
    23
    34bool pswarpThreadTest (psArray *args) {

  • branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout.c

    r12744 r18753  
    11# include "pswarp.h"
    22
    3 bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) {
    4 
     3// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
     4bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config)
     5{
    56    // XXX this implementation currently ignores the use of the region
    67    psImage *region = NULL;
    7     pmCell *cell = NULL;
    88
    9     // select the current recipe
    10     // psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     9    psTimerStart ("warp");
    1110
    12     int outNx = output->image->numCols;
    13     int outNy = output->image->numRows;
     11    // 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");
    1416
    15     psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
    16     psPlane *outPix = psPlaneAlloc();   // Coordinates on the output detector
     17    // load the recipe
     18    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     19    psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
    1720
    18     psPlane *FP = psPlaneAlloc();       // Coordinates on the focal plane
    19     psPlane *TP = psPlaneAlloc();       // Coordinates on the tangent plane
    20     psSphere *sky = psSphereAlloc();    // Coordinates on the sky
     21    // output mask bits
     22    psMaskType maskIn   = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT");
     23    psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config);
     24    psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config);
     25    psAssert (mdok, "MASK.INPUT was not defined");
    2126
    22     cell = input->parent;
    23     pmChip *chipInput = cell->parent;
    24     pmFPA *fpaInput = chipInput->parent;
     27    int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
     28    if (!mdok) nThreads = 0;
    2529
    26     cell = output->parent;
    27     pmChip *chipOutput = cell->parent;
    28     pmFPA *fpaOutput = chipOutput->parent;
     30    // Flux fraction for "poor"
     31    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC");
    2932
    30     psF32 **outData = output->image->data.F32;
    31     psImage *inImage = input->image;
     33    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
     34    // output coordinates to input coordinates
     35    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);
    3236
    33     psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage,
    34                                                                        NULL, NULL, 0, NAN, NAN, 0, 0, 0.0);
     37    // 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);
    3540
    36     // Iterate over the output image pixels
    37     for (int y = 0; y < outNy; y++) {
    38         for (int x = 0; x < outNx; x++) {
    39             // Only transform those pixels requested
    40             if (region && region->data.U8[y][x]) continue;
     41    // 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);
    4143
    42             // XXX double check this 1/2 pixel offset
    43             outPix->x = (double)x + 0.5;
    44             outPix->y = (double)y + 0.5;
     44    if (input->weight && !output->weight) {
     45        output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
     46        psImageInit(output->weight, NAN);
     47    }
     48    if ((input->mask || maskPoor || maskBad) && !output->mask) {
     49        output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
     50        psImageInit(output->mask, maskBad);
     51    }
    4552
    46             psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
    47             psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
    48             psDeproject (sky, TP, fpaOutput->toSky);
     53    // total number of good pixels across all tiles (summed below)
     54    int goodPixels = 0;
    4955
    50             psProject (TP, sky, fpaInput->toSky);
    51             psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
    52             psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
     56    // create jobs and supply them to the threads
     57    for (int gridY = 0; gridY < grid->nYpts; gridY++) {
     58        for (int gridX = 0; gridX < grid->nXpts; gridX++) {
    5359
    54             // XXX get interpolation method from the recipe
    55             double value;
    56             if (!psImageInterpolate(&value, NULL, NULL, inPix->x, inPix->y, interp)) {
    57                 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    58                 psFree(interp);
    59                 psFree(inPix);
    60                 psFree(outPix);
    61                 psFree(FP);
    62                 psFree(TP);
    63                 psFree(sky);
    64                 return false;
    65             }
     60            pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
    6661
    67             outData[y][x] = value;
    68             // modify zero and scale?
     62            // these items are just views to the data; they are not freed with args
     63            args->input = input;
     64            args->output = output;
     65            args->grid = grid;
     66            args->interp = interp;
     67            args->region = region;
     68
     69            args->gridX = gridX;
     70            args->gridY = gridY;
     71            args->goodPixels = 0;
     72
     73            if (nThreads) {
     74                // allocate a job
     75                psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE", 0);
     76
     77                // construct the arguments for this job
     78                // job is pswarpTransformTile (gridX, gridY);
     79                psArrayAdd (job->args, 1, args);
     80                // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);
     81                psThreadJobAddPending (job);
     82            } else {
     83                pswarpTransformTile (args);
     84                goodPixels += args->goodPixels;
     85            }
     86            psFree (args);
     87        }
     88    }
     89
     90    // wait for the threads to finish and manage results
     91    if (nThreads) {
     92        // wait here for the threaded jobs to finish
     93        if (!psThreadPoolWait ()) {
     94            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     95            return false;
     96        }
     97        fprintf (stderr, "success for threaded jobs\n");
     98
     99        // each job records its own goodPixel values; sum them here
     100        // we have only supplied one type of job, so we can assume the types here
     101        psThreadJob *job = NULL;
     102        while ((job = psThreadJobGetDone()) != NULL) {
     103            if (job->args->n < 1) {
     104                fprintf (stderr, "error with job\n");
     105            } else {
     106                pswarpTransformTileArgs *args = job->args->data[0];
     107                // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
     108                goodPixels += args->goodPixels;
     109            }
     110            psFree (job);
     111        }
     112    }
     113    psFree(interp);
     114    psFree(grid);
     115
     116
     117    // Store the variance factor and number of good pixels
     118    if (goodPixels > 0) {
     119        // Variance factor: large factor --> small scale
     120        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp);
     121        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
     122        if (vfItem) {
     123            psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX);
     124            psAssert(goodpixItem, "It should be where we left it!");
     125            psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64,
     126                     "Should be the type we said.");
     127
     128            vfItem->data.F32 += varFactor * goodPixels;
     129            goodpixItem->data.S64 += goodPixels;
     130        } else {
     131            psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0,
     132                             "Variance factor weighted by the good pixels", varFactor * goodPixels);
     133            psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0,
     134                             "Number of good pixels", goodPixels);
    69135        }
    70136    }
    71     psFree(interp);
    72     psFree(inPix);
    73     psFree(outPix);
    74     psFree(FP);
    75     psFree(TP);
    76     psFree(sky);
     137
     138    if (goodPixels > 0) {
     139        if (!pswarpTransformSources (output, input, config)) {
     140            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     141            return false;
     142        }
     143    }
     144
     145    // XXX should we not write anything out if there are no good pixels?
     146    output->data_exists = true;
     147    output->parent->data_exists = true;
     148    output->parent->parent->data_exists = true;
     149
     150    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
    77151
    78152    return true;

  • branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c

    r18752 r18753  
    11# include "pswarp.h"
    22
    3 #define SOURCE_ARRAY_BUFFER 100         // Size to grow the array of sources at a time
     3# define THREADED 1
     4# define SOURCE_ARRAY_BUFFER 100         // Size to grow the array of sources at a time
    45
    56// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
    … …  
    9091            args->goodPixels = 0;
    9192
     93# if (THREADED)     
    9294            // allocate a job
    9395            psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE", 0);
    … …  
    9799            psArrayAdd (job->args, 1, args);
    98100            psThreadJobAddPending (job);
     101# else
     102            pswarpTransformTile (args);
     103            goodPixels += args->goodPixels;
     104# endif
    99105            psFree (args);
    100106        }
    101107    }
     108
     109# if (THREADED)     
    102110    // wait here for the threaded jobs to finish
    103111    if (!psThreadPoolWait ()) {
    … …  
    115123        psFree (job);
    116124    }
     125# endif
    117126    psFree(interp);
    118127    psFree(grid);
     128
    119129
    120130    // Store the variance factor and number of good pixels

  • branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c

    r18752 r18753  
    3737    // int outNrow = args->output->image->numRows;
    3838
    39     // XXX these are wrong (min or max is wrong)
    40     int minX = pswarpMapGridNextGrid_Y (args->grid, args->gridY);
    41     int minY = pswarpMapGridNextGrid_X (args->grid, args->gridX);
    42     int maxX = pswarpMapGridNextGrid_Y (args->grid, args->gridY);
    43     int maxY = pswarpMapGridNextGrid_X (args->grid, args->gridX);
     39    // get the coordinate range for this grid tile
     40    psPlane minPt, maxPt;
     41    pswarpMapGridCoordRange (args->grid, args->gridX, args->gridY, &minPt, &maxPt);
    4442
    45     psF32 **outImageData     = args->output->image->data.F32;
    46     psF32 **outVarData       = args->output->weight->data.F32;
    47     psMaskType **outMaskData = args->output->mask->data.U8;
     43    psF32 **outImageData     = (args->output->image)  ? args->output->image->data.F32  : NULL;
     44    psF32 **outVarData       = (args->output->weight) ? args->output->weight->data.F32 : NULL;
     45    psMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.U8    : NULL;
    4846
    4947    pswarpMap *map = args->grid->maps[args->gridX][args->gridY];
    … …  
    6058    // Iterate over the output image pixels (parent frame)
    6159    long goodPixels = 0;                // Number of input pixels landing on the output image
    62     for (int y = minY; y < maxY; y++) {
    63         for (int x = minX; x < maxX; x++) {
     60    for (int y = minPt.y; y < maxPt.y; y++) {
     61        for (int x = minPt.x; x < maxPt.x; x++) {
    6462
    6563            // Only transform those pixels requested
Note: See TracChangeset for help on using the changeset viewer.