VTK笔记——医学图像配准（Landmark）

发表于 更新于 分类于 Disqus：
本文字数： 6.6k 阅读时长 ≈ 6 分钟

随着现代医学的迅速发展，医学影像学为临床诊断提供了多种模态的医学图像，如X射线，CT，MRI等，当然，他们各自有各自的诊断优势。不过，有时候，医生希望将不同图像信息进行适当的集成。然而不同模态的医学图像成像原理不同，风辨率，成像参数却不相同，因此在图像融合前必须进行图像匹配。

医学图像配准就是通过寻找一种（或一系列）的空间变换（旋转、偏移、缩放和变形等），使两幅图像的对应点达到空间位置和解剖结构上的完全一致。配准结果应使两幅图像上是所有的解剖点（精确匹配），或至少是所有具有诊断意义的点都达到匹配（粗略匹配）。

其中，基于标记点对应关系的配准方法已经得到了广泛的应用，即所谓的特征匹配，通常用于图像分割、特征提取和关键点（landmark）搜寻等场合。

vtkLandmarkTransform

vtkLandmarkTransform是一种比较经典的匹配算法，基于标记点，两个点集在配准后的平均距离最小，要求输入两个点数必须相等，序号一致的点集，做线性变换。它常用于粗略匹配，效率高。
下面简单的示范一下它是如何使用的？

1.源标记点集

1
2
3
4
5
6
7
8
vtkSmartPointer<vtkPoints> sourcePoints =
  vtkSmartPointer<vtkPoints>::New();
double sourcePoint1[3] = {1.0, 0.0, 0.0};
sourcePoints->InsertNextPoint(sourcePoint1);
double sourcePoint2[3] = {0.0, 1.0, 0.0};
sourcePoints->InsertNextPoint(sourcePoint2);
double sourcePoint3[3] = {0.0, 0.0, 1.0};
sourcePoints->InsertNextPoint(sourcePoint3);

2.目标标记点集

1
2
3
4
5
6
7
8
vtkSmartPointer<vtkPoints> targetPoints =
  vtkSmartPointer<vtkPoints>::New();
double targetPoint1[3] = {0.0, 0.0, 1.1};
targetPoints->InsertNextPoint(targetPoint1);
double targetPoint2[3] = {0.0, 1.02, 0.0};
targetPoints->InsertNextPoint(targetPoint2);
double targetPoint3[3] = {-1.11, 0.0, 0.0};
targetPoints->InsertNextPoint(targetPoint3);

3.设置Landmark

1
2
3
4
5
6
vtkSmartPointer<vtkLandmarkTransform> landmarkTransform = 
  vtkSmartPointer<vtkLandmarkTransform>::New();
landmarkTransform->SetSourceLandmarks(sourcePoints);
landmarkTransform->SetTargetLandmarks(targetPoints);
landmarkTransform->SetModeToRigidBody();
landmarkTransform->Update();

4.获取线性变换

1
vtkMatrix4x4* mat = landmarkTransform->GetMatrix();

效果

在这里插入图片描述

LandmarkTransform.cxx

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#include <vtkPoints.h>
#include <vtkSmartPointer.h>
#include <vtkLandmarkTransform.h>
#include <vtkMatrix4x4.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkProperty.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkVertexGlyphFilter.h>

int main(int, char *[])
{
  /*
  This demo creates a coordinate frame (+x, +y, +z) of vectors and a rotated,
  peturbed frame (+z, +y, -x) and aligns the rotated frame to the original as best as possible.
  */
  
  vtkSmartPointer<vtkPoints> sourcePoints =
    vtkSmartPointer<vtkPoints>::New();
  double sourcePoint1[3] = {1.0, 0.0, 0.0};
  sourcePoints->InsertNextPoint(sourcePoint1);
  double sourcePoint2[3] = {0.0, 1.0, 0.0};
  sourcePoints->InsertNextPoint(sourcePoint2);
  double sourcePoint3[3] = {0.0, 0.0, 1.0};
  sourcePoints->InsertNextPoint(sourcePoint3);
  
  vtkSmartPointer<vtkPoints> targetPoints =
    vtkSmartPointer<vtkPoints>::New();
  double targetPoint1[3] = {0.0, 0.0, 1.1};
  targetPoints->InsertNextPoint(targetPoint1);
  double targetPoint2[3] = {0.0, 1.02, 0.0};
  targetPoints->InsertNextPoint(targetPoint2);
  double targetPoint3[3] = {-1.11, 0.0, 0.0};
  targetPoints->InsertNextPoint(targetPoint3);
  
  // Setup the transform
  vtkSmartPointer<vtkLandmarkTransform> landmarkTransform = 
    vtkSmartPointer<vtkLandmarkTransform>::New();
  landmarkTransform->SetSourceLandmarks(sourcePoints);
  landmarkTransform->SetTargetLandmarks(targetPoints);
  landmarkTransform->SetModeToRigidBody();
  landmarkTransform->Update(); //should this be here?

  vtkSmartPointer<vtkPolyData> source =
    vtkSmartPointer<vtkPolyData>::New();
  source->SetPoints(sourcePoints);
  
  vtkSmartPointer<vtkPolyData> target =
    vtkSmartPointer<vtkPolyData>::New();
  target->SetPoints(targetPoints);
  
  vtkSmartPointer<vtkVertexGlyphFilter> sourceGlyphFilter =
    vtkSmartPointer<vtkVertexGlyphFilter>::New();
  sourceGlyphFilter->SetInputData(source);
  sourceGlyphFilter->Update();
  
  vtkSmartPointer<vtkVertexGlyphFilter> targetGlyphFilter =
    vtkSmartPointer<vtkVertexGlyphFilter>::New();
  targetGlyphFilter->SetInputData(target);
  targetGlyphFilter->Update();
    
  vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter =
    vtkSmartPointer<vtkTransformPolyDataFilter>::New();
  transformFilter->SetInputConnection(sourceGlyphFilter->GetOutputPort());
  transformFilter->SetTransform(landmarkTransform);
  transformFilter->Update();
    
  // Display the transformation matrix that was computed
  vtkMatrix4x4* mat = landmarkTransform->GetMatrix();
  std::cout << "Matrix: " << *mat << std::endl;

  // Visualize
  vtkSmartPointer<vtkPolyDataMapper> sourceMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  sourceMapper->SetInputConnection(sourceGlyphFilter->GetOutputPort());

  vtkSmartPointer<vtkActor> sourceActor =
    vtkSmartPointer<vtkActor>::New();
  sourceActor->SetMapper(sourceMapper);
  sourceActor->GetProperty()->SetColor(0,1,0);
  sourceActor->GetProperty()->SetPointSize(4);

  vtkSmartPointer<vtkPolyDataMapper> targetMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  targetMapper->SetInputConnection(targetGlyphFilter->GetOutputPort());

  vtkSmartPointer<vtkActor> targetActor =
    vtkSmartPointer<vtkActor>::New();
  targetActor->SetMapper(targetMapper);
  targetActor->GetProperty()->SetColor(1,0,0);
  targetActor->GetProperty()->SetPointSize(4);
  
  vtkSmartPointer<vtkPolyDataMapper> solutionMapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();
  solutionMapper->SetInputConnection(transformFilter->GetOutputPort());

  vtkSmartPointer<vtkActor> solutionActor =
    vtkSmartPointer<vtkActor>::New();
  solutionActor->SetMapper(solutionMapper);
  solutionActor->GetProperty()->SetColor(0,0,1);
  solutionActor->GetProperty()->SetPointSize(3);
  
  // Create a renderer, render window, and interactor
  vtkSmartPointer<vtkRenderer> renderer =
    vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renderWindow =
    vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);
  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  renderWindowInteractor->SetRenderWindow(renderWindow);

  // Add the actor to the scene
  renderer->AddActor(sourceActor);
  renderer->AddActor(targetActor);
  renderer->AddActor(solutionActor);
  renderer->SetBackground(.3, .6, .3); // Background color green

  // Render and interact
  renderWindow->Render();
  renderWindowInteractor->Start();
  
  return EXIT_SUCCESS;
}

Ref

VTKExamples/Cxx/Filtering/LandmarkTransform
在这里插入图片描述