基於ITK和VTK實現三維體資料的區域生長分割和視覺化
阿新 • • 發佈:2019-01-02
該文由Markdown語法編譯器編輯完成。
1. 前言:
在醫學影像的開源庫中,ITK主要擅長影象分割和影象配准算法的 研究,VTK則擅長三維視覺化的實現。通過結合二者,可以實現基本的影象分割或配准算法的執行和結果顯示。
本文主要介紹ITK中的基於itkConnectedThresholdFilter的區域生長分割演算法,並且利用 VTK將區域生長分割後的體資料進行顯示。
2. 基本步驟:
本文所採用的基本開發環境為:
序號 | 軟體 | 版本號 |
---|---|---|
1 | ITK | 4.10.0 |
2 | VTK | 6.0.0 |
3 | QT | 5.1.0 |
4 | Visual Studio | 2012 |
5 | OS | Windows 7 |
這篇文章是基於當前我正在公司開發的一個專案中抽取的部分需求來撰寫的。其中的需求是輸入患者的一組CT資料,首先基於閾值進行分割,然後在閾值分割的基礎上,再進行區域生長分割,這樣可以快速地獲取使用者需要後處理的感興趣區域。
由於要處理的是醫學體資料,而且是已經進行完閾值分割後的資料,因此它的資料型別設為vtkImageData。
由於VTK和ITK的資料型別不能直接使用,因此首先需要將VTK型別的體資料vtkImageData轉化成itkImage,然後再呼叫ITK中的區域生長的演算法,最後再將分割後的資料再轉換回vtkImageData,進行後續的三維後處理等。
好在,ITK中已經提供了轉化資料型別的filter,只需根據需要選擇使用即可。
以下是主要的實現程式碼:
//.h檔案:
static const int kDimension=3;
typedef itk::Image< int, kDimension > ImageType;
typedef ImageType::IndexType PixelIndexType;
ImageType::Pointer m_itkImage;
ImageType::IndexType m_pixelIndex;
//.cpp檔案:
void ConstructITKImage(const vtkImageData* originalImageData)
{
double imageSpacing[3] = {0.0};
double imageOrigin[3] = {0.0};
int imageDimension[3] = {0};
originalImageData->GetSpacing(imageSpacing);
originalImageData->GetOrigin(imageOrigin);
originalImageData->GetDimensions(imageDimension);
const ImageType::SizeType size = {{imageDimension[0], imageDimension[1], imageDimension[2]}}; //Size along {X,Y,Z}
const ImageType::IndexType start = {{0,0,0}}; // First index on {X,Y,Z}
ImageType::RegionType region;
region.SetSize( size );
region.SetIndex( start );
//Types for converting between ITK and VTK
typedef itk::VTKImageToImageFilter<ImageType> VTKImageToImageType;
//Converting to ITK Image Format
VTKImageToImageType::Pointer vtkImageToImageFilter = VTKImageToImageType::New();
vtkImageToImageFilter->SetInput(originalImageData);
vtkImageToImageFilter->UpdateLargestPossibleRegion();
//將vtk的影象轉化為itk的影象,以便利用itk的分割演算法進行分割.
m_itkImage->SetRegions( region );
m_itkImage->Allocate();
m_itkImage = const_cast<itk::Image<int, kDimension>*>(vtkImageToImageFilter->GetImporter()->GetOutput());
ImageType::SpacingType spacing;
spacing[0] = imageSpacing[0]; // spacing along X
spacing[1] = imageSpacing[1]; // spacing along Y
spacing[2] = imageSpacing[2]; // spacing along Z
m_itkImage->SetSpacing( spacing );
ImageType::PointType newOrigin;
newOrigin[0] = imageOrigin[0];
newOrigin[1] = imageOrigin[1];
newOrigin[2] = imageOrigin[2];
m_itkImage->SetOrigin( newOrigin );
ImageType::DirectionType direction;
direction.SetIdentity();
m_itkImage->SetDirection( direction );
m_itkImage->UpdateOutputInformation();
};
bool CheckSeedPointInsideImage( double worldSeedPoint[3] )
{
typedef itk::Point< double, ImageType::ImageDimension > PointType;
PointType point;
point[0] = worldSeedPoint[0]; // x coordinate
point[1] = worldSeedPoint[1]; // y coordinate
point[2] = worldSeedPoint[2]; // z coordinate
//根據選取的種子點的世界座標,獲取這個種子點在itk影象中的畫素索引值.
bool isInside = m_itkImage->TransformPhysicalPointToIndex(point, m_pixelIndex);
return isInside;
}
int RegionGrowingByConnectedThredholdFilter(const vtkImageData* pImage, double worldSeedPoint[3])
{
ConstructITKImage(originalImageData);
bool isSeedPointInsideImage = CheckSeedPointInsideImage(worldSeedPoint);
if (!isSeedPointInsideImage)
{
return;
}
typedef itk::Image< int, kDimension > InternalImageType;
typedef itk::ConnectedThresholdImageFilter< InternalImageType, InternalImageType > ConnectedFilterType;
ConnectedFilterType::Pointer connectedThreshold = ConnectedFilterType::New( );
connectedThreshold->SetInput(m_itkImage);
//由於區域生長是在閾值分割後進行的,此時的vtkImageData是二值的,灰度只有0和1。因此基於閾值的區域生長的灰度上下限都是1.
//區域生長分割後的影象,要經過vtkLookupTable進行顏色對映.vtkLookupTable也只有0和1兩個預設值.
connectedThreshold->SetLower(1);
connectedThreshold->SetUpper(1);
connectedThreshold->SetReplaceValue(1);
connectedThreshold->SetSeed(m_pixelIndex);
typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
//Converting Back from ITK to VTK Image for Visualization.
ConnectorType::Pointer connector = ConnectorType::New();
connector->SetInput(connectedThreshold->GetOutput());
connector->Update();
}
;
3. 區域生長分割結果及顯示:
完。