CN115619890A - 基于并行随机迭代求解线性方程组的断层成像方法及系统 - Google Patents

基于并行随机迭代求解线性方程组的断层成像方法及系统 Download PDF

Info

Publication number
CN115619890A
CN115619890A CN202211545390.1A CN202211545390A CN115619890A CN 115619890 A CN115619890 A CN 115619890A CN 202211545390 A CN202211545390 A CN 202211545390A CN 115619890 A CN115619890 A CN 115619890A
Authority
CN
China
Prior art keywords
solution
projection
processor
equation set
approximate solution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211545390.1A
Other languages
English (en)
Other versions
CN115619890B (zh
Inventor
田敏
王俊洁
杜伟
王英龙
杨美红
潘景山
张赞军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Computer Science Center National Super Computing Center in Jinan
Original Assignee
Shandong Computer Science Center National Super Computing Center in Jinan
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong Computer Science Center National Super Computing Center in Jinan filed Critical Shandong Computer Science Center National Super Computing Center in Jinan
Priority to CN202211545390.1A priority Critical patent/CN115619890B/zh
Publication of CN115619890A publication Critical patent/CN115619890A/zh
Application granted granted Critical
Publication of CN115619890B publication Critical patent/CN115619890B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明提出了基于并行随机迭代求解线性方程组的断层成像方法及系统，涉及计算机断层成像技术领域,建立以像素点吸收系数为未知数的线性方程组；将线性方程组的求解任务按行划分到各处理器上；迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；本发明使用多个处理器并行计算，有效地解决现有方法存储和计算的瓶颈，高效求解计算机断层成像中高精度扫描所产生的高维方程组，减少求解所消耗的内存空间和求解时间，使得计算机断层成像可以产生更高精度的图像。

Description

基于并行随机迭代求解线性方程组的断层成像方法及系统
技术领域
本发明属于计算机断层成像技术领域，尤其涉及基于并行随机迭代求解线性方程组的断层成像方法及系统。
背景技术
本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。
计算机断层成像是一种使用X射线获取物体内部截面图像的成像技术，被广泛应用于物理学、生物学、考古学、地球物理学和其他科学领域；首先使用扫描仪从不同角度对物体内部的截面进行多次X射线扫描，使用探测器拾取扫描信号并传输到计算机，经处理产生线性方程组，然后通过迭代方法求解该线性方程组，从而产生截面图像。
计算机断层成像扫描产生的线性方程组具有方程数量远远大于未知数数量的特点，因此，通常使用随机Kaczmarz方法求解该线性方程组；但是，随着扫描精度的提高，产生的线性方程组维数越来越大，高达千万维甚至以上；高维度线性方程组占用的存储空间巨大，使用单个处理器求解消耗的时间长，甚至需要花费数天时间，难以求解高精度和超高精度扫描产生的巨大方程组，限制了计算机断层成像的发展。
现有的随机Kaczmarz方法使用单个处理器进行存储和计算，单个处理器的内存难以完整存储整个线性方程组，且单个处理器的计算速度十分有限，使得求解所需要的时间成本非常高，限制了计算机断层成像技术的精度，迫切需要借助并行计算技术来缩短计算时间。
对于并行计算的可行性，随机Kaczmarz方法具有收敛速度慢、单次迭代计算量小、两次迭代间依赖性强的特点，这使得该方法迭代次数非常多、单次迭代计算成本非常小、每次迭代必须依赖前一次迭代的计算结果，该方法大量的迭代次数和非常小的单次迭代成本使得对每一次迭代步骤进行并行计算因计算量太小而无法有效提高计算效率，反而因并行计算所带来的通信开销使得计算效率大大降低；而每次迭代必须依赖前一次迭代的计算结果又使得同时计算不同的迭代步骤不可行，导致该方法难以利用多个处理器共同存储、同时并行计算，几乎没有并行化的机会，无法有效利用计算机强大的计算能力。
也就是说，随机Kaczmarz方法具有收敛速度慢、单次迭代计算量小、两次迭代间依赖性强的特点，使用并行计算，导致非常高的通信量和非常低的浮点运算密度，产生非常低的计算效率和高额时间成本，从而使得在多个处理器上并行求解方程组具有非常大的困难。
因此，如何借助多个处理器实现高效高精度快速断层成像，是一项亟待解决的问题。
发明内容
为克服上述现有技术的不足，本发明提供了基于并行随机迭代求解线性方程组的断层成像方法及系统，使用多个处理器并行计算，有效地解决现有方法存储和计算的瓶颈，高效求解计算机断层成像中高精度扫描所产生的高维方程组，减少求解所消耗的内存空间和求解时间，使得计算机断层成像可以产生更高精度的图像。
为实现上述目的，本发明的一个或多个实施例提供了如下技术方案：
本发明第一方面提供了基于并行随机迭代求解线性方程组的断层成像方法；
基于并行随机迭代求解线性方程组的断层成像方法，包括：
使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
将线性方程组的求解任务按行划分到各处理器上；
以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
其中，各处理器通过贪心采样方法计算各自对当前近似解的最优投影；对各处理器的最优投影进行平均，得到新的近似解。
进一步的，所述线性方程组为
Figure 17359DEST_PATH_IMAGE001
， A是射线穿过断层各个像素点的路径长 度，为系数矩阵，
Figure 301710DEST_PATH_IMAGE002
是断层各个像素点的吸收系数，为待求解的向量， b 是射线的出射强 度，为常数项向量。
进一步的，所述按行划分是对系数矩阵和常数项向量进行按行划分，将划分得到的系数矩阵分量及对应的常数项分量，组成子方程组，分配到不同的处理器上。
进一步的，所述贪心采样方法，具体为：
（1）设置初始投影为当前近似解，投影次数为0；
（2）基于子方程组和当前投影，计算投影距离范数最大的超平面；
（3）将当前投影投影到计算出来的超平面上，得到新的投影；
（4）投影次数加1，返回步骤（2），计算新的超平面和投影，直到投影次数满足要求。
进一步的，所述投影次数满足要求，是投影次数等于子方程组的行数。
进一步的，所述计算投影距离范数最大的超平面，具体为：
从子方程组中随机选择多个超平面，计算每个超平面的投影距离范数，比较投影距离范数的大小，从中计算投影距离范数最大的超平面。
进一步的，所述对各处理器的最优投影进行平均，是对各处理器产生的最优投影的各个元素分别取算数平均，产生新的近似解。
本发明第二方面提供了基于并行随机迭代求解线性方程组的断层成像系统。
断层成像中线性方程组的并行随机迭代求解系统，包括方程建立模块、任务划分模块、并行求解模块和图像绘制模块：
方程建立模块，被配置为：使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
任务划分模块，被配置为：将线性方程组的求解任务按行划分到各处理器上；
并行求解模块，被配置为：以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
图像绘制模块，被配置为：将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
其中，各处理器通过贪心采样方法计算各自对当前近似解的最优投影；对各处理器的最优投影进行平均，得到新的近似解。
本发明第三方面提供了计算机可读存储介质，其上存储有程序，该程序被处理器执行时实现如本发明第一方面所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
本发明第四方面提供了电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本发明第一方面所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
以上一个或多个技术方案存在以下有益效果：
本发明提供的断层成像方法使得计算机断层扫描可以利用多个处理器共同存储、并行计算求解线性方程组，突破了高精度扫描产生的高维方程组无法存储、无法求解的难题，提高了扫描精度和速度。
本发明采用并行存储策略减少了单个处理器内存空间使用量，从而可以求解更高维度的方程组；若使用P个处理器并行存储，则求解精度与单处理器相比可以提高P倍。
本发明采用并行求解策略减少了线性方程组的求解时间，从而提高了图像生成的速度；若使用P个处理器并行求解，则求解速度与单处理器相比有明显提高，理想条件下可以提高P倍，在实际应用中，随着处理器个数的增加，该倍数逐渐减少。
各处理器并行计算近似解的多次投影，每个处理器在一次迭代中进行多次投影，与现有方法一次迭代只进行一次投影不同，解决了现有方法因大量的迭代次数和非常小的单次迭代成本使得对每一次迭代步骤进行并行计算因计算量太小而无法有效提高计算效率，反而因并行计算所带来的通信开销使得计算效率大大降低的问题，从而使求解过程可以利用计算机多处理器强大的计算能力。
本发明的贪心采样过程结合了随机选择策略和贪心选择策略；随机选择策略可能会产生糟糕的投影，从而使收敛速度过慢；而贪心选择策略虽然会产生最优的投影，但代价高昂；本发明的贪心采样策略构造了一个低维随机子空间，并在随机子空间中产生最优投影，贪心选择策略的成本在低维空间中非常小；因此，可以以很少的计算成本产生相对最优的投影，解决了现有方法收敛速度慢的问题。
本发明通过取各处理器计算的投影的平均值的方式来给出新的近似解，与现有方法在投影之后立即进行迭代，把投影作为新的近似解不同，通过把现有方法中投影和迭代的过程分离，消除了每次迭代必须依赖前一次迭代的计算结果的问题。
本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1为第一个实施例的方法流程图。
图2为第一个实施例迭代求解线性方程组的流程图。
图3为第一个实施例中使用100000×100维度的线性方程组进行测试得到的加速比示意图。
图4为第一个实施例中使用500000×500维度的线性方程组进行测试得到的加速比示意图。
图5为第一个实施例中使用250000×100维度的线性方程组进行测试得到的加速比示意图。
图6为第一个实施例中使用1250000×500维度的线性方程组进行测试得到的加速比示意图。
图7为第二个实施例的系统结构图。
具体实施方式
下面结合附图与实施例对本发明作进一步说明。
应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明；除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式；如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
实施例一
本实施例公开了基于并行随机迭代求解线性方程组的断层成像方法，如图1所示，包括：
步骤S1、使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
具体的，使用扫描仪对物体内部断层进行X射线扫描，断层是一个截面，X射线通过 物体内部断层时会发生吸收，不同角度入射的X射线在物体内部断层中穿过的路径不同，吸 收的强弱程度也不同；通过入射角度，确定X射线在物体内部断层中穿过的路径，进而得到 路径长度；使用探测器拾取X射线的出射强度，根据朗伯-比尔定律，通过物体内部断层的X 射线的吸收过程，用线性方程组
Figure 651920DEST_PATH_IMAGE001
表示。
其中，A是射线穿过断层各个像素点的路径长度，为系数矩阵，
Figure 414339DEST_PATH_IMAGE002
是断层各个像素 点的吸收系数，为待求解的向量， b 是射线的出射强度，为常数项向量。系数矩阵A的列，代 表X射线穿过物体内部断层中一个像素点的路径长度。
步骤S2、将线性方程组的求解任务按行划分到各处理器上；
具体的，求解计算机断层成像扫描产生的线性方程组
Figure 32140DEST_PATH_IMAGE001
，其中，系数矩阵
Figure 487392DEST_PATH_IMAGE003
是
Figure 387215DEST_PATH_IMAGE004
的矩阵，常数项向量
Figure 891009DEST_PATH_IMAGE005
是
Figure 927098DEST_PATH_IMAGE006
的向量，系数矩阵
Figure 287672DEST_PATH_IMAGE003
和常数项向量
Figure 674791DEST_PATH_IMAGE005
的行数相同， 用
Figure 919959DEST_PATH_IMAGE007
表示，
Figure 76134DEST_PATH_IMAGE008
是系数矩阵的列数，即物体内部断层中像素点的个数，使用
Figure 332844DEST_PATH_IMAGE009
个处理器的线性 方程组并行随机迭代求解方法，
Figure 941680DEST_PATH_IMAGE009
是处理器个数，为正数；
如图2所示，对系数矩阵
Figure 52855DEST_PATH_IMAGE003
和常数项向量
Figure 63537DEST_PATH_IMAGE005
进行按行分块，分为
Figure 969176DEST_PATH_IMAGE009
块，划分成
Figure 65308DEST_PATH_IMAGE010
，共
Figure 449016DEST_PATH_IMAGE009
个子方程组,
Figure 579783DEST_PATH_IMAGE011
是分配到第
Figure 718640DEST_PATH_IMAGE012
个处理上的
Figure 738287DEST_PATH_IMAGE013
的系数矩阵分量，
Figure 988002DEST_PATH_IMAGE014
是分配到第
Figure 910959DEST_PATH_IMAGE012
个处理上的
Figure 220718DEST_PATH_IMAGE015
的常数项分量，系数矩阵分量
Figure 229125DEST_PATH_IMAGE011
和常数项分量
Figure 282531DEST_PATH_IMAGE014
的行数相同，用
Figure 59995DEST_PATH_IMAGE016
表示，
Figure 540655DEST_PATH_IMAGE008
是系数矩阵分量
Figure 269314DEST_PATH_IMAGE011
的列数，
Figure 126411DEST_PATH_IMAGE012
是处理器的序 号，
Figure 758381DEST_PATH_IMAGE017
；各个处理器存储经过划分的线性方程组，每个处理器的内存仅存储其 中的一块，即存储子方程组
Figure 409942DEST_PATH_IMAGE018
。
步骤S3、以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
具体的，如图2所示，各处理器将初始解
Figure 455259DEST_PATH_IMAGE019
设置为经验解，若没有经验解则设置为 0 向量，迭代次数
Figure 788151DEST_PATH_IMAGE020
设置为 0，
Figure 336944DEST_PATH_IMAGE020
为非负整数，准备开始迭代求解。
迭代求解的具体步骤为：
步骤S301：设置当前近似解
Figure 159407DEST_PATH_IMAGE021
为初始解
Figure 426440DEST_PATH_IMAGE019
，迭代次数k为0；
步骤S302：并行计算当前近似解
Figure 999241DEST_PATH_IMAGE021
的最优投影
Figure 668120DEST_PATH_IMAGE022
；
各处理器通过贪心采样方法，并行计算各自对当前近似解
Figure 864746DEST_PATH_IMAGE021
的最优投影
Figure 353496DEST_PATH_IMAGE022
；贪 心采样过程结合随机选择策略和贪心选择策略，随机选择策略可能会产生糟糕的投影，从 而使收敛速度过慢；而贪心选择策略虽然会产生最优的投影，但代价高昂；本实施例的贪心 采样策略构造了一个低维随机子空间，并在随机子空间中产生最优投影，贪心选择策略的 成本在低维空间中非常小，因此，它可以以很少的计算成本产生相对最优的投影，解决了现 有方法收敛速度慢的问题，具体为：
（1）设置初始投影为当前近似解，投影次数为0；
各处理器将初始投影
Figure 293771DEST_PATH_IMAGE023
设置为当前近似解
Figure 817156DEST_PATH_IMAGE021
，投影次数
Figure 919104DEST_PATH_IMAGE024
设置为 0，
Figure 160729DEST_PATH_IMAGE024
为非负整 数，准备开始投影，其中，
Figure 137650DEST_PATH_IMAGE023
为当前近似解
Figure 515542DEST_PATH_IMAGE021
的第0次投影。
（2）基于子方程组和当前投影，计算投影距离范数最大的超平面；
从子方程组
Figure 116288DEST_PATH_IMAGE018
中随机选择多个超平面，每个超平面
Figure 517313DEST_PATH_IMAGE025
被选中 的概率
Figure 861707DEST_PATH_IMAGE026
与超平面的范数成正比，j为超平面的序号，第j个超平面被选中的概率
Figure 297368DEST_PATH_IMAGE026
使用以下 公式计算：
Figure 803435DEST_PATH_IMAGE027
(1)
其中，
Figure 957336DEST_PATH_IMAGE028
是第j个超平面的向量，
Figure 839841DEST_PATH_IMAGE029
是第j个超平面的右端项常数，
Figure 628544DEST_PATH_IMAGE030
表示 “
Figure 305513DEST_PATH_IMAGE031
” 的欧几里得范数，对于每一个选中的超平面，计算投影距离范数
Figure 681130DEST_PATH_IMAGE032
：
Figure 367326DEST_PATH_IMAGE033
(2)
其中，
Figure 512000DEST_PATH_IMAGE034
代表近似解
Figure 625450DEST_PATH_IMAGE021
的第
Figure 488363DEST_PATH_IMAGE024
次投影，比较所有选中超平面的投影距离范数
Figure 712671DEST_PATH_IMAGE032
的大小，并且挑选出具有最大的投影距离范数
Figure 272703DEST_PATH_IMAGE032
值的超平面
Figure 557054DEST_PATH_IMAGE035
，
Figure 844947DEST_PATH_IMAGE036
分别是 具有最大的投影距离范数的超平面的向量和右端项常数，J是具有最大的投影距离范数的 超平面的序号。
（3）将当前投影
Figure 872946DEST_PATH_IMAGE034
投影到计算出来的超平面上，得到新的投影
Figure 992212DEST_PATH_IMAGE037
；
把当前投影
Figure 447464DEST_PATH_IMAGE034
投影到选中的超平面
Figure 347287DEST_PATH_IMAGE035
上，得到新的投影
Figure 851080DEST_PATH_IMAGE037
：
Figure 887170DEST_PATH_IMAGE038
(3)
其中，
Figure 949541DEST_PATH_IMAGE034
代表近似解
Figure 71081DEST_PATH_IMAGE021
的第i次投影。
（4）投影次数i加1，返回步骤（2），计算新的超平面和投影，直到投影次数i等于子 方程组
Figure 378566DEST_PATH_IMAGE018
的系数矩阵分量
Figure 534740DEST_PATH_IMAGE011
的行数d，得到最优投影
Figure 269478DEST_PATH_IMAGE022
。
各处理器并行计算了近似解的多次投影，每个处理器在一次迭代中进行d次投影，与现有方法一次迭代只进行一次投影不同，解决了现有方法因大量的迭代次数和非常小的单次迭代成本使得对每一次迭代步骤进行并行计算因计算量太小而无法有效提高计算效率，反而因并行计算所带来的通信开销使得计算效率大大降低的问题，从而使求解过程可以利用计算机多处理器强大的计算能力。
步骤S303：基于最优投影
Figure 878314DEST_PATH_IMAGE022
，计算新的近似解
Figure 723910DEST_PATH_IMAGE039
；
经过
Figure 734592DEST_PATH_IMAGE016
次投影，各处理器都产生了各自对近似解
Figure 138766DEST_PATH_IMAGE021
的最优投影
Figure 234898DEST_PATH_IMAGE022
；把各处理 器产生的投影
Figure 884185DEST_PATH_IMAGE022
中的各个元素分别取算数平均，产生新的近似解
Figure 14952DEST_PATH_IMAGE039
，记第
Figure 91493DEST_PATH_IMAGE012
个处理器 产生的投影为
Figure 409341DEST_PATH_IMAGE040
，算数平均为：
Figure 862319DEST_PATH_IMAGE041
(4)
其中，
Figure 847593DEST_PATH_IMAGE009
代表求解使用的处理器的个数。
通过取各处理器计算的投影的平均值的方式来给出新的近似解，与现有方法在投影之后立即进行迭代，把投影作为新的近似解不同，通过把现有方法中投影和迭代的过程分离，消除了每次迭代必须依赖前一次迭代的计算结果的问题。
步骤S304：迭代次数 k加1，若本次迭代产生的近似解
Figure 390308DEST_PATH_IMAGE021
不满足产生图像所需要的 精度，则转到步骤S302，进行新一轮的迭代，否则，近似解
Figure 195453DEST_PATH_IMAGE021
即为线性方程组
Figure 389805DEST_PATH_IMAGE001
的解。
步骤S4、将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
具体的，从不同角度对断层进行多次扫描，通过截面的X射线束的吸收过程可以用线性方程组表示；通过迭代方法求解该线性方程组，可以得到截面每个像素点的吸收系数，把吸收系数的大小作为灰度图像的像素点亮度大小，从而绘制断层图像。
选取四组测试数据对本实施例方法的效果进行验证：
（1）使用100000×100维度的线性方程组进行测试，测试规模从1个处理器到8个处理器，分别记录求解时间，计算加速比，加速比如图3所示。
（2）使用500000×500维度的线性方程组进行测试，测试规模从1个处理器到8个处理器，分别记录求解时间，计算加速比，加速比如图4所示。
（3）使用250000×100维度的线性方程组进行测试，测试规模从1个处理器到8个处理器，分别记录求解时间，计算加速比,加速比如图5所示。
（4）使用1250000×500维度的线性方程组进行测试，测试规模从1个处理器到8个处理器，分别记录求解时间，计算加速比,加速比如图6所示。
基于以上四组测试数据的结果，可以看到本实施例方法具有近线性加速比，将处理器的数量加倍，计算速度也会加倍，这表明本实施例方法充分利用了每个参与计算的处理器的性能，各处理器的计算效率接近单处理器求解的理想计算效率，解决了并行求解高精度扫描产生的高维方程组计算效率低、无法有效利用处理器性能的难题，提高了扫描精度和速度。
实施例二
本实施例公开了基于并行随机迭代求解线性方程组的断层成像系统；
如图7所示，基于并行随机迭代求解线性方程组的断层成像系统，包括方程建立模块、任务划分模块、并行求解模块和图像绘制模块：
方程建立模块，被配置为：使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
任务划分模块，被配置为：将线性方程组的求解任务按行划分到各处理器上；
并行求解模块，被配置为：以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
图像绘制模块，被配置为：将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
其中，各处理器通过贪心采样方法计算各自对当前近似解的最优投影；对各处理器的最优投影进行平均，得到新的近似解。
实施例三
本实施例的目的是提供计算机可读存储介质。
计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如本公开实施例一所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
实施例四
本实施例的目的是提供电子设备。
电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本公开实施例一所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，包括：
使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
将线性方程组的求解任务按行划分到各处理器上；
以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
其中，各处理器通过贪心采样方法计算各自对当前近似解的最优投影；对各处理器的最优投影进行平均，得到新的近似解。
2.如权利要求1所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在 于，所述线性方程组为
Figure 902620DEST_PATH_IMAGE001
，A是射线穿过断层各个像素点的路径长度，为系数矩阵，
Figure 126928DEST_PATH_IMAGE002
是断层各个像素点的吸收系数，为待求解的向量， b 是射线的出射强度，为常数项向量。
3.如权利要求1所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，所述按行划分，是对系数矩阵和常数项向量进行按行划分，将划分得到的系数矩阵分量及对应的常数项分量，组成子方程组，分配到不同的处理器上。
4.如权利要求3所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，所述贪心采样方法，具体为：
（1）设置初始投影为当前近似解，投影次数为0；
（2）基于子方程组和当前投影，计算投影距离范数最大的超平面；
（3）将当前投影投影到计算出来的超平面上，得到新的投影；
（4）投影次数加1，返回步骤（2），计算新的超平面和投影，直到投影次数满足要求。
5.如权利要求4所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，所述投影次数满足要求，是投影次数等于子方程组的行数。
6.如权利要求4所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，所述计算投影距离范数最大的超平面，具体为：
从子方程组中随机选择多个超平面，计算每个超平面的投影距离范数，比较投影距离范数的大小，从中计算投影距离范数最大的超平面。
7.如权利要求1所述的基于并行随机迭代求解线性方程组的断层成像方法，其特征在于，所述对各处理器的最优投影进行平均，是对各处理器产生的最优投影的各个元素分别取算数平均，产生新的近似解。
8.基于并行随机迭代求解线性方程组的断层成像系统，其特征在于，包括方程建立模块、任务划分模块、并行求解模块和图像绘制模块：
方程建立模块，被配置为：使用扫描仪对断层进行射线扫描，依据获取的扫描信号，建立以像素点吸收系数为未知数的线性方程组；
任务划分模块，被配置为：将线性方程组的求解任务按行划分到各处理器上；
并行求解模块，被配置为：以设置的初始解为当前近似解，迭代并行计算各处理器对当前近似解的最优投影和新的近似解，直到新的近似解满足设置的求解精度要求，新的近似解为线性方程组的最终解，即为断层每个像素点的吸收系数；
图像绘制模块，被配置为：将吸收系数的大小作为灰度图像的像素点亮度大小，绘制断层图像；
其中，各处理器通过贪心采样方法计算各自对当前近似解的最优投影；对各处理器的最优投影进行平均，得到新的近似解。
9.计算机可读存储介质，其上存储有程序，其特征在于，该程序被处理器执行时实现如权利要求1-7任一项所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
10.电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-7任一项所述的基于并行随机迭代求解线性方程组的断层成像方法中的步骤。
CN202211545390.1A 2022-12-05 2022-12-05 基于并行随机迭代求解线性方程组的断层成像方法及系统 Active CN115619890B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211545390.1A CN115619890B (zh) 2022-12-05 2022-12-05 基于并行随机迭代求解线性方程组的断层成像方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211545390.1A CN115619890B (zh) 2022-12-05 2022-12-05 基于并行随机迭代求解线性方程组的断层成像方法及系统

Publications (2)

Publication Number Publication Date
CN115619890A true CN115619890A (zh) 2023-01-17
CN115619890B CN115619890B (zh) 2023-04-07

Family

ID=84880188

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211545390.1A Active CN115619890B (zh) 2022-12-05 2022-12-05 基于并行随机迭代求解线性方程组的断层成像方法及系统

Country Status (1)

Country Link
CN (1) CN115619890B (zh)

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140307934A1 (en) * 2011-10-31 2014-10-16 Universiteit Antwerpen Dynamic Tomography Angle Selection
US20140369581A1 (en) * 2013-06-14 2014-12-18 The Regents Of The University Of Michigan Iterative reconstruction in image formation
US20150125059A1 (en) * 2013-11-01 2015-05-07 Lickenbrock Technologies, LLC Fast iterative algorithm for superresolving computed tomography with missing data
US20150161077A1 (en) * 2013-12-05 2015-06-11 Canon Kabushiki Kaisha Solution method and solution apparatus for underdetermined system of linear equations
CN105590332A (zh) * 2015-12-24 2016-05-18 电子科技大学 一种应用于ct成像的快速代数重建方法
CN105608719A (zh) * 2015-12-28 2016-05-25 电子科技大学 一种基于两阶段投影调整的快速ct图像重建方法
CN105894562A (zh) * 2016-04-01 2016-08-24 西安电子科技大学 一种光学投影断层成像中的吸收和散射系数重建方法
CN106407561A (zh) * 2016-09-19 2017-02-15 复旦大学 一种并行gpdt算法在多核soc上的划分方法
CN108280859A (zh) * 2017-12-25 2018-07-13 华南理工大学 一种采样角度受限下的ct稀疏投影图像重建方法及装置
CN108765509A (zh) * 2018-05-22 2018-11-06 西南交通大学 一种用于线性成像系统的快速图像重建方法
CN109145255A (zh) * 2018-06-11 2019-01-04 山东省计算中心（国家超级计算济南中心） 一种稀疏矩阵lu分解行更新的异构并行计算方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140307934A1 (en) * 2011-10-31 2014-10-16 Universiteit Antwerpen Dynamic Tomography Angle Selection
US20140369581A1 (en) * 2013-06-14 2014-12-18 The Regents Of The University Of Michigan Iterative reconstruction in image formation
US20150125059A1 (en) * 2013-11-01 2015-05-07 Lickenbrock Technologies, LLC Fast iterative algorithm for superresolving computed tomography with missing data
US20150161077A1 (en) * 2013-12-05 2015-06-11 Canon Kabushiki Kaisha Solution method and solution apparatus for underdetermined system of linear equations
CN105590332A (zh) * 2015-12-24 2016-05-18 电子科技大学 一种应用于ct成像的快速代数重建方法
CN105608719A (zh) * 2015-12-28 2016-05-25 电子科技大学 一种基于两阶段投影调整的快速ct图像重建方法
CN105894562A (zh) * 2016-04-01 2016-08-24 西安电子科技大学 一种光学投影断层成像中的吸收和散射系数重建方法
CN106407561A (zh) * 2016-09-19 2017-02-15 复旦大学 一种并行gpdt算法在多核soc上的划分方法
CN108280859A (zh) * 2017-12-25 2018-07-13 华南理工大学 一种采样角度受限下的ct稀疏投影图像重建方法及装置
CN108765509A (zh) * 2018-05-22 2018-11-06 西南交通大学 一种用于线性成像系统的快速图像重建方法
CN109145255A (zh) * 2018-06-11 2019-01-04 山东省计算中心（国家超级计算济南中心） 一种稀疏矩阵lu分解行更新的异构并行计算方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
F. NOO: ""General reconstruction theory for multislice X-ray computed tomography with a gantry tilt"" *
周伟;董鹏;孙业全;李耀武;: "基于GPU加速和矩阵优化的医学图像重建" *
徐浩;: "求解稀疏线性方程组的预处理共轭梯度并行算法" *
杜亦疏;殷俊锋;张科;: "求解大型稀疏线性方程组的贪婪距离随机Kaczmarz方法" *
黄杰星: ""数字乳腺断层成像的快速迭代重建算法研究"" *

Also Published As

Publication number Publication date
CN115619890B (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
CN104778688B (zh) 点云数据的配准方法及装置
US20220058396A1 (en) Video Classification Model Construction Method and Apparatus, Video Classification Method and Apparatus, Device, and Medium
CN111028327B (zh) 一种三维点云的处理方法、装置及设备
CN110490937B (zh) 一种加速高光谱视频重建的方法及其装置
JP5349407B2 (ja) 平均値シフト手順を使用してサンプルをクラスタリングするプログラム
Jiang et al. Training binary neural network without batch normalization for image super-resolution
US8897544B2 (en) System and method for segmentation of three-dimensional image data
Hazimeh et al. Learning hierarchical interactions at scale: A convex optimization approach
CN109978888A (zh) 一种图像分割方法、装置及计算机可读存储介质
Pan et al. Wide angular sweeping of dynamic electromagnetic responses from large targets by MPI parallel skeletonization
CN114138231B (zh) 执行矩阵乘法运算的方法、电路及soc
CN110365404B (zh) 无波前传感自适应系统及利用该系统提高收敛速度的方法
CN115619890B (zh) 基于并行随机迭代求解线性方程组的断层成像方法及系统
JP6014738B2 (ja) 三次元画像の投影方法
EP4293623A1 (en) Image depth prediction method and electronic device
JP2020041950A (ja) 測量装置、測量方法、及びプログラム
Song et al. DiffusionBlend: Learning 3D Image Prior through Position-aware Diffusion Score Blending for 3D Computed Tomography Reconstruction
CN116989798B (zh) 无人机航迹规划方法
CN105143865A (zh) 用于相关噪声去除的扩展场迭代重构技术（efirt）
KR20160030871A (ko) 그래픽 프로세싱 유닛을 사용하여 히스토그램 계산을 위한 시스템 및 방법
CN116989797B (zh) 无人机航迹优化方法、装置、电子设备及存储介质
CN111210507B (zh) 一种面向多视图三维重建的初始视图选取方法
CN109345465B (zh) 基于gpu的高分辨率图像实时增强方法
Li et al. A deep neural network framework for dynamic multi-valued mapping estimation and its applications
Valencia-Pérez Parallel MLEM algorithm using GPU

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant