IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 8, 2008, 8:10:29 AM (18 years ago)
Author:
Paul Price
Message:

Adding threads to pmSubtractionCalculateEquation and pmSubtractionOrder. More work needs to be done in threading the final convolution.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r18287 r18962  
    1010#include "pmSubtractionKernels.h"
    1111#include "pmSubtractionStamps.h"
     12#include "pmSubtractionThreads.h"
    1213
    1314#include "pmSubtractionEquation.h"
    … …  
    479480//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    480481
    481 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
     482bool pmSubtractionCalculateEquationThread(const psThreadJob *job)
     483{
     484    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     485
     486    pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps
     487    const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     488    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
     489
     490    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
     491}
     492
     493bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     494                                         int index)
    482495{
    483496    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    484497    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     498    PS_ASSERT_INT_NONNEGATIVE(index, false);
     499    PS_ASSERT_INT_LESS_THAN(index, stamps->num, false);
    485500
    486501    int footprint = stamps->footprint;  // Half-size of stamps
    … …  
    494509    int numParams = numKernels * numSpatial + 1 + numBackground;
    495510
    496     psImage *polyValues = NULL;         // Polynomial terms
     511    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
     512    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
     513
     514    // Generate convolutions
     515    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
     516        psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index);
     517        return NULL;
     518    }
     519
     520#ifdef TESTING
     521    for (int j = 0; j < numKernels; j++) {
     522        if (stamp->convolutions1) {
     523            psString convName = NULL;
     524            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
     525            psFits *fits = psFitsOpen(convName, "w");
     526            psFree(convName);
     527            psKernel *conv = stamp->convolutions1->data[j];
     528            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
     529            psFitsClose(fits);
     530        }
     531
     532        if (stamp->convolutions2) {
     533            psString convName = NULL;
     534            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
     535            psFits *fits = psFitsOpen(convName, "w");
     536            psFree(convName);
     537            psKernel *conv = stamp->convolutions2->data[j];
     538            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
     539            psFitsClose(fits);
     540        }
     541    }
     542#endif
     543
     544    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
     545                                                    stamp->xNorm, stamp->yNorm); // Polynomial terms
     546
     547    bool new = stamp->vector1 ? false : true; // Is this a new run?
     548    if (new) {
     549        stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     550        stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
     551    }
     552#ifdef TESTING
     553    psImageInit(stamp->matrix1, NAN);
     554    psVectorInit(stamp->vector1, NAN);
     555#endif
     556
     557    bool status;                    // Status of least-squares matrix/vector calculation
     558    switch (kernels->mode) {
     559      case PM_SUBTRACTION_MODE_1:
     560        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
     561                                 stamp->weight, polyValues, footprint, true);
     562        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
     563                                  stamp->image2, stamp->weight, polyValues, footprint, true);
     564        break;
     565      case PM_SUBTRACTION_MODE_2:
     566        status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
     567                                 stamp->weight, polyValues, footprint, true);
     568        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
     569                                  stamp->image1, stamp->weight, polyValues, footprint, true);
     570        break;
     571      case PM_SUBTRACTION_MODE_DUAL:
     572        if (new) {
     573            stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
     574            stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
     575            stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
     576        }
     577#ifdef TESTING
     578        psImageInit(stamp->matrix2, NAN);
     579        psImageInit(stamp->matrixX, NAN);
     580        psVectorInit(stamp->vector2, NAN);
     581#endif
     582        status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
     583                                  stamp->weight, polyValues, footprint, true);
     584        status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
     585                                  stamp->weight, polyValues, footprint, false);
     586        status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
     587                                       stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
     588                                       footprint);
     589        status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
     590                                  stamp->image2, stamp->weight, polyValues, footprint, true);
     591        status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
     592                                  stamp->image2, stamp->weight, polyValues, footprint, false);
     593        break;
     594      default:
     595        psAbort("Unsupported subtraction mode: %x", kernels->mode);
     596    }
     597
     598    if (!status) {
     599        stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     600        psWarning("Rejecting stamp %d (%d,%d) because of bad equation",
     601                  index, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     602    } else {
     603        stamp->status = PM_SUBTRACTION_STAMP_USED;
     604    }
     605
     606#ifdef TESTING
     607    if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
     608        psString matrixName = NULL;
     609        psStringAppend(&matrixName, "matrix1_%d.fits", index);
     610        psFits *matrixFile = psFitsOpen(matrixName, "w");
     611        psFree(matrixName);
     612        psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
     613        psFitsClose(matrixFile);
     614
     615        matrixName = NULL;
     616        psStringAppend(&matrixName, "vector1_%d.fits", index);
     617        psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
     618        memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
     619               PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
     620        matrixFile = psFitsOpen(matrixName, "w");
     621        psFree(matrixName);
     622        psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
     623        psFree(dummy);
     624        psFitsClose(matrixFile);
     625
     626        if (stamp->vector2) {
     627            matrixName = NULL;
     628            psStringAppend(&matrixName, "vector2_%d.fits", index);
     629            dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
     630            memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
     631                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
     632            matrixFile = psFitsOpen(matrixName, "w");
     633            psFree(matrixName);
     634            psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
     635            psFree(dummy);
     636            psFitsClose(matrixFile);
     637        }
     638
     639        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     640            matrixName = NULL;
     641            psStringAppend(&matrixName, "matrix2_%d.fits", index);
     642            matrixFile = psFitsOpen(matrixName, "w");
     643            psFree(matrixName);
     644            psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
     645            psFitsClose(matrixFile);
     646
     647            matrixName = NULL;
     648            psStringAppend(&matrixName, "matrixX_%d.fits", index);
     649            matrixFile = psFitsOpen(matrixName, "w");
     650            psFree(matrixName);
     651            psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
     652            psFitsClose(matrixFile);
     653        }
     654    }
     655#endif
     656
     657    psFree(polyValues);
     658
     659    return true;
     660}
     661
     662bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
     663{
     664    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     665    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    497666
    498667    // We iterate over each stamp, allocate the matrix and vectors if
    … …  
    504673        }
    505674
    506         // Generate convolutions
    507         if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
    508             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
    509             psFree(polyValues);
    510             return NULL;
    511         }
    512 
    513 #ifdef TESTING
    514         for (int j = 0; j < numKernels; j++) {
    515             if (stamp->convolutions1) {
    516                 psString convName = NULL;
    517                 psStringAppend(&convName, "conv1_%03d_%03d.fits", i, j);
    518                 psFits *fits = psFitsOpen(convName, "w");
    519                 psFree(convName);
    520                 psKernel *conv = stamp->convolutions1->data[j];
    521                 psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
    522                 psFitsClose(fits);
    523             }
    524 
    525             if (stamp->convolutions2) {
    526                 psString convName = NULL;
    527                 psStringAppend(&convName, "conv2_%03d_%03d.fits", i, j);
    528                 psFits *fits = psFitsOpen(convName, "w");
    529                 psFree(convName);
    530                 psKernel *conv = stamp->convolutions2->data[j];
    531                 psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
    532                 psFitsClose(fits);
    533             }
    534         }
    535 #endif
    536 
    537         polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm);
    538 
    539 
    540             stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    541             stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
    542 #ifdef TESTING
    543             psImageInit(stamp->matrix1, NAN);
    544             psVectorInit(stamp->vector1, NAN);
    545 #endif
    546 
    547         bool status;                    // Status of least-squares matrix/vector calculation
    548         switch (kernels->mode) {
    549           case PM_SUBTRACTION_MODE_1:
    550             status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
    551                                      stamp->weight, polyValues, footprint, true);
    552             status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
    553                                       stamp->image2, stamp->weight, polyValues, footprint, true);
    554             break;
    555           case PM_SUBTRACTION_MODE_2:
    556             status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
    557                                      stamp->weight, polyValues, footprint, true);
    558             status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
    559                                       stamp->image1, stamp->weight, polyValues, footprint, true);
    560             break;
    561           case PM_SUBTRACTION_MODE_DUAL:
    562             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
    563             stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
    564             stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
    565 #ifdef TESTING
    566             psImageInit(stamp->matrix2, NAN);
    567             psImageInit(stamp->matrixX, NAN);
    568             psVectorInit(stamp->vector2, NAN);
    569 #endif
    570             status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
    571                                       stamp->weight, polyValues, footprint, true);
    572             status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
    573                                       stamp->weight, polyValues, footprint, false);
    574             status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
    575                                            stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
    576                                            footprint);
    577             status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
    578                                       stamp->image2, stamp->weight, polyValues, footprint, true);
    579             status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
    580                                       stamp->image2, stamp->weight, polyValues, footprint, false);
    581             break;
    582           default:
    583             psAbort("Unsupported subtraction mode: %x", kernels->mode);
    584         }
    585 
    586         if (!status) {
    587             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    588             psWarning("Rejecting stamp %d (%d,%d) because of bad equation",
    589                       i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     675        if (pmSubtractionThreaded()) {
     676            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION");
     677            psArrayAdd(job->args, 1, stamps);
     678            psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
     679            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
     680            if (!psThreadJobAddPending(job)) {
     681                psFree(job);
     682                return false;
     683            }
     684            psFree(job);
    590685        } else {
    591             stamp->status = PM_SUBTRACTION_STAMP_USED;
    592         }
    593 
    594 #ifdef TESTING
    595         if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
    596             psString matrixName = NULL;
    597             psStringAppend(&matrixName, "matrix1_%d.fits", i);
    598             psFits *matrixFile = psFitsOpen(matrixName, "w");
    599             psFree(matrixName);
    600             psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
    601             psFitsClose(matrixFile);
    602 
    603             matrixName = NULL;
    604             psStringAppend(&matrixName, "vector1_%d.fits", i);
    605             psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
    606             memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
    607                    PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
    608             matrixFile = psFitsOpen(matrixName, "w");
    609             psFree(matrixName);
    610             psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
    611             psFree(dummy);
    612             psFitsClose(matrixFile);
    613 
    614             if (stamp->vector2) {
    615                 matrixName = NULL;
    616                 psStringAppend(&matrixName, "vector2_%d.fits", i);
    617                 dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
    618                 memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
    619                        PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
    620                 matrixFile = psFitsOpen(matrixName, "w");
    621                 psFree(matrixName);
    622                 psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
    623                 psFree(dummy);
    624                 psFitsClose(matrixFile);
    625             }
    626 
    627             if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    628                 matrixName = NULL;
    629                 psStringAppend(&matrixName, "matrix2_%d.fits", i);
    630                 matrixFile = psFitsOpen(matrixName, "w");
    631                 psFree(matrixName);
    632                 psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
    633                 psFitsClose(matrixFile);
    634 
    635                 matrixName = NULL;
    636                 psStringAppend(&matrixName, "matrixX_%d.fits", i);
    637                 matrixFile = psFitsOpen(matrixName, "w");
    638                 psFree(matrixName);
    639                 psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
    640                 psFitsClose(matrixFile);
    641             }
    642         }
    643 #endif
    644 
    645     }
    646     psFree(polyValues);
     686            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
     687        }
     688    }
     689
     690    if (!psThreadPoolWait(true)) {
     691        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     692        return false;
     693    }
    647694
    648695    return true;
Note: See TracChangeset for help on using the changeset viewer.