使用ITK和VTK的CT脑血管分割

Question

我正在做自己的项目，需要分割CT dicom图像的大脑血管。我有2套DICOM图像。

1)CT dicom images of head with DSA (888 slices)
2)CT dicom images of head without DSA(same 888 slices),

所以我用vtkDICOMImageReader分别读取了这两个系列，并用vtkImageMathematics减去了这两个系列，并将结果存储在变量mat中，然后将其复制到vtkImageReslice * reslice;然后，我使用reslice将此VKImageToImageFilter *filter_toitkimage导入到ITK中，并使用以下代码：

using HessianFilterType = itk::HessianRecursiveGaussianImageFilter<ImageType>;
HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
hessianFilter->SetInput(filter_toitkimage->GetOutput());
hessianFilter->SetSigma(.6);

using VesselnessMeasureFilterType = itk::Hessian3DToVesselnessMeasureImageFilter<PixelType>;
VesselnessMeasureFilterType::Pointer vesselnessFilter = VesselnessMeasureFilterType::New();
vesselnessFilter->SetInput(hessianFilter->GetOutput());
vesselnessFilter->SetAlpha1(1);
vesselnessFilter->SetAlpha2(2); 

然后我使用ImageToVTKImageFilter将其转换回VTK，并使用了此：

       vtkSmartPointer<vtkGPUVolumeRayCastMapper> volumeMapper = vtkSmartPointer<vtkGPUVolumeRayCastMapper>::New();
volumeMapper->CroppingOff();
volumeMapper->RemoveAllInputs();
volumeMapper->SetInputData(conv->GetOutput());   //conv is the converted output   
volumeMapper->SetBlendModeToMaximumIntensity();
volumeMapper->AutoAdjustSampleDistancesOff();
volumeMapper->SetSampleDistance(.1);

volumeMapper->Update();

//-----volume color
vtkSmartPointer<vtkColorTransferFunction> volumeColor = vtkSmartPointer<vtkColorTransferFunction>::New();
volumeColor->RemoveAllPoints();

////-----volume color

//---air
volumeColor->AddRGBPoint(0, 0.0, 0.0, 0.0);
volumeColor->AddRGBPoint(220, 1, 1, 1);
//---skin
volumeColor->AddRGBPoint(250, 0, 0, 0);
volumeColor->AddRGBPoint(300, 0, 0, 0);
//---muscle
volumeColor->AddRGBPoint(800, 1, 1, 1);
volumeColor->AddRGBPoint(1000, 1, 1, 1);
volumeColor->AddRGBPoint(1320, 1, 1, 1);
//---soft bone
volumeColor->AddRGBPoint(1350, 0, 0, 0);
volumeColor->AddRGBPoint(1500, 0, 0, 0);
//---hard bone
volumeColor->AddRGBPoint(1550, 0, 0, 0);
volumeColor->AddRGBPoint(2000, 0, 0, 0);
//---hardest bone
volumeColor->AddRGBPoint(2050, 0, 0, 0);
volumeColor->AddRGBPoint(4096, 0, 0, 0);

vtkSmartPointer<vtkPiecewiseFunction> volumeScalarOpacity =
    vtkSmartPointer<vtkPiecewiseFunction>::New();
volumeScalarOpacity->AddPoint(0, 0.00);
volumeScalarOpacity->AddPoint(500, 0.15);
volumeScalarOpacity->AddPoint(1000, 0.15);
volumeScalarOpacity->AddPoint(1150, 0.85);

vtkSmartPointer<vtkPiecewiseFunction> volumeGradientOpacity =
    vtkSmartPointer<vtkPiecewiseFunction>::New();
volumeGradientOpacity->AddPoint(0, 0.0);
volumeGradientOpacity->AddPoint(90, 0.5);
volumeGradientOpacity->AddPoint(100, 1.0);

vtkSmartPointer<vtkVolumeProperty> volumeProperty = vtkSmartPointer<vtkVolumeProperty>::New();
volumeProperty->SetColor(volumeColor);
    volumeProperty->SetScalarOpacity(volumeScalarOpacity);
    volumeProperty->SetGradientOpacity(volumeGradientOpacity);
    volumeProperty->SetInterpolationTypeToLinear();
    volumeProperty->ShadeOn();
    volumeProperty->SetAmbient(0.5);
    volumeProperty->SetDiffuse(0.3);
    volumeProperty->SetSpecular(0.3);

    vtkSmartPointer<vtkVolume> volume =vtkSmartPointer<vtkVolume>::New();
    volume->SetMapper(volumeMapper);
    volume->SetProperty(volumeProperty);
    renderer->AddVolume(volume);

我正在获取肺部血管，但我没有清晰的脑血管图像。有人可以提出建议吗？

Answer 1

您的传递函数假设为CT强度，但是血管过滤器的输出却与此不同。可能将其标准化为0.0-1.0范围，0-255范围或其他值。您应该检查一下。

Answer 2

DSA已经表示“数字减影血管造影”。因此，您不应再执行该操作，它已经由扫描仪的软件完成。您正在寻找的可能是差异的一些中间步骤，但是我除非您想进行新的CT重建，否则我建议您直接使用DSA，而不要再做一次差异。请注意，在PACS中，它可能被称为血管造影的“ XA”模态（您不应该寻找“ CT”模态）。

[每个系统都不同，这可能很复杂，但是我也在研究大脑的血管造影数据，起初这让我感到困惑，所以我认为值得一提。

关于您的血管检测的原型，我建议第一步使用3D Slicer。有一个名为VMTK的库的插件，该插件实现了容器过滤器，但是您也可以添加其他插件，尤其是获得所有ITK过滤器的插件。我认为以它为起点，在实施之前尝试一些东西会更容易。

如果您真的想自己计算DSA，在这种情况下，您需要有2个带有和不带有造影剂的CT（而不是DSA），那么在进行差异之前，必须确保两个体积均已完美记录。如果发现未对齐，我建议您使用简单的elastix框架注册图像。

最后，容器性滤波器的临界点实际上是您使用sigma选择的频率。这基本上定义了容器的大小。在VTMK实施中，您可以使用不同的频率，因此您可以设置开始和结束sigma。在Slicer中，如果您在血管周围输入一些基准，甚至可以找到合适的sigma插件。因此，如果您看到一些血管，但看不到其他血管，我建议将sigma = 0.6降低一些（假设血管比肺血管小）

我希望这会有所帮助！