itk下FFTW的FFT和IFFT

来源：互联网发布：dns域名服务器简答题编辑：程序博客网时间：2024/05/16 07:19

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

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

[cpp] view plain copy

#include <tiffio.h>
#include "fftw3.h"
#include <stdio.h>
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkWrapPadImageFilter.h"
#include "itkForwardFFTImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkComplexToModulusImageFilter.h"
#include "itkIntensityWindowingImageFilter.h"
#include "itkFFTShiftImageFilter.h"
#pragma comment(lib, "libfftw3-3.lib")
#pragma comment(lib, "libfftw3f-3.lib")
#pragma comment(lib, "libfftw3l-3.lib")
#define N 8
typedef itk:: Image<float , 2> FloatImageType;
typedef itk ::Image< unsigned char , 2> UnsignedCharImageType;
typedef float MatrixType;
void CopyMatrix(FloatImageType:: Pointer image , MatrixType **mask);
void CopyImage(FloatImageType:: Pointer image , fftw_complex *mask);
void CopyRealImage (FloatImageType:: Pointer image , float ** mask);//仅仅拷贝实部
void MyCastFloatToUnsichar(FloatImageType::Pointer image, UnsignedCharImageType::Pointer imageNew );//将float转为unsigned char型
int main()
{
//读取图像
FloatImageType::Pointer image;
const char * inputFileName = "F:\\DIPcode\\SIMC\\imagesrc\\2.jpg" ; //argv[1];
const char * outputFileName = "F:\\DIPcode\\SIMC\\2_new.tif" ; //argv[2];
typedef itk ::ImageFileReader< FloatImageType> ReaderType ;
ReaderType::Pointer reader = ReaderType::New ();
reader->SetFileName (inputFileName);
reader->Update ();
image = reader ->GetOutput();
//遍历图像像素，赋值到一个矩阵
MatrixType **pixelColors = new MatrixType*[512]();
for (int i=0;i<512;i++)
{
pixelColors[i] = new MatrixType[512]; //为每个指针都分配8个char元素空间。
}
CopyMatrix(image , pixelColors);
int width =512;
int height=512;
fftw_plan planR;
fftw_complex *inR, *outR, *resultR;
//Allocate arrays.
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
//int currentIndex = ((y * width) + (x)) * 3;
int currentIndex = (y * width) + x;
inR[y * width + x][0] = (double)pixelColors[y][x];
inR[y * width + x][1] = 0.0;
}
}
//Forward plans.
planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE);
//Forward FFT.
fftw_execute(planR);
//Backward plans.
fftw_plan planRR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE);
//Backward fft
fftw_execute(planRR);
//// Have to scale the output values to get back to the original.
int size=height*width;
for(int i = 0; i < size; i++) {
resultR[i][0] = resultR[i][0] / size;
resultR[i][1] = resultR[i][1] / size;
}
//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
//int currentIndex = ((y * width) + (x)) * 3;
pixelColors[y][x] =resultR[y * width + x][0];//
//std::cout<<pixelColors[y][x]<<'\t';
}
}
//将图像写出来
CopyRealImage(image , pixelColors);//将pixelColors中的像素值复制到image中
//调整亮度
typedef itk ::RescaleIntensityImageFilter< FloatImageType , FloatImageType > RescaleFilterType;
RescaleFilterType::Pointer imaginaryRescaleFilter = RescaleFilterType ::New();
imaginaryRescaleFilter->SetInput (image);
imaginaryRescaleFilter->SetOutputMinimum (0);
imaginaryRescaleFilter->SetOutputMaximum (255);
imaginaryRescaleFilter->Update ();
//转为8bit
//类型转换，float转为8bit的。
typedef itk ::CastImageFilter< FloatImageType, UnsignedCharImageType > CastFilterType;
CastFilterType::Pointer castFilter = CastFilterType::New ();
castFilter->SetInput (imaginaryRescaleFilter->GetOutput());
castFilter->Update ();
//将图像写出来
typedef itk ::ImageFileWriter< UnsignedCharImageType > WriterType ;
WriterType::Pointer writer = WriterType::New ();
writer->SetFileName ( outputFileName );
writer->SetInput (castFilter->GetOutput());
try
{
writer->Update ();
}
catch( itk ::ExceptionObject & error )
{
std::cerr << "Error: " << error << std ::endl;
return EXIT_FAILURE ;
}
delete []pixelColors;
fftw_free(inR);
fftw_free(outR);
fftw_free(resultR);
return 0;
}
void CopyMatrix(FloatImageType:: Pointer image , MatrixType **mask)
{
FloatImageType::RegionType region = image->GetLargestPossibleRegion ();
FloatImageType::SizeType regionSize = region.GetSize ();
for (int x=0; x<regionSize [0]; x++)
for(int y=0; y<regionSize [1]; y++)
{ //3、定义像素索引
FloatImageType::IndexType index;
index[0] = x ;
index[1] = y ;
mask[x][y] =(MatrixType) image->GetPixel(index);//获取该位置的像素值
}
}
void CopyImage(FloatImageType:: Pointer image , fftw_complex *mask)
{
FloatImageType::RegionType region = image->GetLargestPossibleRegion ();
FloatImageType::SizeType regionSize = region.GetSize ();
for (int x=0; x<regionSize [0]; x++)
for(int y=0; y<regionSize [1]; y++)
{
double re = mask[x*regionSize[1]+y][0];
double im = mask[x*regionSize[1]+y][1];
double mag = sqrt (re*re + im*im);
//3、定义像素索引
FloatImageType::IndexType index;
index[0] = x ;
index[1] = y ;
image->SetPixel(index, mag);//获取该位置的像素值
}
}
//仅仅将mask中的值，赋给image中！
void CopyRealImage (FloatImageType:: Pointer image , MatrixType ** mask)
{
FloatImageType::RegionType region = image->GetLargestPossibleRegion ();
FloatImageType::SizeType regionSize = region.GetSize ();
int Imgsize = regionSize[0]*regionSize[1];
for (int x=0; x<regionSize [0]; x++)
for(int y=0; y<regionSize [1]; y++)
{ //3、定义像素索引
FloatImageType::IndexType index;
index[0] = x ;
index[1] = y ;
image->SetPixel (index, mask[x ][y]);
//获取该位置的像素值
//std::cout<<image->GetPixel(index);
}
}
void MyCastFloatToUnsichar(FloatImageType::Pointer image, UnsignedCharImageType::Pointer imageNew )//将float转为unsigned char型
{
FloatImageType::RegionType region = image->GetLargestPossibleRegion ();
FloatImageType::SizeType regionSize = region.GetSize ();
for (int x=0; x<regionSize [0]; x++)
for(int y=0; y<regionSize [1]; y++)
{ //3、定义像素索引
FloatImageType::IndexType index;
index[0] = x ;
index[1] = y ;
//std::cout<<image->GetPixel(index)<<'\t';
float temp = static_cast<unsigned char>(image->GetPixel(index));
//std::cout<<temp<<'\t';
imageNew->SetPixel (index, temp);
int tt=imageNew->GetPixel(index);
////获取该位置的像素值
//std::cout<<tt<<'\t';
}
}

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

0 0