结合多篇论文,得到了常规流程:

在使用分割算法对二维图像序列进行分割后,将分割后的二维图像按照分割前一样的Marching cube算法就能重建得到需要的三维面数据。

#include <iostream>
#include <string>
#include <vector>
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageToVTKImageFilter.h"
#include "itkVTKImageToImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkImageFileWriter.h"
#include "vtkSmartPointer.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkStructuredPoints.h"
#include "vtkStructuredPointsWriter.h"

typedef unsigned short PixelType;
const int Dimension = 3;
typedef itk::Image<PixelType, Dimension> ImageType;
typedef itk::ImageSeriesReader<ImageType> ReaderType;
typedef itk::ImageSeriesWriter<ImageType, itk::ImageFileWriter<ImageType>> WriterType;
typedef itk::GDCMImageIO ImageIOType;
typedef itk::GDCMSeriesFileNames NamesGeneratorType;
typedef itk::ImageToVTKImageFilter<ImageType> ImageToVTKFilterType;
typedef itk::VTKImageToImageFilter<ImageType> VTKToImageFilterType;
typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> ThresholdFilterType;
typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectedComponentFilterType;
typedef itk::RelabelComponentImageFilter<ImageType, ImageType> RelabelFilterType;
typedef itk::RescaleIntensityImageFilter<ImageType, ImageType> RescaleFilterType;

int main(int argc, char *argv[])
{
    if (argc < 2) {
        std::cerr << "Usage: " << argv[0] << " dicom_directory" << std::endl;
        return EXIT_FAILURE;
    }

    // Step 1: Read DICOM series
    std::string directory = argv[1];
    ImageIOType::Pointer gdcmIO = ImageIOType::New();
    NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
    namesGenerator->SetInputDirectory(directory);
    const ReaderType::FileNamesContainer &fileNames = namesGenerator->GetInputFileNames();
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetImageIO(gdcmIO);
    reader->SetFileNames(fileNames);
    try {
        reader->Update();
    }
    catch (itk::ExceptionObject &ex) {
        std::cerr << "Exception caught!" << std::endl;
        std::cerr << ex << std::endl;
        return EXIT_FAILURE;
    }

    // Step 2: Thresholding
    ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New();
    thresholdFilter->SetInput(reader->GetOutput());
    thresholdFilter->SetLowerThreshold(200);
    thresholdFilter->SetUpperThreshold(255);
    thresholdFilter->SetInsideValue(1);
    thresholdFilter->SetOutsideValue(0);
    thresholdFilter->Update();

    // Step 3: Connected component labeling
    ConnectedComponentFilterType::Pointer connectedComponentFilter = ConnectedComponentFilterType::New();
    connectedComponentFilter->SetInput(thresholdFilter->GetOutput());
    connectedComponentFilter->Update();

    // Step 4: Relabeling
    RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
    relabelFilter->SetInput(connectedComponentFilter->GetOutput());
    relabelFilter->Update();

    // Step 5: Rescale intensity
    RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
    rescaleFilter->SetInput(relabelFilter->GetOutput());
    rescaleFilter->SetOutputMinimum(0);
    rescaleFilter->SetOutputMaximum(255);
    rescaleFilter->Update();

    // Step 6: Write output series
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput(rescaleFilter->GetOutput());
    writer->SetImageIO(gdcmIO);
    writer->SetFileNames(namesGenerator->GetOutputFileNames());
    writer->Update();

    // Step 7: Convert to VTK image
    ImageToVTKFilterType::Pointer imageToVTKFilter = ImageToVTKFilterType::New();
    imageToVTKFilter->SetInput(rescaleFilter->GetOutput());
    imageToVTKFilter->Update();

    // Step 8: Convert to VTK polydata
    vtkSmartPointer<vtkStructuredPoints> structuredPoints = vtkSmartPointer<vtkStructuredPoints>::New();
    structuredPoints->SetDimensions(imageToVTKFilter->GetOutput()->GetDimensions());
    structuredPoints->SetSpacing(imageToVTKFilter->GetOutput()->GetSpacing());
    structuredPoints->SetOrigin(imageToVTKFilter->GetOutput()->GetOrigin());
    structuredPoints->SetScalarTypeToUnsignedChar();
    structuredPoints->GetPointData()->SetScalars(imageToVTKFilter->GetOutput()->GetPointData()->GetScalars());

    vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkStructuredPointsWriter> writer2 = vtkSmartPointer<vtkStructuredPointsWriter>::New();
    writer2->SetInputData(structuredPoints);
    writer2->SetFileName("blood_vessels.vtk");
    writer2->Write();

    return EXIT_SUCCESS;
}

该程序主要分为以下几个步骤:

1.使用ITK库读取DICOM序列图像。

2.使用二值阈值滤波器将图像二值化。

3.使用连通成分标记滤波器标记和提取感兴趣区域。

4.使用重新标记滤波器对连通成分进行重新标记。

5.使用强度重缩放滤波器将像素值范围从[0,1]重缩放到[0,255]。

6.使用ITK库写入DICOM序列图像。

7.使用ITK-VTK转换器将ITK图像转换为VTK图像。

8.将VTK图像转换为VTK多边形数据,并使用VTK库将数据写入VTK文件。

 

Logo

魔乐社区(Modelers.cn) 是一个中立、公益的人工智能社区,提供人工智能工具、模型、数据的托管、展示与应用协同服务,为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作,由全产业链共同建设、共同运营、共同享有,推动国产AI生态繁荣发展。

更多推荐