本发明公开了一种基于交叉爬山memetic量子进化计算的医学图像配准方法，主要解决现有技术配准效果差的问题。其实现步骤为：读入图像，初始相关参数；用chaotic方法产生初始种群；计算种群中每个个体对应的配准参数；按照配准参数对浮动图像进行图像变换，得到变换后的图像；计算变换后的图像与参考图像之间的相似性，并寻找本代最优配准参数，按照最优配准参数对种群中所有个体进行量子旋转门更新；对最优配准参数进行交叉爬山memetic局部学习；判断最优个体是否满足遗忘条件，对最优个体执行遗忘算子；控制循环条件，若循环结束，则输出配准结果。本发明能获得更好的配准结果，可用于医学图像的配准。
1.一种基于交叉爬山memetic量子进化计算的医学图像配准方法，包括如下步骤：(1)读入两幅待配准图像，将其中一幅固定不变的图像作为参考图像IR，另一幅进行形变的图像作为浮动图像IF，并按间距L＝35个像素，对浮动图像均匀采样，得到由m个元素组成的均匀网格；
(2)用chaotic方法产生一组个体，作为种群，其中每个个体是一个1×m的向量；
(3)将种群中每个个体转化为候选解，再与均匀网格叠加，得到与每个个体对应的控制网格，作为配准参数；
所述将种群中每个个体转化为候选解，是将个体p＝[θ1,θ2,...,θm]映射到配准参数解空间，得到具有m个分量的解S＝[S1,S2,...,Sm]，各分量表达式如下：Si＝ai·θi+bi·(1-θi)，i＝1,2,…m，
其中，Si是S的第i个分量，ai和bi分别是第i个分量的取值下限和上限；
(4)按照配准参数对浮动图像做三次B样条函数变换，并计算变换后的浮动图像与参考图像之间的归一化互信息值，找出本次迭代中最大的归一化互信息值及其对应的配准参数和个体，并将该个体存储为本代最优个体pbest；
(5)依据本代最优个体pbest，利用量子旋转门算子对本代种群中所有个体进行量子旋转门更新，使图像配准参数朝相似性更大的方向移动：(5a)按照迭代次数g构造旋转角度的大小Δθ0：
Δθ0＝0.015π×(g/max_g)+0.15π×exp(-mod(g,100)/10)，其中，max_g是设定的最大迭代次数，mod(g,100)是取g除以100的余数，exp(-mod(g,
100)/10)是取e的-mod(g,100)/10次方；
(5b)取出第g代种群中个体p的第i个分量 判断该分量是否越过最优个体的第i个分量，若越过，则旋转角度 反之，旋转角度 则第g+1代的个体p的第i个分量 满足：
(5c)重复步骤(5b)直至个体p的所有分量全部更新；
(6)选出除最优个体之外的任意一个个体，并与最优个体一起作为交叉爬山memetic局部学习算子的初始点，进行局部学习，进一步优化本次迭代的最优个体pbest；
(7)判断最优个体pbest是否发生改变，若经连续F次迭代后pbest未发生改变，则最优个体满足遗忘条件，并采用遗忘算子将最优个体遗忘，避免配准参数陷入局部极值，否则，执行步骤(8)，F为设定的大于2的正整数；
(8)更新迭代次数g，如迭代次数达到最大值，则结束配准，最优个体对应的控制网格即为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准图像，否则，返回步骤(3)进行下一次迭代。
2.根据权利要求书1所述的基于交叉爬山memetic量子进化计算的医学图像配准方法，其中步骤(2)所述的用chaotic方法产生一组个体，按如下步骤进行：(2a)先生成一个随机值θ1，0<θ1<1，θ1≠0.5，再通过θ1计算整个chaotic个体p＝[θ1,θ2,...,θm]，其中：θi+1＝4·θi·(1-θi)，
其中，i＝1,2,…m-1，m是均匀网格的元素个数；
(2b)重复步骤(2a)产生一组个体P＝[p1,p2,...,pN]，组成个数为N的种群。
3.根据权利要求书1所述的基于交叉爬山memetic量子进化计算的医学图像配准方法，其中步骤(6)所述的选出除最优个体之外的任意一个个体，并与最优个体一起作为交叉爬山memetic局部学习算子的初始点，进行局部学习，是按如下步骤进行：(6a)选出种群中除最优个体之外的任意一个个体pr，设定局部学习次数X＝10次；
(6b)用BLX交叉算子生成所述个体pr与最优个体pbest的一组子代个体off＝[off1,off2,...,offR]，R是子代个体总数，其中第k个子代个体offk为：offk＝(1-ξk)·pr+ξk·pbest，
其中，k＝1,2,…R，ξk＝U[-α,1+α]，U为均匀分布，α是可调参数，0<α<1；
(6c)按照步骤(3)与步骤(4)所述的方法计算子代off＝[off1,off2,...,offR]对应的归一化互信息值 找出NMIoff中的最大值NMIoffmax，以及NMIoffmax对应的最优子代个体offbest；
(6d)将最优子代个体offbest对应的归一化互信息值NMIoffmax与所述个体pr对应的归一化互信息值NMIr比较，若NMIoffmax>NMIr，则令所述个体pr是最优子代个体offbest，即pr＝offbest，并令所述个体pr对应的归一化互信息值NMIr是最优子代个体offbest对应的归一化互信息值NMIoffmax，即NMIr＝NMIoffmax，执行步骤(6f)，否则，执行步骤(6e)；
(6e)将最优子代个体offbest对应的归一化互信息值NMIoffmax与最优个体pbest对应的归一化互信息值NMIbest比较，若NMIoffbest>NMIbest，则令最优个体pbest是最优子代个体offbest，即pbest＝offbest，并令最优个体pbest对应的归一化互信息值NMIbest是最优子代个体offbest对应的归一化互信息值NMIoffmax，即NMIbest＝NMIoffmax，否则，执行步骤(6f)；
(6f)判断局部学习次数是否达到X＝10次，若达到，则停止局部学习，否则，返回步骤(6b)。
基于交叉爬山memetic量子进化计算的医学图像配准方法
技术领域
[0001] 本发明属于图像处理技术领域，涉及医学图像配准，具体的说是一种新的基于交叉爬山memetic量子进化计算的医学图像配准方法，可用于后续的医学图像理解和处理。
背景技术
[0002] 医学图像配准是一种通过优化算法寻找空间变换参数的方法，用该参数对两幅不同的医学图像做空间变换，可使这两幅医学图像在空间位置上达到一一对应，或者是在具有手术意义的点上达到匹配。参与配准的两幅图像可以是不同成像方式的、不同患者的或者是同一患者在不同时间拍摄的医学图像。
[0003] 医学图像配准方法一般包含特征空间、搜索空间、相似性测度和搜索算法。其中：
[0004] 特征空间是图像配准方法所使用的图像特征，包括灰度、拐点、边界线等；搜索空间是图像配准方法中对图像进行形变的方式和范围，包括刚性变换、非刚性变换等；相似性测度是度量配准图像与参考图像直接相似性的一个量化值，包括相关系数、互信息、归一化互信息等；搜索算法是图像配准方法中寻找配准参数所使用的优化算法，包括最速下降法、遗传算法等。按照医学图像配准方法使用的搜索算法的不同，可以将其分为基于传统优化的配准方法和基于智能优化的配准方法。不同的优化算法对于整个配准方法的结果有至关重要的影响。现有的配准方法主要使用传统优化算法，这类方法利用目标函数的数学特性构造优化方法，虽然速度较快，但配准程度不够高，极易陷入局部极值，特别是对于优化参数较多的非刚性配准。近年也出现了基于遗传算法的医学图像配准，但遗传算法编码长度较长，且仍然会收敛到局部极值，影响配准效果。
发明内容
[0005] 本发明的目的在于针对已有配准方法的不足，提出一种新的基于交叉爬山memetic量子进化计算的医学图像配准方法，以减少收敛到局部极值，提高配准效果。
[0006] 为实现上述目的，本发明技术方案包括如下步骤：
[0007] (1)读入两幅待配准图像，将其中一幅固定不变的图像作为参考图像IR，另一幅进行形变的图像作为浮动图像IF，并按间距L＝35个像素，对浮动图像均匀采样，得到由m个元素组成的均匀网格；
[0008] (2)用chaotic方法产生一组个体，作为种群，其中每个个体是一个1×m的向量；
[0009] (3)将种群中每个个体转化为候选解，再与均匀网格叠加，得到与每个个体对应的控制网格，作为配准参数；
[0010] (4)按照配准参数对浮动图像做三次B样条函数变换，并计算变换后的浮动图像与参考图像之间的归一化互信息值，找出本次迭代中最大的归一化互信息值及其对应的配准参数和个体，并将该个体存储为本代最优个体pbest；
[0011] (5)依据本代最优个体pbest，利用量子旋转门算子对本代种群中所有个体进行量子旋转门更新，使图像配准参数朝相似性更大的方向移动；
[0012] (6)选出除最优个体之外的任意一个个体，并与最优个体一起作为交叉爬山memetic局部学习算子的初始点，进行局部学习，进一步优化本次迭代的最优个体pbest；
[0013] (7)判断最优个体pbest是否发生改变，若经连续F次迭代后pbest未发生改变，则最优个体满足遗忘条件，并采用遗忘算子将最优个体遗忘，避免配准参数陷入局部极值，否则，执行步骤(8)，F为设定的大于2的正整数；
[0014] (8)更新迭代次数g，如迭代次数达到最大值，则结束配准，最优个体对应的控制网格即为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准图像，否则，返回步骤(3)进行下一次迭代。
[0015] 本发明与现有的技术相比具有以下优点：
[0016] 1.本发明加入了遗忘算子，可有效避免搜索算法陷入局部极值，提高配准图像与参考图像的相似性，提高配准效果；
[0017] 2.本发明在量子进化计算后加入局部学习阶段，在迭代过程中对量子进化计算得到的最优个体再进行一次局部优化，可很好地将全局搜索与局部搜索结合起来，使搜索算法更有效，提高配准效果。
附图说明
[0018] 图1是本发明的总流程图；
[0019] 图2是本发明仿真实验的参考图像；
[0020] 图3是本发明仿真实验的浮动图像。
具体实施方式
[0021] 参照图1，本发明的具体实现步骤如下：
[0022] 步骤1，读入图像，初始相关参数。
[0023] 读入两幅待配准图像，一幅图像固定不变，另一幅为要变换的图像；将其中固定不变的图像作为参考图像IR，将另一幅进行变换的图像作为浮动图像IF，并按间距L＝35个像素的个数，对浮动图像进行均匀采样，得到由m个元素组成的均匀网格；
[0024] 设定最大迭代次数max_g＝100，设定局部学习次数X＝10次，并将当前迭代次数g置零。
[0025] 步骤2，用chaotic方法产生初始种群。
[0026] (2a)生成一个随机值θ1，0＜θ1＜1，θ1≠0.5，通过θ1计算整个chaotic个体p＝[θ1,θ2,...,θm]，其中：
[0027] θi+1＝4·θi·(1-θi)，
[0028] 其中，i＝1,2,…m-1，m是均匀网格的元素个数；
[0029] (2b)重复步骤(2a)产生一组个体P＝[p1,p2,...,pN]，组成个数为N的种群。
[0030] 步骤3，计算种群中每个个体对应的控制网格。
[0031] (3a)选择种群P＝[p1,p2,...,pN]中的一个个体pn，n＝1,2,…N；
[0032] (3b)将个体pn映射到配准参数解空间，得到具有m个分量的解Sn＝[Sn1,Sn2,...,Snm]，各分量表达式如下：
[0033] Sni＝ai·θni+bi·(1-θni)，i＝1,2,…m，
[0034] 其中，Sni是Sn的第i个分量，ai和bi分别是第i个分量的取值下限和上限；
[0035] (3c)将上述Sn整合为与均匀网格相同大小的矩阵，并与均匀网格叠加得到个体pn对应的控制网格Φn；
[0036] (3d)重复步骤(3a)-(3c)，得到种群P＝[p1,p2,...,pN]对应的控制网格Φ＝[Φ1,Φ2,...,ΦN]。
[0037] 步骤4，按照控制网格对浮动图像做三次B样条函数变换，得到变换后的图像，并计算变换后的图像与浮动图像的相似性。
[0038] (4a)取控制网格Φ＝[Φ1,Φ2,...,ΦN]中的一个网格Φn，n＝1,2,…N；
[0039] (4b)取浮动图像一个像素点，记为(x,y)，按照控制网格Φn计算该像素点的位移值D(x,y)，即：
[0040]
[0041] 其中， 和 分别表示对x和y取整， 是控制网格Φn中该像素点(x,y)周围4×4控制网格点的值，Bq(e)和Bl(f)是三次B样条基函数Bd(t)中的一种，d＝0,1,2,3，0≤t＜1，三次B样条基函数Bd(t)表示如下：
[0042]
[0043] (4c)重复步骤(4b)计算浮动图像每个像素点的位移值，得到Φn对应的变换后的图像In；
[0044] (4d)计算参考图像IR与变换后的图像In之间的归一化互信息值NMIn为：
[0045]
[0046] 其中，H(IR)是参考图像IR的熵，H(In)是变换后的图像In的熵，H(IR,In)是参考图像IR与变换后的浮动图像In的联合熵；
[0047] (4e)重复步骤(4a)-(4d)得到种群P＝[p1,p2,...,pN]对应的归一化互信息值NMI＝[NMI1,NMI2,...,NMIN]；
[0048] (4f)找出NMI＝[NMI1,NMI2,...,NMIN]中的最大值及其对应的个体，并将该个体存储为本次迭代的最优个体pbest，将pbest对应的控制网格储存为本次迭代的最优配准参数Φbest，将pbest对应的归一化互信息值存储为最大归一化互信息NMIbest。
[0049] 步骤5，依据本代最优个体，对本代所有个体进行量子旋转门更新。
[0050] (5a)按照迭代次数g构造旋转角度的大小Δθ0：
[0051] Δθ0＝0.015π×(g/max_g)+0.15π×exp(-mod(g,100)/10)，
[0052] 其中，max_g是设定的最大迭代次数，mod(g,100)是取g除以100的余数，exp(-mod(g,100)/10)是取e的-mod(g,100)/10次方；
[0053] (5b)取出第g代种群中个体p；
[0054] (5c)取出个体p的第i个分量 如果该分量越过了最优个体的第i个分量，则反之， 第g+1代的个体p的第i个分量 满足：
[0055]
[0056] (5d)重复步骤(5c)直至个体p的所有分量全部更新；
[0057] (5e)重复步骤(5b)-(5d)直至种群中所有个体全部更新。
[0058] 步骤6，选出本代最优个体与除最优个体之外的任一个体，进行交叉爬山memetic局部学习。
[0059] (6a)选出种群中除最优个体之外的任意一个个体pr；
[0060] (6b)用BLX交叉算子生成所述个体pr与最优个体pbest的一组子代个体off＝[off1,off2,...,offR]，R是子代个体总数，其中第k个子代个体offk为：
[0061] offk＝(1-ξk)·pr+ξk·pbest，
[0062] 其中，k＝1,2,…R，ξk＝U[-α,1+α]，U为均匀分布，α是可调参数，0＜α＜1；
[0063] (6c)按照步骤(3)与步骤(4)所述的方法计算子代off＝[off1,off2,...,offR]对应的归一化互信息值NMIoff＝[NMIoff1,NMIoff2,...,NMIoffR]，找出NMIoff中的最大值NMIoffmax，以及NMIoffmax对应的最优子代个体offbest；
[0064] (6d)将最优子代个体offbest对应的归一化互信息值NMIoffmax与所述个体pr对应的归一化互信息值NMIr比较，若NMIoffmax＞NMIr，则令所述个体pr是最优子代个体offbest，即pr＝offbest，并令所述个体pr对应的归一化互信息值NMIr是最优子代个体offbest对应的归一化互信息值NMIoffmax，即NMIr＝NMIoffmax，执行步骤(6f)，否则，执行步骤(6e)；
[0065] (6e)将最优子代个体offbest对应的归一化互信息值NMIoffmax与最优个体pbest对应的归一化互信息值NMIbest比较，若NMIoffbest＞NMIbest，则令最优个体pbest是最优子代个体offbest，即pbest＝offbest，并令最优个体pbest对应的归一化互信息值NMIbest是最优子代个体offbest对应的归一化互信息值NMIoffmax，即NMIbest＝NMIoffmax，否则，执行步骤(6f)；
[0066] (6f)判断局部学习次数是否达到X=10次，若达到，则停止局部学习，否则，返回步骤(6b)。
[0067] 步骤7，判断是否满足遗忘条件，执行遗忘算子。
[0068] (7a)设定一个大于2的正整数F,本实例选F＝10；
[0069] (7b)判断最优个体pbest是否发生改变，若经连续F次迭代后pbest未发生改变，则最优个体pbest满足遗忘条件，执行步骤(7c)，否则，执行步骤(8)；
[0070] (7c)删除最优个体pbest，并将pbest对应的最大归一化互信息NMIbest置零。
[0071] 步骤8，更新迭代次数，判断是否达到结束条件。
[0072] 更新迭代次数g＝g+1，若迭代次数达到最大值max_g=100，则结束配准，将最优个体pbest对应的控制网格Φbest作为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准后的图像，输出最大归一化互信息值NMIbest，否则，返回步骤(3)进行下一次迭代。
[0073] 本发明的效果可以通过下面的实验仿真进一步说明：
[0074] 1、实验条件与方法
[0075] 硬件平台为：Intel(R)Core(TM)i5-2450M@2.50GHz、3.91GBRAM.；
[0076] 软件平台为：MATLAB R2012b；
[0077] 实验方法：分别用基于最速下降法的方法、基于遗传算法的方法以及本发明仿真实现医学图像配准。
[0078] 2、仿真内容
[0079] 仿真实验所用图像如图2和图3所示，图2是正常人的脑部MRI图像，图3是病人的脑部MRI图像。
[0080] 为了验证本发明的有效性，分别用基于最速下降法的方法、基于遗传算法的方法以及本发明对图2和图3所示的实验图像进行配准。将每种方法独立仿真100次，得到归一化互信息值的平均值和归一化互信息值的最大值两个结果。
[0081] 3、仿真结果
[0082] 仿真得到的两个结果如表1所示，其中基于最速下降法的方法是方法1，基于遗传算法的方法是方法2。
[0083] 表1三种种配准方法结果对比表
[0084] 方法1 方法2 本发明
平均归一化互信息 1.1422 1.1386 1.1490
最大归一化互信息 1.1422 1.1429 1.1544
[0085] 由表1可见，本发明方法相较于基于最速下降法的方法和遗传算法的方法在平均归一化互信息和最大归一化互信息两个方面都更加优秀，提高了配准效果，算法性能更优，是一种有效的医学图像配准方法。
