ITK笔记：SetRequestedRegion设置ROI无效，滤波器仍然处理了全图

在医学图像处理中，有时仅希望处理图像中的感兴趣区域(ROI)。ITK有三种区域：(1) LargestPossibleRegion：指的是图像最大可能区域，通常是图像整体；(2) BufferedRegion：指的是实际存在于内存中的图像区域；(3) RequestedRegion：要求被处理的图像区域。三种区域可分别通过Set方法进行设置。

按照ITK用户手册[1]中的描述，通过设置RequestedRegion，就可以实现图像指定区域处理。ITK网站提供的示例代码[2]如下：

#include "itkImage.h"#include "itkRandomImageSource.h"#include "itkDerivativeImageFilter.h"int main(int, char *[]){typedef itk::Image<float, 2> ImageType;itk::Size<2> smallSize;smallSize.Fill(10);itk::Index<2> index;index.Fill(0);itk::ImageRegion<2> region(index, smallSize);itk::Size<2> bigSize;bigSize.Fill(10000);itk::RandomImageSource<ImageType>::Pointer randomImageSource = itk::RandomImageSource<ImageType>::New();randomImageSource->SetNumberOfThreads(1); // to produce non-random resultsrandomImageSource->SetSize(bigSize);randomImageSource->GetOutput()->SetRequestedRegion(smallSize);randomImageSource->Update();std::cout << "Created random image." << std::endl;typedef itk::DerivativeImageFilter<ImageType, ImageType > DerivativeImageFilterType;DerivativeImageFilterType::Pointer derivativeFilter = DerivativeImageFilterType::New();derivativeFilter->SetInput( randomImageSource->GetOutput() );derivativeFilter->SetDirection(0); // "x" axisderivativeFilter->GetOutput()->SetRequestedRegion(smallSize);derivativeFilter->Update();std::cout << "Computed derivative." << std::endl;return EXIT_SUCCESS;}

参照这份示例代码，我写了一个简单的阈值分割demo，并且加了一个转换到VTK图像并显示的步骤：

#include "itkImage.h"#include "itkImageSeriesReader.h"#include "itkImageToVTKImageFilter.h"#include "itkGDCMImageIO.h"#include "itkGDCMSeriesFileNames.h"#include "itkRescaleIntensityImageFilter.h"#include "itkCastImageFilter.h"#include "itkBinaryThresholdImageFilter.h"#include "vtkSmartPointer.h"#include "vtkImageViewer2.h"#include "vtkImageData.h"#include "vtkRenderWindow.h"#include "vtkRenderWindowInteractor.h"#include "vtkRenderer.h"#include "vtkInteractorStyleImage.h"typedef float                                 InternalPixelType;typedef std::string                           FilenameType;typedef itk::Image<InternalPixelType, 3>      InternalImageType;typedef itk::Size<3>                          SizeType;typedef itk::Index<3>                         IndexType;typedef itk::ImageRegion<3>                   RegionType;InternalImageType::Pointer LoadDICOM(const FilenameType & dir);//自己写的读取DICOM文件的函数int main(){        // image size: [640,640,100]        InternalImageType::Pointer image = LoadDICOM("D:\\Workspace\\001");                typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> ThresholdFilterType;        ThresholdFilterType::Pointer thrdFilter = ThresholdFilterType::New();        IndexType index{ 0,0,0 };        SizeType size{ 100,100,50 };        RegionType requestedRegion(index, size);                thrdFilter->SetInput(image);        thrdFilter->SetUpperThreshold(100);        thrdFilter->SetInsideValue(255);        thrdFilter->SetOutsideValue(0);        thrdFilter->GetOutput()->SetRequestedRegion(requestedRegion);        thrdFilter->Update();        typedef itk::ImageToVTKImageFilter<InternalImageType> ConnectType;        ConnectType::Pointer connector = ConnectType::New();        connector->SetInput(thrdFilter->GetOutput());        connector->Update();        vtkSmartPointer<vtkImageViewer2> ImageViewer = vtkSmartPointer<vtkImageViewer2>::New();        ImageViewer->SetInputData(connector->GetOutput());        ImageViewer->SetSlice(50);        ImageViewer->Render();        vtkSmartPointer<vtkRenderWindowInteractor> iren = vtkSmartPointer<vtkRenderWindowInteractor>::New();        iren->SetRenderWindow(ImageViewer->GetRenderWindow());        iren->Initialize();        iren->Start();        return EXIT_SUCCESS;}

然后就愉快地运行了。结果，滤波器仍然处理了整个图像范围，根本没有理会我设置的RequestedRegion：

为了弄清原因，我在Update方法处设置断点，单步执行。一直到阈值分割结束，RequestedRegion都是processRegion，然后继续执行，输出结果就是整个图像都被处理。我被这个奇怪的bug折腾了一整天，终于在kitware论坛上搜到了一个类似的话题讨论[3]。其中一个回复作者提到，输出连接到writer，就会导致整个图像范围被处理，而断开管线就可以解决。

问题的关键在后一级的滤波器将RequestedRegion传递到了前级，导致前级滤波器信息更新。按照ITK管线的机制，前级滤波器又一次执行滤波操作。现在回到用户手册，8.3.2节提到了管线执行过程的细节：DataObject::Update()方法依次调用了以下三个方法：

(1) DataObject::UpdateOutputInformation()

(2) DataObject::PropagateRequestedRegion()

(3) DataObject::UpdateOutputData()

其中，第(2)个方法将当前滤波器的RequestedRegion传递到了前级，并且覆盖了原先的RequestedRegion：

为了证实这一点，我也试验了一下，加了一个指针probe，代码如下： 

观察probe的RequestedRegion：

第一个断点：

第二个断点：

果然变了...

本来是打算用RegionOfInterestFilter的，发现这个滤波器把ROI从原图上“挖”出来了，可我又想保留原图尺寸，只在ROI内运行我现在要跑的算法。这么作实在是因为原图太大了，跑起来电脑要跪。看到可以直接设置RequestedRegion还高兴了好一会，结果竟是这样。

目前的解决方案是用RegionOfInterestFilter把ROI提取出来，处理完毕再用PasteImageFilter给贴回一个和原图一样大的空白图上。暂时也没想到更好的。

Reference:

[1] The ITK Software Guide Book 1: Introduction and Development Guidelines Fourth Edition Updated for ITK version 4.11, chap 8.1-8.3

[2] https://itk.org/Doxygen412/html/WikiExamples_2SimpleOperations_2RequestedRegion_8cxx-example.html

[3] https://public.kitware.com/pipermail/insight-users/2010-February/035139.html