[Из песочницы] Фурье-обработка цифровых изображений

         /// 
        ///     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;
                }
            }
        }


© Habrahabr.ru