IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18839


Ignore:
Timestamp:
Jul 31, 2008, 2:13:59 PM (18 years ago)
Author:
eugene
Message:

adding multithread capability

Location:
trunk
Files:
8 added
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppMerge/src/Makefile.am

    r18758 r18839  
    1010        ppMergeFiles.c          \
    1111        ppMergeScaleZero.c      \
    12         ppMergeLoop.c           \
    1312        ppMergeFileGroup.c      \
     13        ppMergeReadChunk.c      \
     14        ppMergeLoop_Threaded.c  \
     15        ppMergeSetThreads.c     \
    1416        ppMergeMask.c
    1517
    16 #       ppMergeLoop_Threaded.c 
    17 #       ppMergeThreadLauncher.c
    18 
     18#       ppMergeLoop.c           
    1919
    2020noinst_HEADERS =                \

  • trunk/ppMerge/src/ppMerge.c

    r18757 r18839  
    1010{
    1111    psLibInit(NULL);
    12     psMemSetThreadSafety(false);
    1312    psTimerStart(TIMERNAME);
    1413

  • trunk/ppMerge/src/ppMerge.h

    r18758 r18839  
    1515#define TIMERNAME "ppMerge"             // Name for timer
    1616#define PPMERGE_RECIPE "PPMERGE"        // Recipe name
    17 #define THREADED 0
     17#define THREADED 1
    1818
    1919// Type of frame to merge
    … …  
    3939    bool read;
    4040    bool busy;
     41    int firstScan;
     42    int lastScan;
    4143} ppMergeFileGroup;
    4244
    … …  
    100102
    101103ppMergeFileGroup *ppMergeFileGroupAlloc();
    102 ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk);
     104ppMergeFileGroup *ppMergeReadChunk (bool *status, psArray *fileGroups, pmConfig *config, int numChunk);
    103105void *ppMergeThreadLauncher (void *data);
    104106
     107bool ppMergeSetThreads ();
     108
    105109#endif

  • trunk/ppMerge/src/ppMergeArguments.c

    r18758 r18839  
    171171    }
    172172
    173 # if (THREADED)
    174173    // Number of threads
    175174    if ((argnum = psArgumentGet(argc, argv, "-threads"))) {
    … …  
    181180        // create the thread pool with number of desired threads, supplying our thread launcher function
    182181        // XXX need to determine the number of threads from the config data
    183         psThreadPoolInit (nThreads, &ppMergeThreadLauncher);
    184     }
    185 # endif
     182        psThreadPoolInit (nThreads);
     183    }
     184    ppMergeSetThreads();
    186185
    187186    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {

  • trunk/ppMerge/src/ppMergeLoop_Threaded.c

    r18757 r18839  
    3636
    3737    // General combination parameters
    38     int rows = psMetadataLookupS32(NULL, arguments, "ROWS"); // Number of rows to read per chunk
    3938    int iter = psMetadataLookupS32(NULL, arguments, "ITER"); // Number of rejection iterations
    4039    float rej = psMetadataLookupF32(NULL, arguments, "REJ"); // Rejection level
    … …  
    175174                fileGroup->read = false;
    176175                fileGroup->busy = false;
     176                fileGroup->lastScan = 0;
     177                fileGroup->firstScan = 0;
    177178                fileGroups->data[i] = fileGroup;
    178179            }
    179180
     181            // call the init functions
     182            switch (type) {
     183              case PPMERGE_TYPE_BIAS:
     184              case PPMERGE_TYPE_FLAT:
     185              case PPMERGE_TYPE_FRINGE:
     186                psAssert (fileGroups->n > 0, "no valid file groups defined");
     187                ppMergeFileGroup *fileGroup = fileGroups->data[0];
     188                if (!pmReadoutCombinePrepare(outRO, fileGroup->readouts, combination)) {
     189                    goto ERROR;
     190                }
     191                break;
     192
     193              default:
     194                fprintf (stderr, "not yet ready");
     195                goto ERROR;
     196            }
     197
    180198            // Read input data by chunks
     199            // psTimerStart ("ppMergeLoop");
    181200            for (int numChunk = 0; true; numChunk++) {
    182                
     201
     202                bool status = false;
    183203                ppMergeFileGroup *fileGroup = ppMergeReadChunk (&status, fileGroups, config, numChunk);
    184204                if (!status) goto ERROR;
    185205                if (!fileGroup) break;
    186206
     207                psThreadJob *job = NULL;
     208
    187209                switch (type) {
    188210                  case PPMERGE_TYPE_SHUTTER:
    189                     if (nThreads) {
    190                         // allocate a job
    191                         psThreadJob *job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION", 0);
    192 
    193                         // construct the arguments for this job
    194                         psArrayAdd (job->args, 1, outRO);
    195                         psArrayAdd (job->args, 1, fileGroup);
    196                         psArrayAdd (job->args, 1, psScalarAlloc(shutterRef, PS_TYPE_F32))
    197                         psArrayAdd (job->args, 1, shutters->data[cellNum]);
    198                         psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32));
    199                         psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32));
    200                         psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8));
    201 
    202                         psThreadJobAddPending (job);
    203                     } else {
    204                         if (!pmShutterCorrectionGenerate(outRO, NULL, fileGroup->readouts, shutterRef, shutters->data[cellNum], iter, rej, maskVal)) {
    205                             goto ERROR;
    206                         }
    207                         fileGroup->busy = false;
     211                    // allocate a job
     212                    job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION");
     213
     214                    // construct the arguments for this job
     215                    psArrayAdd (job->args, 1, outRO);
     216                    psArrayAdd (job->args, 1, fileGroup);
     217                    psArrayAdd (job->args, 1, psScalarAlloc(shutterRef, PS_TYPE_F32));
     218                    psArrayAdd (job->args, 1, shutters->data[cellNum]);
     219                    psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32));
     220                    psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32));
     221                    psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8));
     222
     223                    // call: pmShutterCorrectionGenerate(outRO, NULL, fileGroup->readouts, shutterRef, shutters->data[cellNum], iter, rej, maskVal)
     224                    if (!psThreadJobAddPending (job)) {
     225                        goto ERROR;
    208226                    }
    209227                    break;
    210228                  case PPMERGE_TYPE_DARK:
    211                     if (nThreads) {
    212                         // allocate a job
    213                         psThreadJob *job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE", 0);
    214 
    215                         // construct the arguments for this job
    216                         psArrayAdd (job->args, 1, outCell);
    217                         psArrayAdd (job->args, 1, fileGroup);
    218                         psArrayAdd (job->args, 1, darkOrdinates);
    219                         psArrayAdd (job->args, 1, darkNorm);
    220                         psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32));
    221                         psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32));
    222                         psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8));
    223 
    224                         psThreadJobAddPending (job);
    225                     } else {
    226                         if (!pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter, rej, maskVal)) {
    227                             goto ERROR;
    228                         }
    229                         fileGroup->busy = false;
     229                    // allocate a job
     230                    job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE");
     231
     232                    // construct the arguments for this job
     233                    psArrayAdd (job->args, 1, outCell);
     234                    psArrayAdd (job->args, 1, fileGroup);
     235                    psArrayAdd (job->args, 1, darkOrdinates);
     236                    psArrayAdd (job->args, 1, darkNorm);
     237                    psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32));
     238                    psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32));
     239                    psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8));
     240
     241                    // call: pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter, rej, maskVal);
     242                    if (!psThreadJobAddPending (job)) {                 
     243                        goto ERROR;
    230244                    }
    231245                    break;
    … …  
    233247                  case PPMERGE_TYPE_FLAT:
    234248                  case PPMERGE_TYPE_FRINGE:
    235                     if (nThreads) {
    236                         // allocate a job
    237                         psThreadJob *job = psThreadJobAlloc ("PPMERGE_READOUT_COMBINE", 0);
    238 
    239                         // construct the arguments for this job
    240                         psArrayAdd (job->args, 1, outRO);
    241                         psArrayAdd (job->args, 1, fileGroup);
    242                         psArrayAdd (job->args, 1, zeros);
    243                         psArrayAdd (job->args, 1, scales);
    244                         psArrayAdd (job->args, 1, combination);
    245 
    246                         psThreadJobAddPending (job);
    247                     } else {
    248                         if (!pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination)) {
    249                             goto ERROR;
    250                         }
    251                         fileGroup->busy = false;
     249                    // allocate a job
     250                    job = psThreadJobAlloc ("PPMERGE_READOUT_COMBINE");
     251
     252                    // construct the arguments for this job
     253                    psArrayAdd (job->args, 1, outRO);
     254                    psArrayAdd (job->args, 1, fileGroup);
     255                    psArrayAdd (job->args, 1, zeros);
     256                    psArrayAdd (job->args, 1, scales);
     257                    psArrayAdd (job->args, 1, combination);
     258
     259                    // call: pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination);
     260                    if (!psThreadJobAddPending (job)) {
     261                        goto ERROR;
    252262                    }
    253263                    break;
    … …  
    258268
    259269            // wait for the threads to finish and manage results
    260             if (nThreads) {
    261                 // wait here for the threaded jobs to finish
    262                 if (!psThreadPoolWait ()) {
    263                     psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    264                     return false;
    265                 }
    266                 fprintf (stderr, "success for threaded jobs\n");
    267 
    268                 // we don't care about the results, just dump the done queue jobs
    269                 psThreadJob *job = NULL;
    270                 while ((job = psThreadJobGetDone()) != NULL) {
    271                     psFree (job);
    272                 }
     270            if (!psThreadPoolWait ()) {
     271                psError(PS_ERR_UNKNOWN, false, "Unable to combine images.");
     272                return false;
    273273            }
     274
     275            // we don't care about the results, just dump the done queue jobs
     276            psThreadJob *job = NULL;
     277            while ((job = psThreadJobGetDone()) != NULL) {
     278                psFree (job);
     279            }
     280
     281            psFree(fileGroups);
    274282
    275283            // Get list of cells for concepts averaging
    276284            psList *inCells = psListAlloc(NULL); // List of cells
    277285            for (int i = 0; i < numFiles; i++) {
    278                 pmReadout *readout = readouts->data[i]; // Readout of interest
    279                 psListAdd(inCells, PS_LIST_TAIL, readout->parent);
     286                pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i);
     287                pmCell *inCell = pmFPAviewThisCell(view, input->fpa); // Input cell
     288                psListAdd(inCells, PS_LIST_TAIL, inCell);
    280289            }
    281290            if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) {
    … …  
    286295            }
    287296            psFree(inCells);
    288 
    289             psFree(fileGroups);
     297            // fprintf (stdout, "done ppMergeLoop for cell : %f\n", psTimerMark ("ppMergeLoop"));
    290298
    291299            // Plug supplementary images into their own FPAs

  • trunk/ppMerge/src/ppMergeReadChunk.c

    r18757 r18839  
    11# include "ppMerge.h"
    22
    3 ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk) {
     3ppMergeFileGroup *ppMergeReadChunk (bool *status, psArray *fileGroups, pmConfig *config, int numChunk) {
     4
     5    *status = true;
     6
     7    bool mdok;
     8    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
     9    bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights?
     10    int rows = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of rows to read per chunk
    411
    512    // select an available fileGroup
    6  
    7     bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks?
    8     bool haveWeights = psMetadataLookupBool(&mdok, arguments, "INPUTS.WEIGHTS"); // Do we have weights?
    9 
    1013    while (1) {
    1114        // check for any fileGroups which can read data
    … …  
    1417            if (fileGroup->read) continue;
    1518
     19            // find max last scan so far
     20            int lastScan = 0;
     21            for (int i = 0; i < fileGroups->n; i++) {
     22                ppMergeFileGroup *fileGroup = fileGroups->data[i];
     23                lastScan = PS_MAX (fileGroup->lastScan, lastScan);
     24            }
     25            fileGroup->firstScan = lastScan;
     26            fileGroup->lastScan = lastScan + rows;
     27
    1628            psArray *readouts = fileGroup->readouts;
     29
     30            psTimerStart ("ppMergeReadChunk");
    1731
    1832            psTrace("ppStack", 2, "Reading data for chunk %d into fileGroup %d....n", numChunk, j);
    … …  
    2034                pmReadout *inRO = readouts->data[i]; // Input readout
    2135
     36                // override the recorded last scan
     37                inRO->thisImageScan  = fileGroup->firstScan;
     38                inRO->thisWeightScan = fileGroup->firstScan;
     39                inRO->thisMaskScan   = fileGroup->firstScan;
     40
    2241                // Read a chunk from a file
    2342                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i);
    24                 if (!pmReadoutReadChunk(inRO, file->fits, 0, rows, 0, config)) {
    25                     psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT %d", numChunk, i);
    26                     return NULL;
    27                 }                                                       
    2843
    29                 if (haveMasks) {
    30                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i);
    31                     if (!pmReadoutReadChunkMask(inRO, file->fits, 0, rows, 0, config)) {
    32                         psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.MASK %d", numChunk, NAME, i);
     44                bool keepReading = false;
     45                if (pmReadoutMore(inRO, file->fits, 0, rows, config)) {
     46                    keepReading = true;
     47                    if (!pmReadoutReadChunk(inRO, file->fits, 0, rows, 0, config)) {
     48                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT %d", numChunk, i);
     49                        *status = false;
    3350                        return NULL;
    3451                    }                                                   
    3552                }
    3653
    37                 if (haveWeights) {
    38                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i);
    39                     if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) {
    40                         psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d", numChunk, NAME, i);
     54                if (haveMasks && pmReadoutMoreMask(inRO, file->fits, 0, rows, config)) {
     55                    keepReading = true;
     56                    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i);
     57                    if (!pmReadoutReadChunkMask(inRO, file->fits, 0, rows, 0, config)) {
     58                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.MASK %d", numChunk, i);
     59                        *status = false;
    4160                        return NULL;
    4261                    }                                                   
    4362                }
     63
     64                if (haveWeights && pmReadoutMoreWeight(inRO, file->fits, 0, rows, config)) {
     65                    keepReading = true;
     66                    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i);
     67                    if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) {
     68                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d", numChunk, i);
     69                        *status = false;
     70                        return NULL;
     71                    }                                                   
     72                }
     73                if (!keepReading) {
     74                    return NULL;
     75                }
    4476            }
     77
    4578            fileGroup->read = fileGroup->busy = true;
    4679            return fileGroup;
    … …  
    4881
    4982        // check for any fileGroups which are done processing
    50         bool wait = true;
    51         bool more = true;
     83        bool wait = false;
    5284        for (int j = 0; j < fileGroups->n; j++) {
    5385            ppMergeFileGroup *fileGroup = fileGroups->data[j];
    5486            if (!fileGroup->read || fileGroup->busy) continue;
    55            
    56             wait = false;
    57             psArray *readouts = fileGroup->readouts;
    58             // any more data to be read?
    59             for (int i = 0; i < readouts->n && more; i++) {
    60                 pmReadout *inRO = readouts->data[i];
    61 
    62                 pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i);
    63                 more &= pmReadoutMore(inRO, file->fits, 0, rows, config);
    64 
    65                 if (haveMasks) {
    66                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i);
    67                     more &= pmReadoutMoreMask(inRO, file->fits, 0, rows, config);
    68                 }
    69                 if (haveWeights) {
    70                     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i);
    71                     more &= pmReadoutMoreWeight(inRO, file->fits, 0, rows, config);
    72                 }
    73             }
    7487            fileGroup->read = false;
     88            wait = true;
    7589        }
    76         if (!more) return NULL;
    77 
    7890        if (wait) usleep (10000);
    7991    }

  • trunk/ppMerge/src/ppMergeThreadLauncher.c

    r18757 r18839  
    4141                self->fault = true;
    4242            }
     43
    4344            // we do not have to lock here because this transition is not tied to the job queue
    4445            fileGroup->busy = false;

  • trunk/pswarp/src/Makefile.am

    r18558 r18839  
    1616        pswarpPixelFraction.c           \
    1717        pswarpSetMaskBits.c             \
    18         pswarpTransformReadout_Opt.c    \
     18        pswarpSetThreads.c              \
     19        pswarpTransformReadout.c        \
     20        pswarpTransformSources.c \
     21        pswarpTransformTile.c           \
    1922        pswarpVersion.c           
    2023

  • trunk/pswarp/src/pswarp.h

    r18558 r18839  
    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;
     39
     40typedef struct {
     41    // values which are common to all tiles
     42    pmReadout *input;
     43    pmReadout *output;
     44    pswarpMapGrid *grid;
     45    psImageInterpolateOptions *interp;
     46    psImage *region;
     47
     48    // input values for this tile
     49    int gridX;
     50    int gridY;
     51
     52    // output values for this tile
     53    long goodPixels;
     54} pswarpTransformTileArgs;
     55
     56pswarpTransformTileArgs *pswarpTransformTileArgsAlloc();
     57bool pswarpTransformTile (pswarpTransformTileArgs *args);
    3758
    3859pmConfig *pswarpArguments (int argc, char **argv);
    … …  
    4465void pswarpCleanup (pmConfig *config);
    4566bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);
    46 bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config);
     67bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config);
    4768
    4869bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src);
    … …  
    5374pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
    5475bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
     76bool pswarpMapGridCoordRange (pswarpMapGrid *grid, int gridX, int gridY, psPlane *min, psPlane *max);
    5577int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX);
    5678int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY);
    … …  
    6890    );
    6991
     92// define threads for this program
     93bool pswarpSetThreads ();

  • trunk/pswarp/src/pswarpArguments.c

    r18558 r18839  
    4040        psArgumentRemove(N, &argc, argv);
    4141    }
     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);
     53    }
     54    pswarpSetThreads ();
    4255
    4356    // PSF determination?

  • trunk/pswarp/src/pswarpCleanup.c

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

  • trunk/pswarp/src/pswarpLoop.c

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

  • trunk/pswarp/src/pswarpMapGrid.c

    r12523 r18839  
    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}

  • trunk/pswarp/src/pswarpTransformReadout.c

    r12744 r18839  
    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            // allocate a job
     74            psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
     75
     76            // construct the arguments for this job
     77            // job is pswarpTransformTile (gridX, gridY);
     78            psArrayAdd (job->args, 1, args);
     79            // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);
     80
     81            // call: pswarpTransformTile (args);
     82            if (!psThreadJobAddPending (job)) {
     83                psError(PS_ERR_UNKNOWN, false, "Unable to warp image.");
     84                return false;
     85            }
     86            psFree (args);
     87        }
     88    }
     89
     90    // wait for the threads to finish and manage results
     91    // wait here for the threaded jobs to finish
     92    if (!psThreadPoolWait ()) {
     93        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     94        return false;
     95    }
     96    fprintf (stderr, "success for threaded jobs\n");
     97
     98    // each job records its own goodPixel values; sum them here
     99    // we have only supplied one type of job, so we can assume the types here
     100    psThreadJob *job = NULL;
     101    while ((job = psThreadJobGetDone()) != NULL) {
     102        if (job->args->n < 1) {
     103            fprintf (stderr, "error with job\n");
     104        } else {
     105            pswarpTransformTileArgs *args = job->args->data[0];
     106            // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
     107            goodPixels += args->goodPixels;
     108        }
     109        psFree (job);
     110    }
     111    psFree(interp);
     112    psFree(grid);
     113
     114
     115    // Store the variance factor and number of good pixels
     116    if (goodPixels > 0) {
     117        // Variance factor: large factor --> small scale
     118        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp);
     119        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
     120        if (vfItem) {
     121            psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX);
     122            psAssert(goodpixItem, "It should be where we left it!");
     123            psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64,
     124                     "Should be the type we said.");
     125
     126            vfItem->data.F32 += varFactor * goodPixels;
     127            goodpixItem->data.S64 += goodPixels;
     128        } else {
     129            psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0,
     130                             "Variance factor weighted by the good pixels", varFactor * goodPixels);
     131            psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0,
     132                             "Number of good pixels", goodPixels);
    69133        }
    70134    }
    71     psFree(interp);
    72     psFree(inPix);
    73     psFree(outPix);
    74     psFree(FP);
    75     psFree(TP);
    76     psFree(sky);
     135
     136    if (goodPixels > 0) {
     137        if (!pswarpTransformSources (output, input, config)) {
     138            psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
     139            return false;
     140        }
     141    }
     142
     143    // XXX should we not write anything out if there are no good pixels?
     144    output->data_exists = true;
     145    output->parent->data_exists = true;
     146    output->parent->parent->data_exists = true;
     147
     148    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));
    77149
    78150    return true;
Note: See TracChangeset for help on using the changeset viewer.