CT Brain vessel segmentation using ITK and VTK

Question

I am doing my own project in which i need to segment brain vessels of CT dicom images. I have 2 sets of DICOM images .

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

So i read these two series seperately using vtkDICOMImageReaderand subtracted these two sets using vtkImageMathematics and stored the result in a variable matand then copied it to vtkImageReslice * reslice; Then I imported this reslice to ITK using VKImageToImageFilter *filter_toitkimage and used this code:

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); 

Then i converted it back to VTK using ImageToVTKImageFilter and used this :

       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);

I am getting the vessels in the lungs But i didn't get the clear image of brain vessels. Can anybody give suggestions?

Answer 1

Your transfer function assumes CT intensities, but the output of vesselness filter is nothing like it. It might be normalized to 0.0-1.0 range, or 0-255 range, or something else. You should examine that.

Answer 2

DSA already means "Digital subtraction angiography". Therefore you shouldn't perform the difference again, it's already done by the software of the scanner. What you are looking at are probably some intermediate steps of the difference, but I unless you would like to work on a new CT reconstruction, I would advise to use the DSA directly and not to do again the difference. Note that in the PACS, it might be referred as modality "XA" for angiography (you shouldn't look for the "CT" modality).

Every system is different, this can be complex but I am also working angiographic data of the brain and this was confusing for me at first, so I thought it worth to say that.

For the prototyping of your vessels detection I would advise to use 3D Slicer in a first step. There is a plugin of a library called VMTK that implements a vesselness filter, but you can also add other plugins, especially one to get all ITK filters. I think it's easier to use that as a starting point, to try out some stuff before implementing.

If you really want to calculate the DSA yourself, in which case you need to have 2 CTs with and without contrast agent (but not DSAs), you have to ensure that both volume are perfectly registered before doing the difference. If you see a misalignment, I would advice to use the simple elastix framework to register your images.

Finally, the critical point of the vesselness filter is actually the frequency that you select with the sigma. This basically defined the vessel size. In the VTMK implementation you can have different frequency so you set a start and a stop sigma. In Slicer you even a plugin to find the right sigma if you input some fiducials around the vessels. So if you can see some vessels but not some others, I would advise to lower this sigma=0.6 to something lower (assuming vessels to show smaller than the lung vessels)

I hope this helps!