CN105096349A - 一种图像对称结构的提取方法 - Google Patents

Abstract

本发明提供一种医学图像对称结构的提取方法，包括：提供医学图像及初始对称，根据所述初始对称将所述医学图像划分为位于所述初始对称两侧的两部分；基于所述两侧的像素信息建立基于对称参数的目标函数；迭代优化计算对称参数，进而提取所述对称参数对应的所述医学图像对称结构。本发明借鉴类似图像配准的迭代优化方法，利用对称坐标方程将图像分成位于对称线或对称面两侧的两部分，基于两边像素点的目标函数作为代价函数进行优化迭代，上述优化方法具有更加准确、快速的优点，并且鲁棒性更好。

Description

一种图像对称结构的提取方法

技术领域

本发明涉及医学图像处理领域，尤其涉及一种医学图像对称结构的提取方法。

背景技术

在医学图像处理中，图像的线性对称结构(如二维图像的对称线，三位图像的对称面)的提取具有广泛的应用。

美国专利申请文件US5533143提出了一种检测图像是否对称的方法，判断准则主要依靠图像两边像素值之差是否为0。但是所述专利申请中，假设前提是图像的对称轴是和通过图像中心的一条垂直的中心线一致的，这显然应用是很狭隘的。而且两边在差值为0也是很苛刻的条件，很容易受噪声、病变等影响。

美国专利申请文件US6263097提出以图像两边像素灰度值的均方差为判断准则，不断优化对称线的角度来得到最后的对称线。另一篇美国专利申请文件US7409085首先计算出图像质心(xc，yc)，选定一个搜索区域(-15mm<＝[x-xc]<＝15mm，-2mm<＝[y-yc]<＝2mm，-20。<＝β<＝20。)，计算此区域内可能性最大的一条对称线(x0，y0，β0)，再重定义搜索区域(-3mm<＝[x-x0]<＝3mm，-3mm<＝[y-y0]<＝3mm，-3。<＝[β-β0]<＝3。)，重复计算可能性最大的对称线。显然上述方法是基于遍历寻找最优解的方法，与基于迭代的最优化方法相比，在速度和准确性等方面都有着明显的不足。

文献BabakA.Ardekani，etc，‘AutomaticDetectionoftheMid-SagittalPlanein3-DBrainImages’，IEEEtransactionSonmedicalimaging，vol.，16，N).6，December1997.采用单纯形下降法优化求解了3D头部图像的对称面，代价函数选取互相关函数。考虑到医学图像的复杂性，并且代价函数一般为对称面方程参数的非线性函数，单纯形法并不是非常适合的。

美国专利申请文件US20070276219用优化方法寻找相互正交的三个面，其中第一个平面为所求的对称面。代价函数选取三个面与各自相平行的临近数个面对应点的灰度熵值，采用基于单纯形的Nelder-Mead优化方法。美国专利US20117986823目的是为了保证从多次扫描检测出的对称面的一致性(ConSistent)，把医生(MR设备operator)标定的MSP作为金标准，然后以互信息为代价函数衡量自动检测的MSP同金标准的差异。实现上需要专家检测MSP，如果不同病人的MSP第一次扫描(localizer)均需人工检测时，操作上会非常麻烦。

发明内容

本发明解决的问题是提供一种医学图像的对称结构的提取方法，用以对医学图像的对称结构进行提取。

为了解决上述问题，本发明提供了一种医学图像对称结构的提取方法，包括：提供医学图像及初始对称，根据所述初始对称将所述医学图像划分为位于所述初始对称两侧的两部分；基于所述两侧的像素信息建立基于对称参数的目标函数；迭代优化计算对称参数，进而提取所述对称参数对应的所述医学图像对称结构。

可选的，所述初始对称包括初始中心线和初始中心面：若所述医学图像为二维图像，则所述初始中心线为对称线方程x+k*y+c＝0，其中，k，c为待定的对称参数；若所述医学图像为三维图形，则所述初始中心面为对称面方程x＝k₁y+k₂z+b，其中，K1，K2，b为待定的对称参数。

可选的，所述迭代优化包括：利用lucas-kanade算法、梯度下降法或者BFGS方法对所述目标函数优化迭代。

可选的，建立所述目标函数包括：基于所述两侧像素的均方差、互相关或者互信息建立所述目标函数。

可选的，基于所述目标函数计算对称参数还包括：对所述目标函数进行泰勒展开，以获取增量目标函数。

可选的，包括：通过迭代方法求解所述对称参数，每次迭代中，获取令所述增量目标函数值最小对应的增量，并将所述增量后的对称参数替换当前对称参数，进行迭代优化。

可选的，所述目标函数为C＝∑_Pi∈L[I(R(p_i|Q))-I(p_i)]²，其中，其中，I(p_i)为一侧的点p_i的像素值，L为遍历所述点p_i所在一侧的像素值，I(R(p_i|Q))为与p_i对称的点R(p_i|Q)的像素值，所述对称为以当前对称参数Q的对称。

可选的，通过迭代方法求解对称参数Q，及每次迭代使得所述目标函数最小的增量ΔQ，并令增量后的对称参数Q+ΔQ替换当前对称参数Q进行迭代优化。

可选的，还包括采用泰勒展开对所述目标函数进行展开：若所述医学图像为二维图像，展开的目标函数表达式为 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中， $\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y})$ 为医学图像在点R(p_i|Q)的偏导数；若所述医学图像为三维图像，则所述展开的目标函数为 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中 $\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y},\frac{\&PartialD;I}{\&PartialD;z})$ 为图像I在点R(p_i|Q)的偏导数。

可选的，所述目标函数在ΔQ处具有极小值，即所述目标函数对ΔQ的偏导为0，若所述医学图像为二维图像，所述目标函数对ΔQ的偏导的表达式如下：

$2{\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(R\left(p_i\right|k,c)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right)]=0,$ 获取ΔQ＝H^-1S，其中 $H={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}],S={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(p_i\right)-I\left(R\left(p_i\right|Q)\right)],$ $\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}=(\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;k}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{\&PartialD;k},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;c}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{\&PartialD;c});$

若所述医学图像为三维图像，则所述目标函数对ΔQ的偏导的表达式如下 $2{\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(R\left(p_i\right|k,c)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right)]=0,$ 获取ΔQ＝H^-1S，其中 $H={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}],S={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(p_i\right)-I\left(R\left(p_i\right|Q)\right)],$ $\begin{array}{c}\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}=(\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{{\&PartialD;k}_1}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;z}\frac{\&PartialD;R_z}{{\&PartialD;k}_1},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;z}\frac{{\&PartialD;R}_z}{{\&PartialD;k}_2},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;b}+\\ \frac{\&PartialD;I}{\&PartialD;y}\frac{\&PartialD;R_y}{\&PartialD;b}+\frac{\&PartialD;I}{\&PartialD;z}\frac{{\&PartialD;R}_z}{\&PartialD;b}.\end{array}$

可选的，若所述医学图像为二维图像，设像素沿对称线x+k*y+c＝0的对称点的变换为R(R_x，R_y)，对于点p_i(x，y)则其对称变换为： $R_x=x-\frac{2(x+\mathrm{ky}+c)}{k*k+1},R_y=y-\frac{2k(x+\mathrm{ky}+1)}{k*k+1},$ 对应求得

$\frac{\&PartialD;R_x}{\&PartialD;k}=\frac{4k(x+\mathrm{ky}+c)-2y(k^2+1)}{{(k^2+1)}^2},\frac{\&PartialD;R_y}{\&PartialD;k}=\frac{{4k}^2(x+\mathrm{ky}+c)-(k^2+1)(2x+4\mathrm{ky}+\mathrm{zc})}{{\left(k^2\right|1)}^2},$ $\frac{{\&PartialD;R}_x}{\&PartialD;c}=-\frac{2}{k^2+1},\frac{{\&PartialD;R}_y}{\&PartialD;c}=-\frac{2k}{k^2+1};$ 若所述医学图像为三维图像，则设像素沿对称面x-k₁*y-k₂*z-b＝0的对称点的变换为R(R_x，R_y，R_z)，对于点p_i(x，y，y，)则其对称变换 $R_x=x-\frac{2*F}{\mathrm{KF}},R_y=y+\frac{2*k_1*F}{\mathrm{KF}},R_z=z+\frac{2*k_2*F}{\mathrm{KF}},$ 其中F＝x-k₁*y-k₂*z-b，KF＝1+k₁ ²+k₂ ²，并可求得对称变换R_x，R_y，R_z对参数k₁，k₂，b的偏导数。

可选的，所述迭代的终止条件为：当所述对称参数的增量的模小于最小迭代步长时或者迭代次数不小于最大迭代次数时，所述迭代终止。

本发明对二维图像选取对称线方程x+k*y+c＝0，对于三维图像选取对称面方程x＝k₁y+k₂z+b，借鉴类似图像配准的迭代优化方法，利用对称坐标方程将图像分成位于对称线或对称面两侧的左、右两部分，基于两边像素点的目标函数作为代价函数进行优化迭代，上述优化方法具有更加准确、快速的优点，并且鲁棒性更好。

附图说明

图1所示为本发明一个实施例的医学图像的对称结构提取方法的方法流程示意图。

具体实施方式

在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施的限制。

其次，本发明利用示意图进行详细描述，在详述本发明实施例时，为便于说明，所述示意图只是实例，其在此不应限制本发明保护的范围。

为了解决上述问题，本发明提供了一种医学图像对称结构的提取方法，包括：提供医学图像及初始对称，根据所述初始对称将所述医学图像划分为位于所述初始对称两侧的两部分；基于所述两侧的像素信息建立基于对称参数的目标函数；迭代优化计算对称参数，进而提取所述对称参数对应的所述医学图像对称结构。建立所述目标函数包括：基于所述两侧像素的均方差、互相关或者互信息建立所述目标函数。所述迭代优化包括利用lucas-kanade算法、梯度下降法或者BFGS方法对所述目标函数优化迭代。

进一步地，通过迭代方法求解所述对称参数，每次迭代中，获取令所述增量目标函数的值最小对应的增量，并将所述增量后的对称参数替换当前对称参数，进行迭代优化。基于所述目标函数计算对称参数包括：对所述目标函数进行泰勒展开，以获取增量目标函数。

其中，所述初始对称包括初始中心线和初始中心面：若所述医学图像为二维图像，则所述初始中心线为对称线方程x+k*y+c＝0，其中，k，c为待定的对称参数；若所述医学图像为三维图形，则所述初始中心面为对称面方程x＝k₁y+k₂z+b，其中，K1，K2，b为待定的对称参数。

所述目标函数为C＝∑_Pi∈L[I(R(p_i|Q))-I(p_i)]²，其中，其中，I(p_i)为一侧的点p_i的像素值，L为遍历所述两侧中的其中一侧的像素值，I(R(p_i|Q))为与p_i对称的点R(p_i|Q)的像素值，所述对称为以当前对称参数Q的对称。

具体地，通过迭代方法求解对称参数Q包括：每次迭代获取使得C＝∑_Pi∈L[I(R(p_i|Q+ΔQ))-I(p_i)]²最小的增量ΔQ，通过求解获取的增量后的对称参数Q+ΔQ替换当前对称参数Q，进行迭代优化。

采用泰勒展开对所述目标函数进行优化：若所述医学图像为二维图像，优化的目标函数表达式为 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中，为医学图像在点R(p_i1Q)的偏导数；若所述医学图像为三维图像，则所述优化的目标函数为 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中 $\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y},\frac{\&PartialD;I}{\&PartialD;z})$ 为图像I在点R(p_i|Q)的偏导数。

如图1所示为本发明一个实施例的医学图像的对称结构提取方法，包括：执行步骤S1，输入图像，计算灰度质心并初始化参数；执行步骤S2，判断像素点是否位于左侧。所述左侧为位于当前对称参数对应的对称线或者对称面的左侧。

上述判断步骤S2中若为否，则执行步骤S3，输入下一点像素点，并继续后续的判断和方法流程，即执行步骤S2，判断像素点是否位于左侧。

若所述判断步骤S2为是，则执行步骤S4，计算原坐标点与对称点灰度差，及图像在对称点梯度，并将对称变换对系数求偏导，累加S及H矩阵，具体地，S，H矩阵将下述具体实施例详述。

接着执行步骤S5，判断像素是否遍历完毕。若判断为是，则执行步骤S6，求H逆矩阵，并更新系数。若判断为否，则执行步骤S3，输入下一点像素点。并继续执行步骤S3的后续步骤，即执行步骤S2，判断像素点是否位于左侧。

接着，执行步骤S6的更新系数后，执行步骤S7，判断步长是否小于设定。若是，则执行步骤S9，输出系数。

若所述步骤S7中的判断为否，则执行步骤S8，输入下一遍历中的像素点。并继续执行步骤S8中的后续步骤，即执行步骤S2，像素点是否位于左侧。

下面分别对二维医学图像和三维医学图像的具体实施方法是进行详细说明，若所述医学图像为二维医学图像，则对应的对称提取包括如下步骤：

(1)输入2D图像，计算灰度质心坐标(x₀，y₀)，初始迭代次数为0；

(2)设定被优化参数的初始值k＝0。并根据中心线经过质心的条件计算出c＝-x₀-k*y₀，从而初始中心线为f(x，y|k，c)＝x+k*y+c＝0；

(3)对当前系数进行迭代。

具体地，目标函数为C＝Σ_Pi∈L[I(R(p_i|Q))-I(p_i)]²，而系数Q为二维向量(k，c)，通过迭代方法求解Q，每次迭代求使C＝∑_Pi∈L[I(R(p_i|Q+ΔQ))-I(p_i)]²最小的ΔQ，通过Q←Q+ΔQ，进行迭代。对改进的目标函数公式进行一阶泰勒展开，可得 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中 $\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y})$ 为图像I在点R(p_i|Q)的偏导数。

设像素沿直线x+k*y+c＝0的对称点的变换为R(R_x，R_y)，对于点p_i(x，y)则其对称变换 $R_x=x-\frac{2(x+\mathrm{ky}+c)}{k*k+1},R_y=y-\frac{2k(x+\mathrm{ky}+1)}{k*k+1},$ 可求得 $\frac{\&PartialD;R_x}{\&PartialD;k}=\frac{4k(x+\mathrm{ky}+c)-2y(k^2+1)}{{(k^2+1)}^2},\frac{\&PartialD;R_y}{\&PartialD;k}=\frac{{4k}^2(x+\mathrm{ky}+c)-(k^2+1)(2x+4\mathrm{ky}+\mathrm{zc})}{{\left(k^2\right|1)}^2},$ $\frac{{\&PartialD;R}_x}{\&PartialD;c}=-\frac{2}{k^2+1},\frac{{\&PartialD;R}_y}{\&PartialD;c}=-\frac{2k}{k^2+1}.$

因为目标函数在ΔQ处取极小值，所以目标函数对ΔQ的偏导为0，即令使泰勒展开后的目标函数对ΔQ的偏导为0，得到公式 $2{\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(R\left(p_i\right|k,c)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right)]=0,$ 可转化为ΔQ＝H^-1S，其中 $H={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}],S={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(p_i\right)-I\left(R\left(p_i\right|Q)\right)],$ $\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}=(\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;k}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{\&PartialD;k},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;c}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{\&PartialD;c})$

(4)遍历全部像素，判断像素p_i是否位于当前中心线左侧，如果是进行下一步，否则遍历下一像素。

(5)计算像素p_i右侧镜像对称点的坐标，在图像中插值，计算原像素点与对称点的灰度之差I(p_i)-I(R(p_i|Q))。

(6)计算像素p_i点对应的并对S进行累加，同时对H矩阵进行累加。

(7)当按步骤(4)，(5)，(6)对全部像素遍历完毕，对H求逆矩阵H^-1，并利用公式ΔQ＝H^-1S得到当前迭代的ΔQ，更新迭代系数。迭代次数加一。

(8)当向量ΔQ的模小于最小迭代步长时或者迭代次数不小于最大迭代次数时，迭代终止，输出参数，否则转到步骤(3)。

对应地，若所述医学图像为三维医学图像，则对应的对称提取包括如下步骤：

(1)输入3D图像，计算灰度质心坐标(x₀，y₀，z₀)，初始迭代次数为0；

(2)设定被优化参数的初始值k＝0。并根据d对称面经过质心的条件计算出b＝-x₀-k₁*y₀-k₂*z₀，从而初始中心线为f(x，y，z|k₁，k₂，b)＝x-k₁*y-k₂*z-b＝0;

(3)对当前系数进行迭代。目标函数为C＝Σ_Pi∈L[I(R(p_i|Q))-I(p_i)]²而系数Q为二维向量(k₁，k₂，b)，通过迭代方法求解Q，每次迭代求使C＝∑_Pi∈L[I(R(p_i|Q+ΔQ))-I(p_i)]²最小的ΔQ，通过Q←Q+ΔQ，进行迭代。对改进的目标函数公式进行一阶泰勒展开，可得 $C={\mathrm{\&Sigma;}}_{p_i}[I\left(R\left(P_i\right|Q)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right){]}^2,$ 其中 $\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y},\frac{\&PartialD;I}{\&PartialD;z})$ 为图像I在点R(p_i|Q)的偏导数。

设像素沿对称面x-k₁*y-k₂*z-b＝0的对称点的变换为R(R_x，R_y，R_z)，对于点p_i(x，y，z)则其对称变换 $R_x=x-\frac{2*F}{\mathrm{KF}},R_y=y+\frac{2*k_1*F}{\mathrm{KF}},R_z=z+\frac{2*k_2*F}{\mathrm{KF}},$ 其中F＝x-k₁*y-k₂*z-b，KF＝1+k₁ ²+k₂ ²，并可求得对称变换R_x，R_y，R_z对参数k₁，k₂，b的偏导数。使泰勒展开后的目标函数对ΔQ的偏导为0，得到公式 $2{\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(R\left(p_i\right|k,c)\right)+\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}\mathrm{\&Delta;Q}-I\left(p_i\right)]=0,$ 可转化为ΔQ＝H-¹S，其中 $H={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}],S={\mathrm{\&Sigma;}}_{\mathrm{Pi}\&Element;L}{[\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}]}^T[I\left(p_i\right)-I\left(R\left(p_i\right|Q)\right)],$ $\begin{array}{c}\&dtri;I\frac{\&PartialD;R}{\&PartialD;Q}=(\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{{\&PartialD;k}_1}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;z}\frac{\&PartialD;R_z}{{\&PartialD;k}_1},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;y}\frac{{\&PartialD;R}_y}{{\&PartialD;k}_2}+\frac{\&PartialD;I}{\&PartialD;z}\frac{{\&PartialD;R}_z}{{\&PartialD;k}_2},\frac{\&PartialD;I}{\&PartialD;x}\frac{{\&PartialD;R}_x}{\&PartialD;b}+\\ \frac{\&PartialD;I}{\&PartialD;y}\frac{\&PartialD;R_y}{\&PartialD;b}+\frac{\&PartialD;I}{\&PartialD;z}\frac{{\&PartialD;R}_z}{\&PartialD;b}).\end{array}$

(4)遍历全部像素，判断像素p_i是否位于当前对称面左侧，如果是进行下一步，否则遍历下一像素。

(5)计算像素p_i右侧镜像对称点的坐标，在图像中插值，计算原像素点与对称点的灰度之差I(p_i)-I(R(p_i|Q))。

(6)计算像素p_i点对应的并对S进行累加，同时对H矩阵进行累加。

(7)当按步骤(4)，(5)，(6)对全部像素遍历完毕，对H求逆矩阵H^-1，并利用公式ΔQ＝H^-1S得到当前迭代的ΔQ，更新迭代系数。迭代次数加一。

(8)当向量ΔQ的模小于最小迭代步长时或者迭代次数不小于最大迭代次数时，迭代终止，输出参数，否则转到步骤(3)。

本发明对二维图像选取对称线方程x+k*y+c＝0，对于三维图像选取对称面方程x＝k₁y+k₂z+b，借鉴类似图像配准的迭代优化方法，利用对称坐标方程将图像分成位于所述对称线或对称面两侧的左、右两部分，并分别看作图像配准中的“参考图像”和“浮动图像”；并通过选取两边像素点的均方差作为代价函数；进一步地，根据代价函数，利用lucas-kanade算法进行优化迭代，上述优化方法具有更加准确、快速的优点，并且鲁棒性更好。

在图像处理中，图像的线性对称结构，如二维图像的对称线，三维图像的对称面的提取具有广泛的应用。本文用优化方法求解了图像对称坐标方程的参数。以对称坐标方程左右两侧的像素灰度均方差作为代价函数，优化参数得到收敛，算法具有一定的实用性。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

Claims (11)

1.一种医学图像对称结构的提取方法，其特征在于，包括：提供医学图像及初始对称，根据所述初始对称将所述医学图像划分为位于所述初始对称两侧的两部分；基于所述两侧的像素信息建立基于对称参数的目标函数；迭代优化计算对称参数，进而提取所述对称参数对应的所述医学图像对称结构。

2.如权利要求1所述的医学图像对称结构的提取方法，其特征在于，所述初始对称包括初始中心线和初始中心面：若所述医学图像为二维图像，则所述初始中心线为对称线方程x+k*y+c＝0，其中，k，c为待定的对称参数；若所述医学图像为三维图形，则所述初始中心面为对称面方程x＝k₁y+k₂z+b，其中，k₁，k₂，b为待定的对称参数。

3.如权利要求1所述的医学图像对称结构的提取方法，其特征在于，所述迭代优化包括：：利用lucas-kanade算法、梯度下降法或者BFGS方法对所述目标函数优化迭代。

4.如权利要求1所述的医学图像对称结构的提取方法，其特征在于，建立所述目标函数包括：基于所述两侧像素的均方差、互相关或者互信息建立所述目标函数。

5.如权利要求1所述的医学图像对称结构的提取方法，其特征在于，基于所述目标函数计算对称参数包括：对所述目标函数进行泰勒展开，以获取增量目标函数。

6.如权利要求5所述的医学图像对称结构的提取方法，其特征在于，包括：通过迭代方法求解所述对称参数，每次迭代中，获取令所述增量目标函数值最小对应的增量，并将所述增量后的对称参数替换当前对称参数，进行迭代优化。

7.如权利要求1所述的医学图像对称结构的提取方法，其特征在于，所述目标函数为C＝Σ_Pi∈L[I(R(p_i|Q))-I(p_i)]²，其中，其中，I(p_i)为一侧的点p_i的像素值，L为遍历所述点p_i所在一侧的像素值，I(R(p_i|Q))为与p_i对称的点R(p_i|Q)的像素值，所述对称为以当前对称参数Q的对称。

8.如权利要求7所述的医学图像对称结构的提取方法，其特征在于，通过迭代方法求解对称参数Q及每次迭代使得所述目标函数最小的增量ΔQ，并令增量后的对称参数Q+ΔQ替换当前对称参数Q进行迭代优化。

9.如权利要求7所述的医学图像对称结构的提取方法，其特征在于，还包括采用泰勒展开对所述目标函数进行展开：若所述医学图像为二维图像，展开的目标函数表达式为其中，为医学图像在点R(p_i|Q)的偏导数；若所述医学图像为三维图像，则展开的目标函数为其中为图像I在点R(p_i|Q)的偏导数。

10.如权利要求7所述的医学图像对称结构的提取方法，其特征在于，还包括：所述目标函数在ΔQ处具有极小值，即所述目标函数对ΔQ的偏导为0，若所述医学图像为二维图像，所述目标函数对ΔQ的偏导的表达式如下：

获取ΔQ＝H^-1S，其中 若所述医学图像为三维图像，所述目标函数对ΔQ的偏导的表达式如下：

获取ΔQ＝H-1S，其中

11.如权利要求7所述的医学图像对称结构的提取方法，其特征在于，包括：若所述医学图像为二维图像，设像素沿对称线x+k*y+c＝0的对称点的变换为R(R_x，R_y)，对于点p_i(x，y)则其对称变换为： 对应求得

若所述医学图像为三维图像，则设像素沿对称面x-k₁*y-k₂*z-b＝0的对称点的变换为R(R_x，R_y，R_z)，对于点p_i(x，y，z)则其对称变换为：其中

F＝x-k₁*y-k₂*z-b，KF＝1+k₁ ²+k₂ ²，并可求得对称变换R_x，R_y，R_z对参数k₁，k₂，b的偏导数。

如权利要求1所述的医学图像对称结构的提取方法，其特征在于，所述迭代的终止条件为：当所述对称参数的增量的模小于最小迭代步长时或者迭代次数不小于最大迭代次数时，所述迭代终止。