Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 21324) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 21335) @@ -7,6 +7,6 @@ /// @author Eugene Magnier, IfA /// -/// @version $Revision: 1.82 $ $Name: not supported by cvs2svn $ -/// @date $Date: 2009-01-27 06:39:37 $ +/// @version $Revision: 1.83 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2009-02-05 22:36:19 $ /// /// Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -131,9 +131,6 @@ -psKernel *psKernelGenerate(const psVector *tShifts, - const psVector *xShifts, - const psVector *yShifts, - float totalTime, - bool xyRelative) +psKernel *psKernelGenerate(const psVector *tShifts, const psVector *xShifts, const psVector *yShifts, + float totalTime, bool xyRelative) { PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL); @@ -232,7 +229,5 @@ } -psImage *psImageConvolveDirect(psImage *out, - const psImage *in, - const psKernel *kernel) +psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel) { PS_ASSERT_IMAGE_NON_NULL(in, NULL); @@ -466,7 +461,39 @@ } -bool psImageSmooth(psImage *image, - double sigma, - double Nsigma) + +// Generate normalised Gaussian vector +#define IMAGE_SMOOTH_GAUSS(TARGET, SIZE, SIGMA, TYPE) \ + psVector *TARGET = psVectorAlloc(2 * SIZE + 1, PS_TYPE_##TYPE); /* Gaussian */ \ + double sum = 0.0; /* Sum of Gaussian, for normalisation */ \ + double factor = -0.5/PS_SQR(SIGMA); /* Multiplier for exponential */ \ + for (int i = -SIZE, j = 0; i <= SIZE; i++, j++) { \ + sum += TARGET->data.TYPE[j] = exp(factor * PS_SQR(i)); \ + } \ + for (int i = 0; i <= 2 * SIZE + 1; i++) { \ + TARGET->data.TYPE[i] /= sum; \ + } + +psKernel *psImageSmoothKernel(float sigma, float nSigma) +{ + PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(nSigma, 0.0, NULL); + + int size = sigma * nSigma + 0.5; // Half-size of kernel + psKernel *kernel = psKernelAlloc(-size, size, -size, size); // Kernel to return + + IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32); + psF32 *gauss = &gaussNorm->data.F32[size]; // Convenient kernel-like indexing for Gaussian + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + kernel->kernel[y][x] = gauss[y] * gauss[x]; + } + } + psFree(gaussNorm); + + return kernel; +} + + +bool psImageSmooth(psImage *image, double sigma, double Nsigma) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -474,23 +501,11 @@ // relevant terms int Nrange = sigma*Nsigma + 0.5; // Number of pixels either side for convolution kernel - int Npixel = 2*Nrange + 1; // Total number of pixels in convolution kernel int Nx = image->numCols; // Number of columns int Ny = image->numRows; // Number of rows - #define IMAGESMOOTH_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ +#define IMAGESMOOTH_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ /* generate normalized gaussian */ \ - psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \ - { \ - double sum = 0.0; \ - double factor = -0.5/(sigma*sigma); \ - for (int i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \ - sum += gaussnorm->data.TYPE[i+Nrange]; \ - } \ - for (int i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] /= sum; \ - } \ - } \ + IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, TYPE) \ ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \ \ @@ -711,11 +726,6 @@ } -psImage *psImageSmoothMask(psImage *output, - const psImage *image, - const psImage *mask, - psImageMaskType maskVal, - float sigma, - float numSigma, - float minGauss) +psImage *psImageSmoothMask(psImage *output, const psImage *image, const psImage *mask, + psImageMaskType maskVal, float sigma, float numSigma, float minGauss) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -726,24 +736,8 @@ int numCols = image->numCols, numRows = image->numRows; // Size of image int xLast = numCols - 1, yLast = numRows - 1; // Last index - - // relevant terms int size = sigma * numSigma + 0.5; // Half-size of kernel - int fullSize = 2*size + 1; // Total number of pixels in convolution kernel // Generate normalized gaussian - psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian - { - double sum = 0.0; - double norm = -0.5/(sigma*sigma); - for (int i = 0, u = -size; i <= fullSize; i++, u++) { - float value = expf(norm*PS_SQR(u)); - gaussNorm->data.F32[i] = value; - sum += value; - } - sum = 1.0 / sum; - for (int i = 0; i <= fullSize; i++) { - gaussNorm->data.F32[i] *= sum; - } - } + IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32); const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel @@ -828,14 +822,7 @@ } -bool psImageSmoothMask_ScanRows_F32 (psImage *calculation, - psImage *calcMask, - const psImage *image, - const psImage *mask, - psImageMaskType maskVal, - psVector *gaussNorm, - float minGauss, - int size, - int rowStart, - int rowStop) +bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image, + const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm, + float minGauss, int size, int rowStart, int rowStop) { const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel @@ -873,13 +860,7 @@ } -bool psImageSmoothMask_ScanCols_F32 (psImage *output, - psImage *calculation, - psImage *calcMask, - psImageMaskType maskVal, - psVector *gaussNorm, - float minGauss, - int size, - int colStart, - int colStop) +bool psImageSmoothMask_ScanCols_F32(psImage *output, psImage *calculation, psImage *calcMask, + psImageMaskType maskVal, psVector *gaussNorm, float minGauss, + int size, int colStart, int colStop) { const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel @@ -953,11 +934,6 @@ } -psImage *psImageSmoothMask_Threaded(psImage *output, - const psImage *image, - const psImage *mask, - psImageMaskType maskVal, - float sigma, - float numSigma, - float minGauss) +psImage *psImageSmoothMask_Threaded(psImage *output, const psImage *image, const psImage *mask, + psImageMaskType maskVal, float sigma, float numSigma, float minGauss) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -971,24 +947,8 @@ int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols; int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows; - - // relevant terms int size = sigma * numSigma + 0.5; // Half-size of kernel - int fullSize = 2*size + 1; // Total number of pixels in convolution kernel // Generate normalized gaussian - psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian - { - double sum = 0.0; - double norm = -0.5/(sigma*sigma); - for (int i = 0, u = -size; i <= fullSize; i++, u++) { - float value = expf(norm*PS_SQR(u)); - gaussNorm->data.F32[i] = value; - sum += value; - } - sum = 1.0 / sum; - for (int i = 0; i <= fullSize; i++) { - gaussNorm->data.F32[i] *= sum; - } - } + IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32); switch (image->type.type) { @@ -1089,9 +1049,6 @@ -bool psImageSmoothMaskF32 (psImage *image, - psImage *mask, - psImageMaskType maskVal, - double sigma, - double Nsigma) +bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal, + double sigma, double Nsigma) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -1100,21 +1057,9 @@ // relevant terms int Nrange = sigma*Nsigma + 0.5; // Number of pixels either side for convolution kernel - int Npixel = 2*Nrange + 1; // Total number of pixels in convolution kernel int Nx = image->numCols; // Number of columns int Ny = image->numRows; // Number of rows /* generate normalized gaussian */ - psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_F32); - { - double sum = 0.0; - double factor = -0.5/(sigma*sigma); - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] = exp(factor*i*i); - sum += gaussnorm->data.F32[i+Nrange]; - } - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] /= sum; - } - } + IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, F32); psF32 *gauss = &gaussnorm->data.F32[Nrange]; @@ -1148,10 +1093,4 @@ psFree(calculation); - // XXX test - if (0) { - psFree(gaussnorm); - return true; - } - /** Smooth in Y direction **/ // allocate and save Nrange extra row vectors for storage. these will be Index: trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.h (revision 21324) +++ trunk/psLib/src/imageops/psImageConvolve.h (revision 21335) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-28 22:08:14 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 22:36:19 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -195,4 +195,10 @@ ); +/// Return the kernel used for smoothing +psKernel *psImageSmoothKernel( + float sigma, ///< Width of the smoothing kernel, pixels + float nSigma ///< Size of the smoothing box, sigma + ); + /// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for masked pixels ///