There are several possibilities: firstly, perhaps you didn't make your filter's frequency-domain definition Hermitian-symmetric. It is straightforward to show that the inverse DFT of a sequence is only purely real if it has Hermitian symmetry.

To understand why this is the case, have a look at the definition of the DFT:

x[n] = 1/N sum (k = 0 to N - 1) X(k) exp (j * 2 * pi * n * k / N)

If you make X(k) Hermitian symmetric, the imaginary parts in the summation all cancel out, leaving only real-valued x[n]. If you're not totally clear on this, try playing with the DFT of some real-valued sequences in Matlab and have a look at the real and imaginary components; also try it by hand for some short sequences.

If you have done that and there are still non-zero imaginary components which are all much smaller than the real component (say by a factor of 10^6), then, the problem is due to roundoff errors and you can just use real (h) (for a filter point spread function h).

The best way to do 2D frequency-domain filter design is shown here:

http://www.mathworks.com/help/images/designing-linear-filters-in-the-frequency-domain.html