IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18752


Ignore:
Timestamp:
Jul 27, 2008, 8:00:23 AM (18 years ago)
Author:
eugene
Message:

threaded and unthreaded equivalent versions now compile

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

Legend:

Unmodified
Added
Removed

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

    r18748 r18752  
    1616        pswarpPixelFraction.c           \
    1717        pswarpSetMaskBits.c             \
     18        pswarpThreadLauncher.c          \
    1819        pswarpTransformReadout_Opt.c    \
    19         pswarpThreadLauncher.c          \
     20        pswarpTransformReadout_Threaded.c   \
     21        pswarpTransformReadout_Unthreaded.c \
     22        pswarpTransformTile.c           \
    2023        pswarpVersion.c           
    2124

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

    r18748 r18752  
    1818
    1919    // 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
    2021    psThreadPoolInit (4, &pswarpThreadLauncher);
    2122   
    22     srand48(0);
    23 
    24     for (int i = 0; i < 10; i++) {
    25         // allocate a job which takes two arguments
    26         // XXX probably don't need to pre-specify the number of args
    27         psThreadJob *job = psThreadJobAlloc ("THREAD_TEST", 2);
    28 
    29         // construct the arguments for this job
    30         char *arg1 = NULL;
    31         char *arg2 = NULL;
    32         psStringAppend (&arg1, "arg1.%d", i);
    33         psStringAppend (&arg2, "arg2.%d", i);
    34         psArrayAdd (job->args, 1, arg1);
    35         psArrayAdd (job->args, 1, arg2);
    36 
    37         psThreadJobAddToQueue (job);
    38     }
    39     if (!psThreadPoolWait ()) {
    40         fprintf (stderr, "error in a threaded job\n");
    41         exit (1);
    42     }
    43     fprintf (stderr, "success for threaded jobs\n");
    44     exit (0);
    45 
    4623    // model inits are needed in pmSourceIO
    4724    // models defined in psphot/src/models are not available in psastro

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

    r18748 r18752  
    3636} pswarpMapGrid;
    3737
     38typedef struct {
     39    // values which are common to all tiles
     40    pmReadout *input;
     41    pmReadout *output;
     42    pswarpMapGrid *grid;
     43    psImageInterpolateOptions *interp;
     44    psImage *region;
     45
     46    // input values for this tile
     47    int gridX;
     48    int gridY;
     49
     50    // output values for this tile
     51    long goodPixels;
     52} pswarpTransformTileArgs;
     53
     54bool pswarpTransformReadout_Unthreaded(pmReadout *output, pmReadout *input, pmConfig *config);
     55bool pswarpTransformReadout_Threaded(pmReadout *output, pmReadout *input, pmConfig *config);
     56
     57pswarpTransformTileArgs *pswarpTransformTileArgsAlloc();
     58bool pswarpTransformTile (pswarpTransformTileArgs *args);
     59
    3860pmConfig *pswarpArguments (int argc, char **argv);
    3961bool pswarpOptions(pmConfig *config);

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

    r18747 r18752  
    3333        // we have to lock here so the job queue cannot be empty yet no threads busy
    3434        psThreadLock();
    35         while ((job = psThreadJobGet ()) == NULL) {
     35        while ((job = psThreadJobGetPending ()) == NULL) {
    3636            psThreadUnlock();
    3737            usleep (10000);
    … …  
    5050            continue;
    5151        }
     52
     53        // list all allowed job types here
     54        if (!strcmp (job->type, "PSWARP_TRANSFORM_TILE")) {
     55            pswarpTransformTileArgs *args = job->args->data[0];
     56            bool status = pswarpTransformTile (args);
     57            if (!status) {
     58                self->fault = true;
     59            }
     60            // we do not have to lock here because this transition is not tied to the job queue
     61            self->busy = false; 
     62            continue;
     63        }
    5264    } 
    5365}

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

    r18666 r18752  
    1515    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX");
    1616    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY");
    17     psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
    18                                                                    "INTERPOLATION.MODE");
     17    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE");
     18
    1919    // load the recipe
    2020    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
    21     if (!recipe) {
    22         psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
    23         return false;
    24     }
     21    psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);
    2522
    2623    // output mask bits
    … …  
    3027    psAssert (mdok, "MASK.INPUT was not defined");
    3128
    32     float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor"
    33 
    34     // we need to apply the offset to convert parent coordinates to child coordinates for
    35     // psImagePixelInterpolate below
    36     int inCol0 = input->image->col0;
    37     int inRow0 = input->image->row0;
    38     int outCol0 = output->image->col0;
    39     int outRow0 = output->image->row0;
    40 
    41     // we might want to do the rectangular regions outside of the selection independently
    42     // psImageInit (output->image, NAN);
     29    // Flux fraction for "poor"
     30    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC");
    4331
    4432    // find the output pixel range
    … …  
    5442    psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);
    5543
    56     int gridX, gridY, nextGridX, nextGridY;
     44    int gridX, gridY;
    5745    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY);
    58     nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
    59     pswarpMap *map = grid->maps[gridX][gridY];
    60 
     46# if (0)
     47    // XXX these asserts probably belong in the function (don't have outCol0,Row0 anyway)
    6148    assert ((int)(minX - outCol0) >= 0);
    6249    assert ((int)(maxX - outCol0) <= output->image->numCols);
    6350    assert ((int)(minY - outRow0) >= 0);
    6451    assert ((int)(maxY - outRow0) <= output->image->numRows);
    65 
    66     int gridXo = gridX;
    67     int nextGridXo = pswarpMapGridNextGrid_X (grid, gridX);
     52# endif
    6853
    6954    psImage *inImage = input->image;    // Input image
    … …  
    7661                                                                       maskBad, maskPoor, poorFrac);
    7762
    78     psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
    79     psF32 **outImageData = output->image->data.F32; // Output image pixels
    80     psF32 **outVarData = NULL;          // Output variance pixels
    81     psMaskType **outMaskData = NULL;    // Output mask pixels
    82     if (inVar) {
    83         if (!output->weight) {
    84             output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
    85             psImageInit(output->weight, NAN);
    86         }
    87         outVarData = output->weight->data.F32;
    88     }
    89     if (inMask || maskPoor || maskBad) {
    90         if (!output->mask) {
    91             output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
    92             psImageInit(output->mask, maskBad);
    93         }
    94         outMaskData = output->mask->data.PS_TYPE_MASK_DATA;
    95     }
    96 
    97     // create N threads
    98     for (gridY = gridYmin; gridY < gridYmax; gridY++) {
    99         for (gridX = gridXmin; gridX < gridXmax; gridX++) {
    100             pswarpTransformTile (gridX, gridY);
    101             // send function and arguments to thread manager
    102             // threadFunction (pswarpTransformTile)
    103             // while (threadFull()) {
    104             // usleep (50000);
    105             // }
     63    if (inVar && !output->weight) {
     64        output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
     65        psImageInit(output->weight, NAN);
     66    }
     67    if ((inMask || maskPoor || maskBad) && !output->mask) {
     68        output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
     69        psImageInit(output->mask, maskBad);
     70    }
     71
     72    // total number of good pixels across all tiles (summed below)
     73    int goodPixels = 0;
     74
     75    // create jobs and supply them to the threads
     76    for (gridY = 0; gridY < nGridX; gridY++) {
     77        for (gridX = 0; gridX < nGridY; gridX++) {
     78
     79            pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
     80
     81            // these items are just views to the data; they are not freed with args
     82            args->input = input;
     83            args->output = output;
     84            args->grid = grid;
     85            args->interp = interp;
     86            args->region = region;
     87
     88            args->gridX = gridX;
     89            args->gridY = gridY;
     90            args->goodPixels = 0;
     91
     92            // allocate a job
     93            psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE", 0);
     94
     95            // construct the arguments for this job
     96            // job is pswarpTransformTile (gridX, gridY);
     97            psArrayAdd (job->args, 1, args);
     98            psThreadJobAddPending (job);
     99            psFree (args);
    106100        }
    107101    }
    108 
    109 
    110     // Iterate over the output image pixels (parent frame)
    111     long goodPixels = 0;                // Number of input pixels landing on the output image
    112     for (int y = minY; y < maxY; y++) {
    113         if (y >= nextGridY) {
    114             gridY ++;
    115             nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
    116             map = grid->maps[gridX][gridY];
    117         }
    118 
    119         gridX = gridXo;
    120         nextGridX = nextGridXo;
    121         map = grid->maps[gridX][gridY];
    122         int yOut = y - outRow0;         // Position on image
    123         for (int x = minX; x < maxX; x++) {
    124             if (x >= nextGridX) {
    125                 gridX ++;
    126                 nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
    127                 map = grid->maps[gridX][gridY];
    128             }
    129 
    130             // Only transform those pixels requested
    131             if (region && region->data.U8[y][x]) continue;
    132 
    133             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
    134             // both are in the parent frames of the input and output images.
    135             pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
    136 
    137             if (inPix->x - inCol0 < 0) continue;
    138             if (inPix->x - inCol0 >= inImage->numCols) continue;
    139             if (inPix->y - inRow0 < 0) continue;
    140             if (inPix->y - inRow0 >= inImage->numRows) continue;
    141 
    142             goodPixels++;
    143 
    144             // XXX include mask
    145             // XXX apply scale and offset?
    146             // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
    147             double imageValue, varValue;
    148             psMaskType maskValue;
    149             if (!psImageInterpolate(&imageValue, &varValue, &maskValue,
    150                                     inPix->x - inCol0, inPix->y - inRow0, interp)) {
    151                 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    152                 psFree(interp);
    153                 psFree(inPix);
    154                 psFree(grid);
    155                 return false;
    156             }
    157             int xOut = x - outCol0;     // Position on image
    158             outImageData[yOut][xOut] = imageValue;
    159             if (inVar) {
    160                 outVarData[yOut][xOut] = varValue;
    161             }
    162             if (outMaskData) {
    163                 outMaskData[yOut][xOut] = maskValue;
    164             }
    165         }
     102    // wait here for the threaded jobs to finish
     103    if (!psThreadPoolWait ()) {
     104        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     105        return false;
     106    }
     107    fprintf (stderr, "success for threaded jobs\n");
     108
     109    // each job records its own goodPixel values; sum them here
     110    // we have only supplied one type of job, so we can assume the types here
     111    psThreadJob *job = NULL;
     112    while ((job = psThreadJobGetDone()) != NULL) {
     113        pswarpTransformTileArgs *args = job->args->data[0];
     114        goodPixels += args->goodPixels;
     115        psFree (job);
    166116    }
    167117    psFree(interp);
    168     psFree(inPix);
    169118    psFree(grid);
    170119
    171120    // Store the variance factor and number of good pixels
    172121    if (goodPixels > 0) {
    173         float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inCol0,
    174                                                            inImage->numRows / 2.0 + inRow0,
     122        float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inImage->col0,
     123                                                           inImage->numRows / 2.0 + inImage->row0,
    175124                                                           interp); // Variance factor: large --> small scale
    176125        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
    … …  
    208157            pmModel *model = source->modelPSF; // Model for this source
    209158            float xIn, yIn;             // Coordinates of source
    210             xIn = model->params->data.F32[PM_PAR_XPOS] - inCol0;
    211             yIn = model->params->data.F32[PM_PAR_YPOS] - inRow0;
     159            xIn = model->params->data.F32[PM_PAR_XPOS] - inImage->col0;
     160            yIn = model->params->data.F32[PM_PAR_YPOS] - inImage->row0;
    212161            int xGrid, yGrid;           // Grid coordinates for local map
    213162            if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) {
    … …  
    232181                return false;
    233182            }
    234             xOut += outCol0 - 0.5;
    235             yOut += outRow0 - 0.5;
     183            xOut += output->image->col0 - 0.5;
     184            yOut += output->image->row0 - 0.5;
    236185            if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) {
    237186                // It's not in the output image

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

    r18666 r18752  
     1# include "pswarp.h"
    12
    2 // XXX need to determine all of the needed inputs and put them in a structure
     3void pswarpTransformTileArgsFree (pswarpTransformTileArgs *args) {
     4    if (!args) return;
     5    return;
     6}
    37
    4 typedef struct {
    5     int gridX;
    6     int gridY;
    7 } pswarpTransformTileOpts;
     8pswarpTransformTileArgs *pswarpTransformTileArgsAlloc() {
    89
    9 pswarpTransformTile (int gridX, int gridY) {
     10    pswarpTransformTileArgs *args = (pswarpTransformTileArgs *)psAlloc(sizeof(pswarpTransformTileArgs));
     11    psMemSetDeallocator(args, (psFreeFunc)pswarpTransformTileArgsFree);
    1012
    11     minX = pswarpMapGridNextGrid_Y (grid, gridY);
    12     minY = pswarpMapGridNextGrid_X (grid, gridX);
     13    args->input = NULL;
     14    args->output = NULL;
     15    args->grid = NULL;
     16    args->interp = NULL;
     17    args->region = NULL;
    1318
    14     maxX = pswarpMapGridNextGrid_Y (grid, gridY);
    15     maxY = pswarpMapGridNextGrid_X (grid, gridX);
     19    args->gridX = 0;
     20    args->gridY = 0;
    1621
    17     map = grid->maps[gridX][gridY];
     22    args->goodPixels = 0;
    1823
    19     psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
     24    return args;
     25}
     26
     27bool pswarpTransformTile (pswarpTransformTileArgs *args) {
     28
     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;
     33
     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;
     38
     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);
     44
     45    psF32 **outImageData     = args->output->image->data.F32;
     46    psF32 **outVarData       = args->output->weight->data.F32;
     47    psMaskType **outMaskData = args->output->mask->data.U8;
     48
     49    pswarpMap *map = args->grid->maps[args->gridX][args->gridY];
     50
     51    psImage *region = args->region;
     52
     53    double xInRaw, yInRaw;
     54
     55    // output values for this pixel
     56    double imageValue;
     57    double varValue;
     58    psMaskType maskValue;
    2059
    2160    // Iterate over the output image pixels (parent frame)
    2261    long goodPixels = 0;                // Number of input pixels landing on the output image
    2362    for (int y = minY; y < maxY; y++) {
    24         int yOut = y - outRow0;         // Position on image
    2563        for (int x = minX; x < maxX; x++) {
    2664
    … …  
    3068            // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
    3169            // both are in the parent frames of the input and output images.
    32             pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
     70            pswarpMapApply (&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
    3371
    34             if (inPix->x - inCol0 < 0) continue;
    35             if (inPix->x - inCol0 >= inImage->numCols) continue;
    36             if (inPix->y - inRow0 < 0) continue;
    37             if (inPix->y - inRow0 >= inImage->numRows) continue;
     72            double xIn = xInRaw - outCol0;      // Position on input image
     73            double yIn = yInRaw - outRow0;      // Position on input image
     74
     75            if (xIn < 0) continue;
     76            if (yIn >= inNcol) continue;
     77            if (yIn < 0) continue;
     78            if (yIn >= inNrow) continue;
    3879
    3980            goodPixels++;
    4081
    41             // XXX include mask
    42             // XXX apply scale and offset?
    4382            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
    44             double imageValue, varValue;
    45             psMaskType maskValue;
    46             if (!psImageInterpolate(&imageValue, &varValue, &maskValue,
    47                                     inPix->x - inCol0, inPix->y - inRow0, interp)) {
     83            if (!psImageInterpolate(&imageValue, &varValue, &maskValue, xIn, yIn, args->interp)) {
    4884                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    49                 psFree(interp);
    50                 psFree(inPix);
    51                 psFree(grid);
    5285                return false;
    5386            }
    54             int xOut = x - outCol0;     // Position on image
    55             outImageData[yOut][xOut] = imageValue;
    56             if (inVar) {
     87
     88            int xOut = x - outCol0;     // Position on output image
     89            int yOut = y - outRow0;     // Position on output image
     90
     91            // not all images need be transformed
     92            if (outImageData) {
     93                outImageData[yOut][xOut] = imageValue;
     94            }
     95            if (outVarData) {
    5796                outVarData[yOut][xOut] = varValue;
    5897            }
    … …  
    62101        }
    63102    }
    64 
    65     psFree(inPix);
     103    args->goodPixels = goodPixels;
     104    return true;
    66105}
Note: See TracChangeset for help on using the changeset viewer.