Image Convolution in Frequency Domain (MATLAB)

I'm writing a c++ program for a class to do convolution in the frequency domain, and I noticed the final result had error along the corner. So I tried it out in MATLAB and got the exact same results. e.g.

Using cameraman from http://engronline.ee.memphis.edu/e开发者_JS百科ece7214/images/Downlodable.htm

I did

a = imread('cameraman.pgm');
h = ones(25,25)/25/25;
a(512,512) = 0;
h(512,512) = 0;
c = ifft2(fft2(a).*fft2(h))/256;
c = c(1:256, 1:256);
c = real(c);
imwrite(c,'test2.png')

I took a peek at c before extracting the top-left corner, and I found it was the same answer as imfilter(a, h) except it was translated a tad from the corner. Digital Image Processing by Gonzalez says nothing about this, and Google has left my eyes bleeding from searching with no help at all (everyone repeats the same instructions Gonzalez gives to extract the top left corner).

Unrelated to the main question, I'd also like to know why I had to divide by 256 in this MATLAB code. In my C++ code, I didn't have to scale the result, and I got the same answer as this MATLAB code.

edit: I did a little playing around with 1-dimensional vectors (doing conv and ifft(fft*fft)), and I think the 'error' is from the output showing the 'full' convolution in the top left corner instead of the 'same' convolution. But even if that were the case, I'm not sure how to deterministically code "extract this portion only to get the 'same' instead of the top left 256x256 portion of the 'full'"

edit: More googling has resulted in a possible solution via http://jeremy.fix.free.fr/IMG/pdf/fftconvolution.pdf. It has a lot of math symbols I've never seen before, but from what I can gather, if you are convolving an nxn and a mxm, extract m:(m+n-1) to get the 'same' convolution from the fft approximation. I'd still like to hear from someone who is more expert than I am, so don't choose not to comment based on this update!

---

MATLAB indexes arrays from 1 to N.
Canonical C usage (and most C matrix math libraries) indexes array from 0 to N-1. That may be the source of an off-by-one difference.

To be an identity, the bare ifft(fft(x)) algorithm requires a scale factor of 1/N somewhere. Some C libraries put that 1/N inside the fft(). Many build that scale factor into the ifft(). Some put a 1/sqrt(N) in both. Some add no built-in scale factor, requiring the user to scale as needed to get an identity.

---

You need to shift it by half of the window size. In this case you should do: c = c(13:256+12,13:256+12); instead of: c = c(1:256, 1:256);

You are padding the signal with zeros but not getting any use of that since you take the signal from 1 to 256 - you still have the thrash on the edges and the valuable signal on the right corners which you lose.