COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

### 2D Convolution using FFTW

• Email
• Follow

```Hello everyone,

I am implementing a fast convolution program using FFTW library.  I
understand the constraints of the Discrete Convolution Theorem, and
the issues of wrap around error and periodicity.

Hence, I am taking the following approach at the moment.

/*-------------------*/

Image: image[sizeX][sizeY]
(Symmetrical)Filter: kernel[2k + 1][2k + 1]

Step(1) Define padX = sizeX + k; padY = sizeY + k;

Step (2) Pad image on two sides with zeros so that its size is padX *
padY.  I pad zeros on the right side and the bottom side of the image
buffer.

Step (3) Store kernel in a wrap around fashion with the size padX *

So, if Symmetrical (eg. Gaussian) kernel = [a, b, c
d, e, f
g, h, i]

newKernel = [a, b, 0, 0, 0, 0, 0, c
d, e, 0, 0, 0, 0, 0, f
0, 0, 0, 0, 0, 0, 0, 0
....
....
g,h, 0, 0, 0, 0, 0, i]

Step (4) Take FFT of newKernel and the image padded with zeros
Step (5) Do pixel by pixel complex multiplication of FFT output.
Step (6) Take inverse FFT of the output
Step (7) Read in the output

/*---------*/

Well, this only works for small kernel sizes, because when the kernel
size is bigger I get a shift in the blurred image.

I tried various other ways but I haven't been successful in resolving

The C++ code for the FFT functions is as follows:
(1)

/*-------Building Gaussian FFT--------------------*/

/*----Create Gaussian Kernel---------*/
inline void ConvolutionOld::CreateGaussianKernel()
{
float d;
int counter;
float normalise;
int i;
counter = 0;
int y,x;
normalise = 0;

///////////////////////////////
//-2 -1 0 1 2
//-2  x
//-1  x
// 0  x
// 1  x
// 2  x
///////////////////////////////
//e ^ -(x*x + y*y)/2 * sigma * sigma
////////////////////////////////

//We only consider cases with k.kernelSide as odd

for(y = -((k.kernelSide - 1)/2); y <= (k.kernelSide - 1)/2; y++)
{
for(x = -((k.kernelSide - 1)/2); x <= (k.kernelSide - 1)/2; x++)
{
d = (float)(x*x + y*y)/(float)(2.0f * k.sigma * k.sigma);
kernel[counter] = pow(e,-d);
normalise += kernel[counter];
counter++;
}
}

for(i = 0; i < k.kernelSide*k.kernelSide; i++)
kernel[i] = kernel[i]/normalise;

}

inline void ConvolutionOld::ComputeKernelFFT()
{
float *kernelRow;
float *kernelColumn;
int y,x;

kernelRow = new float[padX * k.kernelSide];

/*---------Perform padding with zeros in wrap around
fashion-------*/

for(y = 0; y < padY; y++)
{
for(x = 0; x < padX; x++)
{
kernel_FFT[y][x].im = 0;
}
}

//Compute the FFT of the Gaussian Kernel
ComputeFFT(kernel_FFT, -1);

//Free the kernels that are no longer required
delete [] kernelRow;
delete [] kernelColumn;

}

/*--------------Wrap the kernel along rows-----------------*/
*kernelRow)
{
int y,x;

for(y = 0; y < k.kernelSide; y++)
{
for(x = 0; x <= (k.kernelSide - 1)/2; x++)
kernelRow[x + y*padX] = kernel[x + y*k.kernelSide];

for(x = (k.kernelSide - 1)/2 + 1; x < padX - (k.kernelSide - 1)/2;
x++)
kernelRow[x + y * padX] = 0;

for(x = 1; x <= (k.kernelSide - 1)/2 ; x++)
y*k.kernelSide];
}
}

/*---------------Wrap the kernelRow along columns--------------*/
*kernelColumn)
{
int y,x;

for(y = 0; y <= (k.kernelSide - 1)/2; y++)
{
for(x = 0; x < padX; x++)
}

for(y = (k.kernelSide - 1)/2 + 1; y < padY - (k.kernelSide - 1)/2;
y++)
{
for(x = 0; x < padX; x++)
}

for(y = 1; y <= (k.kernelSide - 1)/2; y++)
{
for(x = 0; x < padX; x++)
kernelColumn[x + (padY - y) * padX] = kernel[x + (k.kernelSide -
}
}

(2) /*-----Function to compute FFT using FFTW
library----------------*/

inline void ConvolutionOld::ComputeFFT(zomplex** array, int inv)
{
fftw_plan p;
int y,x;
fftw_complex *carray;

carray = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) *

for(y = 0; y < padY; y++)
{
for(x = 0; x < padX; x++)
{
carray[x + y * padX][0] = array[y][x].re;
carray[x + y * padX][1] = array[y][x].im;
}
}

fftw_execute(p);

for(y = 0; y < padY; y++)
{
for(x = 0; x < padX; x++)
{
array[y][x].re = carray[x + y * padX][0];
array[y][x].im = carray[x + y * padX][1];
}
}

fftw_destroy_plan(p);
fftw_free(carray);
}

(3)/*----------Function to compute image FFT-------------*/

/*----------Pad the image with zeros of the right and
bottom--------------*/
{
int y,x;

for(y = 0; y < sizeY; y++)
{
for(x = 0; x < sizeX; x++)
{
image_FFT[y][x].re = image[x + y*sizeX];
image_FFT[y][x].im = 0;
}

for(x = sizeX; x < padX; x++)
{
image_FFT[y][x].re = 0.0f;
image_FFT[y][x].im = 0.0f;
}
}

for(y = sizeY; y < padY; y++)
{

for(x = 0; x < padX; x++)
{
image_FFT[y][x].re = 0.0f;
image_FFT[y][x].im = 0.0f;
}
}

ComputeFFT(image_FFT, -1);

}

(4) /*------Function to perform convolution and take inverse FFT----*/
inline void ConvolutionOld::Convolve(float *output)
{

int i,y,x;

//Do per pixel multiplication of the gaussian kernel and the image
matrix
for(y = 0; y < padY; y++)
{
for(x = 0; x < padX; x++)
{
convolution_FFT[y][x].re = (image_FFT[y][x].re *
kernel_FFT[y][x].re) - (image_FFT[y][x].im * kernel_FFT[y][x].im);
convolution_FFT[y][x].im = (image_FFT[y][x].re *
kernel_FFT[y][x].im) + (image_FFT[y][x].im * kernel_FFT[y][x].re);
}

}

ComputeFFT(convolution_FFT, 1); //Find inverse FFT

for(y = 0; y < sizeY; y++)
{
for(x = 0; x < sizeX; x++)
output[x + y*sizeX] = convolution_FFT[y][x].re;

}

}

This is a lot of code to post, but I am unable to understand the cause
of the error.  If someone has implemented 2D convolution using FFTW
me.

-Swati
```
 0

See related articles to this posting

```>(4) /*------Function to perform convolution and take inverse FFT----*/
>inline void ConvolutionOld::Convolve(float *output)
>{
>
>	int i,y,x;
>
>	//Do per pixel multiplication of the gaussian kernel and the image
>matrix
>	for(y = 0; y < padY; y++)
>	{
>		for(x = 0; x < padX; x++)
>		{
>			convolution_FFT[y][x].re = (image_FFT[y][x].re *
>kernel_FFT[y][x].re) - (image_FFT[y][x].im * kernel_FFT[y][x].im);
>			convolution_FFT[y][x].im = (image_FFT[y][x].re *
>kernel_FFT[y][x].im) + (image_FFT[y][x].im * kernel_FFT[y][x].re);
>		}
>
>	}
>
>	ComputeFFT(convolution_FFT, 1); //Find inverse FFT
>
>
>	for(y = 0; y < sizeY; y++)
>	{
>		for(x = 0; x < sizeX; x++)
>			output[x + y*sizeX] = convolution_FFT[y][x].re;
>
>	}
>
>}
>
>
>
>
>This is a lot of code to post, but I am unable to understand the cause
>of the error.  If someone has implemented 2D convolution using FFTW
>me.
>
>-Swati
>

You need to scale in your Convolution function -

scale = 1.0 / (sizeX * sizeY);
convolution_FFT[y][x].re = (image_FFT[y][x].re * kernel_FFT[y][x].re)
(image_FFT[y][x].im * kernel_FFT[y][x].im) * scale;
convolution_FFT[y][x].im = (image_FFT[y][x].re * kernel_FFT[y][x].im)
(image_FFT[y][x].im * kernel_FFT[y][x].re) * scale;

Otherwise this looks good. Sorry no one got to this sooner, thought I'
reply in case someone found it who was having a similar issue. More o
fftw:

http://www.fftw.org/fftw2_doc/fftw_2.html

- Pat

```
 0
Reply patorjk (1) 11/27/2005 2:17:30 PM

1 Replies
600 Views

Similar Articles

12/6/2013 11:12:15 AM
page loaded in 61796 ms. (0)

Similar Artilces:

using SQLalchemy
We have a very chaotic database (on MySql) at the moment, with for example table names used as keys to query other tables (but that's just an example). We are going to redesign it but first I still have to replace the perl+vbscript system with only one Python program, so I still have to deal with the current state. I'm trying to use SQLalchemy and it looks absolutely great, but in general as a policy we don't use external dependencies.. To try to do an exception in this case: - are there any problems with SQLalchemy on Windows? - are there any possibly drawbacks of using SQLalchemy instead of the MySqlDB interface? For the second point I guess that we might have a bit less fine tuning, but the amount of data is not so much and speed is also not a bit issue (also because all the queries are extremely inefficient now). Any other possible issue? Thanks, Andrea

Using the mouse in a display file
keyword? Should the display file be > compiled with some specific parameters? > > Ewout > I know about the RTNCSRLOC keyword. However, as I said in my previous message, nothing seems to happen when I double-click with the mouse. Control does not return to the program. Is it a setting in the emulator configuration? We are using iSeries Access for Windows (or whatever the name might be at the moment). Ewout You need to enable hotspots. In iSeries Access V5R3, it's under the "Edit" Menu -> "Preferences" -> "Hotspots..." All the options on the dialog that comes up is pretty self evident. -- Bob Comer "Ewout" <ewout.boter@planet.nl> wrote in message news:1132776793.016456.296790@g44g2000cwa.googlegroups.com... >I know about the RTNCSRLOC keyword. However, as I said in my previous > message, nothing seems to happen when I double-click with the mouse. > Control does not return to the program. Is it a setting in the emulator > configuration? We are using iSeries Access for Windows (or whatever the > name might be at the moment). > > Ewout >

quantum visualizations using IDL
I recently finished documenting some IDL procedures and functions that calculate solutions to the Schrodinger equation of quantum mechanics for a number of different situations (e.g. for different potentials, in one and two dimensions, time-independent, time-dependent). For those of you in attendance at the IDL Users' Group meeting in Boulder a few weeks ago, I demonstrated a number of these programs during my presentation. Since a few of you at the meeting expressed interest in the code and documentation, you can find it at the following URL: http://www.ncnr.nist.gov/staff/dimeo/qvis/s

Using xml.xpath question.
Hi, I am trying to get the value of child from xmlstr = """<p:root xmlns:p="http://tempuri.org/string"><p:child DataType="String">Hellpppp</p:child></p:root>""" using doc=parseString(xmlstr) nodeList = xml.xpath.Evaluate("/p:root/p:child/text()", doc) and am getting the following exception: xml.xpath.RuntimeException: Undefined namespace prefix: "p". I am using python 2.2.2 with PyXML 0.8 I think my usage of the Evaluate is wrong or incomplete. I have tried to google for information but to no avail. Can anyone please shed light on this. Thanks Bipin bvelkur@yahoo.com (0wl) wrote in message news:<de5c7cc0.0307101347.2c707bf8@posting.google.com>... > Hi, > > I am trying to get the value of child from > > xmlstr = """<p:root xmlns:p="http://tempuri.org/string"><p:child > DataType="String">Hellpppp</p:child></p:root>""" > > using > doc=parseString(xmlstr) > nodeList = xml.xpath.Evaluate("/p:root/p:child/text()", doc) > > and am getting the following

using glasspane to disable screen
Hi, I have an application that I would like to disable for user input while showing a wait cursor. It should be disabled every time an action event is processed from a button and I have tried using a glasspane for this purpose. The main component for the application is a JWindow. In the action event I have code for disabling the screen and showing the wait cursor, executing the event of the clicked button and finally enabling the screen again and changing back to the default cursor. The cursor changes all right to the wait cursor and everything looks fine, the problem is that it is still...() == blueButton) { blueButton_actionPerformed(); } else if(e.getSource() == pinkButton) { pinkButton_actionPerformed(); } } } On 11/24/2004 at 5:41:39 AM, Kristina Engdahl wrote: > I have an application that I would like to disable for user input > while showing a wait cursor. It should be disabled every time an > action event is processed from a button and I have tried using a > glasspane for this purpose. There is really no need to explicitly disable anything for actions that take place entirely in the event handler. That is because the UI

using S-Functions #2
I am having a problem using S-Function, I want to use S-function, and i have a perfect error free file for my S-function but i don't know how to attach this file with S-function Can any one tell me how to attach files with S-Functions Thanks Regards Perkash

Big problem by using OpenMP
Dear all, I'm using OpenMP to parallelize a part of an old Fortran modul to increase the performance (to make it quicker). Within the Fortran modul I create 4 threads by using the OpenMP direktive "SECTIONS". The 4 sections are independent of each other. But after activating (switching on) OpenMP in the development environment the Fortran modul don't provide the same results like without OpenMP. I get this strange behaviour with and without the "SECTIONS" direktive, it's only necessary to enable (switching on) OpenMP. If I disable (switching off) OpenMP.....@googlemail.com> wrote: > Dear all, > > I'm using OpenMP to parallelize a part of an old Fortran modul to > increase the performance (to make it quicker). > Within the Fortran modul I create 4 threads by using the OpenMP > direktive "SECTIONS". > The 4 sections are independent of each other. > But after activating (switching on) OpenMP in the development > environment the Fortran modul don't provide the same results like > without OpenMP. > I get this strange behaviour with and without the "SECTIONS" > direktive, it's only

Using Compound Aritmetic function
will be true, so if you don't care how many are below 450 (what I'm calling the TRUE state), the output will be what you want, just use a case statement with your task in the TRUE case. Dave. Sorry if it was a bit muddled, I think I know what to do now. Use AND instead of OR. Then use the false statement of the case for my task when all the sensors readings are not more than 450. Here are two functions you might apply to your project. One is the Array Max & Min Function. Using this function use the min value and compare to your min set point (450). The other function you might use is the In Range And Coerce. Using this method you can set a min and max value 450 - 500 (do this first, then wire your array to the X input.) The output will be an array of Booleans all will be true unless they fall out of the range you specify. A second output will be an array of the inputs except any that fall out of the range you specify will be coerced to the min or max value. Example: 465;480;420;475;505 will return T;T;F;T;F and output 465;480;450;475;500. There are many variations where these two functions may be used. They may be found in the comparison or the array

image comparison using histogram
hi,,, i m doing a project on electronic voting matchine with image comparison technology. using MATLAB. for that image comparison i have to compare a image of face of a person with other images which are stored in database.so would u give me algorithm for that? i preffred a histogram matching technique for image comparison. rushi patel: For obvious reasons, that won't work, as I'm sure you've realized immediately upon clicking the "Send" button. It won't be nearly robust enough. Don't you think that the histograms from all people from a certain race will have pretty much the same histogram since they have pretty much the same skin color? Histograms may be one small part of face recognition but is by no means sufficient. I suggest you start here: http://iris.usc.edu/Vision-Notes/bibliography/contentspeople.html#Face%20Recognition,%20Detection,%20Tracking,%20Gesture%20Recognition,%20Fingerprints,%20Biometrics That will show you how published, robust, sophisticated, successful, and accurate systems do it. Pick one that seems reasonable and give it a try. -ImageAnalyst ImageAnalyst <imageanalyst@mailinator.com> wrote