2021年认证杯SPSSPRO杯数学建模A题(第一阶段)医学图像的配准全过程文档及程序

2021年认证杯SPSSPRO杯数学建模

A题 医学图像的配准

原题再现：

  图像的配准是图像处理领域中的一个典型问题和技术难点，其目的在于比较或融合同一对象在不同条件下获取的图像。例如为了更好地综合多种信息来辨识不同组织或病变，医生可能使用多种仪器对患者的同一部位进行成像。在综合多幅图像时，首先需要将它们严格对齐，使得图上同一个坐标的位置对应的是真实对象的同一个点，这个过程称之为配准。现在的许多医学成像技术，包括 CT、MRI、PET 等，最终生成的是人体的断层影像。在这里，我们主要关心的是断层成像的配准问题。
  我们考虑对一个患者的腹部进行断层成像。由于人体组织是柔软的，所以即使使用同一台成像设备，两次成像的结果也并不完全一致。最终输出时还会对图像进行自动放缩，所以输出图片的大小也并不完全相同。想要精确配准，需要将其中一次的成像结果进行某种仿射变换（或非线性变换），以尽可能地匹配另一次的结果（或将两次结果都映射到同一个标准模板中）。求得合适的变换就是图像配准的核心任务。
  第一阶段问题： 对同一个患者进行两次 CT 成像，间隔长达数周乃至更长。两次扫描的断层位置相同，但由于占位性病变的发展，成像的结果在某些区域会有区别。请你设计一个有效的方法，使我们能够对这样两张图像进行配准，以使计算机辅助医疗系统能够通过比较来自动识别出病变的位置，并定量地评估出病变的发展情况。

整体求解过程概述(摘要)

  数字化医学影像技术的高速发展推动了图像配准技术的不断进步，目前刚性组织的医学图像配准技术已经较为成熟，而非刚性部位（如腹部）由于形状和位置会随生理运动而发生变化，因此其图像配准面临巨大的挑战。
  本文是对同一患者腹部的两次 CT 图像进行配准分析，重点考虑占位性病变的影响因素，首先采用 Kullback-Leiber 距离对传统的 Demons 算法模型进行改进，以基于Kullback-Leiber 距离的互信息作为图像配准的精度量度，通过迭代次数以及互信息变化量的阈值来控制配准进程，从而实现收敛速度较快的非刚性医学图像配准；再通过对配准前后图像的相似性精度进行仿真计算，以及对模型性能的评估，验证了配准模型的准确度和有效性。对于配准后的图像，采用 CT 兴趣体积法测定 CT 图像中的病变区域的体积，同时通过单因素 ANOVA 分析法对病变区域的体积和临床严重程度的关系进行定量分析。
  本文最后对于基于腹部（非刚性部位）的改进 Demons 形变模型进行评估，针对模型一些局限性，引入拓扑校正算法，以进一步提高配准的准确性，可以有效避免浮动图像产生拓扑结构无法保持的现象，从而使计算机辅助医疗系统更好地识别出病变的位置。

问题分析：

  通过对题目背景的分析，我们知道由于同一个患者占位性病变的发展，并且腹部本身具有非刚性的一些特征，对图像的配准提出了更高的要求。对此，我们基于 Kullback-Leiber 距离对传统的 Demons 算法模型进行改进，实现收敛速度较快的腹部（非刚性）医学图像配准。在本模型中，我们拟用 Kullback-Leiber 距离定义法定义的互信息作为图像配准的测度，通过迭代次数和互信息变化量的阈值要求控制配准过程的结束。之后我们通过对十个案例的仿真，采用配准前后图像相似性测度（互信息与归一化互信息）变化情况，通过对模型性能进行评估，验证了模型的实用性。最后，在优化模块我们通过将拓扑校正算法引入配准模型，以提高配准的准确性，保证变形后的浮动图像不产生拓扑结构不能保持的现象，以便于计算机辅助医疗系统更好地识别出病变的位置。针对于问题的另一方面，即在医学图像配准后如何定量地评估出病变的发展情况的问题，我们选择采用 CT 兴趣体积法定量测定 CT 图像中病变区域体积，并采用单因素ANOVA 分析比较病变区域体积与临床严重程度关系。

模型假设：

  （1）假设两次扫描使用的仪器都是 CT 机；
  （2）假设两张 CT 图像来自于同一患者；
  （3）假设该患者两次扫描用的是同一台 CT 机；
  （4）假设该 CT 机性能良好，可以正常工作；
  （5）假设两次 CT 扫描成像的断层位置相同；
  （6）假设该患者腹腔内占位性病变的影响占主要因素；
  （7）假设两次扫描时间间隔足够长，使得占位性病变对 CT 图像的影响足够显著。

模型的建立与求解

  在 Demons 形变模型中，有关像素点偏移量的计算思想实际上是来自于光流场模型[3]，在配准的过程中，由于这两幅医学图像都是通过同一台 CT 机扫描得到的，因此二者灰度分布相等，即所对应的灰度点可以近似为在单位时间内参考图像里面的像素点的平面运动。Demons 模型的优点就是可以在一定程度上有效、精准地校正图像的非刚性变形，因此该模型就是基于 Kullback-Leibler 距离对 Demons 多模态图像配准算法[4]进行改进。

基于 Kullback-Leibler 距离的互信息梯度的计算模型

  在多模态配准的过程中，互信息常被用于医学图像配准的测度，而互信息越大则可以说明图像匹配效果越好。当将两图像进行配准时，所包含的信息总量达到最大，与此同时，互信息也会达到最大。而互信息梯度（MIG）代表的含义是互信息变大的方向，如果我们在 Demons 形变力的基础上加上 MIG 这个分量，那么不难发现 Demons 形变图像中的各个像素点都会向着互信息变大的方向发生形变，从而实现多模态医学图像的配准。
  在概率论和信息论中，互信息是两个随机变量间相互依赖性的量度，即统计相关性，一般可以用熵来表示，其物理含义是代表系统的复杂性和不确定性。假设该患者的两幅CT 图像分别是 M 和 N，两图像中任意像素点的灰度值分别用𝛼和𝛽来表示，因此根据互信息的定义，我们可以得到：
                MI(M,N)=H(N)-H(N|M)， （1）

基于 Kullback-Leiber 距离的改进 Demons 非刚性图像配准模型

  由于两个医学图像来源于同一患者，且用同一 CT 机扫描得到，因此二者灰度分布相等，当两幅图像的空间位置实现完美匹配之后，就可以做到两图像间的互信息最大。因此以下模型就是基于 Kullback-Leiber 距离的改进 Demons 非刚性图像配准模型，流程图如图 1 所示。

论文缩略图：

程序代码：

1. loadMnistDataScript; %加载数据
2. ntrain = size(training_data_label,2);
3. mini_batch_size = 100;
4. cnn.ntrain = ntrain;
5. cnn.eta = 1; %学习速率
6. cnn.lambda = 5; %正则化惩罚系数
8. cnn.layer = {
9. % input layer: 'input', mini_size, [height,width] of image
10. {'input',mini_batch_size,[28,28]};
11. % convlution layer: 'conv', kernel_number, [height,width] of kernel
12. {'conv',20,[9,9]};
13. % pooling layer: 'pool', pooling_type, [height,width] of pooling area
14. {'pool','average',[2,2]};
15. % flatten layer: 'flat', a layer for pre-allocated memory
16. {'flat'};
17. % full connect layer: 'full', neuron number
18. {'full',100};
19. {'full',100};
20. % output layer: 'output', neuron number
21. {'output',10};
22. };
24. function cnn = cnn_initialize(cnn)
25. %CNN_INIT initialize the weights and biases, and other parameters
26. %
27. index = 0;
28. num_layer = numel(cnn.layer);
29. for in = 1:num_layer
30. switch cnn.layer{in}{1}
31. case 'input'
32. index = index + 1;
33. height = cnn.layer{in}{3}(1);
34. width = cnn.layer{in}{3}(2);
35. mini_size = cnn.layer{in}{2};
36. cnn.weights{index} = [];
37. cnn.biases{index} = [];
38. cnn.nabla_w{index} = [];
39. cnn.nabla_b{index} = [];
40. %n*n*m
41. cnn.a{index} = [];
42. cnn.z{index} = [];
43. cnn.delta{index} = [];
44. cnn.mini_size = mini_size;
45. case 'conv'
46. index = index + 1;
47. %kernel height, width, number
48. ker_height = cnn.layer{in}{3}(1);
49. ker_width = cnn.layer{in}{3}(2);
50. ker_num = cnn.layer{in}{2};
51. cnn.weights{index} = grand(ker_height,ker_width,ker_num) - 0.5;
52. cnn.biases{index} = grand(1,ker_num) - 0.5;
53. cnn.nabla_w{index} = zeros(ker_height,ker_width,ker_num);
54. cnn.nabla_b{index} = zeros(1,ker_num);
55. height = height - ker_height + 1;
56. width = width - ker_width + 1;
57. cnn.a{index} = zeros(height,width,mini_size,ker_num);
58. cnn.z{index} = zeros(height,width,mini_size,ker_num);
59. cnn.delta{index} = zeros(height,width,mini_size,ker_num);
60. case 'pool'
61. index = index + 1;
62. %kernel height, width, number
63. ker_height = cnn.layer{in}{3}(1);
64. ker_width = cnn.layer{in}{3}(2);
65. cnn.weights{index} = [];
66. cnn.biases{index} = [];
67. cnn.nabla_w{index} = [];
68. cnn.nabla_b{index} = [];
69. height = height / ker_height;
70. width = width / ker_width;
71. cnn.a{index} = zeros(height,width,mini_size,ker_num);
72. cnn.z{index} = [];
73. cnn.delta{index} = zeros(height,width,mini_size,ker_num);
74. case 'flat'
75. index = index + 1;
76. cnn.weights{index} = [];
77. cnn.biases{index} = [];
78. cnn.nabla_w{index} = [];
79. cnn.nabla_b{index} = [];
80.
81. cnn.a{index} = zeros(height*width*ker_num,mini_size);
82. cnn.z{index} = [];
83. cnn.delta{index} = zeros(height*width*ker_num,mini_size);
84. case 'full'
85. index = index + 1;
86. %kernel height, width, number
87. neuron_num = cnn.layer{in}{2};
88. neuron_num0 = size(cnn.a{in-1},1);
89.
90. cnn.weights{index} = grand(neuron_num,neuron_num0) - 0.5;
91. cnn.biases{index} = grand(neuron_num,1) - 0.5;
92. cnn.nabla_w{index} = zeros(neuron_num,neuron_num0);
93. cnn.nabla_b{index} = zeros(neuron_num,1);
95. cnn.a{index} = zeros(neuron_num,mini_size);
96. cnn.z{index} = zeros(neuron_num,mini_size);
97. cnn.delta{index} = zeros(neuron_num,mini_size);
98.
99. case 'output'
100. index = index + 1;
101. %kernel height, width, number
102. neuron_num = cnn.layer{in}{2};
103. neuron_num0 = size(cnn.a{in-1},1);
104.
105. cnn.weights{index} = grand(neuron_num,neuron_num0) - 0.5;
106. cnn.biases{index} = grand(neuron_num,1);
107. cnn.nabla_w{index} = zeros(neuron_num,neuron_num0);
108. cnn.nabla_b{index} = zeros(neuron_num,1);
109.
110. cnn.a{index} = zeros(neuron_num,mini_size);
111. cnn.z{index} = zeros(neuron_num,mini_size);
112. cnn.delta{index} = zeros(neuron_num,mini_size);
113. otherwise
114.
115. end
116. end
117. end
118. function cnn = cnn_feedforward(cnn,x)
119. %CNN_FEEDFORWARD CNN feedforward
120. %
121. num = numel(cnn.layer);
122. for in = 1:num
123.
124. switch cnn.layer{in}{1}
125. case 'input'
126. cnn.a{in} = x;
127. case 'conv'
128. kernel_num = cnn.layer{in}{2};
129. for ik = 1:kernel_num
130. cnn.z{in}(:,:,:,ik) = convn(cnn.a{in-1},...
131. cnn.weights{in}(:,:,ik),'valid')+cnn.biases{in}(ik);
132. end
133. cnn.a{in} = relu(cnn.z{in});
134.
135. case 'pool'
136.
137. ker_h = cnn.layer{in}{3}(1);
138. ker_w = cnn.layer{in}{3}(2);
139. kernel = ones(ker_h,ker_w)/ker_h/ker_w;
140.
141. tmp = convn(cnn.a{in-1},kernel,'valid');
142. cnn.a{in} = tmp(1:ker_h:end,1:ker_w:end,:,:);
143.
144. case 'flat'
145. [height,width,mini_size,kernel_num] = size(cnn.a{in-1});
146. for ik = 1:mini_size
147. cnn.a{in}(:,ik) = reshape(cnn.a{in1}(:,:,ik,:),[height*width*kernel_num,1]);
148. end
149. case 'full'
150. cnn.z{in}= bsxfun(@plus,cnn.weights{in}*cnn.a{in-1},cnn.biases{in});
151. cnn.a{in} = sigmoid(cnn.z{in});
152. case 'output'
153. cnn.z{in}= bsxfun(@plus,cnn.weights{in}*cnn.a{in-1},cnn.biases{in});
154. cnn.a{in} = softmax(cnn.z{in});
155. end
156. end
157. end
158. function cnn = cnn_backpropagation(cnn,y)
159. %CNN_BP CNN backpropagation
160.
161. num = numel(cnn.layer);
162.
163. for in = num:-1:2
164.
165. switch cnn.layer{in}{1}
166. case 'conv'
167.
168. ker_h = cnn.layer{in+1}{3}(1);
169. ker_w = cnn.layer{in+1}{3}(2);
170. kernel = ones(ker_h,ker_w)/ker_h/ker_w;
171. [mini_size,kernel_num] = size(cnn.delta{in+1});
172. cnn.nabla_w{in}(:) = 0;
173. cnn.nabla_b{in}(:) = 0;
174. for ik = 1:kernel_num
175. for im = 1:mini_size
176. cnn.delta{in}(:,:,im,ik) = kron(cnn.delta{in+1}(:,:,im,ik),kernel).
*relu_prime(cnn.z{in}(:,:,im,ik));
177. cnn.nabla_w{in}(:,:,ik) = cnn.nabla_w{in}(:,:,ik) +...
178. conv2(rot90(cnn.a{in1}(:,:,im),2),cnn.delta{in}(:,:,im,ik),'valid');
179. cnn.nabla_b{in}(ik) = cnn.nabla_b{in}(ik) + mean(mean(cnn.delta{in}
(:,:,im,ik)));
180. end
181. cnn.nabla_w{in}(:,:,ik) = cnn.nabla_w{in}(:,:,ik)/mini_size;
182. cnn.nabla_b{in}(ik) = cnn.nabla_b{in}(ik)/mini_size;
183. end
184. case 'pool'
185. [height,width,mini_size,kernel_num] = size(cnn.a{in});
186. for ik = 1:mini_size
187. cnn.delta{in}(:,:,ik,:) = reshape(cnn.delta{in+1}(:,ik),[height,width,k
ernel_num]);
188. end
189. case 'flat'
190. cnn.delta{in} = cnn.weights{in+1}'*cnn.delta{in+1};
191. case 'full'
192. cnn.delta{in}= cnn.weights{in+1}'*cnn.delta{in+1}.*sigmoid_prime(cnn.z{in})
;
193. cnn.nabla_w{in} = cnn.delta{in}*(cnn.a{in-1})'/cnn.mini_size;
194. cnn.nabla_b{in} = mean(cnn.delta{in},2);
195. case 'output'
196. cnn.delta{in}= (cnn.a{in} - y);
197. cnn.nabla_w{in} = cnn.delta{in}*(cnn.a{in-1})'/cnn.mini_size;
198. cnn.nabla_b{in} = mean(cnn.delta{in},2);
199. otherwise
200.
201. end
202.
203. end
204.
205. eta = cnn.eta;
206. lambda = cnn.lambda;
207. ntrain = cnn.ntrain;
208. % update models
209. for in = 1:num
210. cnn.weights{in} = (1-
eta*lambda/ntrain)*cnn.weights{in} - eta*cnn.nabla_w{in};
211. cnn.biases{in} = (1-eta*lambda/ntrain)*cnn.biases{in} - eta*cnn.nabla_b{in};
212. end
213.
214. end
215. cnn = cnn_initialize(cnn);
216. max_iter = 50000;
217. for in = 1:max_iter
218. pos = randi(ntrain-mini_batch_size);
219. x = training_data(:,:,pos+1:pos+mini_batch_size);
220. y = training_data_label(:,pos+1:pos+mini_batch_size);
221. cnn = cnn_feedforward(cnn,x);
222. cnn = cnn_backpropagation(cnn,y);
223. if mod(in,100) == 0
224. disp(in);
225. end
226. if mod(in,5000) == 0
227. disp(['validtion accuracy: ',num2str(...
228. cnn_evaluate(cnn,validation_data,validation_data_label)*100), '%']);
229. end
230. end
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可