I'm doing something relatively simple, like this:
// image and kernel are both 512 x 512
Matrix<Complex> image = GetImage();
Matrix<Complex> kernel = GetKernel();
// Apply FFT to both image and kernel
Fourier.Forward2D(image, FourierOptions.NoScaling);
Fourier.Forward2D(kernel, FourierOptions.NoScaling);
// Multiply matrices to do convolution in frequency domain
Matrix<Complex> convolved = image.PointwiseMultiply(kernel);
// Apply inverse FFT to get image in spatial domain
Fourier.Inverse2D(convolved, FourierOptions.NoScaling);
However, the resulting image has its quadrants swapped. If you were to divide the original image into quadrants like this:
-------
| 1 | 2 |
|-------|
| 3 | 4 |
-------
The final image has the quadrants like this:
-------
| 4 | 3 |
|-------|
| 2 | 1 |
-------
The kernel is just an image with all zeros, except the central pixel (at 256, 256) with a value of 1.0.
Your problem is one of encountering a sign convention in the way that FFTs are computed. By putting your single pixel at (256,256) you have in effect said that you wanted to shift the whole image by +256,+256.
If you want your image stay put then you should centre your PSF on pixel (0,0) with its other parts lurking in the corners of the image. The alternative is to shift the resulting convolved image after doing the transform.
You have to remember that there are implicit periodic boundary conditions in the FFT. For this reason you might want a guard band of zeroes around your image that is just a little wider than half the width of the convolving PSF. Otherwise edge discontinuities will leak from top to bottom, left to right and vice versa.