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

news2024/11/25 16:28:49

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

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/13783.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

5年自动化测试,终于进字节跳动了,年薪30w其实也并非触不可及

一些碎碎念 什么都做了,和什么都没做其实是一样的,走出“瞎忙活”的安乐窝,才是避开弯路的最佳路径。希望我的经历能帮助到有需要的朋友。 在测试行业已经混了5个年头了,以前经常听到开发对我说,天天的点点点有意思没…

java计算机毕业设计springboot+vue+elementUI永加乡精准扶贫信息管理系统

项目介绍 系统设计的主要意义在于,一方面,对于网站来讲,系统上线后可以带来很大的便利性,精准扶贫网站管理属于非常细致的管理模式,要求数据量大,计算机管理可以提高精确性,更为便利的就是信息…

NF-κB 信号通路调节细胞因子转录

NF-κB 大家族哺乳动物 NF-κB 家族由五种成员组成:RelA/p65、c-Rel、RelB、p50 (NF-κB1) 和 p52 (NF-κB2),它们可以形成各种异源二聚体或者同源二聚体 (如常见 p50/RelA 异源二聚体),并通过与启动子的 κB 位点结合来激活大量基因。所有 N…

Mysql常用函数

Mysql常用函数 字段拼接(concat) CONCAT() 函数用于将多个字符串连接成一个字符串 格式: select CONCAT(str1,str2,…) from table_name; #查询商品表,返回一列:商品名称(价格)。 SELECT concat(prod_name,(,prod…

【论文阅读】Weakly Supervised Semantic Segmentation using Out-of-Distribution Data

一篇弱监督分割领域的论文,发表在CVPR2022上: 论文标题: Weakly Supervised Semantic Segmentation using Out-of-Distribution Data 作者信息: 代码地址: https://github.com/naver-ai/w-ood Abstract 作者认为…

专精特新小巨人的申报条件

专精特新企业分为市级专精特新、省级专精特新和国家级专精特新。 在2018年,开展了国家第一批专精特新“小巨人” 企业申报工作。为了引导中小企业积极走“专精特新”发展之路,加快新旧动能转 换步伐,提升自主创新能力、加快转型升级&#xf…

软考的网络工程师对就业有用吗?

考证只是一个结果,学会技能才是最重要的。 视工作而言,软考中级网络工程师的性价比还是非常高的,对于从事同类的技术人员,基础扎实一般可以裸考通过。 含金量嘛,计算机专业可以以考代凭,毕竟证书是人社部和…

安装谷歌服务框架2022最新版本22.45.15失败

在这里(谷歌play服务框架下载安装安卓版-谷歌服务框架2022最新版本(Google Play 服务)下载22.45.15官方手机版-蜻蜓手游网 (qt6.com)http://www.qt6.com/XiaZai/155507.html)下载了谷歌服务框架(Google Play 服务),其应用信息为: 包名:com.g…

Mutated 源代码解析 client (一)

Mutated , a C project https://github.com/scslab/mutated usage Main function in the client directory, mutated_synthetic.cc Line 14 parse the user arguments, such as “-h, -w, -c” parse_synthetic is implemented in client\opts_synthetic.cc Here, use th…

Dive into TensorFlow系列(3)- 揭开Tensor的神秘面纱

TensorFlow计算图是由op和tensor组成,那么tensor一般都用来代表什么呢?显然,像模型的输入数据、网络权重、输入数据经op处理后的输出结果都需要用张量或特殊张量进行表达。既然tensor在TensorFlow体系架构中如此重要,因此本文将带…

Redis通用命令和key的层级结构

目录 1 Redis数据结构介绍 2 Redis 通用命令 3 Redis命令-Key的层级结构 1 Redis数据结构介绍 Redis是一个key-value的数据库,key一般是String类型,不过value的类型多种多样: value的数据类型共有8种,前面5中为基本数据类型&a…

5000立方米球罐设计

目 录 摘 要 I Abstract II 1 文献综述 1 1.1 课题研究的工程背景及理论、实际意义 1 1.2 球罐用钢 1 1.2.1 球罐用钢基本要求分析 1 1.2.2 国内外球罐的常用钢种 2 1.2.3 几种典型球罐用钢的优劣对比 2 1.3 球罐设计 3 1.3.1 球罐设计的执行标准及法规 3 1.3.2 球壳结构 4 1.3…

通过PLC网关如何实现三菱FX3U的远程上下载程序?

FX3U是三菱推出的高性能PLC品牌。基本性能大幅提升,晶体管输出型的基本单元内置了3轴独立最高100kHz的定位功能,并且增加了新的定位指令,从而使得定位控制功能更加强大,使用更为方便,受到企业的青睐。因此,…

PyQt5 QLabel标签

PyQt5 QLabel标签标签显示标签快捷键功能标签显示 QLabel背景色设置: palette QPalette() # 创建调色板 palette.setColor(QPalette.Window, Qt.green) # 设置调色板属性 label.setPalette(palette) # 标签设置Palette label.setAutoFillBackground(True) # 设为T…

【安全测试学习】数据库基础

以上来自学习极客时间《Web 安全攻防实战》课程内容,汇总整理思维导图。

记录多次安装mysql失败问题

首先说明一下,本人电脑已经安装过mysql,不过想从5.7版本升级到8.0 首先是卸载电脑上的mysql5.7版本,卸载过程如下: 进入控制面板,直接卸载所有mysql相关进入安装目录下删除mysql相关文件夹,通常在C:\Prog…

”互联网行业还在等金三银四或是金九银十?“,我劝你还是早做打算

对于今年的IT行业来说,可能真的根本没有所谓的“金三银四”或是“金九银十”。各大招聘网站或者软件上不管是大厂还是中小公司看似挂个招聘需求,但实际上这些公司并不非常缺人也不急着招人;我想今年程序员听的最多的一个词就是”互联网寒冬“…

二十九、Java 数据结构

Java 数据结构 Java工具包提供了强大的数据结构。在Java中的数据结构主要包括以下几种接口和类: 枚举(Enumeration)位集合(BitSet)向量(Vector)栈(Stack)字典&#xff…

谈谈制定数据治理战略路线图的方法

对于商业世界最具前瞻性思维能力的发展来说,如数据分析、机器学习和人工智能,高质量的数据是一个关键的成功因素。因此,当涉及到数字化转型时,数据发挥着至关重要的作用。 然而,如果没有适当的数据治理,组织最终可能会构建腐败的模型,做出低效的决策,甚至违反法律。商…

C++11 thread

目录 线程thread 主要成员函数 简单线程的创建 线程封装 zero_thread.h zero_thread.cpp main.cpp C/CLinux服务器开发/后台架构师【零声教育】-学习视频教程-腾讯课堂 线程thread std::thread 在 #include 头文件中声明,因此使用 std::thread 时需要包含…