本发明涉及图像处理技术，尤其涉及一种海量脑图像数据集的三维非线性配准方法。
背景技术：
图像配准技术作为图像处理领域的一个重要研究方向，目前广泛地应用于遥感、自动导航、模式识别、计算机视觉、生物医学等领域。随着大数据时代的来临，海量三维图像数据集不断增大，往往能达到tb级。如何对这些海量三维数据集进行非线性配准是一个非常大的挑战。
在中国发明专利说明书cn104361590中，提出了一种控制点自适应分布的高分辨率遥感影像配准方法，该方法基于分块策略的自适应控制点提取方法，实现了多尺度j-image影像中的控制点提取，进而采用归一化的互信息测度nmi对控制点进行多尺度进行匹配。但是该方法所处理的数据为二维的遥感数据，适用于特定的应用领域。而且，该方法所描述的高分辨也仅仅是空间分辨率10m，尺寸为512×512像素，数据大小仅为mb级别，无法对tb级三维图像数据进行非线性配准。
在图像配准研究中广泛使用的工具也实现了三维数据集的非线性配准，例如insightsegmentationandregistrationtoolkit(itk)和advancednormalizationtools(ants)。虽然这些工具都提供了丰富的三维非线性配准的方法，但是在配准过程中它们将三维图像视为一体，一次性导入计算机中进行运算，这急剧加大了对内存的消耗。对于几十mb的三维图像数据，有些工具甚至会消耗数百gb的内存，更加难以对tb级的三维图像数据进行非线性配准。
综合上述专利技术和工具，目前三维图像非线性配准存在很多方法，也能解决各种领域的问题。但是随着大数据时代的来临，海量三维图像数据集不断增大，还没有一种能够处理tb级三维图像数据的非线性配准方法。
技术实现要素：
本发明的目的在于，解决上述技术问题，提供一种海量脑图像数据集的三维非线性配准方法及系统，对海量高分辨三维图像数据集进行快速变换，实现海量高分辨三维图像数据集的非线性配准。
本发明的目的通过如下技术方案来实现。
一种海量脑图像数据集的三维非线性配准方法，包括以下步骤：
s0.降采样步骤，对预设高分辨参考图像数据和高分辨待配准图像数据进行降采样，得到低分辨参考图像数据和低分辨待配准图像数据；
s1.低分辨图像配准步骤，其中包括：
s11.线性配准步骤，对低分辨参考图像数据和低分辨待配准图像数据进行线性配准，得到低分辨线性配准的变换矩阵，并计算其逆矩阵低分辨配准的逆矩阵，以及低分辨线性配准后结果；
s12.非线性配准步骤，对低分辨线性配准后结果和低分辨参考图像数据进行非线性配准，得到低分辨非线性配准的位移场，以及低分辨非线性配准后结果；
s2.高分辨图像快速变换步骤，其中包括：
s21.高分辨配准的变换矩阵计算步骤，基于低分辨线性配准的变换矩阵计算出高分辨配准的变换矩阵，以及其逆矩阵高分辨配准的逆矩阵；
s22.高分辨配准后空间范围计算步骤，利用步骤s21所得高分辨配准的逆矩阵，计算出高分辨配准后空间范围；
s23.分块步骤，对步骤s22所得高分辨配准后空间范围进行分块；
s24.高分辨非线性变换步骤，对步骤s23所得每一个分块进行并行计算，计算每一个像素点在原高分辨待配准图像数据的对应位置，最终得到高分辨图像非线性变换结果。
进一步，所述步骤s11中所述的线性配准，将变换映射描述为变换矩阵。
进一步，所述步骤s12中所述的非线性配准，将变换映射描述为像素偏移向量的位移场。
进一步，所述步骤s23中所述的分块，基于计算资源进行分块，并行计算分块中像素点的变换，所述计算资源最大能支持n个线程并行，以长度n进行分块。
进一步，所述步骤s24中所述的高分辨非线性变换，包括以下子步骤：
s241.估算分块对应变换前原高分辨待配准图像数据的范围；
s242.遍历分块中每个像素点，基于低分辨非线性配准的位移场，插值计算出该像素点对应的偏移向量，基于偏移向量和高分辨配准的逆矩阵，计算出该像素点在原高分辨待配准图像数据的对应位置，将其灰度值作为该像素点灰度值。
本发明还揭示了一种海量脑图像数据集的三维非线性配准系统，包括以下模块：
降采样模块，对预设高分辨参考图像数据和高分辨待配准图像数据进行降采样，得到低分辨参考图像数据和低分辨待配准图像数据；
低分辨图像配准模块，其中包括：
线性配准单元，对低分辨参考图像数据和低分辨待配准图像数据进行线性配准，得到低分辨线性配准的变换矩阵，并计算其逆矩阵低分辨配准的逆矩阵，以及低分辨线性配准后结果；
非线性配准单元，对低分辨线性配准后结果和低分辨参考图像数据进行非线性配准，得到低分辨非线性配准的位移场，以及低分辨非线性配准后结果；
高分辨图像快速变换模块，其中包括：
高分辨配准的变换矩阵计算单元，基于低分辨线性配准的变换矩阵计算出高分辨配准的变换矩阵，以及其逆矩阵高分辨配准的逆矩阵；
高分辨配准后空间范围计算单元，利用高分辨配准的变换矩阵计算单元所得高分辨配准的逆矩阵，计算出高分辨配准后空间范围；
分块单元，对高分辨配准后空间范围计算单元所得高分辨配准后空间范围进行分块；
高分辨非线性变换单元，对分块单元所得每一个分块进行并行计算，计算每一个像素点在原高分辨待配准图像数据的对应位置，最终得到高分辨图像非线性变换结果。
进一步，所述线性配准单元中所述的线性配准，将变换映射描述为变换矩阵。
进一步，所述非线性配准单元中所述的非线性配准，将变换映射描述为像素偏移向量的位移场。
进一步，所述分块单元中所述的分块，基于计算资源进行分块，并行计算分块中像素点的变换，所述计算资源最大能支持n个线程并行，以长度n进行分块。
进一步，所述高分辨非线性变换单元中所述的高分辨非线性变换，包括以下单元：
估算单元，估算分块对应变换前原高分辨待配准图像数据的范围；
变换单元，遍历分块中每个像素点，基于低分辨非线性配准的位移场，插值计算出该像素点对应的偏移向量，基于偏移向量和高分辨配准的逆矩阵，计算出该像素点在原高分辨待配准图像数据的对应位置，将其灰度值作为该像素点灰度值。
本发明具有的优点为：利用低分辨图像数据配准信息，对海量高分辨三维图像数据集进行变换，实现海量高分辨三维图像数据集的非线性配准。本发明适用于基于变换矩阵的所有线性配准和基于位移场的所有非线性配准，能够快速对tb级数据集进行非线性配准。
附图说明
图1是本发明优选实施方式的海量脑图像数据集的三维非线性配准方法的流程示意图；
图2是本发明优选实施方式中高分辨图像快速变换的示意图；
图3是本发明优选实施方式的海量脑图像数据集的三维非线性配准系统的结构示意图。
具体实施方式
下面结合附图对本发明作进一步的详细描述。图1示出了本发明优选实施方式提供的海量脑图像数据集的三维非线性配准方法的流程示意图。
一种海量脑图像数据集的三维非线性配准方法，包括以下步骤：
s0.降采样步骤，对预设高分辨参考图像数据和高分辨待配准图像数据进行降采样，得到低分辨参考图像数据和低分辨待配准图像数据。
其中，所述高分辨参考图像数据d1，高分辨待配准图像数据d2。对两个高分辨图像数据进行降采样，得到低分辨参考图像数据d3和低分辨待配准图像数据d4。降采样的过程要保证两点：低分辨图像数据三维各方向分辨率保持一致；两个低分辨图像数据d3、d4整体分辨率一致。
s1.低分辨图像配准步骤，其中包括：
s11.线性配准步骤，对低分辨参考图像数据和低分辨待配准图像数据进行线性配准，得到低分辨线性配准的变换矩阵，并计算其逆矩阵低分辨配准的逆矩阵，以及低分辨线性配准后结果；
其中，利用基于仿射变换的线性配准，对低分辨图像数据d3、d4进行线性配准，得到低分辨线性配准的变换矩阵m1，并计算其逆矩阵低分辨线性配准的逆矩阵m2，以及低分辨线性配准后结果d5。
其中，所述步骤s11中所述的线性配准，为能将变换映射描述为变换矩阵的所有线性配准方法，例如：平移、旋转、缩放等线性配准。
s12.非线性配准步骤，对低分辨线性配准后结果和低分辨参考图像数据进行非线性配准，得到低分辨非线性配准的位移场，以及低分辨非线性配准后结果；
其中，对所述低分辨线性配准后结果d5和低分辨参考图像数据d3进行基于位移场的非线性配准，得到低分辨非线性配准的位移场t1，以及低分辨非线性配准后结果。
其中，所述步骤s12中所述的非线性配准，为能将变换映射描述为像素偏移向量的位移场的所有非线性配准方法，例如：弹性模型、流体模型变换的非线性配准
s2.高分辨图像快速变换步骤，其中包括：
s21.高分辨配准的变换矩阵计算步骤，基于低分辨线性配准的变换矩阵计算出高分辨配准的变换矩阵，以及其逆矩阵高分辨配准的逆矩阵；
其中，基于低分辨线性配准的变换矩阵m1计算出高分辨配准的变换矩阵m3，以及其逆矩阵高分辨配准的逆矩阵m4，公式如下：
m3＝scalem×m3×scalem′
m4＝m3-1
其中scalem和scalem′分别是缩小矩阵和放大矩阵，缩小和放大的倍数与高分辨与低分辨倍数一致。
s22.高分辨配准后空间范围计算步骤，利用步骤s21所得高分辨配准的逆矩阵，计算出高分辨配准后空间范围；
其中，利用高分辨逆矩阵m4对高分辨待配准图像数据d2的8个顶点进行变换，取他们x、y、z方向的最大最小值作为高分辨配准后图像数据d6的高分辨配准后空间范围。
s23.分块步骤，对步骤s22所得高分辨配准后空间范围进行分块；
其中，基于计算资源进行分块，并行计算分块中像素点的变换。根据计算环境，对所述高分辨配准后空间范围进行z方向分块。若计算环境最大能支持n个线程并行，则以长度n进行分块，这使得能够并行计算同一段数据。若每一分块需载入内存数据超过内存，则分块长度减少一半，直到满足内存需求。
s24.高分辨非线性变换步骤，对步骤s23所得每一个分块进行并行计算，计算每一个像素点在原高分辨待配准图像数据的对应位置，最终得到高分辨图像非线性变换结果。
所述步骤s24中所述的高分辨非线性变换，包括以下子步骤：
s241.估算分块对应变换前原高分辨待配准图像数据的范围；
s242.遍历分块中每个像素点，基于低分辨非线性配准的位移场，插值计算出该像素点对应的偏移向量，基于偏移向量和高分辨配准的逆矩阵，计算出该像素点在原高分辨待配准图像数据的对应位置，将其灰度值作为该像素点灰度值。
其中，对于每一个分块首先计算其对应原高分辨待配准图像数据d2所在z方向范围。计算分块对应低分辨待配准图像数据d4的分块,对分块中每一个像素的位置，先加上低分辨非线性配准的位移场t1相应位置三个方向的偏移向量，然后再乘以低分辨线性配准的逆矩阵m2，得到其对应原低分辨待配准图像数据d4所在位置。求出分块所有像素点对应在原低分辨待配准图像数据d4的位置后，计算三个方向的最大最小值，得到原低分辨待配准图像数据d4的z方向范围，最后再计算出对应原高分辨待配准图像数据d2的z方向范围。
如图2所示。对于像素点p0(x0,y0,z0)，计算其在低分辨待配准图像数据d4中的坐标像素p1(x1,y1,z1),则其对应低分辨非线性配准的位移场t1的坐标为(x1,y1,z1)，由于x1,y1,z1不一定为整数，所以需要对低分辨非线性配准的位移场t1的对应位置进行三线性插值计算出像素点p1对应的偏移向量d1。根据计算得到的p1点偏移向量d1计算出p0点对应的偏移向量d0。将像素点p0坐标值与偏移向量相加d0，将得到的坐标p2再乘以高分辨配准的逆矩阵m4，最后得到坐标p3为p0逆变换的坐标，即像素点p0逆变换对应的原高分辨待配准图像数据d2中的点，将这个点的灰度值作为p0的灰度值。
所述的计算出分块中每个像素对应的灰度值，每段计算完成后以二维图像的形式写入磁盘。循环计算每一个分块，直到整个高分辨待配准图像数据计算完成。
由以上内容可以看出，本发明利用低分辨图像数据配准信息，对海量高分辨三维图像数据集进行变换，实现海量高分辨三维图像数据集的非线性配准。本发明适用于基于变换矩阵的所有线性配准和基于位移场的所有非线性配准，能够快速对tb级数据集进行非线性配准。
请参见图3所示，本发明一实施方式提供的海量脑图像数据集的三维非线性配准系统的结构示意图。
一种海量脑图像数据集的三维非线性配准系统，包括以下模块：
降采样模块，对预设高分辨参考图像数据和高分辨待配准图像数据进行降采样，得到低分辨参考图像数据和低分辨待配准图像数据；
低分辨图像配准模块，其中包括：
线性配准单元，对低分辨参考图像数据和低分辨待配准图像数据进行线性配准，得到低分辨线性配准的变换矩阵，并计算其逆矩阵低分辨配准的逆矩阵，以及低分辨线性配准后结果；
非线性配准单元，对低分辨线性配准后结果和低分辨参考图像数据进行非线性配准，得到低分辨非线性配准的位移场，以及低分辨非线性配准后结果；
高分辨图像快速变换模块，其中包括：
高分辨配准的变换矩阵计算单元，基于低分辨线性配准的变换矩阵计算出高分辨配准的变换矩阵，以及其逆矩阵高分辨配准的逆矩阵；
高分辨配准后空间范围计算单元，利用高分辨配准的变换矩阵计算单元所得高分辨配准的逆矩阵，计算出高分辨配准后空间范围；
分块单元，对高分辨配准后空间范围计算单元所得高分辨配准后空间范围进行分块；
高分辨非线性变换单元，对分块单元所得每一个分块进行并行计算，计算每一个像素点在原高分辨待配准图像数据的对应位置，最终得到高分辨图像非线性变换结果。
所述线性配准单元中所述的线性配准，将变换映射描述为变换矩阵。所述非线性配准单元中所述的非线性配准，将变换映射描述为像素偏移向量的位移场。
所述分块单元中，所述的分块，基于计算资源进行分块，并行计算分块中像素点的变换，所述计算资源最大能支持n个线程并行，以长度n进行分块。
所述高分辨非线性变换单元中所述的高分辨非线性变换，包括以下单元：估算单元，估算分块对应变换前原高分辨待配准图像数据的范围；变换单元，遍历分块中每个像素点，基于低分辨非线性配准的位移场，插值计算出该像素点对应的偏移向量，基于偏移向量和高分辨配准的逆矩阵，计算出该像素点在原高分辨待配准图像数据的对应位置，将其灰度值作为该像素点灰度值。
由以上内容可以看出，本发明利用低分辨图像数据配准信息，对海量高分辨三维图像数据集进行变换，实现海量高分辨三维图像数据集的非线性配准。本发明适用于基于变换矩阵的所有线性配准和基于位移场的所有非线性配准，能够快速对tb级数据集进行非线性配准。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。