一种三维数字体积图像变形测量方法技术

本发明专利技术涉及一种三维数字体积图像变形测量方法，其包括步骤：（１）在变形前数字体积图像上选取离散采样点；（２）确定参考子体积及搜索区域尺寸；（３）基于三维快速傅立叶变换方法对参考子体积及搜索区域进行三维图像相关运算；（４）建立图像灰度三维求和表和图像能量三维求和表；（５）利用步骤（３）运算结果及三维求和表求解三维零均值归一化互相关系数矩阵；（６）运用基于梯度的三维亚体素位移定位算法计算采样点的亚体素位移值；（７）根据步骤（３）～步骤（６）计算所有采样点处的准确位移值，得到整个变形后体积图像相对于变形前体积图像的三维位移场；（８）根据连续介质力学中拉格朗日应变张量理论计算数字体积图像的三维应变场。本发明专利技术能精确、高效地测量数字体积图像的三维变形。

【技术实现步骤摘要】

本专利技术涉及一种测量方法，特别是关于一种用于定量分析三维数字体积图像变形 测量方法。
技术介绍
现有的二维数字图像相关技术的基本原理在于对比变形前的数字图像(参考图) 与变形后的图像(变形图)，即将参考图与变形图进行对比。再通过做一系列的交叉相关 运算，来获得变形后的图像相对于变形前图像的位移场和应变场。目前，传统二维数字图 像相关方法已经成为一种经典的实验力学技术在材料力学行为测试、微尺度变形测量、微 电子机械系统(MEMS)结构动态表征以及生物组织力学性能评估等诸多领域都得到了广泛 应用。在二维数字图像相关方法的基础上，Bay与其合作者于1999年提出了数字体积相 关(Digital Volume Correlation, DVC)的概念(Bay B. K.，Smith Τ. S.，Fyhrie D. P.， and Saad Μ. , Digital νοlumecorreIation Three-dimensional strain mapping using x-ray tomography, Experimental Mechanics, 1999，39 (3)，217-226 ；贝尔，史密斯，费锐， 塞得，数字体积图像相关使用X射线断层摄影技术进行三维应变表征，实验力学，1999， 39 (3)，217-226)，并首先应用于疏松骨质三维X射线微结构的变形测量。数字体积相关方 法的基本原理与数字图像相关方法类似，其核心思想在于把变形前体图像中特定的参考子 体积(Reference Subvolume)与变形后体图像中相应搜索区域内所有可能的目标子体积 (Target Subvolume)之间作三维图像互相关运算，通过相关函数的极值条件来确定参考体 子体积在变形后体图像中最可能的位置，从而得到数字体积图像的三维位移场和应变场。近年来随着微CT (Micro-computed Tomography, μ CT)，核磁共振成像(Magnetic Resonance Imaging, MRI)等三维微结构形貌探测技术的发展，人们开始借助于数字体积 相关方法定量分析多孔蜂窝材料(例如骨骼，木质材料，粘土质材料以及泡沫结构材料 等)的三维变形。随着激光共聚焦显微方法(laser-scanning confocal microscopy, LSCM)的快速发展，人们开始探索用数字体积相关方法研究含有随机标记物且透光性 较好软材料(例如含有随机分布蛋白纤维的胶原材料)的三维变形(Roeder B. A.， Kokini K. , Robinson J. P. , andVoytik-Harbiη S. L. , Local, three-dimensional strain measurements withinlargely deformed extracellular matrix constructs,Journal of BiomechanicalEngineering, 2004，126 (6)，699-708 ；饶德，考克尼，罗宾逊，瓦提克-哈斌， 细胞外基质结构局部三维大变形应变测量，生物力学工程杂志，2004，126(6)，699-708)。然 而，与发展较为成熟的数字图像相关方法对比，数字体积相关方法在求解精度和计算效率 方面还亟待发展完善，其应用领域仍需进一步拓展。目前，数字体积相关方法在应用过程中还只是主要考虑参考子体积与目标子体积 之间的整体素(Voxel，也称立体像素)相关匹配运算，且通常在匹配过程中至多考虑子体 积之间的刚体平移和旋转效应，数字体积相关方法的亚体素(Subvoxel)定位算法框架还 没有系统建立起来，特别是基于激光共聚焦成像的三维图像亚体素定位算法研究还是空白。另一方面，由数字体积相关方法的基本原理可知，对于每一个选定的参考子体积，都需 要与相应搜索区域内所有可能的目标子体积逐个进行三维体积相关运算，因此，不难发现 其相关运算的计算量要比类似情况下数字图像相关方法的搜索计算量大得多。例如，当子 体积尺寸为32 X 32 X 32体素，搜索区域为62 X 62 X 62体素时，参考子体积与目标子体积之 间至少需要进行30 X 30 X 30 = 27000次三维体积相关运算，而每一次体积相关计算又需要 至少进行三个三维求和运算。经粗略估计表明，此种情况下数字体积相关的计算量要比类 似情况下(当子区尺寸为32X32像素，搜索区尺寸为62X62像素时)数字图像相关的计 算量增加近似三个量级。以上传统数字体积相关方法在计算上的复杂性已经成为制约其进 一步广泛应用的主要障碍，特别是在需要进行大样本统计分析的情况下，传统数字体积相 关方法的计算效率根本不能满足要求。
技术实现思路
针对上述问题，本专利技术的目的是提供一种计算效率较高、具有亚体素测量精度的 三维数字体积图像变形测量方法。为实现上述目的，本专利技术采取以下技术方案一种三维数字体积图像变形测量方 法，其包括如下步骤(1)将未发生变形的三维数字体积图像作为参考体积图像F(X，y，z)， 将发生变形的三维数字体积图像作为变形体积图像G(x，y，ζ)；在参考体积图像F(x，y， ζ)上选定一系列离散采样点，设其中任意一个采样点的坐标为(α，β，Α) ； (2)分别以所 述步骤(1)中的采样点为中心，选定包含LxXLyXLz个体素的长方体作为参考子体积f(x， y,z)；在变形体积图像G(x，y，z)中，分别以采样点坐标(α，β，Α)为中心，各自选定包含 MxXMyXMz个体素的长方体作为搜索区域g(x，y，ζ)，且Mx = My = Mz > Lx = Ly = Lz ； (3) 对于每一个采样点(α，β，￥)，分别对其对应的参考子体积￡0^，7，2)和相应的搜索区域 g(x，y，z)进行三维图像相关运算，其运算结果P为？ =巧'(g)]，其中，*代表复共轭；FT3[ ·]和分别表示三维快速傅立叶正变换和逆变换；f = [/(AhzMvvii和g = !^^:^)]^^分别代表参考子体积f(x，y，z)图像灰度三维矩阵和搜索区域g(x，y，ζ)体积图像灰度三维矩阵；(4)根据每一个采样点(α，β，Y)所对应的搜索区域g(x， y，ζ)的图像灰度值，利用快速递推法分别建立变形后数字体积图像灰度三维求和表\及 图像能量三维求和表； (5)通过所述步骤(4)所建立的三维求和表Sg和&2进行快速查 表运算，并根据所述步骤(3)的相关运算结果，计算采样点(α，β，Y)所对应的参考子体 积f(x，y，z)和相应的搜索区域g(x，y，z)之间的三维零均值归一化互相关系数矩阵[C(u， V, w)](a,0,Y) ； (6)采用基于梯度的三维亚体素位移定位算法，在零均值归一化互相关系数 矩阵的最大峰值所在位置附近进行亚体素插值运算，求得参考数字体积图像中采样点(a， β , Y)在相应变形后体积图像中的准确位置坐标(U，V，W)，其中，U = U+Δ u，V = ν+Δ ν， W = w+Aw, (u,v,w)代表由所述三维零均值归一化相关系数矩阵[C(U，V，w)](a,e,Y)中最 大元素所确定的整体素位移值，(本文档来自技高网...

【技术保护点】
一种三维数字体积图像变形测量方法，其包括如下步骤：（１）将未发生变形的三维数字体积图像作为参考体积图像Ｆ（ｘ，ｙ，ｚ），将发生变形的三维数字体积图像作为变形体积图像Ｇ（ｘ，ｙ，ｚ）；在参考体积图像Ｆ（ｘ，ｙ，ｚ）上选定一系列离散采样点，设其中任意一个采样点的坐标为（α，β，γ）；（２）分别以所述步骤（１）中的采样点为中心，选定包含Ｌ↓［ｘ］×Ｌ↓［ｙ］×Ｌ↓［ｚ］个体素的长方体作为参考子体积ｆ（ｘ，ｙ，ｚ）；在变形体积图像Ｇ（ｘ，ｙ，ｚ）中，分别以采样点坐标（α，β，γ）为中心，各自选定包含Ｍ↓［ｘ］×Ｍ↓［ｙ］×Ｍ↓［ｚ］个体素的长方体作为搜索区域ｇ（ｘ，ｙ，ｚ），且Ｍ↓［ｘ］＝Ｍ↓［ｙ］＝Ｍ↓［ｚ］＞Ｌ↓［ｘ］＝Ｌ↓［ｙ］＝Ｌ↓［ｚ］；（３）对于每一个采样点（α，β，γ），分别对其对应的参考子体积ｆ（ｘ，ｙ，ｚ）和相应的搜索区域ｇ（ｘ，ｙ，ｚ）进行三维图像相关运算，其运算结果Ｐ为：Ｐ＝ＦＴ↓［３］↑［－１］［ＦＴ↓［３］（ｆ）ＦＴ↓［３］↑［＊］（ｇ）］，其中，＊代表复共轭；ＦＴ↓［３］［.］和ＦＴ↓［３］↑［－１］［.］分别表示三维快速傅立叶正变换和逆变换；ｆ＝［ｆ（ｘ，ｙ，ｚ）］↓［Ｌ↓［ｘ］×Ｌ↓［ｙ］×Ｌ↓［ｚ］］和ｇ＝［ｇ（ｘ，ｙ，ｚ）］↓［Ｍ↓［ｘ］×Ｍ↓［ｙ］×Ｍ↓［ｚ］］分别代表参考子体积ｆ（ｘ，ｙ，ｚ）图像灰度三维矩阵和搜索区域ｇ（ｘ，ｙ，ｚ）体积图像灰度三维矩阵；（４）根据每一个采样点（α，β，γ）所对应的搜索区域ｇ（ｘ，ｙ，ｚ）的图像灰度值，利用快速递推法分别建立变形后数字体积图像灰度三维求和表Ｓ↓［ｇ］及图像能量三维求和表Ｓ↓［ｇ↑［２］］；（５）通过所述步骤（４）所建立的三维求和表Ｓ↓［ｇ］和Ｓ↓［ｇ↑［２］］进行快速查表运算，并根据所述步骤（３）的相关运算结果，计算采样点（α，β，γ）所对应的参考子体积ｆ（ｘ，ｙ，ｚ）和相应的搜索区域ｇ（ｘ，ｙ，ｚ）之间的三维零均值归一化互相关系数矩阵［Ｃ（ｕ，ｖ，ｗ）］↓［（α，β，γ）］；（６）采用基于梯度的三维亚体素位移定位算法，在零均值归一化互相关系数矩阵的最大峰值所在位置附近进行亚体素插值运算，求得参考数字体积图像中采样点（α，β，γ）在相应变形后体积图像中的准确位置坐标（Ｕ，Ｖ，Ｗ），其中，Ｕ＝ｕ＋Δｕ，Ｖ＝ｖ＋Δｖ，Ｗ＝ｗ＋Δｗ，（ｕ，ｖ，ｗ）代表由所述三维零均值归一化相关系数矩阵［Ｃ（ｕ，ｖ，ｗ）］↓［（α，β，γ）］中最大元素所确...

【技术特征摘要】
一种三维数字体积图像变形测量方法，其包括如下步骤(1)将未发生变形的三维数字体积图像作为参考体积图像F(x，y，z)，将发生变形的三维数字体积图像作为变形体积图像G(x，y，z)；在参考体积图像F(x，y，z)上选定一系列离散采样点，设其中任意一个采样点的坐标为(α，β，γ)；(2)分别以所述步骤(1)中的采样点为中心，选定包含Lx×Ly×Lz个体素的长方体作为参考子体积f(x，y，z)；在变形体积图像G(x，y，z)中，分别以采样点坐标(α，β，γ)为中心，各自选定包含Mx×My×Mz个体素的长方体作为搜索区域g(x，y，z)，且Mx＝My＝Mz＞Lx＝Ly＝Lz；(3)对于每一个采样点(α，β，γ)，分别对其对应的参考子体积f(x，y，z)和相应的搜索区域g(x，y，z)进行三维图像相关运算，其运算结果P为 <mrow><mi>P</mi><mo>=</mo><msubsup> <mi>FT</mi> <mn>3</mn> <mrow><mo>-</mo><mn>1</mn> </mrow></msubsup><mo>[</mo><msub> <mi>FT</mi> <mn>3</mn></msub><mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo></mrow><msubsup> <mi>FT</mi> <mn>3</mn> <mo>*</mo></msubsup><mrow> <mo>(</mo> <mi>g</mi> <mo>)</mo></mrow><mo>]</mo><mo>,</mo> </mrow>其中，*代表复共轭；FT3[·]和分别表示三维快速傅立叶正变换和逆变换；和分别代表参考子体积f(x，y，z)图像灰度三维矩阵和搜索区域g(x，y，z)体积图像灰度三维矩阵；(4)根据每一个采样点(α，β，γ)所对应的搜索区域g(x，y，z)的图像灰度值，利用快速递推法分别建立变形后数字体积图像灰度三维求和表Sg及图像能量三维求和表(5)通过所述步骤(4)所建立的三维求和表Sg和进行快速查表运算，并根据所述步骤(3)的相关运算结果，计算采样点(α，β，γ)所对应的参考子体积f(x，y，z)和相应的搜索区域g(x，y，z)之间的三维零均值归一化互相关系数矩阵[C(u，v，w)](α，β，γ)；(6)采用基于梯度的三维亚体素位移定位算法，在零均值归一化互相关系数矩阵的最大峰值所在位置附近进行亚体素插值运算，求得参考数字体积图像中采样点(α，β，γ)在相应变形后体积图像中的准确位置坐标(U，V，W)，其中，U＝u+Δu，V＝v+Δv，W＝w+Δw，(u，v，w)代表由所述三维零均值归一化相关系数矩阵[C(u，v，w)](α，β，γ)中最大元素所确定的整体素位移值，(Δu，Δv，Δw)是由基于梯度的三维亚体素位移定位算法计算所得的亚体素位移值；(7)重复所述步骤(3)～步骤(6)，计算所有参考体积图像中的采样点在变形后数字体积图像中的准确位置，进而得到整个变形后体积图像相对于参考体积图像的三维位移场；(8)根据连续介质力学中拉格朗日应变张量理论计算数字体积图像的三维应变场 <mrow><msub> <mi>&epsiv;</mi> <mi>ij</mi></msub><mo>=</mo><mfrac> <mn>1</mn> <mn>2</mn></mfrac><mo>[</mo><mfrac> <mrow><mo>&PartialD;</mo><msub> <mi>U</mi> <mi>i</mi></msub> </mrow> <mrow><mo>&PartialD;</mo><msub> <mi>x</mi> <mi>j</mi></msub> </mrow></mfrac><mo>+</mo><mfrac> <mrow><mo>&PartialD;</mo><msub> <mi>U</mi> <mi>j</mi></msub> </mrow> <mrow><mo>&PartialD;</mo><msub> <mi>x</mi> <mi>i</mi></msub> </mrow></mfrac><mo>+</mo><munderover> <mi>&Sigma;</mi> <mrow><...

【专利技术属性】
技术研发人员：黄建永，潘晓畅，李姗姗，彭小玲，熊春阳，方竞，
申请(专利权)人：北京大学，
类型：发明
国别省市：11[]