Lecture
The classical technique for comparing a current image with a reference image is based on treating images as two-dimensional brightness functions (discrete two-dimensional intensity matrices). Here, either the distance between the images or a measure of their similarity is computed.
As a rule, the following formula is used to compute the distance between images, being the sum of the absolute values or squares of the intensity differences:
If, in addition to simply comparing two images, one needs to solve the problem of detecting the position of a fragment of one image within another, the classical "entry-level" method, which consists of iterating over all coordinates and computing the distance using the formula above, as a rule fails in practical use because of the large amount of computation required.
One of the methods that allows the amount of computation to be significantly reduced is the use of Fourier transforms and discrete Fourier transforms to calculate a measure of match between two images at various offsets relative to each other. In this case, the computation takes place simultaneously for different combinations of image shifts relative to one another.
For example, find:
By definition, the correlation <F,G> of two functions F and G is the quantity:
This quantity is well known from the course on linear spaces in mathematics and geometry, where it is called the dot product. We will use the following formula as a measure between images:
or
This quantity is obtained from the dot product of vectors (treating images as vectors in a multidimensional space). Moreover, this same formula also represents the standard statistical test formula for the hypothesis that two probability distributions coincide.
Note:
When computing the correlation between image fragments, if one image is smaller than the other, we will divide only by the value of the norms of the overlapping parts.
By definition, the convolution of two functions F and G is the function F×G:
It is also obvious that F×G'(t) is equal to the correlation obtained as a result of shifting one vector relative to the other by a step t (this is easy to verify by direct substitution of values into the correlation formula).
The Fourier transform (ℱ) is an operation that maps one function of a real variable to another function, also of a real variable. This new function describes the coefficients ("amplitudes") in the expansion of the original function into elementary components — harmonic oscillations with different frequencies.
The Fourier transform of a function f of a real variable is an integral transform given by the following formula:
Different sources may give definitions differing from the one above in the choice of the coefficient in front of the integral, as well as in the sign "−" in the exponent. But all the properties remain the same, although the form of some formulas may change.
In addition, there are various generalizations of this concept.
The Fourier transform of functions defined on the space ℝ^n is given by the formula:
The inverse transform in this case is given by the formula:
As before, in different sources the definitions of the multidimensional Fourier transform may differ in the choice of the constant in front of the integral.
The discrete Fourier transform (DFT in the English-language literature) is one of the Fourier transforms widely used in digital signal processing algorithms (its modifications are applied in MP3 audio compression, JPEG image compression, etc.), as well as in other fields related to frequency analysis of a discrete (for example, digitized analog) signal. The discrete Fourier transform requires a discrete function as input. Such functions are often created by sampling (taking values from continuous functions). Discrete Fourier transforms help solve partial differential equations and perform operations such as convolutions. Discrete Fourier transforms are also actively used in statistics, in time series analysis. Multidimensional discrete Fourier transforms also exist.
Forward transform:
Inverse transform:
The discrete Fourier transform is a linear transform that maps a vector of time samples to a vector of spectral samples of the same length. Thus the transform can be implemented as the multiplication of a symmetric square matrix by a vector:
One of the remarkable properties of Fourier transforms is the ability to quickly compute the correlation of two functions defined either on a real argument (using the classical formula) or on a finite ring (using discrete transforms).
And although similar properties are inherent to many linear transforms, for practical application, to compute the convolution operation, according to the definition we have given, the following formula is used
Where
Verifying the correctness of the equality is fairly easy – by explicitly substituting into the Fourier transform formulas and simplifying the resulting formulas
Let <F,G>(t) be equal to the correlation obtained as a result of shifting one vector relative to the other by a step t
Then, as already shown earlier, the following holds
If implementations of the Fourier transform algorithm using complex numbers are used, then such transforms have yet another remarkable property:
Where CONJUGATE ( FFT(G) ) – is the matrix made up of the conjugate elements of the matrix FFT(G)
Thus, we obtain
Using the formula for estimating the distance between images under a shift (i,j) relative to one another
we obtain that
Where
When solving the problem of searching for a single sample, additional normalization of the sample is unnecessary, and also the computation of the norm of the common part can be replaced by the sum of pixel brightness values in that common part, or by the sum of squares of brightness values in that common part
Using the formula for estimating the distance between images under a shift (i,j) relative to one another
we obtain that
Where
Note:
When using the discrete Fourier transform, matrix M also contains elements from a cyclic shift of the images relative to each other. Therefore, if it is not necessary to analyze cyclic frame shift, the search for the maximum element in matrix M should be limited to the region (0,0)-(N1-M1, N2-M2).
The implemented algorithms are part of the open-source library FFTTools. Web address: github.com/dprotopopov/FFTTools
Software used
/// <summary>
/// Catch pattern bitmap with the Fastest Fourier Transform
/// </summary>
/// <returns>Matrix of values</returns>
private Matrix<double> Catch(Image<Gray, double> 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<double>(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<double> 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<double> 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;
}
/// <summary>
/// Copy 3D array to 2D array (sizes can be different)
/// Flip copied data
/// Reduce last dimension
/// </summary>
/// <param name="input">Input array</param>
/// <param name="output">Output array</param>
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];
}
/// <summary>
/// Copy 3D array to 2D array (sizes can be different)
/// Replace items copied by value
/// Flip copied data
/// Reduce last dimension
/// </summary>
/// <param name="input">Input array</param>
/// <param name="output">Output array</param>
/// <param name="value">Value to replace copied data</param>
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;
}
/// <summary>
/// Find a maximum element in the matrix
/// </summary>
/// <param name="matrix">Matrix of values</param>
/// <param name="x">Index of maximum element</param>
/// <param name="y">Index of maximum element</param>
/// <param name="value">Value of maximum element</param>
public void Max(Matrix<double> 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;
}
}
}
/// <summary>
/// Catch pattern bitmap with the Fastest Fourier Transform
/// </summary>
/// <returns>Array of values</returns>
public Matrix<double> Catch(Bitmap bitmap)
{
using (var image = new Image<Gray, Byte>(bitmap))
return Catch(image);
}
Comments