【ITK库学习】使用itk库进行图像滤波ImageFilter：频域滤波

1、itkFFTConvolutionImageFilter 快速傅里叶变换计算

该类使用傅里叶域中的乘法将给定图像与任意图像内核进行卷积。

该过滤器产生的输出相当于 ConvolutionImageFilter 的输出，当核较大时，它利用卷积定理来加速卷积计算。由于FFT的特性，它可以大大加速卷积计算，尤其是对于卷积核较大的情况。但是对于某些特殊的卷积核，如具有特定形状或非对称性的核，可能会导致较大的计算误差。

警告
该过滤器忽略内核图像的间距、原点和方向，并将它们视为与输入图像中的相同。

常用的成员函数：

• SetKernelImage()：设置卷积内核

示例代码：

#include "itkImage.h"
#include "itkFFTConvolutionImageFilter.h"

typedef itk::Image<float, 3> FloatImageType;

bool FFTConvolutionImageFilter(FloatImageType* image, FloatImageType* outputImage)
{
    typename FloatImageType::Pointer kernelImage = FloatImageType::New();
    unsigned int width = 10;
	createKernel(kernelImage, width);
	
    typedef itk::FFTConvolutionImageFilter<FloatImageType, FloatImageType> FFTConvolutionFilterType;
	typename FFTConvolutionFilterType::Pointer fftFilter = FFTConvolutionFilterType::New();
	fftFilter->SetInput(image);
	fftFilter->SetKernelImage(kernelImage);
	try
	{
		fftFilter->Update();
	}
	catch (itk::ExceptionObject& ex)
	{
		//读取过程发生错误
		std::cerr << "Error: " << ex << std::endl;
		return false;
	}

	outputImage = fftFilter->GetOutput();
	return true;
}


void CreateKernel(FloatImageType::Pointer kernel, unsigned int width)
{
  FloatImageType::IndexType start;
  start.Fill(0);
 
  FloatImageType::SizeType size;
  size.Fill(width);
 
  FloatImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
 
  kernel->SetRegions(region);
  kernel->Allocate();
 
  itk::ImageRegionIterator<FloatImageType> imageIterator(kernel, region);
 
  while (!imageIterator.IsAtEnd())
  {
    // imageIterator.Set(255);
    imageIterator.Set(1);
 
    ++imageIterator;
  }
}

2、扩展：itkConvolutionImageFilter.h

该类将给定图像与任意图像内核进行卷积。

该滤波器的工作原理是将翻转后的内核集中在图像中的每个像素处，并计算图像中的像素值与内核中的像素值之间的内积。

内核的中心定义为[(2*i+s?1)/2]，其中 i 是索引，s 是内核图像的最大可能区域的大小。

对于所有维度上具有奇数大小的内核，这对应于中心像素，如果核图像的某个维度具有偶数大小，则该维度中核的中心索引将是小于图像中心的连续索引的最大积分索引。

可以选择使用 NormalizeOn() 将内核标准化为总和为 1，默认情况下标准化处于关闭状态。

警告
该过滤器忽略内核图像的间距、原点和方向，并将它们视为与输入图像中的相同。

其常用的成员函数与代码结构与FFTConvolutionImageFilter类似，参见1中内容

3、itkFFTShiftImageFilter 频率转移滤波器

将傅里叶变换的零频率分量移至图像的中心。

傅立叶变换生成的图像中，零频率分量位于图像的角落，使其难以理解，该滤波器将组件移动到图像的中心。

注意：对于奇数尺寸的图像，如果不对两个过滤器中的一个（且仅一个）使用 SetInverse(true)，则应用此过滤器两次将不会生成与原始图像相同的图像。

常用的成员函数：

• Set/GetInverse()：设置/获取过滤器是否必须反转变换，这个选项对于输入图像的尺寸都是偶数的情况没有影响，但是如果至少有一个维度的尺寸是奇数，则需要使用这个选项来恢复原始图像，简而言之，这个选项的作用是针对尺寸为奇数的输入图像进行处理，使其尺寸与原始图像保持一致
• InverseOn/Off()：同上

示例代码：

#include "itkImage.h"
#include "itkFFTShiftImageFilter.h"

typedef itk::Image<float, 3> FloatImageType;

bool FFTShiftImageFilter(FloatImageType* image, FloatImageType* outputImage)
{
    typedef itk::FFTShiftImageFilter<FloatImageType, FloatImageType> FFTShiftFilterType;
	typename FFTShiftFilterType::Pointer fftFilter = FFTShiftFilterType::New();
	fftFilter->SetInput(image);
	try
	{
		fftFilter->Update();
	}
	catch (itk::ExceptionObject& ex)
	{
		//读取过程发生错误
		std::cerr << "Error: " << ex << std::endl;
		return false;
	}

	outputImage = fftFilter->GetOutput();
	return true;
}

4、itkFFTNormalizedCorrelationImageFilter FFT实现的归一化相关滤波器

使用 FFT 计算归一化互相关。

该滤波器使用 FFT 而不是空间相关来计算两个图像的归一化互相关 (NCC)，对于相当大的结构元素，它比空间关联要快得多，该过滤器是更通用的 MaskedFFTNormalizedCorrelationImageFilter 的子类，其操作原理是将该算法中的掩码设置为图像，正如下面的参考文献中详细描述的，利用更通用的掩模算法没有计算消耗，因为计算仍然需要对图像进行 FFT。

输入：需要两个图像作为输入，fixedImage 和 movingImage，在相关性的背景下，输入通常被定义为：“图像”和“模板”，该过滤器能够关联任意两个图像，并且不限于小型移动图像（模板）。

可选参数：RequiredNumberOfOverlappingPixels 使用户能够指定两个图像的多少体素必须重叠； 相关图中由少于此数量的体素产生的任何位置将被设置为零，较大的值会将相关图像周围较大边界上的像素清零，因此，较大的值消除了不太稳定的计算，但也限制了捕获范围，如果“RequiredNumberOfOverlappingPixels”设置为 0（默认值），则不会进行归零。

图像大小：fixedImage 和 movingImage 不必具有相同的大小。 此外，虽然一些算法要求“模板”小于“图像”，因为两者不完全重叠的区域中存在误差，但该过滤器没有这样的限制。

图像间距：由于计算是在像素域中完成的，因此所有输入图像必须具有相同的间距。

输出： 输出是 RealPixelType 的图像，即两幅图像的 NCC，其值范围为 -1.0 到 1.0。 根据定义，该 NCC 图像的大小为 size(fixedImage) + size(movingImage) - 1。

常用的成员函数：

• Set/GetInverse()：设置/获取过滤器是否必须反转变换，这个选项对于输入图像的尺寸都是偶数的情况没有影响，但是如果至少有一个维度的尺寸是奇数，则需要使用这个选项来恢复原始图像，简而言之，这个选项的作用是针对尺寸为奇数的输入图像进行处理，使其尺寸与原始图像保持一致
• InverseOn/Off()：同上

示例代码：

#include "itkImage.h"
#include "itkFFTNormalizedCorrelationImageFilter.h"

typedef itk::Image<float, 3> FloatImageType;

bool FFTNormalizedCorrelationImageFilter(FloatImageType* image, FloatImageType* outputImage)
{
	itk::Index<3> offset;
	offset[0] = 5;
	offset[1] = 6;
	offset[2] = 6;

	typename FloatImageType::Pointer fixedImage = FloatImageType::New();
	itk::Index<3> cornerOfFixedSquare;
	cornerOfFixedSquare[0] = 3;
	cornerOfFixedSquare[1] = 8;
	cornerOfFixedSquare[2] = 6;
	createImage(fixedImage, cornerOfFixedSquare);

	typename FloatImageType::Pointer movingImage = FloatImageType::New();
	itk::Index<3> cornerOfMovingSquare;
	cornerOfMovingSquare[0] = cornerOfFixedSquare[0] + offset[0];
	cornerOfMovingSquare[1] = cornerOfFixedSquare[1] + offset[1];
	cornerOfMovingSquare[2] = cornerOfFixedSquare[2] + offset[2];
	createImage(movingImage, cornerOfMovingSquare);
	
    typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> FFTNormalizedCorrelationFilterType;
	typename FFTNormalizedCorrelationFilterType::Pointer fftFilter = FFTNormalizedCorrelationFilterType::New();
	fftFilter->SetFixedImage(fixedImage);
	fftFilter->SetMovingImage(movingImage);
	try
	{
		fftFilter->Update();
	}
	catch (itk::ExceptionObject& ex)
	{
		//读取过程发生错误
		std::cerr << "Error: " << ex << std::endl;
		return false;
	}

	outputImage = fftFilter->GetOutput();
	return true;
}

void createImage(FloatImageType::Pointer image, const itk::Index<3>& cornerOfSquare)
{
	FloatImageType::IndexType start;
	start.Fill(0);

	FloatImageType::SizeType size;
	size.Fill(51);

	FloatImageType::RegionType region(start, size);

	image->SetRegions(region);
	image->Allocate();
	image->FillBuffer(0);

	itk::ImageRegionIterator<FloatImageType> imageIterator(image, region);

	FloatImageType::IndexValueType squareSize = 8;

	while (!imageIterator.IsAtEnd())
	{
		if (imageIterator.GetIndex()[0] > cornerOfSquare[0] &&
			imageIterator.GetIndex()[0] < cornerOfSquare[0] + squareSize &&
			imageIterator.GetIndex()[1] > cornerOfSquare[1] && imageIterator.GetIndex()[1] < cornerOfSquare[1] + squareSize)
		{
			imageIterator.Set(255);
		}
		++imageIterator;
	}
}