本发明公开了一种基于量子进化计算和B样条变换的医学图像配准方法，主要解决现有技术匹配程度低和配准预处理过程复杂的问题。其实现步骤为：读入图像，初始相关参数；随机产生量子进化计算的种群；计算种群中每个个体对应的配准参数；按照配准参数对浮动图像做三次B样条函数变换，得到变换后的浮动图像；计算变换后的浮动图像与参考图像之间的归一化互信息值，并寻找本代最优配准参数，对本代种群中所有个体进行量子旋转门更新；用更新后的种群进行量子重组操作；控制循环条件，输出配准图像和配准参数。本发明与基于遗传算法的方法相比，不仅能获得更好的配准结果，并且缩短了配准，可用于医学图像的配准。
1.一种基于量子进化计算和B样条变换的医学图像配准方法，包括如下步骤：(1)读入两幅待配准图像，将其中一幅固定不变的图像作为参考图像IR，另一幅进行变换的图像作为浮动图像IF，并按间距L＝40个像素，对浮动图像均匀采样，得到由m个元素组成的均匀网格；
(2)随机产生一组由量子比特位构成的个体，作为量子进化计算的初始迭代点，称为种群；
(3)将种群中每个个体转化为候选解，再与均匀网格叠加，得到每个个体对应的控制网格；
(4)按照每个个体对应的控制网格对浮动图像做三次B样条函数变换，得到变换后的浮动图像；
(5)计算变换后的浮动图像与参考图像之间的归一化互信息值NMI，找出本次迭代中最大的归一化互信息值及其对应的个体，并存储为本代最优解pbest；
(6)依据本代最优解，利用量子旋转门算子对本代种群中所有个体进行量子旋转门更新，使图像配准参数朝相似性更大的方向移动；
(7)利用量子重组算子对任意选择的本代种群中两个不同个体pa和pb进行量子重组操作，增加种群个体多样性，避免图像配准参数陷入局部极值；
(8)更新迭代次数g，如迭代次数达到最大值或变换后的浮动图像与参考图像之间的归一化互信息值不再增加，则结束配准，即本代最优解对应的控制网格为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准后的图像，否则，返回步骤(3)进行下一次迭代。
2.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(2)所述的随机产生一组由量子比特位构成的个体，是按照量子比特位的方式产生，每个个体p的形态表示为：其中，θi∈[0,2π],i＝1,2,…m，m是均匀网格的元素个数，θi是随机生成的量子比特位角度。
3.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(3)所述的将种群中每个个体转化为候选解，按如下步骤进行：(3a)将个体p分解为两个子个体pc和ps：
pc＝[cos(θ1) cos(θ2) ... cos(θm)]ps＝[sin(θ1) sin(θ2) ... sin(θm)]，其中，pc是余弦个体，ps是正弦个体；
(3b)将余弦个体pc和正弦个体ps映射到配准参数解空间，得到具有m个分量的余弦解Sc和正弦解Ss，各分量表达式如下：其中，θi∈[0,2π],i＝1,2,…m，m是均匀网格的元素个数，θi是随机生成的量子比特位角度，Sci和Ssi分别是Sc和Ss的第i个分量，ai和bi分别是第i个分量的取值下限和上限，i＝1,
4.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(4)所述的按照每个个体对应的控制网格对浮动图像做三次B样条函数变换，是通过计算浮动图像中每个像素点的位移值D(x,y)完成，即：其中， 是取整操作，φ(k+q)(m′+l)是
控制网格中该像素点(x,y)周围4×4控制网格点的值，Bq(e)和Bl(f)是三次B样条基函数Bd(t)中的一种，d＝0,1,2,3，0≤t<1，表示如下：
5.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(5)所述的计算变换后的浮动图像I与参考图像IR之间的归一化互信息值NMI，按如下公式进行：其中，H(IR)是参考图像IR的熵，H(I)是变换后的浮动图像I的熵，H(IR,I)是参考图像IR与变换后的浮动图像I的联合熵。
6.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(6)所述的依据本代最优解，利用量子旋转门算子对本代种群中所有个体进行量子旋转门更新，按如下步骤进行：(6a)按照迭代次数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次方；
(6b)取出第g代种群中个体p的第i个分量 如果该分量越过了最优解的第i个分量，则 反之， 该分量对应的旋转矩阵 只由旋转角 度决定，即：
(6c)重复步骤(6b)直至个体p的所有分量全部更新；
其中，θi∈[0,2π],i＝1,2,…m，m是均匀网格的元素个数，θi是随机生成的量子比特位角度。
7.根据权利要求1所述的基于量子进化计算和B样条变换的医学图像配准方法，其中步骤(7)所述的利用量子重组算子对任意选择的本代种群中两个不同个体pa和pb进行量子重组操作，是先随机选择一个比特位h，再按如下公式进行：其中， 是个体pa的第h个分量，
基于量子进化计算和B样条变换的医学图像配准方法
技术领域
[0001] 本发明属于图像处理技术领域，涉及医学图像配准，具体的说是一种新的基于量子进化计算和B样条变换的图像配准方法，可用于医学图像融合。
背景技术
[0002] 医学图像按照其获取原则和途径的不同，可以分为X射线图像，US图像，CT图像，MRI图像，PET图像，fMRI图像，SPECT图像等等。其中，根据其表示信息的不同又主要分为解剖结构图像和功能图像，这种图像分别包含了器官的详细解剖结构和器官的新陈代谢功能。这两方面的信息对于医学诊断、疾病监测和治疗指导等方面同样重要。于是，将同一患者不同时刻的医学图像或不同患者相同器官的医学图像以及患者的与正常的医学图像里的不同信息融合在一起显得尤为重要。而图像配准作为图像融合的基础直接影响着融合效果。如何将两幅不同成像方式、不同患者或者同一患者在不同时刻的医学图像进行空间位置上的匹配，即配准，成为该领域亟待解决的问题之一。
[0003] 医学图像配准算法的发展主要分为三个阶段，基于外部特征的配准方法，以刚性变换为主的基于内部特征的配准方法和以非线性变换为主的配准方法。F.Maes等将归一化互信息作为图像配准的相似性测度，见F.Maes, A.Collignon, D.Vandermeulen, etal. Multimodality image registration by maximization of mutual information. IEEE Trans on Medical Imaging,vol.16,no.2,1997,187:198；归一化互信息的引入使得图像配准目标函数更加复杂，更加依赖于有效的优化方法。Yang Chen等在Yang Chen,Zheng Qin,Mingyue Song.A continuous medical image registration approach based on image segmentation and genetic algorithm.ICIE-CS,2009,1:4中将遗传算法用于医学图像配准中，但遗传算法编码长度较长，导致在空间上消耗相对较多，且仍然会收敛到局部极值，使算法到达全局最优值之前提前结束，影响配准效果。Hongkui Xu等在Hongkui Xu,Mingyan Jiang, Mingqiang Yang. An Image Registration method combining feature constraint with multilevel Strategy中利用图像形状特征实现图像配准，但需要通过复杂的步骤提取图像的形状特征，且由于非刚性组织结构的复杂性，有些分界面不是很明显，不易正确提取形状特征。
发明内容
[0004] 本发明的目的在于针对已有配准方法的不足，提出一种新的基于量子进化计算和B样条变换的医学图像配准方法，以提高配准效果和简化配准过程。
[0005] 为实现上述目的，本发明包含如下步骤：
[0006] (1)读入两幅待配准图像，将其中一幅固定不变的图像作为参考图像IR，另一幅进行变换的图像作为浮动图像IF，并按间距L＝40个像素，对浮动图像均匀采样，得到由m个元素组成的均匀网格；
[0007] (2)随机产生一组由量子比特位构成的个体，作为量子进化计算的初始迭代点，称为种群；
[0008] (3)将种群中每个个体转化为候选解，再与均匀网格叠加，得到每个个体对应的控制网格；
[0009] (4)按照每个个体对应的控制网格对浮动图像做三次B样条函数变换，得到变换后的浮动图像；
[0010] (5)计算变换后的浮动图像与参考图像之间的归一化互信息值NMI，找出本次迭代中最大的归一化互信息值及其对应的个体，并存储为本代最优解pbest；
[0011] (6)依据本代最优解，利用量子旋转门算子对本代种群中所有个体进行量子旋转门更新，使图像配准参数朝相似性更大的方向移动；
[0012] (7)利用量子重组算子对任意选择的本代种群中两个不同个体pa和pb进行量子重组操作，增加种群个体多样性，避免图像配准参数陷入局部极值；
[0013] (8)更新迭代次数g，如迭代次数达到最大值或变换后的浮动图像与参考图像之间的归一化互信息值不再增加，则结束配准，即本代最优解对应的控制网格为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准后的图像，否则，返回步骤(3)进行下一次迭代。
[0014] 本发明与现有的技术相比具有以下优点：
[0015] 1.本发明使用图像灰度信息作为配准特征，省去了基于图像特征的配准方法复杂的提取图像特征步骤，避免了因为提取图像特征算法对配准结果产生的影响，使配准过程更为简洁；
[0016] 2.本发明采用三次B样条函数变换作为图像变换方式，可对图像进行非刚性配准，相比于刚性配准方法，其适用性更为广泛；
[0017] 3.本发明采用量子旋转门算子和量子重组算子作为种群更新算子，对于配准参数较多的三次B样条图像变换，不易陷入局部极值，配准效果更好。
附图说明
[0018] 图1是本发明的总流程图；
[0019] 图2是本发明仿真实验的参考图像；
[0020] 图3是本发明仿真实验的浮动图像；
[0021] 图4是用现有的基于遗传算法对图2和图3进行配准的仿真结果图；
[0022] 图5是用本发明对图2和图3进行配准的仿真结果图。
具体实施方式
[0023] 参照图1，本发明的实现步骤如下：
[0024] 步骤1，读入图像，初始相关参数。
[0025] 读入两幅待配准图像，一幅固定不变，另一幅为要变换的图像；将其中固定不变的图像作为参考图像IR，另一幅进行变换的图像作为浮动图像IF，并按间距L＝40个像素的个数，对浮动图像进行均匀采样，得到由m个元素组成的均匀网格；
[0026] 设定最大迭代次数max_g＝100，并将当前迭代次数g置零。
[0027] 步骤2，随机产生量子进化计算的初始迭代点。
[0028] (2a)按照量子比特位的方式产生一个个体p，其形态表示为：
[0029]
[0030] 其中，θi∈[0,2π],i＝1,2,...m，m是均匀网格的元素个数，θi是随机生成的量子比特位角度；
[0031] (2b)重复步骤(2a)产生一组个体P＝[p1,p2,...pN]，组成个数为N的种群。
[0032] 步骤3，计算种群中每个个体对应的控制网格。
[0033] (3a)选择种群P＝[p1,p2,...pN]中的一个个体pn，n＝1,2,...N；
[0034] (3b)将个体pn分解为两个子个体pnc和pns：
[0035] pnc＝[cos(θn1) cos(θn2)... cos(θnm)]
[0036] pns＝[sin(θn1) sin(θn2)... sin(θnm)]，
[0037] 其中，pnc是余弦个体，pns是正弦个体；
[0038] (3c)将余弦个体pnc和正弦个体pns映射到配准参数解空间，得到具有m个分量的余弦解Snc和正弦解Sns，各分量表达式如下：
[0039]
[0040]
[0041] 其中，Snci和Snsi分别是Snc和Sns的第i个分量，ai和bi分别是第i个分量的取值下限和上限，i＝1,2,...m,由此得到的Snc和Sns是一个1×m的向量；
[0042] (3d)将上述Snc和Sns整合为与均匀网格相同大小的矩阵，并与均匀网格叠加得到Snc对应的控制网格Φnc和Sns对应的控制网格Φns；
[0043] (3e)重复步骤(3a)至(3d)，得到种群P＝[p1,p2,...pN]中N个个体对应的2N个控制网格Φ＝{{Φ1c,Φ1s},{Φ2c,Φ2s},...{ΦNc,ΦNs}}。
[0044] 步骤4，按照控制网格对浮动图像做三次B样条函数变换，得到变换后的浮动图像。
[0045] (4a)取控制网格Φ＝{{Φ1c,Φ1s},{Φ2c,Φ2s},...{ΦNc,ΦNs}}中的一组控制网格{Φnc,Φns}，n＝1,2,...N。
[0046] (4b)取浮动图像一个像素点，记为(x,y)，分别按照控制网格Φnc和Φns计算该像素点的位移值D(x,y)，即：
[0047]
[0048] 其中， 和 分别表示对x和y取整， 是控制网格Φnc和Φns中该像素点(x,y)周围4×4控制网格点的值，Bq(e)和Bl(f)是三次B样条基函数Bd(t)中的一种，d＝0,1,2,3，0≤t＜1，三次B样条基函数Bd(t)表示如下：
[0049]
[0050] (4c)重复步骤(4b)计算浮动图像每个像素点的位移值，得到Φnc对应的变换后的浮动图像Inc和Φns对应的变换后的浮动图像Ins。
[0051] (4d)重复步骤(4a)至(4c)，得到种群P＝[p1,p2,...pN]对应的变换后的浮动图像[0052] I＝{{I1c,I1s},{I2c,I2s},...{INc,INs}}。
[0053] 步骤5，计算变换后的浮动图像与参考图像之间的归一化互信息值，并寻找本代最优解。
[0054] (5a)在变换后的浮动图像I＝{{I1c,I1s},{I2c,I2s},...{INc,INs}}中取一组图像{Inc,Ins}，n＝1,2,...N，分别计算参考图像IR与Inc，Ins之间的归一化互信息值；
[0055] (5a1)计算参考图像IR与变换后的浮动图像Inc之间的归一化互信息值NMInc为：
[0056]
[0057] 其中，H(IR)是参考图像IR的熵，H(Inc)是变换后的浮动图像Inc的熵，H(IR,Inc)是参考图像IR与变换后的浮动图像Inc的联合熵；
[0058] (5a2)计算参考图像IR与变换后的浮动图像Ins之间的归一化互信息值NMIns为：
[0059]
[0060] 其中，H(Ins)是变换后的浮动图像Ins的熵，H(IR,Ins)是参考图像IR与变换后的浮动Ins的联合熵；
[0061] (5b)将NMInc与NMIns比较，如果NMInc≥NMIns，则个体pn对应的控制网格Φn是pn的余弦个体pnc对应的控制网格Φnc，即Φn＝Φnc，个体pn对应的变换后的浮动图像In是 pnc对应的变换后的浮动图像Inc，即In＝Inc，个体pn对应的归一化互信息值NMIn是pnc对应的归一化互信息值NMInc，即NMIn＝NMInc；否则：个体pn对应的控制网格Φn是pn的正弦个体pns对应的控制网格Φns，即Φn＝Φns，个体pn对应的变换后的浮动图像In是pns对应的变换后的浮动图像Ins，即In＝Ins，个体pn对应的归一化互信息值NMIn是pns对应的归一化互信息值NMIns，即NMIn＝NMIns；
[0062] (5c)重复步骤(5a)-(5b)得出种群P＝[p1,p2,...pN]对应的归一化互信息值NMI＝[NMI1,NMI2,...NMIN]；
[0063] (5d)找出NMI＝[NMI1,NMI2,...NMIN]中的最大值及其对应的个体，并将该个体存储为第g代的最优解pbest，该个体对应的控制网格储存为第g代的最优配准参数Φbest，该个体对应的归一化互信息值存储为最大归一化互信息NMIbest。
[0064] 步骤6，对本代种群中所有个体进行量子旋转门更新。
[0065] (6a)按照迭代次数g构造旋转角度的大小Δθ0：
[0066] Δθ0＝0.015π×(g/max_g)+0.15π×exp(-mod(g,100)/10)，
[0067] 其中，max_g是设定的最大迭代次数，mod(g,100)是取g除以100的余数，exp(-mod(g,100)/10)是取e的-mod(g,100)/10次方；
[0068] (6b)取出第g代种群中个体p；
[0069] (6c)取出个体p的第i个分量 ，如果该分量越过了最优解pbest的第i个分量，则 ，反之， 该分量对应的旋转矩阵 只由旋转角度 决定，即：
[0070]
[0071] 则第g+1代的个体p的第i个分量 为：
[0072]
[0073] (6d)重复步骤(6c)直至个体p的所有分量全部更新；
[0074] (6e)重复步骤(6b)-(6d)直至种群中所有个体全部更新。
[0075] 步骤7，对本代种群进行量子重组操作。
[0076] (7a)从本代种群中任意选择两个不同个体pa和pb；
[0077] (7b)随机选择一个比特位h，按如下公式进行重组算子：
[0078]
[0079]
[0080] 其中， 是个体pa的第h个分量，
[0081] 是个体pb的第h个分量，
[0082] 是重组后个体pa的第h个分量，
[0083] 是重组后个体pb的第h个分量。
[0084] 步骤8，更新迭代次数，判断是否达到结束条件。
[0085] 更新迭代次数g＝g+1，如果迭代次数达到最大值max_g或者变换后的浮动图像与参考图像之间的归一化互信息值不再增加，则结束配准，将本代最优解pbest对应的控制网格Φbest作为配准参数，按此参数对浮动图像做三次B样条函数变换，得到配准后的图像，输出配准后的图像，输出最大归一化互信息值NMIbest，否则，返回步骤(3)进行下一次迭代。
[0086] 本发明的效果可以通过下面的仿真实验进一步说明：
[0087] 1、实验条件与方法
[0088] 硬件平台为：Intel(R)Core(TM)i5-2450M@2.50GHz、3.91GBRAM.；
[0089] 软件平台为：MATLABR2012b；
[0090] 实验方法：分别用现有的基于遗传算法的方法和本发明仿真实现医学图像配准。
[0091] 2、仿真内容与结果
[0092] 实验所用图像为肝脏CT图像，如图2和图3所示。
[0093] 仿真一，用基于遗传算法的方法对图2和图3所示的实验图像进行配准，设置最大迭代次数为100，独立仿真50次，50次仿真所用的平均时间为216.0906367s，得到最佳配准图像和归一化互信息值平均值两个结果，其中最佳配准图像如图4所示，50次仿真输出的归一化互信息值的平均值为1.262681361。
[0094] 仿真二，用本发明对图2和图3所示的实验图像进行配准，最大迭代次数为100，独立仿真50次，50次仿真所用的平均时间为176.221776s，得到最佳配准图像和归一化互信息值平均值两个结果，其中最佳配准图像如图5所示，50次仿真输出的归一化互信息值的平均值为1.300652358。
[0095] 由上述仿真结果可知，本发明相较于基于遗传算法的方法可以有效地提高配准程度，不仅能获得了更好的配准结果，并且缩短了配准时间，是一种有效的医学图像配准方法。
