IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 5, 2012, 5:19:48 PM (14 years ago)
Author:
mhuber
Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:
branches/meh_branches/ppstack_test
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/meh_branches/ppstack_test

  • branches/meh_branches/ppstack_test/psLib/src/imageops/psImageConvolve.c

    r30595 r33415  
    569569// Generate normalised Gaussian vector
    570570#define IMAGE_SMOOTH_GAUSS(TARGET, SIZE, SIGMA, TYPE) \
    571     psVector *TARGET = psVectorAlloc(2 * SIZE + 1, PS_TYPE_##TYPE); /* Gaussian */ \
     571    TARGET = psVectorAlloc(2 * SIZE + 1, PS_TYPE_##TYPE); /* Gaussian */ \
    572572    double sum = 0.0; /* Sum of Gaussian, for normalisation */ \
    573573    double factor = -0.5/PS_SQR(SIGMA); /* Multiplier for exponential */ \
    … …  
    587587    psKernel *kernel = psKernelAlloc(-size, size, -size, size); // Kernel to return
    588588
     589    psVector *gaussNorm = NULL;
    589590    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    590591    psF32 *gauss = &gaussNorm->data.F32[size]; // Convenient kernel-like indexing for Gaussian
    … …  
    602603bool psImageSmooth(psImage *image, double  sigma, double  Nsigma)
    603604{
    604     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     605    PS_ASSERT_IMAGE_NON_NULL(image, false);
    605606
    606607    // relevant terms
    … …  
    612613  case PS_TYPE_##TYPE: { \
    613614        /* generate normalized gaussian */ \
     615        psVector *gaussnorm = NULL; \
    614616        IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, TYPE) \
    615617        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
    … …  
    622624                ps##TYPE *vo = calculation->data.TYPE; \
    623625                int xMax = PS_MIN(Nrange, Nx); \
    624                 int convRange = PS_MIN(Nrange + 1, Nx); \
    625626                /* Smooth first Nrange pixels, with renorm */ \
    626627                for (int i = 0; i < xMax; i++, vi++, vo++) { \
     628                    int convRange = PS_MIN(Nrange + 1, Nx - i);  \
    627629                    ps##TYPE *vr = vi - i; \
    628630                    ps##TYPE *vg = gauss - i; \
    … …  
    648650                    } \
    649651                    /* Smooth last Nrange pixels, with renorm */ \
    650                     /* XXX does this miss the last column? */ \
    651                     for (int i = Nx - Nrange; i < Nx; i++, vi++, vo++) { \
     652                    /* if Nx < 2*Nrange, this pass starts at i == Nrange */ \
     653                    int xMin = PS_MAX(Nx - Nrange, Nrange); \
     654                    for (int i = xMin; i < Nx; i++, vi++, vo++) { \
    652655                        ps##TYPE *vr = vi - Nrange; \
    653656                        ps##TYPE *vg = gauss - Nrange; \
    … …  
    670673        /* Smooth the first Nrange pixels, with renorm */ \
    671674        int yMax = PS_MIN(Nrange, Ny); \
    672         int convRange = PS_MIN(Nrange + 1, Ny); \
    673675        for (int j = 0; j < yMax; j++) { \
     676            int convRange = PS_MIN(Nrange + 1, Ny - j);                \
    674677            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
    675678            /* Zero the output row */ \
    … …  
    715718                calculation = save; \
    716719            } \
    717             /* Smooth last Nrange pixels, with renorm */ \
    718             for (int j = Ny - Nrange; j < Ny; j++) { \
     720            /* Smooth last Nrange pixels, with renorm */                \
     721            /* if Ny < 2*Nrange, this pass starts at j == Nrange */ \
     722            int yMin = PS_MAX(Ny - Nrange, Nrange); \
     723            for (int j = yMin; j < Ny; j++) { \
    719724                /* save the Nrange-offset output row, then zero */ \
    720725                memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
    … …  
    765770}
    766771
     772void psImageSmooth_PreAlloc_DataFree (psImageSmooth_PreAlloc_Data *smdata) {
     773    psFree (smdata->resultX);
     774    psFree (smdata->resultY);
     775    psFree (smdata->kernel);
     776}
     777
     778psImageSmooth_PreAlloc_Data *psImageSmooth_PreAlloc_DataAlloc (psImage *image, double sigma, double Nsigma) {
     779
     780    psImageSmooth_PreAlloc_Data *smdata = psAlloc(sizeof(psImageSmooth_PreAlloc_Data));
     781    psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth_PreAlloc_DataFree);
     782
     783    if (!image) {
     784        // relevant terms
     785        smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
     786        smdata->Nx = 0;
     787        smdata->Ny = 0;
     788        smdata->kernel = NULL;
     789        smdata->resultX = NULL;
     790        smdata->resultY = NULL;
     791        return smdata;
     792    }
     793
     794    // relevant terms
     795    smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
     796    smdata->Nx = image->numCols;            // Number of columns
     797    smdata->Ny = image->numRows;            // Number of rows
     798
     799    IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32);
     800       
     801    // use a temp running buffer for X and Y directions.
     802    smdata->resultX = psAlloc(smdata->Nx * sizeof(psF32));
     803    memset (smdata->resultX, 0, smdata->Nx*sizeof(psF32));
     804
     805    smdata->resultY = psAlloc(smdata->Ny * sizeof(psF32));
     806    memset (smdata->resultY, 0, smdata->Ny*sizeof(psF32));
     807
     808    return smdata;
     809}
     810
     811// we can use the same DATA structure on multiple images of the same size
     812bool psImageSmooth_PreAlloc_F32(psImage *image, psImageSmooth_PreAlloc_Data *smdata)
     813{
     814    PS_ASSERT_IMAGE_NON_NULL(image, false);
     815    // assert on data type
     816
     817    // relevant terms
     818    int Nrange = smdata->Nrange;    // Number of pixels either side for convolution kernel
     819    int Nx = smdata->Nx;            // Number of columns
     820    int Ny = smdata->Ny;            // Number of rows
     821
     822    psF32 *gauss = &smdata->kernel->data.F32[Nrange];
     823    psF32 *resultX = smdata->resultX;
     824    psF32 *resultY = smdata->resultY;
     825       
     826    /* Smooth in X direction */
     827    {
     828        for (int j = 0; j < Ny; j++) {
     829            psF32 *vi = image->data.F32[j];
     830            int xMax = PS_MIN(Nrange, Nx);
     831            /* Smooth first Nrange pixels, with renorm */
     832            for (int i = 0; i < xMax; i++, vi++) {
     833                int convRange = PS_MIN(Nrange + 1, Nx - i);
     834                psF32 *vr = vi - i;
     835                psF32 *vg = gauss - i;
     836                double g = 0.0;
     837                double s = 0.0;
     838                for (int n = -i; n < convRange; n++, vr++, vg++) {
     839                    s += *vg * *vr;
     840                    g += *vg;
     841                }
     842                resultX[i] = s / g;
     843            }
     844            /* If that's all the pixels we have, then we're done already */
     845            if (Nx > Nrange) {
     846                /* Smooth middle pixels; if Nx < 2*Nrange, this pass is skipped */
     847                for (int i = Nrange; i < Nx - Nrange; i++, vi++) {
     848                    psF32 *vr = vi - Nrange;
     849                    psF32 *vg = gauss - Nrange;
     850                    double s = 0;
     851                    for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) {
     852                        s += *vg * *vr;
     853                    }
     854                    resultX[i] = s;
     855                }
     856                /* Smooth last Nrange pixels, with renorm */
     857                // if Nx < 2*Nrange, this pass starts at i == Nrange
     858                int xMin = PS_MAX(Nx - Nrange, Nrange);
     859                for (int i = xMin; i < Nx; i++, vi++) {
     860                    psF32 *vr = vi - Nrange;
     861                    psF32 *vg = gauss - Nrange;
     862                    double g = 0.0;
     863                    double s = 0.0;
     864                    for (int n = -Nrange; n < Nx - i; n++, vr++, vg++) {
     865                        s += *vg * *vr;
     866                        g += *vg;
     867                    }
     868                    resultX[i] = s / g;
     869                }
     870            }
     871            memcpy(image->data.F32[j], resultX, Nx*sizeof(psF32));
     872        }
     873    }
     874       
     875    // this section probably hits the cache poorly for large images, but is probably OK for small ones
     876    /* Smooth in Y direction */
     877    {
     878        for (int i = 0; i < Nx; i++) {
     879            int yMax = PS_MIN(Nrange, Ny);
     880            /* Smooth first Nrange pixels, with renorm */
     881            for (int j = 0; j < yMax; j++) {
     882                int convRange = PS_MIN(Nrange + 1, Ny - j);
     883                psF32 *vg = gauss - j;
     884                double g = 0.0;
     885                double s = 0.0;
     886                for (int n = -j; n < convRange; n++, vg++) {
     887                    psF32 vr = image->data.F32[j+n][i];
     888                    s += *vg * vr;
     889                    g += *vg;
     890                }
     891                resultY[j] = s / g;
     892            }
     893            /* If that's all the pixels we have, then we're done already */
     894            if (Ny > Nrange) {
     895                /* Smooth middle pixels */
     896                for (int j = Nrange; j < Ny - Nrange; j++) {
     897                    psF32 *vg = gauss - Nrange;
     898                    double s = 0;
     899                    for (int n = -Nrange; n < Nrange + 1; n++, vg++) {
     900                        psF32 vr = image->data.F32[j+n][i];
     901                        s += *vg * vr;
     902                    }
     903                    resultY[j] = s;
     904                }
     905                /* Smooth last Nrange pixels, with renorm */
     906                // if Ny < 2*Nrange, this pass starts at j == Nrange
     907                int yMin = PS_MAX(Ny - Nrange, Nrange);
     908                for (int j = yMin; j < Ny; j++) {
     909                    psF32 *vg = gauss - Nrange;
     910                    double g = 0.0;
     911                    double s = 0.0;
     912                    for (int n = -Nrange; n < Ny - j; n++, vg++) {
     913                        psF32 vr = image->data.F32[j+n][i];
     914                        s += *vg * vr;
     915                        g += *vg;
     916                    }
     917                    resultY[j] = s / g;
     918                }
     919            }
     920            // loop here
     921            for (int j = 0; j < Ny; j++) {
     922                image->data.F32[j][i] = resultY[j];
     923            }
     924        }
     925    }
     926    return true;
     927}
     928
    767929static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
    768930                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
    … …  
    8491011
    8501012    // Generate normalized gaussian
     1013    psVector *gaussNorm = NULL;
    8511014    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    8521015
    … …  
    8591022
    8601023            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
    861             psArrayAdd(job->args, 1, out);
    862             psArrayAdd(job->args, 1, (psImage*)image);
    863             psArrayAdd(job->args, 1, (psImage*)mask);
     1024            psArrayAdd(job->args, 0, out);
     1025            psArrayAdd(job->args, 0, (psImage*)image);
     1026            psArrayAdd(job->args, 0, (psImage*)mask);
    8641027            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
    865             psArrayAdd(job->args, 1, (psVector*)x);
    866             psArrayAdd(job->args, 1, (psVector*)y);
    867             psArrayAdd(job->args, 1, gaussNorm);
     1028            psArrayAdd(job->args, 0, (psVector*)x);
     1029            psArrayAdd(job->args, 0, (psVector*)y);
     1030            psArrayAdd(job->args, 0, gaussNorm);
    8681031            PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
    8691032            PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32);
    … …  
    8761039            }
    8771040        }
    878         if (!psThreadPoolWait(true)) {
     1041        if (!psThreadPoolWait(true, true)) {
    8791042            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
    8801043            psFree(gaussNorm);
    … …  
    9741137
    9751138    // Generate normalized gaussian
     1139    psVector *gaussNorm = NULL;
    9761140    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    9771141    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
    … …  
    11851349
    11861350    // Generate normalized gaussian
     1351    psVector *gaussNorm = NULL;
    11871352    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    11881353
    … …  
    12031368                  // allocate a job, construct the arguments for this job
    12041369                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
    1205                   psArrayAdd(job->args, 1, calculation);
    1206                   psArrayAdd(job->args, 1, calcMask);
    1207                   psArrayAdd(job->args, 1, (psImage *) image); // cast away const
    1208                   psArrayAdd(job->args, 1, (psImage *) mask); // cast away const
     1370                  psArrayAdd(job->args, 0, calculation);
     1371                  psArrayAdd(job->args, 0, calcMask);
     1372                  psArrayAdd(job->args, 0, (psImage *) image); // cast away const
     1373                  psArrayAdd(job->args, 0, (psImage *) mask); // cast away const
    12091374                  PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    1210                   psArrayAdd(job->args, 1, gaussNorm);
     1375                  psArrayAdd(job->args, 0, gaussNorm);
    12111376                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
    12121377                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
    … …  
    12251390              }
    12261391              // wait here for the threaded jobs to finish (NOP if threading is not active)
    1227               if (!psThreadPoolWait(true)) {
     1392              if (!psThreadPoolWait(true, true)) {
    12281393                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
    12291394                  psFree(calculation);
    … …  
    12501415                  // allocate a job, construct the arguments for this job
    12511416                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
    1252                   psArrayAdd(job->args, 1, output);
    1253                   psArrayAdd(job->args, 1, calculation);
    1254                   psArrayAdd(job->args, 1, calcMask);
     1417                  psArrayAdd(job->args, 0, output);
     1418                  psArrayAdd(job->args, 0, calculation);
     1419                  psArrayAdd(job->args, 0, calcMask);
    12551420                  PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    1256                   psArrayAdd(job->args, 1, gaussNorm);
     1421                  psArrayAdd(job->args, 0, gaussNorm);
    12571422                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
    12581423                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
    … …  
    12721437
    12731438              // wait here for the threaded jobs to finish (NOP if threading is not active)
    1274               if (!psThreadPoolWait(true)) {
     1439              if (!psThreadPoolWait(true, true)) {
    12751440                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
    12761441                  psFree(calculation);
    … …  
    13211486
    13221487    /* generate normalized gaussian */
     1488    psVector *gaussnorm = NULL;
    13231489    IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, F32);
    13241490    psF32 *gauss = &gaussnorm->data.F32[Nrange];
    … …  
    15721738
    15731739            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
    1574             psArrayAdd(job->args, 1, conv);
    1575             psArrayAdd(job->args, 1, (psImage*)mask); // Casting away const to put on arguments
     1740            psArrayAdd(job->args, 0, conv);
     1741            psArrayAdd(job->args, 0, (psImage*)mask); // Casting away const to put on arguments
    15761742            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
    15771743            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
    … …  
    15861752            }
    15871753        }
    1588         if (!psThreadPoolWait(true)) {
     1754        if (!psThreadPoolWait(true, true)) {
    15891755            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
    15901756            psFree(conv);
    … …  
    16071773
    16081774            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
    1609             psArrayAdd(job->args, 1, out);
    1610             psArrayAdd(job->args, 1, conv);
     1775            psArrayAdd(job->args, 0, out);
     1776            psArrayAdd(job->args, 0, conv);
    16111777            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
    16121778            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
    … …  
    16211787            }
    16221788        }
    1623         if (!psThreadPoolWait(true)) {
     1789        if (!psThreadPoolWait(true, true)) {
    16241790            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
    16251791            psFree(conv);
Note: See TracChangeset for help on using the changeset viewer.