肺部CT图像分割及重建系统

悉尼听巴黎叹息

一、 系统概况

我们实现了一个系统，可以从CT图像中将肺部从胸腔中分离出来，并且通过三维重建实现可视化。该系统是基于Visual Studio 2013平台，借助VTK-7.0和Qt5.6开源库通过C++语言实现。
二、 系统设计

肺部CT图像分割及重建系统的实现需要几个方面的工作，一是CT图像的肺部分割；二是CT图像的三维重建；三是可视化界面的设计。
根据工作内容的不同，需要用到不同的开源库。我们通过VTK实现了CT图像的三维重建，包含体绘制和面绘制两种重建算法；通过Qt实现了UI的设计。因此在windows10环境下用VS2013编译支持Qt的VTK库，是前期的准备工作。
在肺部分割这一部分工作上，我们采用以阈值分割为主，其余简单图像处理操作如滤波操作等为辅，进行肺部的分割。
三、系统实现（含源代码）

3.1 前期准备工作——Windows10环境下用VS2013编译支持Qt的VTK库

网上有大量的相关文章，但由于VTK，VS，Qt的版本不同，可能导致Qt与VTK的联动失败。在此，我们参考了两位博主的搭建方式：
(1) http://www.cnblogs.com/tianhu9102/p/7641397.html
(2) https://blog.csdn.net/u011017966/article/details/40984473
这两篇文章对应两种不同的搭建方式，一种是在两种库编译成功后通过直接添加库目录，包含目录的方法实现VTK与Qt联动；另一种则跳过了添加库目录，包含目录等繁琐操作，程序全部用cmake来管理，直接通过cmake进行编译实现一个工程文件。
经过验证，两种方法均可实现VTK + QT的联动。
3.2三维重建的实现

我们采用了两种三维重建方法：体绘制和面绘制。
面绘制是将感兴趣的部分以等值面的形式抽取出来便于利用真实感技术，通过任意旋转和变换光照效果来生成高质量的三维图像，并可以方便的对其进行观察和分析。
面绘制的代码如下：

1. void MainWindow::open(){
2.        
3.        
4.         //抽取等值面为骨头的信息  
5.         vtkSmartPointer< vtkMarchingCubes > boneExtractor =
6.                 vtkSmartPointer< vtkMarchingCubes >::New();
7.         boneExtractor->SetInputConnection(threshould->GetOutputPort());
8.         //boneExtractor->SetValue(0, 500); //设置提取的等值信息
9.         boneExtractor->SetValue(0, 1000);
10.         boneExtractor->SetNumberOfContours(1);
11.         boneExtractor->Update();
13.         //剔除旧的或废除的数据单元，提高绘制速度(可略去这一步)  
14.         vtkSmartPointer< vtkStripper > boneStripper =
15.                 vtkSmartPointer< vtkStripper >::New(); //三角带连接  
16.         boneStripper->SetInputConnection(boneExtractor->GetOutputPort());
17.         boneStripper->Update();
19.         //建立映射  
20.         vtkSmartPointer< vtkPolyDataMapper > boneMapper =
21.                 vtkSmartPointer< vtkPolyDataMapper >::New();
22.         boneMapper->SetInputData(boneStripper->GetOutput());
23.         boneMapper->ScalarVisibilityOff();
25.         //建立角色  
26.         vtkSmartPointer< vtkActor > bone =
27.                 vtkSmartPointer< vtkActor >::New();
28.         bone->SetMapper(boneMapper);
29.         bone->GetProperty()->SetDiffuseColor(.1, .94, .52);
30.         bone->GetProperty()->SetSpecular(.3);
31.         bone->GetProperty()->SetSpecularPower(20);
32.         bone->GetProperty()->SetOpacity(1.0);
33.         bone->GetProperty()->SetColor(1,1,1);
34.         bone->GetProperty()->SetRepresentationToSurface();
35.         //定义绘制器  
36.         vtkSmartPointer< vtkRenderer > aRenderer =
37.                 vtkSmartPointer< vtkRenderer >::New();
38.         //定义绘制窗口  
39.         vtkSmartPointer< vtkRenderWindow > renWin = ui->qvtkWidget_Volume->GetRenderWindow();
40.                 //vtkSmartPointer< vtkRenderWindow >::New();
41.         renWin->AddRenderer(aRenderer);
43.         //定义窗口交互器  
44.         vtkSmartPointer< vtkRenderWindowInteractor > iren =
45.                 vtkSmartPointer< vtkRenderWindowInteractor >::New();
46.         iren->SetRenderWindow(renWin);
48.         //创建一个camera  
49.         vtkSmartPointer< vtkCamera > aCamera =
50.                 vtkSmartPointer< vtkCamera >::New();
51.         aCamera->SetViewUp(0, 0, -1);
52.         aCamera->SetPosition(0, 1, 0);
53.         aCamera->SetFocalPoint(0, 0, 0);
55.         aRenderer->AddActor(bone);
56.         aRenderer->SetActiveCamera(aCamera);
57.         aRenderer->ResetCamera();
58.         aCamera->Dolly(1.5);
59.         aRenderer->SetBackground(0, 0, 0);
60.         aRenderer->ResetCameraClippingRange();
62.         vtkInteractorStyleTrackballCamera *style =       //设置交互方式
63.                 vtkInteractorStyleTrackballCamera::New();
64.         iren->SetInteractorStyle(style);
65.         iren->Initialize();
66.         iren->Start();
67. }

复制代码

体绘制的原理和面绘制完全不相同。面绘制需要生成中间图元，而体绘制则是直接在原图上进行绘制，内容需求较面绘制小。每切换一个视角需要重新对所有的像素点进行颜色和透明度计算，需要时间比面绘制长。VTK中基于体绘制实现三维重建，使用的是光线投射法(Ray-casting)。
体绘制的具体代码如下：

1. void MainWindow::volume(){
2.         VTK_MODULE_INIT(vtkRenderingVolumeOpenGL2);
3.         vtkSmartPointer<vtkFixedPointVolumeRayCastMapper> volumeMapper =
4.                 vtkSmartPointer<vtkFixedPointVolumeRayCastMapper>::New();
5.         volumeMapper->SetInputConnection(threshould->GetOutputPort());
6.         volumeMapper->SetBlendModeToComposite();
7.         volumeMapper->SetSampleDistance(0.3);
8.         volumeMapper->AutoAdjustSampleDistancesOff();
10.         vtkSmartPointer<vtkVolumeProperty> volumeProperty =
11.                 vtkSmartPointer<vtkVolumeProperty>::New();
12.         volumeProperty->SetInterpolationTypeToLinear();
13.         volumeProperty->ShadeOn();  //打开或者关闭阴影测试  
14.         volumeProperty->SetAmbient(0.2);
15.         volumeProperty->SetDiffuse(1.2);  //漫反射  
16.         volumeProperty->SetSpecular(0.1); //镜面反射  
17.         volumeProperty->SetSpecularPower(10);
18.         //设置不透明度  
19.         vtkSmartPointer<vtkPiecewiseFunction> compositeOpacity =
20.                 vtkSmartPointer<vtkPiecewiseFunction>::New();
21.         compositeOpacity->AddPoint(70, 0.00);
22.         compositeOpacity->AddPoint(90, 0.40);
23.         compositeOpacity->AddPoint(180, 0.60);
25.         volumeProperty->SetScalarOpacity(compositeOpacity); //设置不透明度传输函数  
26.         //设置梯度不透明属性  
28.         vtkSmartPointer<vtkPiecewiseFunction> volumeGradientOpacity =
29.                 vtkSmartPointer<vtkPiecewiseFunction>::New();
30.         volumeGradientOpacity->AddPoint(10, 0.0);
31.         volumeGradientOpacity->AddPoint(90, 0.5);
32.         volumeGradientOpacity->AddPoint(100, 1.0);
33.         volumeProperty->SetGradientOpacity(volumeGradientOpacity);//设置梯度不透明度效果对比  
34.         //设置颜色属性  
35.         vtkSmartPointer<vtkColorTransferFunction> color =
36.                 vtkSmartPointer<vtkColorTransferFunction>::New();
37.         color->AddRGBPoint(64.00, 1.00, 0.52, 0.30);
38.         volumeProperty->SetColor(color);
39.         /********************************************************************************/
40.         vtkSmartPointer<vtkVolume> volume =
41.                 vtkSmartPointer<vtkVolume>::New();
42.         volume->SetMapper(volumeMapper);
43.         volume->SetProperty(volumeProperty);
45.         vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
46.         ren->SetBackground(0, 0, 0);
47.         ren->AddVolume(volume);
49.         vtkSmartPointer<vtkRenderWindow> rw =ui->qvtkWidget_Volume_6->GetRenderWindow();
50.         rw->AddRenderer(ren);
51.         rw->Render();
53.         vtkSmartPointer<vtkRenderWindowInteractor> rwi =
54.                 vtkSmartPointer<vtkRenderWindowInteractor>::New();
55.         rwi->SetRenderWindow(rw);
57.         vtkInteractorStyleTrackballCamera *st =       //设置交互方式
58.                 vtkInteractorStyleTrackballCamera::New();
59.         rwi->SetInteractorStyle(st);
60.         /********************************************************************************/
61.         ren->ResetCamera();
62.         rw->Render();
63. }

复制代码

3.3 UI设计

Qt的UI界面设计如下图所示：

                               
登录/注册后可看大图

我们主要用到了Qt的几个常用Widget:
PushButton：用于启动某个功能
horizontalScrollBar: 用于手动调整上界阈值和下界阈值
textEdit：用于输入图像分割的矩形框大小（后面分割会提到）
textBrowser: 用于显示程序运行状态
3.4 分割操作的实现

为了尽可能提高分割的效果，我们采用了多管齐下的方法。
首先，我们设置了以交互的方式，让用户输入上界阈值和下界阈值，则CT图像的CT值在这个阈值内的保留原值不变，小于下界阈值的变黑，大于上界阈值的变白。

1. void MainWindow::segmentBetween(){
2.         int upper, lower;
3.        
4.     // 从滑条里分别获取上界和下界阈值
5.         upper = ui->horizontalScrollBar->value();
6.         lower = ui->horizontalScrollBar_2->value();
8.         threshould->ThresholdBetween(lower, upper);// 利用ThresholdBetween进行分割
9.         threshould->SetInValue(1024);
10.         threshould->SetOutValue(-1024);
11.         threshould->Update(); //算法执行后必须添加的更新消息
12.         slice();
13. }

复制代码

其次，因为CT原图像中有CT机床等干扰，我们在原CT图像上进行截取，只需输入矩形框的两个对角顶点的坐标，就可以轻松截取框内的图像，框外的图像全部置黑。
最后，我们发现通过简单的阈值分割并不能很好地去除一些细节上的噪声和胸腔组织，因此我们写了一个简单的滤波算法，进行滤波操作，根据该像素点的邻域是否具有像素值进行判断，从而极大程度上去除了噪声，取得了不错的分割效果。

1. void MainWindow::Filter(){
2.         int dims[3];
3.         int upperNum, upperNum2, upperNum3, upperNum4;
4.         QString upper, upper2, upper3, upper4;
5.         upper = ui->textEdit->document()->toPlainText();
6.         upperNum = upper.toInt();
7.         upper2 = ui->textEdit_2->document()->toPlainText();
8.         upperNum2 = upper2.toInt();
9.         upper3 = ui->textEdit_3->document()->toPlainText();
10.         upperNum3 = upper3.toInt();
11.         upper4 = ui->textEdit_4->document()->toPlainText();
12.         upperNum4 = upper4.toInt();
13.         reader->GetOutput()->GetDimensions(dims);
14.        
15.         for (int k = 0; k < dims[2] ; k++)
16.         {
17.                 for (int j = 0; j < dims[1] ; j++)
18.                 {
19.                         for (int i = 0; i < dims[0] ; i++)
20.                         {
21.                                
22.                                 short *pixel = (short *)(threshould->GetOutput()->GetScalarPointer(i, j, k));
23.                                 if (j < upperNum2 || j > upperNum4)
24.                                         *pixel = -1024;
25.                                 if (i < upperNum || i > upperNum3)
26.                                         *pixel = -1024;
27.                                
28.                                 if ((j > 5) && (j < dims[1] - 6))
29.                                 {
30.                                         short *pixel2 = (short *)(threshould->GetOutput()->GetScalarPointer(i, j + 5, k));
31.                                         short *pixel3 = (short *)(threshould->GetOutput()->GetScalarPointer(i, j - 5, k));
32.                                         if ((*pixel2 < -500) && (*pixel3  <-500))
33.                                                 *pixel = -1024;
34.                                 }
35.                         }
37.                 }
38.         }
39.         threshould->Update();
40.         slice();
41. }

复制代码

四、结果展示

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

分割：手动设置上界阈值和下界阈值，并在去噪输入框内输入矩形框的对角顶点坐标，则可达到分割的效果。
显示：左上为三维重建的体绘制（左）和面绘制（右），可以通过鼠标拖动，从不同视角观看肺部结构。右边为水平面，矢状面，冠状面的显示，通过拖动进度条可以看不同层的截面显示。
资源已经上传：
https://download.csdn.net/download/h_swhite/10901427

来源：https://blog.csdn.net/H_SWhite/article/details/86006886
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！