Предисловие
Цифровая фотография или иное растровое изображение представляет собой массив чисел, зафиксированных сенсорами уровней яркости, в двумерной плоскости. Зная что с математической точки зрения тонкая линза выполняет преобразование Фурье изображений, размещённых в фокальных плоскостях, можно создать алгоритмы обработки изображений, являющихся аналогами обработки изображений классической оптической системой.
Формула таких алгоритмов будет выглядеть следующим образом:
- Z=FFT(X) – прямое двухмерное преобразование Фурье
- Z′=T(Z) – применение функции или транспаранта к Фурье-образу изображения
- Y=BFT(Z′) – обратное двухмерное преобразование Фурье
Для вычисления преобразований Фурье используются алгоритмы быстрого дискретного преобразования Фурье. Хотя оптическая система линз осуществляет преобразование Фурье на непрерывном диапазоне аргумента и для непрерывного спектра, но при переходе к цифровой обработке данных формулы преобразования Фурье могут быть заменены на формулы дискретного преобразования Фурье.
Примеры реализации
- Алгоритм размытия изображения
- Алгоритм повышения резкости изображения
- Алгоритм масштабирования изображения
Реализованные алгоритмы являются частью библиотеки с открытым исходным кодом FFTTools. Интернет-адрес: github.com/dprotopopov/FFTTools
Алгоритм размытия изображения
В оптических системах диафрагма, размещённая в фокальной плоскости, представляет собой простое отверстие в экране. В результате прохождения светового потока через диафрагму, волны высоких частот (с более короткими длинами волн) проходят через препятствие, а волны низких частот (с более длинными длинами волн) отсекаются экраном. Таким образом повышается резкость получаемого изображения. Если заменить отверстие в экране на препятствие в экране, то в результате будет получено размытое изображение, поскольку оно будет сформировано из частот волн больших длин.
Алгоритм:
- Пусть X(N1,N2) – массив яркостей пикселей изображения.
- Вычислить Px = средняя (среднеквадратичная) яркость пикселей в массиве X
- Вычислить массив Z=FT(X) – прямое двухмерное дискретное преобразование Фурье
- Вычислить массив Z′=T(Z), где T – обнуление строк и столбцов, находящихся в заданных внутренних областях матрицы-аргумента соответствующих высоким 5. частотам (то есть обнуление коэффициентов Фурье-разложения, соответствующих высоким частотам)
- Вычислить массив Y=RFT(Z′) – обратное двухмерное дискретное преобразование Фурье
- Вычислить Py = средняя (среднеквадратичная) яркость пикселей в массиве Y
- Нормировать массив Y(N1,N2) по среднему уровню яркости Px/Py
Алгоритм повышения резкости изображения
В оптических системах диафрагма, размещённая в фокальной плоскости, представляет собой простое отверстие в экране. В результате прохождения светового потока через диафрагму, волны высоких частот (с более короткими длинами волн) проходят через препятствие, а волны низких частот (с более длинными длинами волн) отсекаются экраном. Таким образом повышается резкость получаемого изображения.
Алгоритм:
- Пусть X(N1,N2) – массив яркостей пикселей изображения.
- Вычислить Px = средняя (среднеквадратичная) яркость пикселей в массиве X
- Вычислить массив Z=FT(X) – прямое двухмерное дискретное преобразование Фурье
- Сохранить значение L=Z(0,0) – соответствующее средней яркости пикселей исходного изображения
- Вычислить массив Z′=T(Z), где T – обнуление строк и столбцов, находящихся в заданных внешних областях матрицы-аргумента, соответствующих низким 6. частотам (то есть обнуление коэффициентов Фурье-разложения, соответствующих низким частотам)
- Восстановить значение Z’(0,0)=L – соответствующее средней яркости пикселей исходного изображения
- Вычислить массив Y=RFT(Z′) – обратное двухмерное дискретное преобразование Фурье
- Вычислить Py = средняя (среднеквадратичная) яркость пикселей в массиве Y
- Нормировать массив Y(N1,N2) по среднему уровню яркости Px/Py
Алгоритм масштабирования изображения
В оптических системах световой поток в фокальной плоскости системы представляет собой Фурье-преобразование исходного изображения. Размер получаемого на выходе оптической системы изображения определяется соотношением фокальных расстояний объектива и окуляра.
Алгоритм:
- Пусть X(N1,N2) – массив яркостей пикселей изображения.
- Вычислить Px = средняя (среднеквадратичная) яркость пикселей в массиве X
- Вычислить массив Z=FT(X) – прямое двухмерное дискретное преобразование Фурье
- Вычислить массив Z′=T(Z), где T – либо добавление нулевых строк и столбцов матрицы соответствующих высоким частотам, либо удаление строк и столбцов матрицы соответствующих высоким частотам для получения требуемого размера итогового изображения
- Вычислить массив Y=RFT(Z′) – обратное двухмерное дискретное преобразование Фурье
- Вычислить Py = средняя (среднеквадратичная) яркость пикселей в массиве Y
- Нормировать массив 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;
}
}
}
Скриншоты программ
Размытие изображения
Масштабирование изображения
Литература
- А.Л. Дмитриев. Оптические методы обработки информации. Учебное пособие. СПб. СПюГУИТМО 2005. 46 с.
- А.А.Акаев, С.А.Майоров «Оптические методы обработки информации» М.:1988
- Дж.Гудмен «Введение в Фурье-оптику» М.: Мир 1970