What is the Fastest Way Ever to Do 2D Correlation? Here is what i gathered so far:


Any sliding window classification, image filtering or similar can be fastly done by a FFT (flip the signal and do convolution).
Also, if the template/filter kernel is separable it is possible to do fast correlation by simply separating into multiple kernels and applying then sequentialy.
I think it is possible to do fast correlation by using integral images, but I am not sure how.


The question is: is FFT the fastest way? If now, which is? Consider both GPU and CPU implementations

If it is not: is there algorithm capable of separating a non-separable kernel (in aproximated way, logically)?

Context:

My filter is an Inner Product Detector, which is found by

$$
h = R^{-1}\mu_1
$$
where $ R^{-1} $ is the autocorrelation matrix and $\mu_1$ is the mean of a given class.

The matrix $ R $ is hermitian. And can be calculated by

$$
\sum x^Tx
$$

where $x$ are d-dimensional vectors and $[.]^T$ denotes the transpose.

Also: kernels are squares of size 21,23,25,27 or 35
Images are 128x128
In your case, since you have multiple images while you have a given set of kernels the DFT based Correlation would be the best fit.

Pay attention that the DFT Based Convolution / Correlation Is Equivalent to Convolution / Correlation with Periodic / Circular Boundary Conditions.
It means that if you need different boundary conditions (Like padding with Nearest Neighbor) then it means you need to pad the image large enough so the cyclic effect won't hurt you.

In that case you might be better doing, given enough memory, what's being don in DNN.
Apply the Correlation / Convolution with Matrix Multiplication.  

Matrix Multiplication is so optimized in our days that even if its complexity is higher, in real world it is faster in many cases.
So you build the Convolution Matrix (You can even chain them) and the images matrix and multiply. I can show a MATLAB Example.

By the way, even if you do it with DFT in the GPU, I suggest you do it with a DNN library (Something which utilize CUDNN, ROCm - Chainer, PyTorch, TensorFlow, MATLAB Neural Network Toolbox, Knet, etc...) since this case of multiple images (Batch) with given kernel is the case optimized in those libraries with automatic mechanism to chose between DFT based or Matrix based implementation.

Good reference - Algorithms for Efficient Computation of Convolution by Karas Pavel and Svoboda David.