Фур'є-обробка цифрових зображень

Передмова
Цифрова фотографія або інше растрове зображення являє собою масив чисел, зафіксованих сенсорами рівнів яскравості, у двомірній площині. Знаючи що з математичної точки зору тонка лінза виконує перетворення Фур'є зображень, розміщених у фокальних площинах, можна створити алгоритми обробки зображень, що є аналогами обробки зображень класичної оптичною системою.

Формула таких алгоритмів буде виглядати наступним чином:
  1. Z=FFT(X) – пряме двовимірне перетворення Фур'є
  2. Z'=T(Z) – застосування функції або заголовку до Фур'є-образу зображення
  3. Y=BFT(Z') – зворотне двовимірне перетворення Фур'є
Для обчислення перетворень Фур'є використовуються алгоритми швидкого дискретного перетворення Фур'є. Хоча оптична система лінз здійснює перетворення Фур'є на безперервному діапазоні аргументу і для безперервного спектру, але при переході до цифрової обробки даних формули перетворення Фур'є можуть бути замінені на формули дискретного перетворення Фур'є.

Приклади реалізації
  • Алгоритм розмиття зображення
  • Алгоритм підвищення різкості зображення
  • Алгоритм масштабування зображення
Реалізовані алгоритми є частиною бібліотеки з відкритим вихідним кодом FFTTools. Інтернет-адреса: github.com/dprotopopov/FFTTools

Алгоритм розмиття зображення
В оптичних системах діафрагма, розміщена в фокальній площині, являє собою просте отвір в екрані. В результаті проходження світлового потоку через діафрагму, хвилі високих частот (з більш короткими довжинами хвиль) проходять через перешкоду, а хвилі низьких частот (з більш довгими довжинами хвиль) відсікаються екраном. Таким чином підвищується різкість одержуваного зображення. Якщо замінити отвір в екрані на перешкоду в екрані, то в результаті буде отримано розмите зображення, оскільки воно буде сформовано з частот великих довжин хвиль.

Алгоритм:
  1. Нехай X(N1,N2) – масив яскравостей пікселів зображення.
  2. Обчислити Px = середня (середньоквадратична) яскравість пікселів в масиві X
  3. Обчислити масив Z=FT(X) – пряме двомірне дискретне перетворення Фур'є
  4. Обчислити масив Z'=T(Z), де T – обнулення рядків і стовпців, які перебувають у заданих внутрішніх областях матриці-аргументу відповідають високим 5. частотам (тобто обнулення коефіцієнтів Фур'є-розкладу, що відповідають високим частотам)
  5. Обчислити масив Y=RFT(Z') – зворотне двомірне дискретне перетворення Фур'є
  6. Обчислити Py = середня (середньоквадратична) яскравість пікселів в масиві Y
  7. Нормувати масив Y(N1,N2) за середнім рівнем яскравості Px/Py
Алгоритм підвищення різкості зображення
В оптичних системах діафрагма, розміщена в фокальній площині, являє собою просте отвір в екрані. В результаті проходження світлового потоку через діафрагму, хвилі високих частот (з більш короткими довжинами хвиль) проходять через перешкоду, а хвилі низьких частот (з більш довгими довжинами хвиль) відсікаються екраном. Таким чином підвищується різкість одержуваного зображення.

Алгоритм:
  1. Нехай X(N1,N2) – масив яскравостей пікселів зображення.
  2. Обчислити Px = середня (середньоквадратична) яскравість пікселів в масиві X
  3. Обчислити масив Z=FT(X) – пряме двомірне дискретне перетворення Фур'є
  4. Зберегти значення L=Z(0,0) – відповідне середньої яскравості пікселів вихідного зображення
  5. Обчислити масив Z'=T(Z), де T – обнулення рядків і стовпців, які перебувають у заданих зовнішніх областях матриці-аргумент, відповідних низьким 6. частотам (тобто обнулення коефіцієнтів Фур'є-розкладу, відповідних низьким частотам)
  6. Відновити значення Z'(0,0)=L – відповідне середньої яскравості пікселів вихідного зображення
  7. Обчислити масив Y=RFT(Z') – зворотне двомірне дискретне перетворення Фур'є
  8. Обчислити Py = середня (середньоквадратична) яскравість пікселів в масиві Y
  9. Нормувати масив Y(N1,N2) за середнім рівнем яскравості Px/Py
Алгоритм масштабування зображення
В оптичних системах світловий потік в фокальній площині системи являє собою Фур'є-перетворення вихідного зображення. Розмір одержуваного на виході оптичної системи зображення визначається співвідношенням фокальних відстаней об'єктиву і окуляра.

Алгоритм:
  1. Нехай X(N1,N2) – масив яскравостей пікселів зображення.
  2. Обчислити Px = середня (середньоквадратична) яскравість пікселів в масиві X
  3. Обчислити масив Z=FT(X) – пряме двомірне дискретне перетворення Фур'є
  4. Обчислити масив Z'=T(Z), де T – небудь додавання нульових рядків і стовпців матриці відповідають високим частотам, або видалення рядків і стовпців матриці відповідають високим частотам для отримання необхідного розміру вихідного зображення
  5. Обчислити масив Y=RFT(Z') – зворотне двомірне дискретне перетворення Фур'є
  6. Обчислити Py = середня (середньоквадратична) яскравість пікселів в масиві Y
  7. Нормувати масив Y(M1,M2) за середнім рівнем яскравості Px/Py
програмне забезпечення
  • Microsoft Visual Studio 2013 C# — середовище і мова програмування
  • EmguCV/OpenCV – C++ бібліотека структур і алгоритмів для обробки зображень
  • FFTWSharp/FFTW – C++ бібліотека реалізує алгоритми швидкого дискретного перетворення Фур'є
Алгоритм розмиття зображення
Код алгоритму
/ / / < summary>
/// Clear internal region of array
/ / / < /summary>
/ / / < param name="data">Array of values</param>
/ / / < param name="size">Internal blind region size</param>
private static void Blind(Complex[,,] data, Size size)
{
int n0 = data.GetLength(0);
int n1 = data.GetLength(1);
int n2 = data.GetLength(2);
int s0 = Math.Max(0, (n0 - size.Height)/2);
int s1 = Math.Max(0, (n1 - size.Width)/2);
int e0 = Math.Min((n0 + size.Height)/2, n0);
int e1 = Math.Min((n1 + size.Width)/2, n1);
for (int i = s0; i < e0; i++)
{
Array.Clear(data, i*n1*n2, n1*n2);
}
for (int i = 0; i < s0; i++)
{
Array.Clear(data, i*n1*n2 + s1*n2, (e1 - s1)*n2);
}
for (int i = e0; i < n0; i++)
{
Array.Clear(data, i*n1*n2 + s1*n2, (e1 - s1)*n2);
}
}
/// <summary>
/// Blur bitmap with the Fastest Fourier Transform
/ / / < /summary>
/ / / < returns>Blured bitmap</returns>
public Bitmap Blur(Bitmap bitmap)
{
using (var image = new Image<Bgr, double>(bitmap))
{
int length = image.Data.Length;
int n0 = image.Data.GetLength(0);
int n1 = image.Data.GetLength(1);
int n2 = image.Data.GetLength(2);

var doubles = new double[length];
Buffer.BlockCopy(image.Data, 0, doubles, 0, length*sizeof (double));
double power = Math.Sqrt(doubles.Average(x => x*x));

var input = new fftw_complexarray(doubles.Select(x => new Complex(x, 0)).ToArray());
var output = new fftw_complexarray(length);
fftw_plan.dft_3d(n0, n1, n2, input, output,
fftw_direction.Forward,
fftw_flags.Estimate).Execute();
Complex[] complex = output.GetData_Complex();

var data = new Complex[n0, n1, n2];
var buffer = new double[length*2];

GCHandle complexHandle = GCHandle.Alloc(complex, GCHandleType.Pinned);
GCHandle dataHandle = GCHandle.Alloc(data, GCHandleType.Pinned);
IntPtr complexPtr = complexHandle.AddrOfPinnedObject();
IntPtr dataPtr = dataHandle.AddrOfPinnedObject();

Marshal.Copy(complexPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, dataPtr, buffer.Length);
Blind(data, _blinderSize);
Marshal.Copy(dataPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, complexPtr, buffer.Length);

complexHandle.Free();
dataHandle.Free();

input.SetData(complex);

fftw_plan.dft_3d(n0, n1, n2, input, output,
fftw_direction.Backward,
fftw_flags.Estimate).Execute();
double[] array2 = output.GetData_Complex().Select(x => x.Magnitude).ToArray();
double power2 = Math.Sqrt(array2.Average(x => x*x));
doubles = array2.Select(x => x*power/power2).ToArray();
Buffer.BlockCopy(doubles, 0, image.Data, 0, length*sizeof (double));
return image.Bitmap;
}
}


Алгоритм підвищення різкості зображення
Код алгоритму
/ / / < summary>
/// Clear external region of array
/ / / < /summary>
/ / / < param name="data">Array of values</param>
/ / / < param name="size">External blind region size</param>
private static void Blind(Complex[,,] data, Size size)
{
int n0 = data.GetLength(0);
int n1 = data.GetLength(1);
int n2 = data.GetLength(2);
int s0 = Math.Max(0, (n0 - size.Height)/2);
int s1 = Math.Max(0, (n1 - size.Width)/2);
int e0 = Math.Min((n0 + size.Height)/2, n0);
int e1 = Math.Min((n1 + size.Width)/2, n1);
for (int i = 0; i < s0; i++)
{
Array.Clear(data, i*n1*n2, s1*n2);
Array.Clear(data, i*n1*n2 + e1*n2, (n1 - e1)*n2);
}
for (int i = e0; i < n0; i++)
{
Array.Clear(data, i*n1*n2, s1*n2);
Array.Clear(data, i*n1*n2 + e1*n2, (n1 - e1)*n2);
}
}
/// <summary>
/// Sharp bitmap with the Fastest Fourier Transform
/ / / < /summary>
/ / / < returns>Sharped bitmap</returns>
public Bitmap Sharp(Bitmap bitmap)
{
using (var image = new Image<Bgr, double>(bitmap))
{
int length = image.Data.Length;
int n0 = image.Data.GetLength(0);
int n1 = image.Data.GetLength(1);
int n2 = image.Data.GetLength(2);

var doubles = new double[length];
Buffer.BlockCopy(image.Data, 0, doubles, 0, length*sizeof (double));
double power = Math.Sqrt(doubles.Average(x => x*x));

var input = new fftw_complexarray(doubles.Select(x => new Complex(x, 0)).ToArray());
var output = new fftw_complexarray(length);
fftw_plan.dft_3d(n0, n1, n2, input, output,
fftw_direction.Forward,
fftw_flags.Estimate).Execute();
Complex[] complex = output.GetData_Complex();

Complex level = complex[0];

var data = new Complex[n0, n1, n2];
var buffer = new double[length*2];

GCHandle complexHandle = GCHandle.Alloc(complex, GCHandleType.Pinned);
GCHandle dataHandle = GCHandle.Alloc(data, GCHandleType.Pinned);
IntPtr complexPtr = complexHandle.AddrOfPinnedObject();
IntPtr dataPtr = dataHandle.AddrOfPinnedObject();

Marshal.Copy(complexPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, dataPtr, buffer.Length);
Blind(data, _blinderSize);
Marshal.Copy(dataPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, complexPtr, buffer.Length);

complexHandle.Free();
dataHandle.Free();

complex[0] = level;

input.SetData(complex);

fftw_plan.dft_3d(n0, n1, n2, input, output,
fftw_direction.Backward,
fftw_flags.Estimate).Execute();
double[] array2 = output.GetData_Complex().Select(x => x.Magnitude).ToArray();
double power2 = Math.Sqrt(array2.Average(x => x*x));
doubles = array2.Select(x => x*power/power2).ToArray();
Buffer.BlockCopy(doubles, 0, image.Data, 0, length*sizeof (double));
return image.Bitmap;
}
}


Алгоритм масштабування зображення
Код алгоритму
/ / / < summary>
/// Copy arrays
/ / / < /summary>
/ / / < param name="input">Input array</param>
/ / / < param name="output">Output array</param>
private static void Copy(Complex[,,] input, Complex[,,] output)
{
int n0 = input.GetLength(0);
int n1 = input.GetLength(1);
int n2 = input.GetLength(2);
int m0 = output.GetLength(0);
int m1 = output.GetLength(1);
int m2 = output.GetLength(2);
int ex0 = Math.Min(n0, m0)/2;
int ex1 = Math.Min(n1, m1)/2;
int ex2 = Math.Min(n2, m2);
Debug.Assert(n2 == m2);
for (int k = 0; k < ex2; k++)
{
for (int i = 0; i <= ex0; i++)
{
for (int j = 0; j <= ex1; j++)
{
int ni = n0 - i - 1;
int nj = n1 - j - 1;
int mi = m0 - i - 1;
int mj = m1 - j - 1;
output[i, j, k] = input[i, j, k];
output[mi, j, k] = input[ni, j, k];
output[i, mj, k] = input[i, nj, k];
output[mi, mj, k] = input[ni, nj, k];
}
}
}
}
/// <summary>
/// Resize bitmap with the Fastest Fourier Transform
/ / / < /summary>
/ / / < returns>Resized bitmap</returns>
public Bitmap Stretch(Bitmap bitmap)
{
using (var image = new Image<Bgr, double>(bitmap))
{
int length = image.Data.Length;
int n0 = image.Data.GetLength(0);
int n1 = image.Data.GetLength(1);
int n2 = image.Data.GetLength(2);
var doubles = new double[length];
Buffer.BlockCopy(image.Data, 0, doubles, 0, length*sizeof (double));
double power = Math.Sqrt(doubles.Average(x => x*x));

var input = new fftw_complexarray(doubles.Select(x => new Complex(x, 0)).ToArray());
var output = new fftw_complexarray(length);
fftw_plan.dft_3d(n0, n1, n2, input, output,
fftw_direction.Forward,
fftw_flags.Estimate).Execute();
Complex[] complex = output.GetData_Complex();

using (var image2 = new Image<Bgr, double>(_newSize))
{
int length2 = image2.Data.Length;
int m0 = image2.Data.GetLength(0);
int m1 = image2.Data.GetLength(1);
int m2 = image2.Data.GetLength(2);
var complex2 = new Complex[length2];

var data = new Complex[n0, n1, n2];
var data2 = new Complex[m0, m1, m2];

var buffer = new double[length*2];
GCHandle complexHandle = GCHandle.Alloc(complex, GCHandleType.Pinned);
GCHandle dataHandle = GCHandle.Alloc(data, GCHandleType.Pinned);
IntPtr complexPtr = complexHandle.AddrOfPinnedObject();
IntPtr dataPtr = dataHandle.AddrOfPinnedObject();

Marshal.Copy(complexPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, dataPtr, buffer.Length);

complexHandle.Free();
dataHandle.Free();

Copy(data, data2);

buffer = new double[length2*2];
complexHandle = GCHandle.Alloc(complex2, GCHandleType.Pinned);
dataHandle = GCHandle.Alloc(data2, GCHandleType.Pinned);
complexPtr = complexHandle.AddrOfPinnedObject();
dataPtr = dataHandle.AddrOfPinnedObject();

Marshal.Copy(dataPtr, buffer, 0, buffer.Length);
Marshal.Copy(buffer, 0, complexPtr, buffer.Length);

complexHandle.Free();
dataHandle.Free();

var input2 = new fftw_complexarray(complex2);
var output2 = new fftw_complexarray(length2);
fftw_plan.dft_3d(m0, m1, m2, input2, output2,
fftw_direction.Backward,
fftw_flags.Estimate).Execute();
double[] array2 = output2.GetData_Complex().Select(x => x.Magnitude).ToArray();
double power2 = Math.Sqrt(array2.Average(x => x*x));
double[] doubles2 = array2.Select(x => x*power/power2).ToArray();
Buffer.BlockCopy(doubles2, 0, image2.Data, 0, length2*sizeof (double));
return image2.Bitmap;
}
}
}


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

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

0 коментарів

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