Pull to refresh

Фурье-обработка цифровых изображений

ProgrammingAlgorithmsImage processingMathematics
Sandbox

Предисловие


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

Формула таких алгоритмов будет выглядеть следующим образом:
  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;
                }
            }
        }



Скриншоты программ


Размытие изображения
image


Масштабирование изображения
image


Литература


  1. А.Л. Дмитриев. Оптические методы обработки информации. Учебное пособие. СПб. СПюГУИТМО 2005. 46 с.
  2. А.А.Акаев, С.А.Майоров «Оптические методы обработки информации» М.:1988
  3. Дж.Гудмен «Введение в Фурье-оптику» М.: Мир 1970
Tags:Фурье-оптика
Hubs: Programming Algorithms Image processing Mathematics
Total votes 28: ↑22 and ↓6+16
Views33K

Top of the last 24 hours