Speed up SampledDensityMap by taking advantage of separability

[cgreenberg]

The value of a (diagonal) 3D gaussian at some voxel is the product of the x,y, and z 1D gaussians. That means for an NxNxN grid, you don't have to compute N^3 exponentials, you only have to compute 3N of them. Maybe even fewer if you take into account symmetry.

[cgreenberg]

For the GaussianKernel, we could rewrite the internal_resample loop to calculate the kernel at each voxel plane and just multiply them. But the SphereKernel isn't separable - at the corners you can have X < R, Y < R, but X^2+Y^2>R (e.g. if X=Y=R/sqrt(2)). So we may need to get rid of the whole (elegant but maybe silly) template-based kernel thing.