August 31, 2015 at 11:05Fourier Digital Image Processing
Foreword
A digital photograph or other raster image is an array of numbers recorded by sensors of brightness levels in a two-dimensional plane. Knowing that, from a mathematical point of view, a thin lens performs the Fourier transform of images placed in focal planes, it is possible to create image processing algorithms that are analogues of image processing by the classical optical system.
The formula of such algorithms will look as follows:
- Z = FFT (X) - direct two-dimensional Fourier transform
- Z ′ = T (Z) - application of a function or transparency to the Fourier transform of the image
- Y = BFT (Z ′) is the inverse two-dimensional Fourier transform
To calculate the Fourier transforms, fast discrete Fourier transform algorithms are used. Although the optical system of lenses performs the Fourier transform on a continuous range of the argument for the continuous spectrum, but when switching to digital data processing, the Fourier transform formulas can be replaced by discrete Fourier transform formulas.
Implementation examples
- Image Blur Algorithm
- Image sharpening algorithm
- Image Scaling Algorithm
The implemented algorithms are part of the open source library FFTTools. Internet Address: github.com/dprotopopov/FFTTools
Image Blur Algorithm
In optical systems, the diaphragm located in the focal plane is a simple hole in the screen. As a result of the passage of light through the diaphragm, high-frequency waves (with shorter wavelengths) pass through the obstacle, and low-frequency waves (with longer wavelengths) are cut off by the screen. Thus, the sharpness of the resulting image is increased. If you replace the hole in the screen with an obstacle in the screen, the result will be a blurry image, since it will be formed from the frequencies of long wavelengths.
Algorithm:
- Let X (N1, N2) be an array of brightnesses of image pixels.
- Calculate Px = average (rms) brightness of pixels in array X
- Compute array Z = FT (X) - direct two-dimensional discrete Fourier transform
- Calculate the array Z ′ = T (Z), where T is the zeroing of rows and columns located in the specified internal regions of the argument matrix corresponding to high 5. frequencies (i.e., the zeroing of the Fourier expansion coefficients corresponding to high frequencies)
- Compute array Y = RFT (Z ′) - inverse two-dimensional discrete Fourier transform
- Calculate Py = average (rms) brightness of pixels in array Y
- Normalize array Y (N1, N2) according to average brightness level Px / Py
Image sharpening algorithm
In optical systems, the diaphragm located in the focal plane is a simple hole in the screen. As a result of the passage of light through the diaphragm, high-frequency waves (with shorter wavelengths) pass through the obstacle, and low-frequency waves (with longer wavelengths) are cut off by the screen. Thus, the sharpness of the resulting image is increased.
Algorithm:
- Let X (N1, N2) be an array of brightnesses of image pixels.
- Calculate Px = average (rms) brightness of pixels in array X
- Compute array Z = FT (X) - direct two-dimensional discrete Fourier transform
- Save value L = Z (0,0) - corresponding to the average brightness of the pixels of the original image
- Calculate the array Z ′ = T (Z), where T is the zeroing of rows and columns located in the specified external regions of the argument matrix corresponding to low 6. frequencies (i.e., the zeroing of the Fourier expansion coefficients corresponding to low frequencies)
- Restore the value Z '(0,0) = L - corresponding to the average brightness of the pixels of the original image
- Compute array Y = RFT (Z ′) - inverse two-dimensional discrete Fourier transform
- Calculate Py = average (rms) brightness of pixels in array Y
- Normalize array Y (N1, N2) according to average brightness level Px / Py
Image Scaling Algorithm
In optical systems, the luminous flux in the focal plane of the system is the Fourier transform of the original image. The size of the image obtained at the output of the optical system is determined by the ratio of the focal distances of the lens and eyepiece.
Algorithm:
- Let X (N1, N2) be an array of brightnesses of image pixels.
- Calculate Px = average (rms) brightness of pixels in array X
- Compute array Z = FT (X) - direct two-dimensional discrete Fourier transform
- Calculate the array Z ′ = T (Z), where T is either adding zero rows and columns of the matrix corresponding to high frequencies, or deleting rows and columns of the matrix corresponding to high frequencies to obtain the required size of the final image
- Compute array Y = RFT (Z ′) - inverse two-dimensional discrete Fourier transform
- Calculate Py = average (rms) brightness of pixels in array Y
- Normalize array Y (M1, M2) according to average brightness level Px / Py
Software used
- 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
Image Blur Algorithm
Algorithm code
///
/// Clear internal region of array
///
/// Array of values
/// Internal blind region size
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);
}
}
///
/// Blur bitmap with the Fastest Fourier Transform
///
/// Blured bitmap
public Bitmap Blur(Bitmap bitmap)
{
using (var image = new Image(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;
}
}
Image sharpening algorithm
Algorithm code
///
/// Clear external region of array
///
/// Array of values
/// External blind region size
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);
}
}
///
/// Sharp bitmap with the Fastest Fourier Transform
///
/// Sharped bitmap
public Bitmap Sharp(Bitmap bitmap)
{
using (var image = new Image(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;
}
}
Image Scaling Algorithm
Algorithm code
///
/// Copy arrays
///
/// Input array
/// Output array
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];
}
}
}
}
///
/// Resize bitmap with the Fastest Fourier Transform
///
/// Resized bitmap
public Bitmap Stretch(Bitmap bitmap)
{
using (var image = new Image(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(_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;
}
}
}
Program screenshots
Image blur
Image scaling
Literature
- A.L. Dmitriev. Optical methods of information processing. Tutorial. SPb. SPYUGUITMO 2005.46 p.
- A.A. Akayev, S.A. Mayorov “Optical methods of information processing” M.: 1988
- J. Goodman "Introduction to Fourier Optics" M .: World 1970