itk下FFTW的FFT和IFFT

原文：http://blog.csdn.net/ljp1919/article/details/41890425

1：先进行FFT变换，再将图像进行逆变换，并对比两者的结果。

[cpp] view plain copy

1. #include <tiffio.h>  
2. #include "fftw3.h"  
3. #include <stdio.h>  
4. #include "itkImageFileReader.h"  
5. #include "itkImageFileWriter.h"  
6. #include "itkWrapPadImageFilter.h"  
7. #include "itkForwardFFTImageFilter.h"  
8. #include "itkCastImageFilter.h"  
9. #include "itkRescaleIntensityImageFilter.h"  
10. #include "itkComplexToModulusImageFilter.h"  
11. #include "itkIntensityWindowingImageFilter.h"  
12. #include "itkFFTShiftImageFilter.h"  
13.   
14. #pragma comment(lib, "libfftw3-3.lib")  
15. #pragma comment(lib, "libfftw3f-3.lib")  
16. #pragma comment(lib, "libfftw3l-3.lib")  
17. #define N 8  
18. typedef itk:: Image<float , 2> FloatImageType;  
19. typedef itk ::Image< unsigned char , 2> UnsignedCharImageType;  
20. typedef float MatrixType;  
21.   
22. void CopyMatrix(FloatImageType:: Pointer image , MatrixType **mask);  
23. void CopyImage(FloatImageType:: Pointer image , fftw_complex *mask);  
24. void CopyRealImage (FloatImageType:: Pointer image , float ** mask);//仅仅拷贝实部  
25. void MyCastFloatToUnsichar(FloatImageType::Pointer image, UnsignedCharImageType::Pointer imageNew );//将float转为unsigned char型  
26. int main()  
27. {  
28.     //读取图像  
29.      FloatImageType::Pointer image;  
30.       
31.      const char * inputFileName = "F:\\DIPcode\\SIMC\\imagesrc\\2.jpg" ; //argv[1];  
32.      const char * outputFileName = "F:\\DIPcode\\SIMC\\2_new.tif" ; //argv[2];  
33.   
34.      typedef itk ::ImageFileReader< FloatImageType> ReaderType ;  
35.      ReaderType::Pointer reader = ReaderType::New ();  
36.      reader->SetFileName (inputFileName);  
37.      reader->Update ();  
38.      image = reader ->GetOutput();  
39.   
40.     //遍历图像像素，赋值到一个矩阵  
41.      MatrixType **pixelColors = new MatrixType*[512]();  
42.      for (int i=0;i<512;i++)  
43.      {  
44.          pixelColors[i] = new MatrixType[512]; //为每个指针都分配8个char元素空间。  
45.      }  
46.      CopyMatrix(image ,  pixelColors);  
47.   
48.     int width =512;  
49.     int height=512;  
50.   
51.     fftw_plan planR;  
52.     fftw_complex *inR,  *outR, *resultR;  
53.   
54.     //Allocate arrays.  
55.     inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);  
56.     outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);  
57.     resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);  
58.   
59.     //Fill in arrays with the pixelcolors.  
60.     for (int y = 0; y < height; y++) {  
61.         for (int x = 0; x < width; x++) {  
62.             //int currentIndex = ((y * width) + (x)) * 3;  
63.             int currentIndex = (y * width) + x;  
64.             inR[y * width + x][0] = (double)pixelColors[y][x];  
65.             inR[y * width + x][1] = 0.0;      
66.         }  
67.     }  
68.   
69.     //Forward plans.  
70.     planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE);  
71.       
72.     //Forward FFT.  
73.     fftw_execute(planR);  
74.       
75.     //Backward plans.  
76.     fftw_plan planRR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE);  
77.       
78.     //Backward fft  
79.     fftw_execute(planRR);  
80.   
81.      Have to scale the output values to get back to the original.  
82.     int size=height*width;  
83.     for(int i = 0; i < size; i++) {  
84.         resultR[i][0] = resultR[i][0] / size;  
85.         resultR[i][1] = resultR[i][1] / size;  
86.     }  
87.   
88.     //Overwrite the pixelcolors with the result.  
89.     for (int y = 0; y < height; y++)   
90.     {  
91.         for (int x = 0; x < width; x++)  
92.         {  
93.             //int currentIndex = ((y * width) + (x)) * 3;  
94.             pixelColors[y][x] =resultR[y * width + x][0];//  
95.             //std::cout<<pixelColors[y][x]<<'\t';  
96.         }  
97.     }  
98.       
99.     //将图像写出来  
100.     CopyRealImage(image , pixelColors);//将pixelColors中的像素值复制到image中  
101.     //调整亮度  
102.     typedef itk ::RescaleIntensityImageFilter< FloatImageType , FloatImageType > RescaleFilterType;  
103.     RescaleFilterType::Pointer imaginaryRescaleFilter = RescaleFilterType ::New();  
104.     imaginaryRescaleFilter->SetInput (image);  
105.     imaginaryRescaleFilter->SetOutputMinimum (0);  
106.     imaginaryRescaleFilter->SetOutputMaximum (255);  
107.     imaginaryRescaleFilter->Update ();  
108.     //转为8bit  
109.     //类型转换，float转为8bit的。  
110.     typedef itk ::CastImageFilter< FloatImageType, UnsignedCharImageType > CastFilterType;  
111.     CastFilterType::Pointer castFilter = CastFilterType::New ();  
112.     castFilter->SetInput (imaginaryRescaleFilter->GetOutput());  
113.     castFilter->Update ();  
114.     //将图像写出来  
115.     typedef itk ::ImageFileWriter< UnsignedCharImageType > WriterType ;  
116.     WriterType::Pointer writer = WriterType::New ();  
117.     writer->SetFileName ( outputFileName );  
118.     writer->SetInput (castFilter->GetOutput());  
119.     try  
120.     {  
121.         writer->Update ();  
122.     }  
123.     catch( itk ::ExceptionObject & error )  
124.     {  
125.         std::cerr << "Error: " << error << std ::endl;  
126.         return EXIT_FAILURE ;  
127.     }  
128.     delete []pixelColors;  
129.     fftw_free(inR);  
130.     fftw_free(outR);  
131.     fftw_free(resultR);  
132.   
133.     return 0;  
134. }  
135. void CopyMatrix(FloatImageType:: Pointer image , MatrixType **mask)  
136. {  
137.     FloatImageType::RegionType region = image->GetLargestPossibleRegion ();  
138.     FloatImageType::SizeType regionSize = region.GetSize ();  
139.   
140.     for (int x=0; x<regionSize [0]; x++)  
141.         for(int y=0; y<regionSize [1]; y++)  
142.         {   //3、定义像素索引  
143.             FloatImageType::IndexType index;  
144.             index[0] = x ;  
145.             index[1] = y ;  
146.             mask[x][y] =(MatrixType) image->GetPixel(index);//获取该位置的像素值  
147.         }  
148.   
149. }  
150.   
151. void CopyImage(FloatImageType:: Pointer image , fftw_complex *mask)  
152. {  
153.     FloatImageType::RegionType region = image->GetLargestPossibleRegion ();  
154.     FloatImageType::SizeType regionSize = region.GetSize ();  
155.   
156.     for (int x=0; x<regionSize [0]; x++)  
157.         for(int y=0; y<regionSize [1]; y++)  
158.         {  
159.             double re = mask[x*regionSize[1]+y][0];  
160.             double im = mask[x*regionSize[1]+y][1];  
161.             double mag = sqrt (re*re + im*im);  
162.             //3、定义像素索引  
163.             FloatImageType::IndexType index;  
164.             index[0] = x ;  
165.             index[1] = y ;  
166.             image->SetPixel(index, mag);//获取该位置的像素值  
167.         }  
168.   
169. }  
170. //仅仅将mask中的值，赋给image中！  
171. void CopyRealImage (FloatImageType:: Pointer image , MatrixType ** mask)  
172. {  
173.     FloatImageType::RegionType region = image->GetLargestPossibleRegion ();  
174.     FloatImageType::SizeType regionSize = region.GetSize ();  
175.     int Imgsize = regionSize[0]*regionSize[1];  
176.   
177.     for (int x=0; x<regionSize [0]; x++)  
178.         for(int y=0; y<regionSize [1]; y++)  
179.         {   //3、定义像素索引  
180.             FloatImageType::IndexType index;  
181.             index[0] = x ;  
182.             index[1] = y ;  
183.               
184.             image->SetPixel (index, mask[x ][y]);   
185.             //获取该位置的像素值  
186.             //std::cout<<image->GetPixel(index);  
187.         }  
188.   
189. }  
190. void MyCastFloatToUnsichar(FloatImageType::Pointer image, UnsignedCharImageType::Pointer imageNew )//将float转为unsigned char型  
191. {  
192.     FloatImageType::RegionType region = image->GetLargestPossibleRegion ();  
193.     FloatImageType::SizeType regionSize = region.GetSize ();  
194.       
195.     for (int x=0; x<regionSize [0]; x++)  
196.         for(int y=0; y<regionSize [1]; y++)  
197.         {   //3、定义像素索引  
198.             FloatImageType::IndexType index;  
199.             index[0] = x ;  
200.             index[1] = y ;  
201.             //std::cout<<image->GetPixel(index)<<'\t';  
202.             float temp = static_cast<unsigned char>(image->GetPixel(index));  
203.             //std::cout<<temp<<'\t';  
204.             imageNew->SetPixel (index, temp);   
205.             int tt=imageNew->GetPixel(index);  
206.             获取该位置的像素值  
207.             //std::cout<<tt<<'\t';  
208.         }  
209.   
210. }  

起初在不进行cast进行图像类型转换的时候，因为是float所以结果是32位的，但是直接对image进行cast，映射为unsigned char （8bit）之后，结果一片黑色，说明在进行映射过程中发生了数据截断。后来经过调整，先进行rescale亮度。再重新进行映射，则可以重建出其逆变换的结果。但是对比两者结果图还是有些不一样的。