Spent some time thinking about gridding convolution functions (gcf) (see "Synthesis Imaging in Radio Astronomy" chapter 7, section 3.3).
A gcf is a numerical trick to mitigate aliasing of power outside the image back into it when preforming a gridded Fourier Transform (an FFT). By convolving the visibilties onto the grid when preparing for the transform we are multiplying by a tapering function in the image plane. The taper continues outside of the image so the aliased power folded back in is attenuated more than the proper part of the image which we want. We then divide by the FT of the gcf to correct for the taper which we artificially introduced. The result is closer to the direct Fourier Transform which we would like to do, but which is ussually too slow.
gcf's are ussually seperable functions - ie. C(u,v)=C(u)C(v). The 2D FT of a seperable function is the product of the 1D FT's - ie. c(u,v)=c(u)c(v).
My initial make_dmap.m and Nils mapplot.m both assigned each vis to the closest grid point. This is implicitly a "square pillbox" gcf with half side length del_u/2, whose 1D FT is a sinc function. It is a very poor gcf. But we had neglected to make the gcf correction. In this case it actually makes the disagreement worse by scaling up the outer region of the image.
I looked at the Difmap source code. Martin uses a truncated Gaussian gcf over a 5x5 grid and 0.7 HWHM. sigma=HWHM/sqrt(2*log(2))=0.59, w=sigma*sqrt(2)=0.84. This works OK, again diverging from the direct FT result with increasing distance from field center. The reason the images I read in from Difmap get so close to the direct FT is that they are generated on a double size grid, but only the central quarter is output. If I do the same thing with my function I get equally close.
SIIRA generalizes the truncated Gaussian gcf (calling it an Exponential). They recommend w=1 (w=1 where w=sigma*sqrt(2)), and a 6x6 box. I have implemented an odd sized box algorithm like Martin, where the box runs from -ngcf to +ngcf, and used ngcf=3 for most tests.
I notice that I can push the region where the result diverges from the direct FT further and further out by increasing w and ngcf. eg: 2 and 10 agree almost perfectly except for a small rectangular border around the edge. Doing this is exploiting the double precision which Matlab uses to the max. (Going further it breaks.)
SIIRA mentions using a sinc function which works quite well.
SIIRA also describes an "Exponential times sinc" function which works very well. I have put this in make_dmap.m for now. It does as well as exp with w=2 and ngcf=10 but using a much smaller ngcf=3 box (and not pushing precision limits). Due to what seems to be the very high overhead of each Matlab function call the exp-sinc is not actually any faster than the plain exp with the larger box. However, if I get around to implementing the loop over visibilities in C then it presumably it would be much faster...
SIIRA mentions that there are better ways to get the FT of the gcf than the FFT that I am using now in all cases except the pillbox. Martin appears to compute the FT of his truncated Gaussian numerically (he uses a cosine transform). I am not sure if my exp, sinc, and exp-sinc algorithms would be any better with smarter calculations of the gcf FT...
Here is a figure showing DASI field A6 restricted to long baselines to enhance point source detection. The top two frames show the result of a direct FT and the exp-sinc which seems to be the best of the above functions. The lower four panels show the difference between the direct FT and the other gcf's considered all using a +/-3 grid point truncation box.