Фур'є-обчислення для порівняння зображень

Традиційна техніка «початкового рівня», порівняння поточного зображення з еталоном ґрунтується на розгляді зображень як двовимірних функцій яскравості (дискретних двовимірних матриць інтенсивності). При цьому вимірюється або відстань між зображеннями, або міра їх близькості.
Як правило, для обчислення відстаней між зображеннями використовується формула, яка є сумою модулів або квадратів різниць інтенсивності
d(X,Y) = SUM ( X[i,j] — Y[i,j] )^2
Якщо крім простого порівняння двох зображень потрібно вирішити задачу виявлення позиції фрагмента одного зображення в інше, то класичний метод «початкового рівня», що полягає в переборі всіх координат та обчислення відстані за вказаною формулою, як правило, зазнає невдачі практичного використання через необхідного великої кількості обчислень.
Одним з методів, що дозволяють значно скоротити кількість обчислень, є застосування Фур'є перетворень і дискретних перетворень Фур'є для розрахунку міри збігу двох зображень при різних зсувах їх між собою. Обчислення при цьому відбуваються одночасно для різних комбінацій зрушень зображень відносно один одного.
Наявність великої кількості бібліотек, які реалізують Фур'є перетворень (у всіляких варіантах швидких версій), робить реалізацію алгоритмів порівняння зображень не дуже складним завданням для програмування.


Постановка завдання
  • Нехай дано два зображення X і Y – зображення і зразок, розмірів (N1,N2) і (M1,M2) відповідно і Ni > Mi
  • Потрібно знайти координати зразка Y в повному зображенні X і обчислити оціночну величину — міру близькості.
Наприклад, знайти
зразок
зразок
зображення
Чеширський кіт

Кореляція як міра між зображеннями
Згідно з визначенням, кореляцією <F,G> двох функцій F і G називається величина
<F,G> = SUM F(i)*G(i)
Ця величина добре відома з курсу математики і геометрії, присвяченого лінійним простором, де носить назву скалярного твори.
Будемо використовувати в якості запобіжного між зображеннями формулу
m(X,Y) = SUM ( X[i,j] * Y[i,j] ) / ( SQRT ( SUM X[i,j] ^2 ) * SQRT ( SUM Y[i,j] ^2 ) )
або
m(X,Y) = <X,Y> / ( SQRT (<X,X> ) * SQRT (<Y,Y> ) )
Дана величина отримана з операції скалярного добутку векторів (розглядаючи зображення як вектори у багатовимірному просторі). І навіть більше — ця ж формула являє собою стандартну статистичну формулу критерію для гіпотези про збіг двох імовірнісних розподілів.
Примітка
При обчисленні кореляції між фрагментами зображень, якщо одне зображення менше іншого, будемо ділити лише на значення норм у перетинається частин.

Згортка двох функцій
Згідно з визначенням, згорткою двох функцій F і G називається функція FхG
FхG(t) = SUM F(i)*G(i)|i+j=t
Нехай G'(t) = G(-t) і F'(t) = F(-t), тоді, очевидна справедливість рівностей
  • FхF'(0) = SUM F(i)^2 – скалярний добуток вектора F на самого себе
  • GхG'(0) = SUM G(j)^2– скалярний добуток вектора G на самого себе
  • FхG'(0) = SUM F(i)*G(i) – скалярний добуток двох векторів F і G
Так само очевидно, що FхG'(t) дорівнює кореляції отримується в результаті зсуву одного вектора, щодо іншого на крок t (це легко перевірити явною підстановкою значень у формулу кореляції).

Перетворення Фур'є
Перетворення Фур'є (ℱ) — операція, що зіставляє однієї функції дійсної змінної іншу функцію, також речової змінної. Ця нова функція описує коефіцієнти («амплітуди») при розкладанні вихідної функції на елементарні складові — гармонійні коливання з різними частотами.
Перетворення Фур'є функції f речової змінної є інтегральним і визначається наступною формулою:
Перетворення Фур&#39;є
Різні джерела можуть давати визначення, що відрізняються від наведеного вище вибором коефіцієнта перед інтегралом, а також знака «−» в показнику експоненти. Але всі властивості будуть ті ж, хоча вигляд деяких формул може змінитися.
Крім того, існують різноманітні узагальнення даного поняття.

Багатовимірне перетворення Фур'є
Перетворення Фур'є функцій, заданих на просторі ℝ^n, визначається формулою
image
Зворотне перетворення в цьому випадку задається формулою
image
Як і раніше, в різних джерелах визначення багатовимірного перетворення Фур'є можуть відрізнятися вибором константи перед інтегралом.

Дискретне перетворення Фур'є
Дискретне перетворення Фур'є (в англомовній літературі DFT, Discrete Fourier Transform) — це одне з перетворень Фур'є, широко застосовуються в алгоритмах цифрової обробки сигналів (його модифікації застосовуються в стиску звуку в MP3, стиску зображень в JPEG та ін), а також в інших областях, пов'язаних з аналізом частот в дискретному (наприклад, оцифроване аналоговому) сигналі. Дискретне перетворення Фур'є вимагає в якості входу дискретну функцію. Такі функції часто створюються шляхом дискретизації (вибірки значень неперервних функцій). Дискретні перетворення Фур'є допомагають вирішувати диференціальні рівняння в приватних похідних і виконувати такі операції, як згортки. Дискретні перетворення Фур'є також активно використовуються в статистиці, при аналізі тимчасових рядів. Існують багатовимірні дискретні перетворення Фур'є.

Формули дискретних перетворень
Пряме перетворення:
image
Зворотне перетворення:
image
Дискретне перетворення Фур'є є лінійним перетворенням, яке переводить вектор часових відліків вектор спектральних відліків тієї ж довжини. Таким чином перетворення може бути реалізовано як множення симетричної квадратної матриці на вектор:
image

Фур'є-перетворення для обчислення згортки
Одним з чудових властивостей перетворень Фур'є є можливість швидкого обчислення кореляції двох функцій певних, або на дійсному аргументі (при використанні класичної формули), або на кінцевому кільці (при використанні дискретних перетворень).
І хоча подібні властивості притаманні багатьом лінійним перетворенням, для практичного застосування, для обчислення операції згортки, згідно з даного нами визначення, використовується формула
FхG = BFT ( FFT(F)*FFT(G) )
Де
  • FFT – операція прямого перетворення Фур'є
  • BFT – операція зворотного перетворення Фур'є
Перевірити правильність рівності досить легко – явно підставивши у формули Фур'є-перетворень і скоротивши отримані формули

Фур'є-перетворення для рішення задачі
При використанні формули для оцінки відстані між зображеннями при зсуві (i,j) відносно один одного
m(X,Y) (i,j) = <X,Y>(i,j) / ( |X|(i,j) ) * Y(i,j) ),
отримуємо, що
  • <X,Y> = XxY' = BFT ( FFT(X)*FFT(Y') )
  • |X| = XxX'xE' = BFT ( FFT(X)*FFT(X')*FFT(E') ) = BFT ( FFT(X)*FFT(X)'*FFT(E') )
  • |Y| = YxY'xE' = BFT ( FFT(Y)*FFT(Y')*FFT(E') ) = BFT ( FFT(Y)*FFT(Y)'*FFT(E') )
Де
  • <X,Y>(i,j) – скалярний добуток двох зображень, одержуваних при зсуві (i,j) відносно один одного зображення X і Y
  • E – зображення розміру рівному мінімальних розмірів X і Y, і заповнене одиничними значеннями (тобто «кадр» в якому порівнюються X та Y)
  • |X|(i,j) – норма загальної частини зображення при зсуві X (i,j)
  • |Y|(i,j) – норма загальної частини зображення Y при зсуві (i,j)
  • FFT – операція прямого двовимірного дискретного перетворення Фур'є
  • BFT – операція зворотного двовимірного дискретного перетворення Фур'є


Спрощення формул для вирішення поставленого завдання
При рішенні задачі для пошуку одного зразка, додаткове унормування зразка є зайвим, а також обчислення норми у загальній частині може бути замінено на суму яскравостей пікселів в цій загальній частині
При використанні формули для оцінки відстані між зображеннями при зсуві (i,j) відносно один одного
m(X,Y) (i,j) = <X,Y>(i,j) / |X|(i,j),
отримуємо, що
  • <X,Y> = XxY' = BFT ( FFT(X)*FFT(Y') )
  • |X| = XxE' = BFT ( FFT(X)*FFT(E') )
Де
  • <X,Y>(i,j) – скалярний добуток двох зображень, одержуваних при зсуві (i,j) відносно один одного зображення X і Y
  • E – зображення розміру рівному мінімальних розмірів X і Y, і заповнене одиничними значеннями (тобто «кадр» в якому порівнюються X і Y)
  • |X|(i,j) – норма (сума яскравостей пікселів) загальної частини зображення при зсуві X (i,j)
  • FFT – операція прямого двовимірного дискретного перетворення Фур'є
  • BFT – операція зворотного двовимірного дискретного перетворення Фур'є


Алгоритм пошуку фрагмента в повному зображенні
  • Нехай дано два зображення X і Y – зображення і зразок, розмірів (N1,N2) і (M1,M2) відповідно і Ni > Mi
  • Потрібно знайти координати зразка Y в повному зображенні X і обчислити оціночну величину — міру близькості.
  1. Розширити зображення Y до розміру (N1,N2), доповнивши його нулями
  2. Сформувати зображення E з одиниць розміру (M1,M2) і розширити до розміру (N1,N2), доповнивши його нулями
  3. Обчислити матрицю Y'[i,j[=Y[(N1-i) mod N1, (N2-i) mod N2]
  4. Обчислити матрицю E'[i,j]=E[(N1-i) mod N1, (N2-i) mod N2]
  5. Обчислити XxY' = BFT ( FFT(X)*FFT(Y') )
  6. Обчислити XxE' = BFT ( FFT(X)*FFT(E') )
  7. Обчислити M[i,j] = (f + XxY'[i,j])/(f + XxE'[i,j])
  8. В матриці M знайти елемент з максимальним значенням – координати елемента і є шуканої позицією зразка в повному зображенні, а значення дорівнює оцінці заходи порівняння.
Примітка
При використанні дискретного перетворення Фур'є, матриця M містить також елементи від циклічного зсуву зображень між собою. Тому, якщо не потрібно аналізувати циклічний зсув кадрів, пошук максимального елемента в матриці M потрібно обмежити областю (0,0)-(N1-M1, N2-M2)

Приклади реалізації
Реалізовані алгоритми є частиною бібліотеки з відкритим вихідним кодом FFTTools
Інтернет-адреса: github.com/dprotopopov/FFTTools
програмне забезпечення
  • Microsoft Visual Studio 2013 C# — середовище і мова програмування
  • EmguCV/OpenCV – C++ бібліотека структур і алгоритмів для обробки зображень
  • FFTWSharp/FFTW – C++ бібліотека реалізує алгоритми швидкого дискретного перетворення Фур'є


/ / / < 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(complex.Zip(output.GetData_Complex(), (x, y) => x*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*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[(n0 - i)%n0, (n1 - j)%n1] = 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[(n0 - i)%n0, (n1 - j)%n1] += 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 copied to replace 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[(n0 - i)%n0, (n1 - j)%n1] = 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, double value out)
{
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(Image<Bgr, Byte> bitmap)
{
using (Image<Gray, Byte> image = bitmap.Convert<Gray, Byte>())
return Catch(image);
}



Попався, який кусався
image

Література
  1. А. Л. Дмитрієв. Оптичні методи обробки інформації. Навчальний посібник. СПб. СПюГУИТМО 2005. 46 с.
  2. А. А. Акаєв, С. А. Майоров «Оптичні методи обробки інформації» М.:1988
  3. Дж.Гудмен «Введення в Фур'є-оптики» М.: Світ 1970

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.