本发明涉及组织探测领域，特别涉及一种加速图像重建的方法。
背景技术：
荧光分子层析成像技术在生物医学研究领域中占据重要地位，其广泛应用于癌症诊断，药物研发以及基因表达可视化的研究。
大多数生物组织均为高散射的三维浑浊介质，建立一个高精度和高效率的计算方法对于荧光分子层析成像的定量精度有着重要的意义。蒙特卡罗是一种基于随机抽样过程的离散统计理论方法。相比较于其它方法，蒙特卡罗方法可模拟任意几何形状，边界条件和光学参数下的光子输运过程。由于其广泛的适用性，蒙特卡罗方法被认为是模拟光子输运实际物理过程的最直接，最有效和最可信的方法。因而，它成为了评价其它特定应用方法的金标准，如扩散近似方法。然而，高精度的蒙特卡罗模拟是耗时的。光在介质中传输时其强度会随着传输距离增大而呈指数性衰减，介质尺寸越大，为了达到同样的高计算精度，蒙特卡罗仿真所需要的时间也是呈指数性的增长。巨大的计算量往往对计算机造成了较大的负担。
蒙特卡罗方法的基本特征是对随机性问题进行仿真，能有效解决随机性问题，甚至对许多确定性方法所难以解决的随机性问题都能较方便地解决。在生物组织中光子传输问题的研究中，蒙特卡罗也表示出它相对于其它方法(如扩散近似)的优势；严密且灵活，不受组织的光学特性限制，易于使用，其模拟精度高，被该领域研究者奉为金标准。
几个基于历史路径的蒙特卡罗模拟方法已经被提出。例如：一种基于时间分辨荧光的高效蒙特卡罗方法，这种方法通过保存大量光子路径信息，建立了透射光子权重与组织光学参数的解析关系；一种基于计算时间门控的荧光雅克比矩阵方法的高效蒙特卡罗方法，这是一种能够对两个寿命不同的荧光团同时成像的方法，在这种方法中，从激发源到探测器的背景权重矩阵可由光子路径和荧光吸收系数计算得出；此外，还有一种直接计算荧光的解耦合荧光蒙特卡罗模型，通过解耦合激发和发射转化过程，将光子路径及相应组织光学参数与逸出荧光光子联系在一起，此方法适用于任意浑浊介质的任意荧光分布。以上的几种方法都保留了大量的光子路径信息，在这些方法中，光子的路径信息是蒙特卡罗模拟和图像重建过程中非常重要的一部分。
图像重建过程即是寻找最优解使目标函数值最小，图像重建过程会涉及到矩阵的计算，以及一些迭代算法，例如共轭梯度法。随着模拟的光子数目增加，模拟源的数目不断增加，计算机模拟的时间也在成倍增加，产生的数据量越来越大，数据的计算、传输等方面所需的时间也越来越长，现有的基于蒙特卡罗模型的图像重建方法的速度已不能满足需求。
技术实现要素：
本发明的目的是提供一种加速图像重建的方法，以提高图像重建的速度。
为实现上述目的，本发明提供了如下方案：
一种加速图像重建的方法，所述方法包括如下步骤：
预设生物组织的初始尺寸、初始光学参数，并确定激发光光源的初始化信息，并将所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息存储至计算机主节点的显卡中的常数寄存器；
根据所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息建立生物组织模型，并确定所述生物组织模型的体素索引，并将所述体素索引存储至计算机主节点的纹理寄存器；
计算机主节点对所述生物组织模型进行蒙特卡罗模拟，计算每个激发光光源对应的所述生物组织模型出射的荧光强度，并根据每个所述荧光强度计算荧光强度与所述体素索引的荧光系数的雅可比矩阵；
计算机主节点将所述雅克比矩阵分割成多个子雅克比矩阵，并将多个所述子雅克比矩阵分配到多个计算机节点；
每个所述计算机节点根据其对应的子雅克比矩阵，计算出其对应的▽f(x)的分量和Hsk的分量，并将▽f(x)的分量和Hsk的分量传递给计算机主节点；所述计算机主节点对多个▽f(x)的分量进行合成并对多个Hsk的分量进行合成，获得重建后的图像，其中，f(x)是待求解函数，H是Hesse矩阵，sk是第k轮迭代的搜索方向。
可选的，所述将所述体素索引存储至计算机主节点的纹理寄存器，具体包括：
对所述体素索引组成的数组进行三维纹理绑定，获得绑定后的体素索引；
将所述绑定后的体素索引存储至纹理寄存器。
可选的，每个所述计算机节点根据其对应的子雅克比矩阵，计算出其对应的▽f(x)的分量和Hsk的分量，并将▽f(x)的分量和Hsk的分量传递给计算机主节点；所述计算机主节点对多个▽f(x)的分量进行合成并对多个Hsk的分量进行合成，获得重建后的图像，具体包括：
每个计算机节点根据其对应的雅克比矩阵子矩阵和公式▽f(x)＝Hx+b＝(JTJ+λ)x-JTΔD计算出▽f(x)在该计算机节点的一个分量其中，n为计算机节点的编号，f(x)是待求解函数，x代表体素索引按索引顺序组成的向量，H是Hesse矩阵，H＝JTJ+λ，λ是正则化参数，b为探测器的探测值所组成的向量；xk是共轭梯度法的第k轮迭代结果，它代表组织模型中各体素中参数值按索引顺序组成的向量；Jn,(n＝0,1,···,i)是分配到n个计算机节点的子雅克比矩阵，对应Jn,(n＝0,1,···,i)的转置矩阵；ΔD为模拟探测值与实际探测值的差；
所有计算机节点将▽f(x)的一个分量传到主计算机节点并求和，得到▽f(x)；
主计算机节点判断▽f(x)是否收敛，若是，则停止算法，根据▽f(x)获得重建后的图像；若否，则计算共轭梯度法的搜索方向sk，并将结果广播到各计算机节点；
各计算机节点分别根据公式Hsk＝JTJsk计算Hsk的分量；其中，sk是第k轮迭代的搜索方向，是一个向量；
各计算机节点把Hsk的分量发送到主计算机节点并在主计算机节点求和，得到Hsk；
主计算机节点根据Hsk，采用共轭梯度法，求得mk,gk+1,βk；gk为最优判别参数，gk＝▽f(xk)；mk为最优搜索步长，βk为搜索方向更新参数，
主计算机节点根据mk和sk，利用共轭梯度法公式xk+1＝xk+mksk来计算xk+1,并将xk+1传递给各计算机节点，返回每个计算机节点根据其对应的雅克比矩阵子矩阵和公式▽f(x)＝Hx+b＝(JTJ+λ)x-JTΔD计算出▽f(x)的一个分量
可选的，将所述雅克比矩阵分割成多个子雅克比矩阵，并将多个所述子雅克比矩阵分配到多个计算机节点，之前还包括：
计算机主节点的CPU的一个核开辟第一进程，将计算机主节点的显卡计算得到的荧光强度从显存复制到内存，复制完毕后，CPU另一个核开辟第二进程，将荧光强度写入硬盘，计算机主节点的显卡继续计算下一个激发光光源对应的所述生物组织模型出射的荧光强度。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种加速图像重建的方法，所述方法包括：首先，预设生物组织的初始尺寸、初始光学参数，并确定激发光光源的初始化信息，并将所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息存储至GPU中的常数寄存器；建立生物组织模型，并其体素索引，将所述体素索引存储至纹理寄存器；然后，计算机主节点对所述生物组织模型进行蒙特卡罗模拟，计算每个激发光光源对应的所述生物组织模型出射的荧光强度，并计算雅可比矩阵；最后，将雅克比矩阵分割成多个子雅克比矩阵，并将多个子雅克比矩阵分配到多个计算机节点；由每个所述计算机节点根据其对应的子雅克比矩阵，计算出其对应的▽f(x)的分量和Hsk的分量，并将▽f(x)的分量和Hsk的分量传递给计算机主节点；由计算机主节点对多个▽f(x)的分量进行合成并对多个Hsk的分量进行合成，获得重建后的图像。本发明按照用途对数据进行分类，不同用途的数据采取不同的存储方式进行存储，以提高图像重建中数据的读取和拷贝速度，在图像重建过程中,将雅可比矩阵进行分割，并把所有子雅克比矩阵分配到各个计算机节点，每个计算机节点根据接收到的子雅克比矩阵计算出最终结果的一个分量，然后各节点将计算得到的分量传递到主计算机节点合并。从而将大规模矩阵的传递转化成向量的传递，提高了数据的传递效率，进而提高了图像重建的速度。
本发明在计算雅克比矩阵的过程中让主计算机节点的CPU的一个核开辟进程，调用显存完成数据的拷贝和计算，同时，CPU的另一个核也开辟进程，把上一轮循环中已经复制到内存的荧光强度和组织中的路径信息写入硬盘，等到两个CPU进程均完成后，一个模拟周期结束。这种进程时间流的优化实现了数据写入和数据拷贝的并行，提高了整个模拟过程的并行度，从而缩短了模拟时间。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明提供的一种加速图像重建的方法的流程图；
图2为本发明提供的一种加速图像重建的方法的获得重建后的图像的流程图；
图3为本发明提供的一种加速图像重建的方法的获得重建后的图像的进程优化的加速原理图。
图4为本发明提供的一种加速图像重建的方法的具体实施例的蒙特卡罗模型(a)及重建结果图(b)。
具体实施方式
本发明的目的是提供一种加速图像重建的方法，以提高图像重建的速度。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。
如图1所示，一种加速图像重建的方法，所述方法包括如下步骤：
步骤101，预设生物组织的初始尺寸、初始光学参数，并确定激发光光源的初始化信息，并将所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息存储至主计算机节点的GPU(Graphics Processing Unit，显示芯片或显卡)中的常数寄存器。
步骤102，根据所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息建立生物组织模型，并确定所述生物组织模型的体素索引，并将所述体素索引存储至计算机主节点的纹理寄存器；具体的，对所述体素索引组成的数组进行三维纹理绑定，获得绑定后的体素索引；将所述绑定后的体素索引存储至纹理寄存器。
本发明在调用GPU进行雅可比矩阵的计算时，把只需进行读取操作的数据，即所述初始尺寸、所述初始光学参数和所述激发光光源的初始化信息用常数寄存器来存储，从而加快只读数据的读取速度；将生物组织模型中各体素索引的参数组成的数组进行纹理绑定，采用纹理寄存器来存储，相对于将数据存入全局存储器加快组织模型的索引速度。采用页锁定内存保存模拟得到的路径信息，页锁定内存位于CPU(CentralProcessing Unit/Processor，中央处理器)，用映像存储器(mappedmemory)处理页锁定内存中的数据，使这块存储器既有内存地址，又有显存地址，GPU和CPU都可以直接访问这块存储器，加快路径信息的访问速度。
步骤103，计算机主节点对所述生物组织模型进行蒙特卡罗模拟，计算每个激发光光源对应的所述生物组织模型出射的荧光强度，并根据每个所述荧光强度计算荧光强度与所述体素索引的荧光系数的雅可比矩阵；具体的，计算机主节点的CPU的一个核开辟第一进程，将计算机主节点的显卡计算得到的荧光强度从显存复制到内存，复制完毕后，CPU另一个核开辟第二进程，将荧光强度写入硬盘，计算机主节点的显卡继续计算下一个激发光光源对应的所述生物组织模型出射的荧光强度。
本发明在计算雅克比矩阵的过程中让主计算机节点的CPU的一个核开辟进程，调用显存完成数据的拷贝和计算，同时，CPU的另一个核也开辟进程，把上一轮循环中已经复制到内存的荧光强度和组织中的路径信息写入硬盘，等到两个CPU进程均完成后，一个模拟周期结束。这种进程时间流的优化实现了数据写入和数据拷贝的并行，提高了整个模拟过程的并行度，从而缩短了模拟时间。
步骤104，计算机主节点将所述雅克比矩阵分割成多个子雅克比矩阵，并将多个所述子雅克比矩阵分配到多个计算机节点；
步骤105，每个所述计算机节点根据其对应的子雅克比矩阵，计算出其对应的▽f(x)的分量和Hsk的分量，并将▽f(x)的分量和Hsk的分量传递给计算机主节点；所述计算机主节点对多个▽f(x)的分量进行合成并对多个Hsk的分量进行合成，获得重建后的图像，其中，f(x)是待求解函数，H是Hesse矩阵，sk是第k轮迭代的搜索方向。具体包括：
每个计算机节点根据其对应的雅克比矩阵子矩阵和公式▽f(x)＝Hx+b＝(JTJ+λ)x-JTΔD计算出▽f(x)在该计算机节点的一个分量其中，n为计算机节点的编号，f(x)是待求解函数，x代表体素索引按索引顺序组成的向量，H是Hesse矩阵，H＝JTJ+λ，λ是正则化参数，b为探测器的探测值所组成的向量；xk是共轭梯度法的第k轮迭代结果，它代表组织模型中各体素中参数值按索引顺序组成的向量；Jn,(n＝0,1,···,i)是分配到n个计算机节点的子雅克比矩阵，对应Jn,(n＝0,1,···,i)的转置矩阵；ΔD为模拟探测值与实际探测值的差；
所有计算机节点将▽f(x)的一个分量传到主计算机节点并求和，得到▽f(x)；
主计算机节点判断▽f(x)是否收敛，若是，则停止算法，根据▽f(x)获得重建后的图像；若否，则计算共轭梯度法的搜索方向sk，并将结果广播到各计算机节点；
各计算机节点分别根据公式Hsk＝JTJsk计算Hsk的分量；其中，sk是第k轮迭代的搜索方向，是一个向量；
各计算机节点把Hsk的分量发送到主计算机节点并在主计算机节点求和，得到Hsk；
主计算机节点根据Hsk，采用共轭梯度法，求得mk,gk+1,βk；gk为最优判别参数，gk＝▽f(xk)；mk为最优搜索步长，βk为搜索方向更新参数
主计算机节点根据mk和sk，利用共轭梯度法公式xk+1＝xk+mksk来计算xk+1,并将xk+1传递给各计算机节点，返回每个计算机节点根据其对应的雅克比矩阵子矩阵和公式▽f(x)＝Hx+b＝(JTJ+λ)x-JTΔD计算出▽f(x)的一个分量
具体的，如图2所示，图中，Noden(n＝0,1,···,i)代表n号节点，其中，Node 0代表主节点，其他的节点都是子节点；
2.1、根据源的不同对雅克比矩阵进行分割并分配到各个节点，每个节点根据雅克比矩阵的子矩阵和公式计算出▽f(x)的一个分量，即n为节点编号；其中，▽f(x)＝Hx+b＝(JTJ+λ)x-JTΔD，f(x)是待求解函数，x代表组织模型中各体素中参数值按索引顺序组成的向量，H是Hesse矩阵，H＝JTJ+λ，λ是正则化参数，b为探测器的探测值所组成的向量；xk是共轭梯度法的第k轮迭代结果，它代表组织模型中各体素中参数值按索引顺序组成的向量；Jn,(n＝0,1,···,i)是分配到n号节点雅克比矩阵的子矩阵，对应Jn,(n＝0,1,···,i)的转置矩阵；ΔD为模拟探测值与实际探测值的差；
2.2、所有子节点将计算结果传到主节点并求和；
2.3、主节点判断是否达到收敛条件，若达到则停止算法；若未达到则继续(2.4)；
2.4、主节点计算共轭梯度法的搜索方向，并将结果广播到子节点；
2.5、所有节点根据Hsk＝JTJsk计算Hsk的一个分量；sk是第k轮迭代的搜索方向，是一个向量；
2.6、所有子节点把计算结果发送到主节点并在主节点求和，得到Hsk；
2.7、根据Hsk的结果以及共轭梯度法的公式，在主节点求得mk,gk+1,βk；gk＝▽f(xk)；代表最优搜索步长；
2.8、主节点根据共轭梯度法公式xk+1＝xk+mksk来计算xk+1,返回(2.3)；
(3)在图像重建的过程中，本发明让CPU的一个核开辟进程，并调用GPU完成数据的复制和计算，同时，CPU的另一个核也开辟进程，把上一轮循环中已经复制到内存的光子的状态信息和组织中的路径信息写入硬盘。等到两个CPU进程均完成后，一个模拟周期结束，如图3所示。
本发明还提供了一个具体实施例，来阐述本发明的方法：
通过仿真实验来测试实现本发明的方法。所设计的蒙特卡罗模型如图4(a)所示，该蒙特卡罗模型为匀质模型，其中包含三根荧光系数分别为0.2cm-1，0.4cm-1，0.6cm-1的荧光棒。激发光和荧光波段下的光学参数是一致的，量子效率设置为1，在模拟中所使用的其他光学参数如表1所示。该模型的尺寸为2.1cm×2.1cm×3.0cm，单个体素的尺寸为0.1cm，模型共含有13230个体素。绿色的区域为重建区域，位于圆柱中间14mm的范围内，共包含4270个体素。在重建过程中，本实施例模拟120个源，其覆盖360度角范围，处于圆柱中间7.5mm的范围内，均匀分布在6层，每层20个源；我们选取240个探测器，其覆盖360度角范围，位于圆柱中间20mm的范围内，均匀分布在10层，每层24个探测器。每个源所模拟的光子数均为2000000个。
表1用于模拟的蒙特卡罗模型的光学参数
其中，n为折射率，μs为散射系数，μa为吸收系数，g为各向异性因子，μaf为荧光团吸收系数。
根据各节点GPU设备的数目及性能和各节点的硬盘读写性能，三个节点上分配的所需模拟的光源数目分别为25，50，45。图4(b)为三维重建结果，表2显示了在并行架构上图像重建时间。
表2图像重建的时间
根据本发明提供的具体实施例，本发明公开了以下技术效果：
由于图像重建过程需要预设很多只读的数据。而GPU读取不同的存储器上的数据的读取速度是不一样的，本发明把光子的初始化信息(位置和方向)，组织的光学参数，模型的尺寸这些信息量较小的信息存入常数寄存器。GPU的纹理寄存器能够加快索引速度，所以采用每个体素的参数进行纹理绑定，采用纹理寄存器存储，GPU读取常数寄存器和纹理寄存器的读取速度快于全局寄存器。因为CPU计算的数据位于计算机内存，GPU计算的数据位于显存，位于内存(内存条)的数据和位于显存(显卡)的数据并不是公共的，如果CPU需要读取显存的数据，则需要先把数据从显存转移到内存再读取，GPU读取内存的数据也一样。而采用映像存储器处理页锁定内存，可以让这块存储器既有内存地址，又有显存地址，即这块存储器中的数据既可以直接被GPU使用，也可以直接被CPU使用，这样，将既会被GPU使用，也会被CPU使用的数据采用页锁定内存的方式进行存储，加快读取速度。
本发明的方法在所述根据所述实际权重和所述雅克比矩阵，重建所述生物组织的图像的过程中还将雅克比矩阵根据源进行分割，每个源对应一个子矩阵，并把所有子矩阵分配到各个节点，每个节点根据接收到的子矩阵计算出最终结果的一个分量，然后各节点将计算得到的分量传递到主节点合并。先将雅克比矩阵分块，实现并行计算，这显然可以加快计算速度。但是当每个子节点计算完成后需要把结果传递到主节点，在每一轮迭代中，雅克比矩阵会被更新，如果某一节点计算的矩阵块是m×n，那么，计算完毕后，会得到一个新的m×n矩阵块，以及得到一个维度为m的向量，如果传递回主节点的是更新后的矩阵，这需要传递的数据量为m×n，如果传递的是计算得到的向量，那么需要传递的数据量为m，实现通过减小数据传递量的方式进行加速。
如图3所示，本发明还对进程时间进行了优化，由于本算法需要使用GPU进行计算，所以需要有CPU将数据从内存写入显存，当GPU计算完成后，会将计算结果传递到内存，再由CPU开辟进程，将内存中的结果写入硬盘。如果用图3(a)中的时间流，则GPU会出现等待的时间，如果用图3(b)的时间流，使用CPU的另外一个核开辟进程将内存的数据写入硬盘，则GPU不会有等待时间，GPU会把所有的时间花费在计算和输出计算结果上。(由于本发明的计算结果一般有几十个G，所以把数据从内存写入硬盘会花费较长时间，这个时间在整个重建过程中占着相当大的比例)。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。