[Numpy-discussion] shift for optimal superimposition of two 3D matrices according to correlation computed using FFT

[Zachary Pincus]

> I have two 3D density maps (meaning volumetric data, each one
> essentially a IxJxK matrix containing real values that sum to one) and
> want to find translation of between the two that maximises  
> correlation.
> This can be done by computing the correlation between the two
> (correlation theorem -> I use FFT for this) and then shifting one of  
> the
> two maps by the vector that corresponds to the index of the maximum in
> the 3D correlation matrix.
>
> My problem with this method is, however, that the shift wraps around  
> the
> edges of the matrix. e.g: for a 30x30x30 matrix that needs to be
> slightly shifted to fit the seconds 3D matrix, sometimes shift vectors
> like (0, 29, 0) appear, that actually mean a shift by (0, -1, 0). This
> means that the size of the matrix needs to be subtracted from some
> components of the shift vector in special cases.
>
> So far, I have been unable to determine the cutoff value, outside of
> which I need to subtract the matrix size from a component of the shift
> vector, as this cutoff seems to vary between different data sets.

Does it work to use a cutoff of half the size of the input arrays in  
each dimension? This is equivalent to calculating both shifts (the  
positive and negative) and using whichever has a smaller absolute value.

Alternately, you could use numpy.roll to shift the data (one axis at a  
time). Since roll wraps around, you wouldn't need to bother figuring  
out which shift is "correct".

Finally, you could not use FFTs but instead directly optimize a  
transformation between the two, using scipy.ndimage's affine  
transforms and scipy.optimize's numerical optimizers.

Zach