CN111739139B - 一种噪声ct图像三维重建方法 - Google Patents
一种噪声ct图像三维重建方法 Download PDFInfo
- Publication number
- CN111739139B CN111739139B CN202010591066.8A CN202010591066A CN111739139B CN 111739139 B CN111739139 B CN 111739139B CN 202010591066 A CN202010591066 A CN 202010591066A CN 111739139 B CN111739139 B CN 111739139B
- Authority
- CN
- China
- Prior art keywords
- image
- voxel
- isosurface
- gray
- output image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 92
- 238000001914 filtration Methods 0.000 claims abstract description 46
- 230000002146 bilateral effect Effects 0.000 claims abstract description 22
- 230000000877 morphologic effect Effects 0.000 claims abstract description 14
- 238000012545 processing Methods 0.000 claims description 10
- 239000002131 composite material Substances 0.000 claims description 9
- 238000010606 normalization Methods 0.000 claims description 7
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000002591 computed tomography Methods 0.000 description 19
- 238000010586 diagram Methods 0.000 description 16
- 230000014509 gene expression Effects 0.000 description 7
- 238000005070 sampling Methods 0.000 description 5
- 230000003044 adaptive effect Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000007797 corrosion Effects 0.000 description 3
- 238000005260 corrosion Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000012805 post-processing Methods 0.000 description 3
- 238000003745 diagnosis Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 238000012800 visualization Methods 0.000 description 2
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 description 1
- 241000256010 Manduca Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000013170 computed tomography imaging Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 210000004185 liver Anatomy 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000008961 swelling Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T15/00—3D [Three Dimensional] image rendering
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T19/00—Manipulating 3D models or images for computer graphics
- G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10004—Still image; Photographic image
- G06T2207/10012—Stereo images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
- G06T2207/20028—Bilateral filtering
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Graphics (AREA)
- Architecture (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种噪声CT图像三维重建方法，包括步骤：S1、利用自适应形态学滤波算法对灰度图像进行滤波，得到第一输出图像；S2、利用双边滤波器对第一输出图像进行滤波，得到第二输出图像；S3、利用基于区域增长法和三线性插值法的移动立方体算法对第二输出图像进行重建，得到三维重建图像。该三维重建方法不仅能够平滑的对图像进行滤波，而且能够很好地保留图像的边缘信息和几何信息，避免了将新的噪声或伪影引入图像中；同时，利用基于区域增长法和三线性插值法的移动立方体算法进行三维重建，区域增长法可以有效加快重建效率，三线性插值法可以有效提高重建精度，复杂度低、耗时短，从而扩大了临床应用领域。
Description
技术领域
本发明属于图像处理技术领域，具体涉及一种噪声CT图像三维重建方法。
背景技术
随着医学数字成像技术的飞速发展及其在临床诊断中的广泛应用，计算机断层(Computed Tomography，CT)图像的三维重建技术在医学辅助诊疗工作中凸显了其优势，取得了巨大的进展。比起二维断层图像来说三维重建后的图像更加立体真实，在医学图像领域它能更加直观清晰的直达病灶部位进行分析，为临床医学实际操作提供有效途径。与此同时，医学界迫切希望CT三维重建的实时性和精确性能进一步提高，以提高医疗诊断的准确性。
医学图像的三维重建和可视化是科学可视化的重要研究领域,在成像和传输过程中，CT图像会受到许多因素的影响，例如散焦模糊的电子辐射检测效果、检测器的数量和间距、帧的几何形状、重构算法等。噪声直接影响CT图像的密度分辨率和空间分辨率，而且使CT图像的三维重建精度大大降低，对临床分辨病灶带来极大干扰。
目前常用的CT去噪方法可分为投影域方法、迭代重建方法和后处理方法三大类。Balda提出的结构自适应滤波和Manduca提出的双边滤波是两种有效的投影域去噪方法，但是该类算法在对投影数据进行降噪处理的过程中易引起数据不一致、过校正或欠校正现象，经重建后可能会将新的噪声或伪影引入图像中。Liu J提出的一种3D字典学习正则化方法和Ding提出的基于压缩感知的方法是两种典型的迭代重建算法，但这类算法的缺点是计算复杂度高、耗时长，在临床应用中有很大的局限性。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种噪声CT图像三维重建方法。本发明要解决的技术问题通过以下技术方案实现：
本发明实施例提供了一种噪声CT图像三维重建方法，包括步骤：
S1、利用自适应形态学滤波算法对灰度图像进行滤波，得到第一输出图像；
S2、利用双边滤波器对所述第一输出图像进行滤波，得到第二输出图像；
S3、利用基于区域增长法和三线性插值法的移动立方体算法对所述第二输出图像进行重建，得到三维重建图像。
在本发明的一个实施例中，步骤S1包括：
S11、获取所述灰度图像；
S12、设定包括形状序列和尺寸序列的结构元素集Anm，其中，Anm＝{A11,A12,A13,…,A21,A22,…,Anm}，n代表形状序列，m代表尺寸序列，
并利用所述结构元素集中的相同尺寸、不同形状的结构元素腐蚀所述灰度图像，根据所述结构元素可填入所述灰度图像的次数自适应确定不同形状的结构元素的权重值：
...
其中，α1,α2,…,αn为不同形状的结构元素权重值，β1,β2,…,βn为不同形状的结构元素可填入灰度图像的次数；
S13、对同一形状的所述结构元素，按尺寸从小到大的顺序将所述结构元素对所述灰度图像形态学处理的过程进行排序，构成串联滤波器；
S14、结合所述权重值，按所述结构元素的不同形状将所述串联滤波器进行并联，构成复合滤波器，并利用所述复合滤波器对所述灰度图像进行滤波，得到所述第一输出图像：
F(x)＝∑αifi(x)
其中，fi(x)为串联滤波器对灰度图像的滤波结果，αi为权重值；i＝1,2,…,n。
在本发明的一个实施例中，步骤S2包括：
S21、获取基于空间距离的高斯权重区域滤波器：
其中，c(ξ,x)为基于空间距离的高斯权重，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，kd(x)为归一化因子，
S22、获取基于像素间相似度的边缘滤波器：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，s(f(ξ),f(x))为基于像素相似度的高斯权重，kr(x)为归一化因子，
S23、将所述高斯权重区域滤波器与所述边缘滤波器结合，得到所述双边滤波器：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，k(x)为归一化因子，
c(ξ,x)为基于空间距离的高斯权重，
且有d(ξ,x)＝d(ξ-x)＝||ξ-x||，σd为基于距离的正态分布标准差，s(f(ξ),f(x))为基于像素相似度的高斯权重，
且有
σr为基于空间的正态分布标准差；
并利用所述双边滤波器对所述第一输出图像进行滤波，得到所述第二输出图像。
在本发明的一个实施例中，步骤S3包括：
S31、根据所述第二输出图像创建多维数据集，并从所述多维数据集中随机选取种子体素；
S32、以所述种子体素的邻域体素为目标区域，利用所述区域增长法提取所述目标区域中包含等值面的所有体素；
S33、根据所述体素中顶点的第一灰度值与所述等值面的第二灰度值的关系，确定所述等值面在每个所述体素中的剖分方式；
S34、采用所述三线性插值法计算在所述剖分方式下所述等值面与每个所述体素的交点位置；
S35、判断是否遍历所有包含所述等值面的体素；若是，则得到所述三维重建图像；若否，则跳转至步骤S31。
在本发明的一个实施例中，步骤S31包括：
扫描所述第二输出图像中相邻的两个二维图像，根据其中一个所述二维图像上的四个相邻像素与另一个所述二维图像上的四个相邻像素创建所述多维数据集，并利用二次取样方法从所述多维数据集中选取包含所述等值面的所述种子体素。
在本发明的一个实施例中，步骤S32包括：
以所述种子体素的邻域体素为目标区域，利用区域增长法判断所述目标区域中是否包含所述等值面，若是，则跳转步骤S33，若否，则跳转步骤S31。
在本发明的一个实施例中，利用区域增长法判断所述目标区域中是否包含所述等值面，包括：
获取所述种子体素与所述邻域体素的公共面的四个顶点的灰度值；
判断所述四个顶点的灰度值与预设阈值的关系；当所述四个顶点的灰度值同时大于或者小于所述预设阈值时，则所述邻域体素中不包含所述等值面；当所述四个顶点的灰度值不是同时大于或者小于所述预设阈值时，则所述邻域体素中含所述等值面。
在本发明的一个实施例中，步骤S33包括：
判断所述第一灰度值和所述第二灰度值的关系，当所述第一灰度值大于或等于所述第二灰度值时，则所述顶点位于所述等值面之外，当所述第一灰度值小于所述第二灰度值时，则所述顶点位于所述等值面之内；
根据所述体素中所有所述顶点与所述等值面的关系确定所述等值面在每个所述体素中的剖分方式。
在本发明的一个实施例中，步骤S34包括：
根据所述剖分方式确定所述等值面与所述体素相交的体素棱边；
利用所述三线性插值法计算所述等值面与所述体素棱边的交点坐标：
Vxyz＝V000(1-x)(1-y)(1-z)+V100x(1-y)(1-z)+V010(1-x)y(1-z)+V001(1-x)(1-y)z+V101x(1-y)z+V011(1-x)yz+V110xy(1-z)+V111xyz
其中，Vxyz为等值面的灰度值，V000,V100,V010,......,V111为体素中各个顶点的灰度值；
由所述交点坐标得到所述等值面与每个所述体素的交点位置。
与现有技术相比，本发明的有益效果：
本发明的三维重建方法先利用自适应形态学滤波算法进行滤波，然后利用双边滤波器进行滤波，将二者进行结合不仅能够平滑的对图像进行滤波，而且能够很好地保留图像的边缘信息和几何信息，避免了将新的噪声或伪影引入图像中；同时，本发明利用基于区域增长法和三线性插值法的移动立方体算法进行三维重建，区域增长发可以有效加快重建效率，三线性插值法可以有效提高重建精度，复杂度低、耗时短，从而扩大了临床应用领域。
附图说明
图1为本发明实施例提供的一种噪声CT图像三维重建方法的流程示意图；
图2为本发明实施例提供的一种串联滤波器的结构示意图；
图3为本发明实施例提供的一种复合滤波器的结构示意图；
图4为本发明实施例提供的一种体素模型的示意图；
图5为本发明实施例提供的一种相邻体素等值面与棱边交点情况的示意图；
图6为本发明实施例提供的一种等值面与体素中一个面的交点情况的示意图；
图7为本发明实施例提供的一种Marching Cubes算法的15种基本剖分方式的示意图；
图8为本发明实施例提供的一种体素顶点和边的索引示意图；
图9为本发明实施例提供的一种单元立方体坐标示意图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
实施例一
本实施例中采用图像后处理方法对重建图像进行处理，采用数字图像处理技术来抑制噪声；图像后处理方法直接对图像进行操作来降低噪声和抑制伪像，具有处理速度快、可移植性较强的特点。
请参见图1，图1为本发明实施例提供的一种噪声CT图像三维重建方法的流程示意图。该噪声CT图像三维重建方法包括步骤：
S1、利用自适应形态学滤波算法对灰度图像进行滤波，得到第一输出图像。
数学形态学是一门建立在格论和拓扑学基础之上的图像分析学科，以形态结构元素为基础对图像进行分析，其基本运算有4个：膨胀、腐蚀、开启和闭合。基于这些基本运算还可以推导和组合成各种数学形态学实用算法。
具体地，步骤S1包括以下步骤：
S11、获取灰度图像。其中，灰度图像为二维图像。
S12、设定包括形状序列和尺寸序列的结构元素集，并利用结构元素集中的相同尺寸、不同形状的结构元素腐蚀灰度图像，根据结构元素可填入灰度图像的次数自适应确定不同形状的结构元素的权重值。
具体地，为了得到好的滤波效果，适当地选择结构元素的尺寸和形状是非常关键的，在复杂图像的处理过程当中，采用多结构元素应包括尺寸和形状两个方面。设结构元素集为Anm，n代表形状序列，m代表尺寸序列，则，Anm＝{A11,A12,A13,…,A21,A22,…,Anm}，其中，
当利用结构元素集Anm中的结构元素腐蚀一幅图像时，腐蚀的过程相当于对可以填入的结构元素的位置作标记的过程。利用相同尺寸、不同形状的结构元素探测一幅图像时，它们可填入的次数是不同的。如果结构元素的选取较好地适合图像(包括边缘)信息，则可填入的次数多，反之则少。利用这一形态学的结构元素“探针”原理，可在腐蚀过程当中，根据结构元素可填入灰度图像的次数的不同自适应确定权重值。
设不同形状的结构元素权重值分别为：α1,α2,…,αn，腐蚀过程当中不同形状的结构元素可填入灰度图像的次数分别为：β1,β2,…,βn，则自适应确定的权重值为：
其中，α1,α2,…,αn为不同形状的结构元素的权重值，β1,β2,…,βn为不同形状的结构元素可填入灰度图像的次数。
S13、对同一形状的结构元素，按尺寸从小到大的顺序将结构元素对灰度图像形态学处理的过程进行排序，构成串联滤波器。
在形态滤波过程中，采用由小到大的结构元素进行处理以达到几乎全部滤除噪声的目的，为了保持更多细节，采用多结构元素以保持更多有用的几何信息。具体地，对于同一形状的结构元素按尺寸从小到大的顺序将结构元素对灰度图像形态学处理的过程进行排序，构成串联滤波器，请参见图2，图2为本发明实施例提供的一种串联滤波器的结构示意图，其中，CO代表开运算，OC代表闭运算，Ai1,Ai2,...,Aim代表同一形状、不同尺寸的结构元素，Input代表输入图像，output代表输出图像。
S14、结合权重值，按结构元素的不同形状将串联滤波器进行并联，构成复合滤波器，并利用复合滤波器对灰度图像进行滤波，得到第一输出图像。
具体地，将不同形状结构元素构成的串联滤波器加以并联，结合自适应确定的权重值，构成串、并复合滤波器，请参见图3，图3为本发明实施例提供的一种复合滤波器的结构示意图，其中，CO代表开运算，OC代表闭运算，Ai1,Ai2,...,Aim、A21,A22,...,A2m、…、An1,An2,...,Anm分别代表同一形状、不同尺寸的结构元素，Ai1,Ai2,...,Aim、A21,A22,...,A2m、…、An1,An2,...,Anm之间的形状不同，Input代表输入图像，output代表输出图像。
根据图3所示，设输入的灰度图像为f(x)，某种形状结构元素串行滤波结果即串联滤波器对灰度图像的滤波结果为fi(x)，i＝1,2,…,n；输出图像为F(x)，结构元素通过公式(1)确定的权重值为α1,α2,…,αn，则第一输出图像为：
F(x)＝∑αifi(x) (2)
其中，fi(x)为串联滤波器对灰度图像的滤波结果，αi为权重值；i＝1,2,…,n。
S2、利用双边滤波器对第一输出图像进行滤波，得到第二输出图像。
双边滤波器是一种非线性去噪算法，旨在消除噪声的同时保留信号细节。一般的高斯滤波在进行采样时主要考虑了像素间的空间距离关系，但是却并没有考虑像素值之间的相似程度，因此这样得到的模糊结果通常是一团模糊的整张图片。双边滤波的优势就在于在采样时不仅考虑像素在空间距离上的关系，同时加入了像素间的相似程度考虑，在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。对上述像素在空间距离上的关系和像素间的相似程度进行结合即可以得到基于空间距离、相似程度综合考量的双边滤波器。
具体地，步骤S2包括：
S21、获取基于空间距离的高斯权重区域滤波器：
其中，c(ξ,x)为基于空间距离的高斯权重，用于度量邻域中心x与附近点ξ之间的几何接近度，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，kd(x)为归一化因子。
kd(x)的表达式为：
S22、获取基于像素间相似度的边缘滤波器：
高斯滤波在低通滤波算法中有不错的表现，但是其却有另外一个问题，那就是只考虑了像素间的空间位置上的关系，因此滤波的结果会丢失边缘的信息。而双边滤波就是在高斯滤波中加入了另外的一个权重部分来解决这一问题。双边滤波中对于边缘的保持通过下述表达式来实现：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，s(f(ξ),f(x))为基于像素相似度的高斯权重，kr(x)为归一化因子。
kr(x)用来对结果进行单位化，其表达式为：
S23、将高斯权重区域滤波器与边缘滤波器结合，得到双边滤波器，并利用双边滤波器对第一输出图像进行滤波，得到第二输出图像。
具体地，将步骤S21基于空间距离的高斯权重区域滤波器和步骤S22基于像素间相似度的边缘滤波器进行结合，即可以得到基于空间距离、相似程度综合考量的双边滤波器：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，k(x)为归一化因子，c(ξ,x)为基于空间距离的高斯权重，s(f(ξ),f(x))为基于像素相似度的高斯权重。
式(7)中，k(x)综合了两种高斯权重于一起而得到，表示单位化分部，用于对结果进行单位化，其表达式为：
基于空间距离的高斯权重c(ξ,x)的计算可以详细描述如下：
且有d(ξ,x)＝d(ξ-x)＝||ξ-x||，σd为基于距离的正态分布标准差。
基于像素相似度的高斯权重s(f(ξ),f(x))的计算可以详细描述如下：
且有
σr为基于空间的正态分布标准差。
利用式(7)中的公式即可对步骤S1中的第一输出图像进行滤波，进而得到第二输出图像。
进一步的，式(7)给出的表达式均是在空间上的无限积分，因而在使用前需要对其进行离散化；而且也不需要对于每个局部像素从整张图像上进行加权操作，距离超过一定程度的像素实际上对当前的目标像素影响很小，可以忽略；因此，限定局部子区域后的离散化公就可以简化为如下形式：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，k(x)为归一化因子，c(ξ,x)为空间距离的高斯权重，s(f(ξ),f(x))为基于像素相似度的高斯权重。
本实施例中，自适应形态学滤波能够很好地保留图像的边缘信息且能够很好地与其他算法结合，双边滤波器不仅能够平滑的对图像进行滤波，而且可以很好地保留图像边缘信息。本实施例结合两者特点，在自适应形态学对灰度图像进行滤波的基础上再利用双边滤波器进行滤波，不仅能够平滑的对图像进行滤波，而且能够很好地保留图像的边缘信息和几何信息，实现了很好的滤波效果，使接下来的三维重建精度得到有效提升。
S3、利用基于区域增长法和三线性插值法的移动立方体算法对第二输出图像进行重建，得到三维重建图像。
传统的移动立方体(Marching Cubes，MC)算法是逐个处理数据场中的立方体(体素)，分离出与等值面相交的立方体，采用线性插值计算出等值面与立方体边的交点。根据立方体每一顶点与等值面的相对位置，将等值面与立方体边的交点按一定方式连接生成等值面，作为等值面在该立方体内的一个逼近表示。MC算法的前提是假设原始数据是离散的三维空间规则数据场，CT医学图像就属于这种数据类型。
本实施例提出了一种基于区域增长法与三线性插值的MC(Marching CubesAlgorithm Based on region growing and trilinear interpolant，RT-MC)算法，该算法对传统的MC算法进行了两个方面的改进：首先，通过选定一个种子体素，处理它的6邻域体素，筛选出与种子体素相邻的、包含等值面的体素，利用区域增长的思想，把所有包含等值面的体素全部筛选出来再继续进行计算，这样有效的加快了重建效率；其次，利用三线性插值算法计算体素与等值面的交点，充分利用体素8个顶点的灰度值计算交点，有效地提高了重建精度。本实施例的RT-MC算法将上述两者进行结合，在重建效率和精度上都取得显著提升。
步骤S3具体包括步骤：
S31、根据第二输出图像创建多维数据集，并从多维数据集中随机选取种子体素。
具体地，对灰度图像进行自适应形态学滤波和双边滤波后，得到第二输出图像的多个切片，每个切片均为一个二维图像，之后将多个二维图像读入内存。接着，扫描两个相邻的二维图像，并用一个二维图像上的四个相邻像素和下一个二维图像上的四个相邻像素创建一个多维数据集；然后从多维数据集中随机选取一个种子体素。
请参见图4，图4为本发明实施例提供的一种体素模型的示意图。移动立方体(Marching Cubes，MC)将表面定位在由八个像素创建的逻辑立方体中，每四个像素点来自相邻的两个切片(即二维图像)。在离散的三维空间数据场中，数据在断层面上排列紧密，而断层之间排列疏松，因此数据分布在由长方体组成的三维网格点上，网格之间的距离可以不等，即在x、y、z方向上不等，但在同一方向上是相等的，这类数据可看作规则的三维网格结构化数据。在这种数据场中，由上下两层相临的八个网格顶点构成的六面体就称体素，如图4所示。图4中，体素的八个顶点坐标为：0(i,j,k)，1(i+1,j,k)，2(i+1,j,k+1)，3(i,j,k+1)，4(i,j+1,k)，5(i+1,j+1,k)，6(i+1,j+1,k+1)，7(i,j+1,k+1)。本实施例中，种子体素是随机选取的一个体素，用于在区域增长法中作为起始体素。
进一步地，为便于后续区域增长法的进行，本实施例中通过二次取样的方法选取种子体素；同时，为在不影响数据精确度的同时加快遍历速度，本实施例找到一个包含等值面的体素作为种子体素。
其中，等值面是空间中所有具有某个相同值的点的集合，它可以表示成：
{(x,y,z)|f(x,y,z)＝c}，c是常数 (12)
等值面是空间中的一张曲面，在该曲面上函数F(x,y,z)的值等于某一给定值。准确地讲，在某一网格空间中，假若每一结点保存着三变量函数F(x,y,z)，而且网格单元在x、y、z方向上的连续采样值为F(x,y,z)，则对于某一给定值Fi，等值面是由所有满足S＝{(x,y,z)|F(x,y,z)＝Fi}的点组成的一张曲面。按照此严格定义下得到的等值面表达式如下：
F(x,y,z)＝a0+a1x+a2y+a3z+a4xy+a5yz+a6zx+a7xyz (13)
其中，a0,a1,...,a7为体素的八个顶点坐标。
S32、以种子体素的邻域体素为目标区域，利用区域增长法提取目标区域中包含等值面的所有体素。
具体地，以种子体素的6个邻域体素形成一个目标区域，利用区域增长法判断该目标区域中是否包含等值面，若是，则跳转步骤S33，若否，则跳转步骤S31。进一步地，利用区域增长法判断该目标区域中是否包含等值面的方法为：判断四个顶点的灰度值与预设阈值的关系；当四个顶点的灰度值同时大于或者小于预设阈值时，则邻域体素中不包含等值面；当四个顶点的灰度值不是同时大于或者小于预设阈值时，则邻域体素中含等值面。
请参见图5，图5为本发明实施例提供的一种相邻体素等值面与棱边交点情况的示意图。本实施例中区域增长的基本思想是将具有相似性质的体素集合起来构成区域。在改进算法时，问题在于怎样判断种子体素周围的体素是否包含等值面，每个体素都含有8个顶点，6个面，在图5中，四个立方体相邻，AB为公共的棱边，下方的两个立方体(体素)左右相邻，有一个面是公共的面即面ABCD，包含4个顶点A、B、C、D，如果这4个顶点的灰度值不是同时大于或者小于预设阈值(预设阈值是指等值面的灰度值)，表示与种子体素相邻的体素中也有等值面。图5内的线段EF为等值面与左下方种子体素的交线，但EF同时又存在于右下方的相邻体素中，因此右边的体素也包含等值面。
因此，要判断种子体素周围的体素是否包含等值面，只需要判断在该方向上的面中4个顶点的状态就可以，对于该面中4个顶点的状态有如下4种情况，请参见图6，图6为本发明实施例提供的一种等值面与体素中一个面的交点情况的示意图，图6中，(a)为有1个顶点与其他状态不同，(b)、(c)为有2个顶点与其他状态不同，(d)为全部顶点具有相同的状态，因此只要4个顶点处于(a)、(b)、(c)三种状态之一，则与种子体素相邻的体素就含有等值面。
通过上述规则可知，种子体素的选择至关重要，因此，本实施例在试验中通过二次取样的方法进行种子体素的选择，同时在不影响数据精确度的同时加快遍历速度，找到一个包含等值面的体素作为种子体素，确定种子体素后就可以监察相邻的体素了。
S33、根据体素中顶点的第一灰度值与等值面的第二灰度值的关系，确定等值面在每个体素中的剖分方式。
本实施例中，首先判断体素中各个顶点的第一灰度值和等值面的第二灰度值的关系，当第一灰度值大于或等于第二灰度值时，则顶点位于等值面之外，当第一灰度值小于第二灰度值时，则顶点位于等值面之内；之后根据体素中所有顶点与等值面的关系确定等值面在每个体素中的剖分方式。
具体地，体素的每个顶点都存在对应的标量值(即灰度值)。如果体素顶点上的灰度值大于或等于等值面的灰度值，则定义该顶点位于等值面之外，标记为“0”；而如果体素顶点上的值小于等值面值，则定义该顶点位于等值面之内，标记为“1”。由于每个体素单元有8个顶点，那么共存在2^8＝256种情形，请参见图7，图7为本发明实施例提供的一种Marching Cubes算法的15种基本剖分方式的示意图，除图7中的15种基本情形外，其他241种情形可以通过这15种基本情形的旋转、映射等方式实现。
S34、采用三线性插值法计算在剖分方式下等值面与每个体素的交点位置。
本实施例中，首先根据步骤S33确定的剖分方式确定等值面与体素相交的体素棱边；然后利用三线性插值法计算等值面与体素棱边的交点坐标；最后由该交点坐标得到等值面与每个体素的交点位置。
具体地，根据体素的状态表索引，找出与等值面相交的体素棱边，然后采用三线性插值的方法，计算出各个交点的位置坐标。请参见图8，图8为本发明实施例提供的一种体素顶点和边的索引示意图，其中，edge index表示边的索引，vertex表示体素顶点的索引，isosurface facet表示等值面；在图8的(a)中，假如体素下方的顶点3的灰度值小于等值面的灰度值，其他顶点上的值都大于等值面值，那么可以生成一个与体素棱边2、3、11相交的三角面片，如图8的(b)图所示。三角面片顶点的具体位置则需要根据等值面的灰度值和体素顶点的灰度值由三线性插值法计算得到。
三线性插值法是在三维离散采样数据的张量积网格上进行线性插值的方法，这种方法通过网格上数据点在局部的矩形棱柱上线性地近似计算点(x,y,z)的值。对于体单元即体素内的任一点，其值只能从体素的八个顶点的采样值来插值估算。
请参见图9，图9为本发明实施例提供的一种单元立方体坐标示意图。在图9的单元立方体中，其原点是下/左/基顶点，则每个顶点的灰度值将表示为V000,V100,V010,...,V111。
图9中立方体内(x,y,z)的值将被表示为等值面的灰度值Vxyz与每个顶点的灰度值V000,V100,V010,...,V111的关系，并由下式得到：
其中，Vxyz为等值面的灰度值，V000,V100,V010,......,V111为体素中各个顶点的灰度值。
然后，将体素中八个顶点的灰度值和等值面的灰度值带入公式(14)，就可得出等值面与体素在这条棱边上交点的坐标值，进而得到等值面与体素的焦点位置。
S35、判断是否遍历所有包含等值面的体素；若是，则得到三维重建图像，完成绘制；若否，则跳转至步骤S31。
在CT成像的过程中，受系统噪声、量子噪声、重建算法等因素影响会不可避免的产生噪声。针对含噪CT图像的三维重建场景，本实施例包含两个部分：首先，提出一种自适应形态学算法与双边滤波器结合的滤波算法，不仅能够平滑的对图像进行滤波，而且能够很好地保留图像的边缘信息和几何信息，提高了含噪CT图像的质量。其次，本实施例提出了基于区域增长法和三线性插值的移动立方体算法进行CT图像的三维重建；利用区域增长法，把所有包含等值面的体素全部筛选出来再进行计算，这样有效的加快了重建效率；通过三线性插值算法计算体素与等值面的交点，充分利用体素8个顶点的灰度值计算交点，有效地提高了重建精度。结果表明：对于含噪CT图像，通过峰值信噪比与结构相似性指标的结果，本实施例的滤波算法对比单一滤波算法取得了较好的效果。进一步，对于肝脏病人，本实施例的RT-MC算法在重建精度与效率上均得到了提高，其中，重建精度较传统的MC算法、S-MC算法(SHUAI Ren-jun提出的MC算法)、G-MC算法(Masala G L提出的MC算法)提高了14.2％，13.6％，2.2％。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。
Claims (4)
1.一种噪声CT图像三维重建方法，其特征在于，包括步骤：
S1、利用自适应形态学滤波算法对灰度图像进行滤波，得到第一输出图像；包括步骤：
S11、获取所述灰度图像；
S12、设定包括形状序列和尺寸序列的结构元素集Anm，其中，Anm＝{A11,A12,A13,…,A21,A22,…,Anm}，n代表形状序列，m代表尺寸序列，
并利用所述结构元素集中的相同尺寸、不同形状的结构元素腐蚀所述灰度图像，根据所述结构元素可填入所述灰度图像的次数自适应确定不同形状的结构元素的权重值：
...
其中，α1,α2,…,αn为不同形状的结构元素权重值，β1,β2,…,βn为不同形状的结构元素可填入灰度图像的次数；
S13、对同一形状的所述结构元素，按尺寸从小到大的顺序将所述结构元素对所述灰度图像形态学处理的过程进行排序，构成串联滤波器；
S14、结合所述权重值，按所述结构元素的不同形状将所述串联滤波器进行并联，构成复合滤波器，并利用所述复合滤波器对所述灰度图像进行滤波，得到所述第一输出图像：
F(x)＝∑αifi(x)
其中，fi(x)为串联滤波器对灰度图像的滤波结果，αi为权重值；i＝1,2,…,n；
S2、利用双边滤波器对所述第一输出图像进行滤波，得到第二输出图像；包括步骤：
S21、获取基于空间距离的高斯权重区域滤波器：
其中，c(ξ,x)为基于空间距离的高斯权重，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，kd(x)为归一化因子，
S22、获取基于像素间相似度的边缘滤波器：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，s(f(ξ),f(x))为基于像素相似度的高斯权重，kr(x)为归一化因子，
S23、将所述高斯权重区域滤波器与所述边缘滤波器结合，得到所述双边滤波器：
其中，h(x)为输出图像，f(ξ)为输入图像，ξ为邻域中心x的附近点，k(x)为归一化因子，
c(ξ,x)为基于空间距离的高斯权重，
且有d(ξ,x)＝d(ξ-x)＝||ξ-x||，σd为基于距离的正态分布标准差，s(f(ξ),f(x))为基于像素相似度的高斯权重，
且有
σr为基于空间的正态分布标准差；
并利用所述双边滤波器对所述第一输出图像进行滤波，得到所述第二输出图像；
S3、利用基于区域增长法和三线性插值法的移动立方体算法对所述第二输出图像进行重建，得到三维重建图像；包括步骤：
S31、根据所述第二输出图像创建多维数据集，并从所述多维数据集中随机选取种子体素；具体包括：扫描所述第二输出图像中相邻的两个二维图像，根据其中一个所述二维图像上的四个相邻像素与另一个所述二维图像上的四个相邻像素创建所述多维数据集，并利用二次取样方法从所述多维数据集中选取包含等值面的所述种子体素；
S32、以所述种子体素的邻域体素为目标区域，利用所述区域增长法提取所述目标区域中包含等值面的所有体素；具体包括：获取所述种子体素与所述邻域体素的公共面的四个顶点的灰度值；判断所述四个顶点的灰度值与预设阈值的关系；当所述四个顶点的灰度值同时大于或者小于所述预设阈值时，则所述邻域体素中不包含所述等值面；当所述四个顶点的灰度值不是同时大于或者小于所述预设阈值时，则所述邻域体素中含所述等值面；
S33、根据所述体素中顶点的第一灰度值与所述等值面的第二灰度值的关系，确定所述等值面在每个所述体素中的剖分方式；
S34、采用所述三线性插值法计算在所述剖分方式下所述等值面与每个所述体素的交点位置；
S35、判断是否遍历所有包含所述等值面的体素；若是，则得到所述三维重建图像；若否，则跳转至步骤S31。
2.如权利要求1所述的噪声CT图像三维重建方法，其特征在于，步骤S32包括：
以所述种子体素的邻域体素为目标区域，利用区域增长法判断所述目标区域中是否包含所述等值面，若是，则跳转步骤S33，若否，则跳转步骤S31。
3.如权利要求1所述的噪声CT图像三维重建方法，其特征在于，步骤S33包括：
判断所述第一灰度值和所述第二灰度值的关系，当所述第一灰度值大于或等于所述第二灰度值时，则所述顶点位于所述等值面之外，当所述第一灰度值小于所述第二灰度值时，则所述顶点位于所述等值面之内；
根据所述体素中所有所述顶点与所述等值面的关系确定所述等值面在每个所述体素中的剖分方式。
4.如权利要求1所述的噪声CT图像三维重建方法，其特征在于，步骤S34包括：
根据所述剖分方式确定所述等值面与所述体素相交的体素棱边；
利用所述三线性插值法计算所述等值面与所述体素棱边的交点坐标：
Vxyz＝V000(1-x)(1-y)(1-z)+
V100 x(1-y)(1-z)+
V010(1-x)y(1-z)+
V001(1-x)(1-y)z+
V101 x(1-y)z+
V011(1-x)y z+
V110 x y(1-z)+
V111 x y z
其中，Vxyz为等值面的灰度值，V000,V100,V010,......,V111为体素中各个顶点的灰度值；
由所述交点坐标得到所述等值面与每个所述体素的交点位置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010591066.8A CN111739139B (zh) | 2020-06-24 | 2020-06-24 | 一种噪声ct图像三维重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010591066.8A CN111739139B (zh) | 2020-06-24 | 2020-06-24 | 一种噪声ct图像三维重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111739139A CN111739139A (zh) | 2020-10-02 |
| CN111739139B true CN111739139B (zh) | 2023-04-07 |
Family
ID=72651084
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010591066.8A Active CN111739139B (zh) | 2020-06-24 | 2020-06-24 | 一种噪声ct图像三维重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111739139B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116402917B (zh) * | 2023-06-09 | 2023-08-15 | 之江实验室 | 宽谱光散斑自相关成像的待重建图像的确定方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101692286A (zh) * | 2009-09-15 | 2010-04-07 | 上海交通大学 | 医学图像三视图获取方法 |
| CN105913396A (zh) * | 2016-04-11 | 2016-08-31 | 湖南源信光电科技有限公司 | 一种噪声估计的图像边缘保持混合去噪方法 |
| CN106373168A (zh) * | 2016-11-24 | 2017-02-01 | 北京三体高创科技有限公司 | 一种基于医疗图像的分割与三维重建方法、3d打印系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE102010043975B4 (de) * | 2010-11-16 | 2021-07-29 | Siemens Healthcare Gmbh | Verfahren zur Reduktion der verwendeten Strahlendosis im Rahmen einer bildgebenden Röntgenuntersuchung und Computersystem |
-
2020
- 2020-06-24 CN CN202010591066.8A patent/CN111739139B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101692286A (zh) * | 2009-09-15 | 2010-04-07 | 上海交通大学 | 医学图像三视图获取方法 |
| CN105913396A (zh) * | 2016-04-11 | 2016-08-31 | 湖南源信光电科技有限公司 | 一种噪声估计的图像边缘保持混合去噪方法 |
| CN106373168A (zh) * | 2016-11-24 | 2017-02-01 | 北京三体高创科技有限公司 | 一种基于医疗图像的分割与三维重建方法、3d打印系统 |
Non-Patent Citations (3)
| Title |
|---|
| Extraction of lung blood vessels from chest X-ray CT image;G Fukano等;《leice Technical Report》;20011231;全文 * |
| 基于MC算法的高质量脊柱CT图像三维重建;许婉露等;《计算机应用与软件》;20130615(第06期);全文 * |
| 改进的Canny算子在裂缝检测中的应用;赵芳等;《电子测量技术》;20181023(第20期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111739139A (zh) | 2020-10-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Fokina et al. | Microstructure synthesis using style-based generative adversarial networks | |
| Miklos et al. | Discrete scale axis representations for 3D geometry | |
| Chen et al. | Nonlocal prior Bayesian tomographic reconstruction | |
| JP5038642B2 (ja) | ボリュメトリック画像強調システム及び方法 | |
| CN112132959B (zh) | 数字岩心图像处理方法、装置、计算机设备及存储介质 | |
| US20060214932A1 (en) | Fast graph cuts: a weak shape assumption provides a fast exact method for graph cuts segmentation | |
| CN107341765B (zh) | 一种基于卡通纹理分解的图像超分辨率重建方法 | |
| CN112116543B (zh) | 基于检测式生成框架的图像修复方法、系统及装置 | |
| CN109544478B (zh) | 一种基于奇异值分解的非局部均值ct图像降噪方法 | |
| Alpers et al. | Geometric reconstruction methods for electron tomography | |
| CN111724307B (zh) | 一种基于最大后验概率和非局部低秩先验的图像超分辨重建方法，终端及可读存储介质 | |
| CN113744136A (zh) | 基于通道约束多特征融合的图像超分辨率重建方法和系统 | |
| CN108364297A (zh) | 血管图像分割方法、终端、存储介质 | |
| CN111739139B (zh) | 一种噪声ct图像三维重建方法 | |
| Jakhongir et al. | 3D Volume Reconstruction from MRI Slices based on VTK | |
| Hadwiger et al. | Sparse pdf maps for non-linear multi-resolution image operations. | |
| CN117575915A (zh) | 一种图像超分辨率重建方法、终端设备及存储介质 | |
| CN117523084A (zh) | 一种移动扫描点云的自动三维重建方法及电子设备 | |
| Mandal et al. | Employing structural and statistical information to learn dictionary (s) for single image super-resolution in sparse domain | |
| Suryanarayana et al. | Deep Learned Singular Residual Network for Super Resolution Reconstruction. | |
| Zheng et al. | VDVR: Verifiable volume visualization of projection-based data | |
| CN108876734B (zh) | 图像去噪方法、装置、电子设备及存储介质 | |
| Song et al. | Super resolution reconstruction of medical image based on adaptive quad-tree decomposition | |
| WO2024078049A1 (en) | System and method for near real-time and unsupervised coordinate projection network for computed tomography images reconstruction | |
| Kiatpapan et al. | Super-resolution based on back-projection of interpolated image |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |