专利名称：一种三维数字体积图像变形测量方法
技术领域：
本发明涉及一种测量方法，特别是关于一种用于定量分析三维数字体积图像变形 测量方法。
背景技术：
现有的二维数字图像相关技术的基本原理在于对比变形前的数字图像(参考图) 与变形后的图像(变形图)，即将参考图与变形图进行对比。再通过做一系列的交叉相关 运算，来获得变形后的图像相对于变形前图像的位移场和应变场。目前，传统二维数字图 像相关方法已经成为一种经典的实验力学技术在材料力学行为测试、微尺度变形测量、微 电子机械系统(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)中最 大元素所确定的整体素位移值，(Διι，Δ ν, Aw)是由基于梯度的三维亚体素位移定位算法 计算所得的亚体素位移值；(7)重复所述步骤(3) 步骤(6)，计算所有参考体积图像中的 采样点在变形后数字体积图像中的准确位置，进而得到整个变形后体积图像相对于参考体 积图像的三维位移场；(8)根据连续介质力学中拉格朗日应变张量理论计算数字体积图像的三维应变场du. du, ，其中，1 ( i，j，m< 3,U1 = U,U2 = V,U3 =W。所述步骤(2)中，在所述LxXLyXLz个体素的长方体中，Lx = Ly = Lz，且它们取 30 60之间的整数，单位为体素；在所述MxXMyXMz个体素的长方体中，Mx = My = Mz,且它 们取35 100之间的整数，单位为体素，且Mx = My = Mz > Lx = Ly = Lz。所述步骤(4)中，所述三维求和表Sg和分别sg(χ, y, z) = g(x, y, ζ)+sg(χ, y-1, z)+sg(x_l, y, z)+sg(x_l, y-1, z_l)+sg(x, y, z_l)_sg(x, y-1, z_l)_sg(x_l, y, z_l)_sg(x_l, y-1, z)，
权利要求
一种三维数字体积图像变形测量方法，其包括如下步骤(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>ϵ</mi> <mi>ij</mi></msub><mo>=</mo><mfrac> <mn>1</mn> <mn>2</mn></mfrac><mo>[</mo><mfrac> <mrow><mo>∂</mo><msub> <mi>U</mi> <mi>i</mi></msub> </mrow> <mrow><mo>∂</mo><msub> <mi>x</mi> <mi>j</mi></msub> </mrow></mfrac><mo>+</mo><mfrac> <mrow><mo>∂</mo><msub> <mi>U</mi> <mi>j</mi></msub> </mrow> <mrow><mo>∂</mo><msub> <mi>x</mi> <mi>i</mi></msub> </mrow></mfrac><mo>+</mo><munderover> <mi>Σ</mi> <mrow><mi>m</mi><mo>=</mo><mn>1</mn> </mrow> <mn>3</mn></munderover><mfrac> <mrow><mo>∂</mo><msub> <mi>U</mi> <mi>m</mi></msub> </mrow> <mrow><mo>∂</mo><msub> <mi>x</mi> <mi>i</mi></msub> </mrow></mfrac><mfrac> <mrow><mo>∂</mo><msub> <mi>U</mi> <mi>m</mi></msub> </mrow> <mrow><mo>∂</mo><msub> <mi>x</mi> <mi>j</mi></msub> </mrow></mfrac><mo>]</mo><mo>,</mo> </mrow>其中，1≤i，j，m≤3，U1＝U，U2＝V，U3＝W。FSA00000320344800012.tif,FSA00000320344800013.tif,FSA00000320344800014.tif,FSA00000320344800015.tif,FSA00000320344800016.tif
2.如权利要求1所述的一种三维数字体积图像变形测量方法，其特征在于所述步骤 ⑵中，在所述Lx X Ly X Lz个体素的长方体中，Lx = Ly = Lz，且它们取30 60之间的整数， 单位为体素；在所述MxXMyXMz个体素的长方体中，Mx = Mv = Mz，且它们取35 100之间 的整数，单位为体素，且Mx = My = Mz > Lx = Ly = Lz。
3.如权利要求1所述的一种三维数字体积图像变形测量方法，其特征在于所述步骤 (4)中，所述三维求和表Sg和&2分别sg(x, y, ζ) = g(x, y, z)+sg(x, y-1, z)+sg(x_l, y, z)+sg(x_l, y-1, z_l)+sg(x, y, z-1)-sg(x, y-1, z-l)-sg(x-l, y, z-l)-sg(x-l, y-1, ζ), 2(x,_y,z) = [g(x，>^)]2+ 2 (x,y-l,z) + sg2(x-i,y,z) + sg29+Sgl (x,y,z-\)~sg2 (x,y-l,z-\)-sg2 (x-l,y,z-l)-sg2 (x-l,y-l,z)其中，三维求和表Sg和5>分别包含MxXMyXMz个元素；g(x，y，z)代表搜索区内中心坐 标为(χ，ι, ζ)的体素所对应的灰度值；而χ，y，ζ取整数，且当χ，y，ζ彡0时，& = & = 0。
4.如权利要求2所述的一种三维数字体积图像变形测量方法，其特征在于所述步骤 (4)中，所述三维求和表Sg和&2分别Sg(x, y, ζ) = g(x, y, z)+sg(x, y-1, z)+sg(x_l, y, z)+sg(x_l, y-1, z_l)+sg(x, y, z-1)-sg(x, y-1, z_l)_sg(x_l, y, z_l)_sg(x_l, y-1, z)，sg2 (x,y,z) = [_g(x,y,z)^ +sgl {x,y-\,z) +sgl (x-l,y,z) +sg2 {x-\,y-],z-\) +^2 (x,y,z-\)-sg2 (x,y-l,z-\)-sg2 (x-l,y,z-l)-sg2 (x-\,y-\,z)其中，三维求和表Sg和&2分别包含MxXMyXMz个元素；g(x，y，z)代表搜索区内中心坐 标为(χ，ι, ζ)的体素所对应的灰度值；而χ, y，ζ取整数，且当χ，y，ζ彡0时，& = ^2 = 0。
5.如权利要求1或2或3或4所述的一种三维数字体积图像变形测量方法，其特征在 于所述步骤(5)中，所述三维零均值归一化互相关系数矩阵[C(u，V，w)(a,e,Y)为[c(M,v,w)1 P(丫」，L v4Fp(u,v,w)i, h L1Jt=Iz=\Lx LV Lz、Lx Lv Lz一 ‘X、^Z、1λ=1 y=\ z=l^xLyLzX=I 户1Lx Ly L.2 Lx+u Ly+v ι +w^ =1z=l*=!+ >'=1+VZ=1+Wx^y^zLx+U 人ν+ · L2+wΣ Σ Σ 外,少，ζ)λ:=1+μ y=]+v ζ=\·￥\ν Lx+u Ly+ν L.+wΣ Σ Σx=l+ j=l+v r=!+ -Lx Ly Lz.ΣΣΣ/Ο^,ζ)_v=l z=lL、 厶 ν LzLx Ly Lzλχ=\ y=\ ζ=1X=I y=l Z=ILxLyLz_JLx Lv LzJLx Ly LZ且 7=“，“ ΣΣΣ/(Χ，·^) ； ^=+ + + ；Lx^LyX [ζ ^=I y=l z=lj^x ^ j^y ^ htZ ■*=^=IP(u, V, w)由所述步骤(3)中的三维快速傅立叶变换方法快速求解，即 V{u,v,w) = FT；1 [FT3 (f)FT； (g)] ； Q(u, ν, w)和 G(u，ν, w)中的三重求和项通过所述步骤(5)中的快速查表方式求解；F通过算术运算求解。
全文摘要
本发明涉及一种三维数字体积图像变形测量方法，其包括步骤(1)在变形前数字体积图像上选取离散采样点；(2)确定参考子体积及搜索区域尺寸；(3)基于三维快速傅立叶变换方法对参考子体积及搜索区域进行三维图像相关运算；(4)建立图像灰度三维求和表和图像能量三维求和表；(5)利用步骤(3)运算结果及三维求和表求解三维零均值归一化互相关系数矩阵；(6)运用基于梯度的三维亚体素位移定位算法计算采样点的亚体素位移值；(7)根据步骤(3)～步骤(6)计算所有采样点处的准确位移值，得到整个变形后体积图像相对于变形前体积图像的三维位移场；(8)根据连续介质力学中拉格朗日应变张量理论计算数字体积图像的三维应变场。本发明能精确、高效地测量数字体积图像的三维变形。
文档编号G06T17/00GK101980304SQ20101052067
公开日2011年2月23日 申请日期2010年10月20日 优先权日2010年10月20日
发明者彭小玲, 方竞, 李姗姗, 潘晓畅, 熊春阳, 黄建永 申请人:北京大学