Index: /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.c =================================================================== --- /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.c (revision 35634) +++ /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.c (revision 35635) @@ -1231,4 +1231,6 @@ } +/*********************** smooth with mask, threaded version ***********************/ + bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image, const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm, @@ -1484,4 +1486,224 @@ +/*********************** smooth no-mask ***********************/ + +bool psImageSmoothNoMask_ScanRows_F32(psImage *calculation, const psImage *image, + psVector *gaussNorm, float minGauss, + int size, int rowStart, int rowStop) +{ + const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel + + int numCols = image->numCols; + int xLast = numCols - 1; // Last index + + for (int j = rowStart; j < rowStop; j++) { + for (int i = 0; i < numCols; i++) { + int xMin = PS_MAX(i - size, 0); + int xMax = PS_MIN(i + size, xLast); + const psF32 *imageData = &image->data.F32[j][xMin]; + int uMin = - PS_MIN(i, size); /* Minimum kernel index */ + const psF32 *gaussData = &gauss[uMin]; + double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ + for (int x = xMin; x <= xMax; x++, imageData++, gaussData++) { + sumIG += *imageData * *gaussData; + sumG += *gaussData; + } + if (sumG > minGauss) { + /* BW */ + calculation->data.F32[i][j] = sumIG / sumG; + } + } + } + return true; +} + +bool psImageSmoothNoMask_ScanCols_F32(psImage *output, psImage *calculation, + psVector *gaussNorm, float minGauss, + int size, int colStart, int colStop) +{ + const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel + + int numRows = output->numRows; + int yLast = numRows - 1; // Last index + + for (int i = colStart; i < colStop; i++) { + for (int j = 0; j < numRows; j++) { + int yMin = PS_MAX(j - size, 0); + int yMax = PS_MIN(j + size, yLast); + const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */ + int vMin = - PS_MIN(j, size); /* Minimum kernel index */ + const psF32 *gaussData = &gauss[vMin]; + double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ + for (int y = yMin; y <= yMax; y++, imageData++, gaussData++) { + sumIG += *imageData * *gaussData; + sumG += *gaussData; + } + output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN; + } + } + return true; +} + +bool psImageSmoothNoMask_ScanRows_F32_Threaded(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert (job->args, "programming error: job args not alloced?"); + psAssert (job->args->n == 7, "programming error: job Nargs mismatch"); + + psImage *calculation = job->args->data[0]; // calculation image + const psImage *image = job->args->data[1]; // input image + + psVector *gaussNorm = job->args->data[2]; // gauss kernel + float minGauss = PS_SCALAR_VALUE(job->args->data[3],F32); + int size = PS_SCALAR_VALUE(job->args->data[4],S32); + int rowStart = PS_SCALAR_VALUE(job->args->data[5],S32); + int rowStop = PS_SCALAR_VALUE(job->args->data[6],S32); + return psImageSmoothNoMask_ScanRows_F32(calculation, image, + gaussNorm, minGauss, size, + rowStart, rowStop); +} + +bool psImageSmoothNoMask_ScanCols_F32_Threaded(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert (job->args, "programming error: job args not alloced?"); + psAssert (job->args->n == 7, "programming error: job Nargs mismatch"); + + psImage *output = job->args->data[0]; // input image + psImage *calculation = job->args->data[1]; // calculation image + + psVector *gaussNorm = job->args->data[2]; // gauss kernel + float minGauss = PS_SCALAR_VALUE(job->args->data[3],F32); + int size = PS_SCALAR_VALUE(job->args->data[4],S32); + int colStart = PS_SCALAR_VALUE(job->args->data[5],S32); + int colStop = PS_SCALAR_VALUE(job->args->data[6],S32); + return psImageSmoothNoMask_ScanCols_F32(output, calculation, + gaussNorm, minGauss, size, + colStart, colStop); +} + +psImage *psImageSmoothNoMask_Threaded(psImage *output, const psImage *image, float sigma, float numSigma, float minGauss) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + + int numCols = image->numCols, numRows = image->numRows; // Size of image + + int nThreads = psThreadPoolSize(); + int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols; + int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows; + int size = sigma * numSigma + 0.5; // Half-size of kernel + + // Generate normalized gaussian + psVector *gaussNorm = NULL; + IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32); + + switch (image->type.type) { + // Smooth an image with masked pixels + // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or + // [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW" + + case PS_TYPE_F32: { + psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */ + + if (threaded) { + /** Smooth in X direction **/ + for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) { + int rowStop = PS_MIN (rowStart + scanRows, numRows); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS"); + psArrayAdd(job->args, 0, calculation); + psArrayAdd(job->args, 0, (psImage *) image); // cast away const + psArrayAdd(job->args, 0, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStop, PS_TYPE_S32); + // -> psImageSmoothNoMask_ScanRows_F32 (calculation, image, gauss, minGauss, size, rowStart, rowStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } else if (!psImageSmoothNoMask_ScanRows_F32(calculation, image, gaussNorm, minGauss, size, 0, numRows)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + + output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32); + + /** Smooth in Y direction **/ + if (threaded) { + for (int colStart = 0; colStart < numCols; colStart+=scanCols) { + int colStop = PS_MIN (colStart + scanCols, numCols); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS"); + psArrayAdd(job->args, 0, output); + psArrayAdd(job->args, 0, calculation); + psArrayAdd(job->args, 0, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStop, PS_TYPE_S32); + // -> psImageSmoothMask_ScanCols_F32 (output, calculation, gauss, minGauss, size, colStart, colStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } + + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } else if (!psImageSmoothNoMask_ScanCols_F32(output, calculation, gaussNorm, minGauss, size, 0, numCols)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + + psFree(calculation); + psFree(gaussNorm); + + return output; + } + default: + psFree(gaussNorm); + char *typeStr; + PS_TYPE_NAME(typeStr,image->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Specified psImage type, %s, is not supported."), + typeStr); + return false; + } + + psAbort("Should never reach here."); + return false; +} + +/****************************************************************************************/ + bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal, double sigma, double Nsigma) @@ -1900,4 +2122,16 @@ } { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS", 7); + task->function = &psImageSmoothNoMask_ScanRows_F32_Threaded; + psThreadTaskAdd(task); + psFree(task); + } + { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS", 7); + task->function = &psImageSmoothNoMask_ScanCols_F32_Threaded; + psThreadTaskAdd(task); + psFree(task); + } + { psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9); task->function = &psImageSmoothMaskPixelsThread; @@ -1909,4 +2143,6 @@ psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS"); psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS"); + psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS"); + psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS"); psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_PIXELS"); } Index: /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.h =================================================================== --- /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.h (revision 35634) +++ /branches/eam_branches/ipp-20130509/psLib/src/imageops/psImageConvolve.h (revision 35635) @@ -242,5 +242,6 @@ ); -/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for masked pixels +/// Smooth an imageby parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS /// /// Applies a circularly symmetric Gaussian smoothing first in x and then in y @@ -256,5 +257,5 @@ ); -/// Smooth particular pixels on an image, allowing for masked pixels +/// Smooth particular pixels on an image, allowing for MASKED PIXELS /// /// Applies a circularly symmetric Gaussian smoothing first in x and then in y @@ -271,4 +272,9 @@ ); +/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS : THREADED VERSION +/// +/// Applies a circularly symmetric Gaussian smoothing first in x and then in y +/// directions with just a vector. This process is 2N faster than 2D convolutions (in general). psImage *psImageSmoothMask_Threaded(psImage *output, const psImage *image, @@ -279,4 +285,17 @@ float minGauss); +/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS : THREADED VERSION +/// +/// Applies a circularly symmetric Gaussian smoothing first in x and then in y +/// directions with just a vector. This process is 2N faster than 2D convolutions (in general). +psImage *psImageSmoothNoMask_Threaded(psImage *output, + const psImage *image, + float sigma, + float numSigma, + float minGauss); + +/// Smooth an image by parts IN-SITU using 1D Gaussian independently in x and y, allowing +/// for MASKED PIXELS bool psImageSmoothMaskF32( psImage *image, ///< the image to be smoothed