(t) = FxG '(t) = BFT (FFT (F) * FFT (G'))

FFT (G ') = CONJUGATE (FFT (G))

(t) = BFT (FFT (F) * CONJUGATE (FFT (G)))

Fourier transforms to solve the problem

m (X, Y) (i, j) = (i, j) / (| X | (i, j)) * | Y | (i, j)),

• = XxY '= BFT (FFT (X) * CONJUGATE (FFT (Y)))
• | X | ^ 2 = = XxX'xE '= BFT (FFT (X) * CONJUGATE (FFT (X)) * CONJUGATE (FFT (E))) = BFT (SQUAREMAGNITUDE (FFT (X)) * CONJUGATE (FFT (E)))
• | Y | ^ 2 = = YxY'xE '= BFT (FFT (Y) * CONJUGATE (FFT (Y)) * CONJUGATE (FFT (E)))

• (i, j) is the scalar product of two images obtained by shifting (i, j) relative to each other images X and Y
• E is an image of size equal to the minimum sizes of X and Y, and filled with unit values ​​(that is, a “frame” in which X and Y are compared)
• | X | (i, j) is the norm of the general part of the image X with a shift (i, j)
• | Y | (i, j) is the norm of the common part of the image Y when shifting (i, j)
• FFT - operation of direct two-dimensional discrete Fourier transform
• BFT - operation of the inverse two-dimensional discrete Fourier transform
• CONJUGATE - operation of calculating a matrix of conjugate elements
• SQUAREMAGNITUDE– operation of calculating the matrix of squares of element amplitudes

Simplification of formulas to solve the problem

m (X, Y) (i, j) = (i, j) / | X | ^ 2 (i, j),

• = BFT (FFT (X) * CONJUGATE (FFT (Y)))
• = BFT (SQUAREMAGNITUDE (FFT (X)) * CONJUGATE (FFT (E)))

• (i, j) is the scalar product of two images obtained by shifting (i, j) relative to each other images X and Y
• E is an image of size equal to the minimum sizes of X and Y, and filled with unit values ​​(that is, a “frame” in which X and Y are compared)
• (i, j) is the norm (the sum of the brightness of pixels) of the total part of the image X when shifting (i, j)
• FFT - operation of direct two-dimensional discrete Fourier transform
• BFT - operation of the inverse two-dimensional discrete Fourier transform
• CONJUGATE - operation of calculating a matrix of conjugate elements
• SQUAREMAGNITUDE– operation of calculating the matrix of squares of element amplitudes

Fragment search algorithm in full image

• Let two images X and Y be given - an image and a sample, of sizes (N1, N2) and (M1, M2), respectively, and Ni> Mi
• It is required to find the coordinates of sample Y in the full image X and calculate the estimated value - a measure of proximity.

1. Extend image Y to size (N1, N2), padding it with zeros
2. Form an image E from units of size (M1, M2) and expand to size (N1, N2), adding zeros to it
3. Calculate = BFT (FFT (X) * CONJUGATE (FFT (Y)))
4. Calculate = BFT (SQUAREMAGNITUDE (FFT (X)) * CONJUGATE (FFT (E)))
5. Calculate M [i, j] = (f + [i, j]) / (f + [i, j])
6. In the matrix M, find the element with the maximum value - the coordinates of this element are the desired position of the sample in the full image, and the value is equal to the estimate of the measure of comparison.

Implementation examples

• Microsoft Visual Studio 2013 C # - environment and programming language
• EmguCV / OpenCV - C ++ library of structures and algorithms for image processing
• FFTWSharp / FFTW - C ++ library implementing fast discrete Fourier transform algorithms

/// 

        ///     Catch pattern bitmap with the Fastest Fourier Transform
        /// 
        /// Matrix of values
        private Matrix Catch(Image image)
        {
            const double f = 1.0;
            int length = image.Data.Length;
            int n0 = image.Data.GetLength(0);
            int n1 = image.Data.GetLength(1);
            int n2 = image.Data.GetLength(2);
            Debug.Assert(n2 == 1);
            // Allocate FFTW structures
            var input = new fftw_complexarray(length);
            var output = new fftw_complexarray(length);
            fftw_plan forward = fftw_plan.dft_3d(n0, n1, n2, input, output,
                fftw_direction.Forward,
                fftw_flags.Estimate);
            fftw_plan backward = fftw_plan.dft_3d(n0, n1, n2, input, output,
                fftw_direction.Backward,
                fftw_flags.Estimate);
            var matrix = new Matrix(n0, n1);
            double[,,] patternData = _patternImage.Data;
            double[,,] imageData = image.Data;
            double[,] data = matrix.Data;
            var doubles = new double[length];
            // Calculate Divisor
            Copy(patternData, data);
            Buffer.BlockCopy(data, 0, doubles, 0, length*sizeof (double));
            input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray());
            forward.Execute();
            Complex[] complex = output.GetData_Complex();
            Buffer.BlockCopy(imageData, 0, doubles, 0, length*sizeof (double));
            input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray());
            forward.Execute();
            input.SetData(output.GetData_Complex().Zip(complex, (x, y) => x*Complex.Conjugate(y)).ToArray());
            backward.Execute();
            IEnumerable doubles1 = output.GetData_Complex().Select(x => x.Magnitude);
            if (_fastMode)
            {
                // Fast Result
                Buffer.BlockCopy(doubles1.ToArray(), 0, data, 0, length*sizeof (double));
                return matrix;
            }
            // Calculate Divider (aka Power)
            input.SetData(doubles.Select(x => new Complex(x*x, 0)).ToArray());
            forward.Execute();
            complex = output.GetData_Complex();
            CopyAndReplace(_patternImage.Data, data);
            Buffer.BlockCopy(data, 0, doubles, 0, length*sizeof (double));
            input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray());
            forward.Execute();
            input.SetData(complex.Zip(output.GetData_Complex(), (x, y) => x*Complex.Conjugate(y)).ToArray());
            backward.Execute();
            IEnumerable doubles2 = output.GetData_Complex().Select(x => x.Magnitude);
            // Result
            Buffer.BlockCopy(doubles1.Zip(doubles2, (x, y) => (f + x*x)/(f + y)).ToArray(), 0, data, 0,
                length*sizeof (double));
            return matrix;
        }

/// 

        ///     Copy 3D array to 2D array (sizes can be different)
        ///     Flip copied data
        ///     Reduce last dimension
        /// 
        /// Input array
        /// Output array
        private static void Copy(double[,,] input, double[,] output)
        {
            int n0 = output.GetLength(0);
            int n1 = output.GetLength(1);
            int m0 = Math.Min(n0, input.GetLength(0));
            int m1 = Math.Min(n1, input.GetLength(1));
            int m2 = input.GetLength(2);
            for (int i = 0; i < m0; i++)
                for (int j = 0; j < m1; j++)
                    output[i, j] = input[i, j, 0];
            for (int k = 1; k < m2; k++)
                for (int i = 0; i < m0; i++)
                    for (int j = 0; j < m1; j++)
                        output[i, j] += input[i, j, k];
        }
        /// 

        ///     Copy 3D array to 2D array (sizes can be different)
        ///     Replace items copied by value
        ///     Flip copied data
        ///     Reduce last dimension
        /// 
        /// Input array
        /// Output array
        /// Value to replace copied data
        private static void CopyAndReplace(double[,,] input, double[,] output, double value = 1.0)
        {
            int n0 = output.GetLength(0);
            int n1 = output.GetLength(1);
            int m0 = Math.Min(n0, input.GetLength(0));
            int m1 = Math.Min(n1, input.GetLength(1));
            int m2 = input.GetLength(2);
            for (int i = 0; i < m0; i++)
                for (int j = 0; j < m1; j++)
                    output[i, j] = value;
        }
        /// 

        ///     Find a maximum element in the matrix
        /// 
        /// Matrix of values
        /// Index of maximum element
        /// Index of maximum element
        /// Value of maximum element
        public void Max(Matrix matrix, out int x, out int y, out double value)
        {
            double[,] data = matrix.Data;
            int n0 = data.GetLength(0);
            int n1 = data.GetLength(1);
            value = data[0, 0];
            x = y = 0;
            for (int i = 0; i < n0; i++)
            {
                for (int j = 0; j < n1; j++)
                {
                    if (data[i, j] < value) continue;
                    value = data[i, j];
                    x = j;
                    y = i;
                }
            }
        }
        /// 

        ///     Catch pattern bitmap with the Fastest Fourier Transform
        /// 
        /// Array of values
        public Matrix Catch(Bitmap bitmap)
        {
            using (var image = new Image(bitmap))
                return Catch(image);
        }

Caught that biting

Literature

1. A.L. Dmitriev. Optical methods of information processing. Tutorial. SPb. SPYUGUITMO 2005.46 p.
2. A.A. Akayev, S.A. Mayorov “Optical methods of information processing” M.: 1988
3. J. Goodman "Introduction to Fourier Optics" M .: World 1970