精华内容
下载资源
问答
  • 迭代重建算法中投影矩阵的计算
    2021-08-10 10:32:51

            在前面学习的重建算法都是属于解析法,它是以Radon变换为理论基础,首先对投影数据在连续域进行一些处理,然后将其离散化进行重建。接下来学习到的是图像重建的迭代算法,该算法的主要思路是:从一幅假设的初始图象出发,采用逐步逼近的方法,将理论投影值与实际测量投影值进行比较,在某种最优准则下寻找最优解。

            该算法本质上类似于解方程组,但是在实际的图像重建过程中,由于运算量巨大、方程的欠定性以及测量误差、噪声的影响等原因,直接通过求逆矩阵来解方程组是比较困难的,很难将其应用在实际生活中。于是产生了一系列迭代方法来解决这个问题。常用的迭代方法包括代数重建法(ART)、联合代数重建算法(SART)、最大似然期望最大化算法(MLEM)和有序子集期望最大化算法(OSEM),其中最后两种算法主要用于PET和SPECT。接下来主要对ART算法进行解释。

           ART算法是一个‘行运算’算法,它每次考虑一条射线就更新一次图像。在整个ART运算中,投影矩阵是迭代重建的关键因素,是获取投影数据与断层图像矢量相联系的桥梁。如果求解的投影矩阵越精确,那么重建后的图像的质量越高。在我们计算投影矩阵的时候,主要有四种模

    更多相关内容
  • 基于线积分模型的代数迭代算法程序,整个程序是基于matlab开发。压缩包内有两个*.m文件,ite.m是执行程序,medfuncMatirx.m是计算正投影矩阵的程序。正投影矩阵的计算是基于线积分模型,结合西北大学张顺利老师的...
  • 利用迭代重建法SIRT,用户可以通过扫描得到的数据或者读取有关的数据内容,输入到工具箱,然后通过工具箱的计算得到该事物的雏形。 图象重建工具箱的优势在于应用范围很广,无论医疗,建筑,工艺制造等等,都能用...
  • 基于CUDA的GPU加速代数迭代重建算法.pdf
  • 提出了一种基于代数迭代的Mojette变换层析重建算法。在确定最佳投影角度的基础上, 结合传统层析技术中的乘性代数迭代算法进行重建。利用该算法对轴对称火焰进行了层析重建。数值模拟结果表明, 与基于角的重建算法...
  • 代数迭代算法ART

    2016-12-23 19:25:39
    ART代数迭代法 基于vc6.0 C++
  • 广义代数重建算法(E-ART)可以重建得到高质量的双基材料图像,但是该算法收敛速度较慢。针对该问题,提出一种加速收敛算法(AE-ART),其基本思想是通过为基函数添加权重来增大高能和低能投影曲线的夹角或减小高能和低能...
  • 针对传统重建算法对火焰重建精度低、重建速度慢的问题,提出了基于正则先验的全变差代数迭代(ARTTV)算法,以提高对称与非对称火焰的重建精度。同时,为了提高重建速度,建立了基于“ARTTV-粒子群算法(PSO)内核”的极限...
  • 构造了平滑约束矩阵作为先验信息引入到重建迭代过程, 建立了一种平滑约束OSEM(SC-OSEM)迭代重建算法。分别将中值滤波、全变差最小(TVM)方法作为平滑约束条件, 通过数值模拟, 针对不完备理想投影数据、含金属不完备...
  • SIRT-FISTA-TV是一种规则化的迭代重建算法,对嘈杂和模糊的数据非常健壮,并且可以大大减少丢失的楔形伪影。 它包括三个步骤: SIRT更新(也可以使用SART和OS-SART) 电视最小化(使用渐变下降) FISTA技术可加快...
  • 采用傅里叶变换相位展开技术处理莫尔条纹,提取轴向任意截面的投影信息,使用一种新的莫尔层析迭代算法重建温度分布,该迭代算法与已有的代数重建算法不同,它由莫尔偏折基本公式和网格重建技术推导得出,直接使用...
  • ART迭代算法图像重建

    2020-12-06 02:00:13
    采用代数迭代算法进行图像重建的MATLAB程序,具有借鉴的价值,欢迎大家分享,谢谢! 采用代数迭代算法进行图像重建的MATLAB程序,具有借鉴的价值,欢迎大家分享,谢谢!
  • 为研究层析成像算法对温度场二维重建质量的影响,实现了两种典型重建算法:代数迭代重建算法(ART)和模拟退火(SA)算法。在不同的射线分布和吸收谱线数目情况下,使用两种算法对给定单峰温度场和双峰温度场分别...
  • 迭代重建技术(ART)简要介绍

    千次阅读 2020-11-12 20:44:28
    本文是基于matlab实现的最简单的迭代重建算法,当然,基于这个算法可以做很多优化。 1. fomula 迭代方法,其实说来都比较简单,最关键就是怎么从当前的x得到下一步的x。 Algebraic reconstruction technique is all ...

    0 说明

    本文是基于matlab实现的最简单的迭代重建算法,当然,基于这个算法可以做很多优化。

    1. fomula

    迭代方法,其实说来都比较简单,最关键就是怎么从当前的x得到下一步的x。
    Algebraic reconstruction technique is all about a system of equations which is

    Ax = p

    In this formula, x = [x1, x2, …, xn]T , p = [p1, p2, …, pm]T , where x is one-dimensional vector of n elements, each representing the pixels of a certain 2-D object, p is one-dimensional vector of m elements, each representing a measurements taken at a certain detector in a certain projection direction, A is called the projection matrix and each value in this matrix called a projection rate defines the contribution of a specific pixel to a specific projection value.

    in general medical imaging applications, A is so large that the entire matrix cannot be stored in the computer at the same time. In the process of solving the equations, the value of that row of matrix A is generated only when needed and released immediately after used up. Based on this limitation, some algorithms,for example, diagonalizing or triangularizing matrix A, are not applicable. In fact, any algorithm that needs to change matrix A cannot be used.

    Let’s focus on the implementation of the iterative reconstruction algorithm. Look at the following formula, which is basically the same as the ordinary gradient descent algorithm, where lambda is our relaxation factor determining the convergence rate.

    在这里插入图片描述

    2 Situations of different projection data

    The entire iterative process can be seen in the two figure on the right. x0 is an arbitrary initial value. The first step is to project the point x0 vertically to L1 to get the new value x1. The next step is to project the point x1 vertically onto L2 to get the new value x2, and so on. Each step is to project the current estimated point onto the next straight line so that it satisfies the next equation. Eventually, this algorithm will converge to the solution of the system of equations.
    下图是相容方程组,即针对方程组是有一个精确解的。迭代次数足够时,能够得到一个收敛值。
    在这里插入图片描述
    If this system of equations is incompatible, this algorithm will cause the “solution” of the system of equations to jump back and forth without converging. one iteration is defined as all equations are accessed by the algorithm once.
    方程组不相容时,就会出现如下状况,在迭代次数到一定数值后,我们会在中间的三角形区域内一直跳来跳去。这代表三角形区域内任何一个点都是该方程组的数值解,比较好的方法是我们取得几个解,然后做一下平均。
    在这里插入图片描述
    上面两种不一样的投影数据,第一种明显是更好的。这里有一个问题值得注意,实际操作中应该尽可能的让我们的projection data靠近第一种。

    3 implement

    one iteration之后的图片:
    在这里插入图片描述
    这是基于90个angle得到的投影数据,然后通过泊松分布加入噪声后,再进行迭代重建后得到的图片,可以发现,这时图片非常模糊。

    在这里插入图片描述
    ground truth 是我们的原始图片,第二排的两张是将fbp重建的图片和迭代200次后重建的图片进行比对。
    The figure on the left is rendered after one iteration, and it is still blurry compared to the original image. The figure on the right is based on the projection data of 90 projection angles, which are reconstructed using filtered back projection and iterative reconstruction techniques. We can clearly see the difference between the two algorithms. The edge of the image after the filtered back projection reconstruction is clearer and the noise is greater, and the image obtained by the iterative reconstruction technique is smoother and more blurred. We compare the mean square error and structural similarity of the two methods. The numerical display shows that ART is closer to the original image, but it has a lot of shortcomings, as explained above.

    4 contrast

    • 首先,看一下迭代次数改变后,均方误差还有结构相似性的一个变化曲线。
      在这里插入图片描述
      图片有点不是很清晰,需要看仔细一点,可以发现在迭代次数超过一定的数值后,均方误差不减反增。
      Based on my personal understanding, the possible reason is that there are only 90-angle projection data, and the complexity of the projection data is less than the complexity of the image to be reconstructed. after a certain number of iterations, because of the excess expressive ability of the model, it will train some features that can only meet the projection data but not the original image, which will lead to a increase in mse.
    • 再看一下fbp和art的一个比对。整个图像来看可能看的不是很clear, 所以我单独将image matrix的其中一条直线plot出来了。这时对比就很清晰了。
      在这里插入图片描述
      左边是fbp和原始图片比对的曲线,右边是art和原始图片比对的曲线。两种算法个有优劣,可以根据需求选择不同的算法。
      The difference is very obvious. The fbp image has larger fluctuations and ART is more stable. Although fbp fluctuates greatly, there is no loss of energy. Observation of the ART curve can clearly find that although it performs better in a smooth place, there is a big gap between the curve of the original image at the peak and this will cause image details Lost.

    5 写在后面的话

    非常感谢各位jmm和xdm能看到这,本人研究生萌新,方向医学图像重建,欢迎有兴趣的朋友们一起交流学习。

    如果对filtered back project 算法不太了解的话,请查看:
    https://blog.csdn.net/sharon_1995/article/details/109454562

    如果想简单了解x-ray,请查看:
    https://blog.csdn.net/sharon_1995/article/details/109652834

    展开全文
  • 通过本文的研究,对原有的ART迭代重建算法进行了改进,得到一种成像质量更高且重建速度更快的图像重建算法,也适用不完全投影的情况,来降低辐射剂量,缩短扫描时间,进而潜在的减少病人的运动伪影. 展开

    摘要:

    自CT技术诞生以来,它就被认为是二十世纪影响人类发展的十大技术之一.其应用更是涉及了多个领域,包括医学,生物,工业,安全等等.特别是在医学方面,CT更是成为了医学检测的主要手段.而口腔CT(Dental CT),也称为牙科CT,则是20世纪90年代后期才兴起的一项新兴技术,它是利用锥束CT原理来进行口腔三维成像的一种技术.口腔CT作为一种无创伤的检测技术具有很多的优越性,它的出现也受到了对CT技术进行研究的学者以及牙科临床医师的广泛关注.而口腔CT技术的核心就是重建算法.用于CT设备的成像算法主要分为两类:解析重建算法和代数重建算法.目前,口腔CT系统大多都采用的是解析重建算法.解析重建算法在重建的过程中计算量较小,重建的速度也较快,而且在获得完整且均匀的投影数据的情况下可以得到高质量的重建图像.但是,它要求在投影数据采集的过程中必须完全均匀,且采集到的投影数据必须完整,在实际应用的过程中,由于存在一些客观的原因,往往不能够满足其要求,采集的投影数据往往不均匀或者不完整.对于这种情况而言,代数重建方法则可以很好的弥补这些缺点.在代数重建方法中最具代表性的则是ART (AlgebraicReconstruction Techniques)算法.这种方法是把抽象的图像重建的问题转换为具体的求解线性方程组的问题.遇到缺少投影数据的问题的时候,我们不妨将其认为是少了若干个方程,它仍然能够通过迭代的算法来求出最佳的解.这样以来,可以在一定程度上淡化掉缺少投影数据的情况,所以可适用于不完全投影情况下的图像重建. 在本文中对图像重建算法中的ART迭代重建算法进行了深入的分析与讨论,并对其性能进行了改进.在本文中主要研究了下面这些内容:对论文的课题来源,研究背景意义以及国内外的研究现状进行了介绍,对本论文研究的主要内容作了综述;对CT成像的基本原理以及CT研究领域一些主要的基础理论进行了介绍;对口腔CT系统的整体设计结构进行了介绍,分析了一些相关参数的选择,系统的整体设计,关键器件的选型,系统控制以及数据采集和预处理等方面的相关内容;提出了一种对三维图像进行投影的方法,并对投影过程进行仿真模拟,且不同于常用的二维Shepp-Logan模体,本文中还使用Fortran语言对三维仿真模体Shepp-Logan进行编写,并针对实验要求对其参数作出了相应的调整,将其作为验证ART迭代重建算法性能的依据;在实现三维图像的ART重建算法的过程中,指定了对重建图像进行质量评价的标准,对影响ART重建算法重建结果的各个因素进行了分析与实验,并选出最优的参数以得到最佳的成像结果;为了提高重建图像的质量,将POCS约束加入到三维图像重建的ART算法中去.经过仿真实验与误差分析,可以看出加入POCS约束之后的ART算法较原算法来说成像质量更高,空间分辨率更高,误差也更小,特别是对物体周围的伪影基本上全部被消除了.而为了减少辐射剂量,使其在不完全投影的情况下进行重建,也能够达到较好的成像效果;为了进一步提高算法的性能,我们对加入POCS约束的ART算法进行了改进.针对投影矩阵为大的稀疏矩阵,在进行迭代重建的过程中会出现一些不必要的计算而使得计算量增大这一不足之处,在迭代重建的过程中对投影矩阵进行处理,限定一个条件来忽略掉一些无效的投影值来减少迭代过程中的计算量.经过仿真实验与误差分析,可以看出改进之后的算法在改善图像质量的同时,大大加快了迭代重建的速度. 通过本文的研究,对原有的ART迭代重建算法进行了改进,得到一种成像质量更高且重建速度更快的图像重建算法,也适用不完全投影的情况,来降低辐射剂量,缩短扫描时间,进而潜在的减少病人的运动伪影.

    展开

    展开全文
  • 为此提出了一种新的光学层析技术的代数迭代重建算法,在算法中引入了包含先验知识的属性矩阵,并摒弃了通常所采用的对超松弛系数人为的确定取法,采用了变超松弛系数。实验计算结果表明,引入属性矩阵和变超松弛系数的...
  • ART与SART代数重建算法

    千次阅读 2020-12-06 21:27:23
    c: 穿过格子的序号 d2: 穿过的格子的右侧与该方格右下角顶点的距离 """ k = kin # 入射的格子即为穿过的第一个格子 c = 0 # c为穿过格子的序号 d2 = d1 + m * delta # 与方格右侧交点 # 当方格数在1~n^2内迭代计算 ...

    1.ART方法:

    1.代码实现:

    import numpy as np
    import matplotlib.pyplot as plt
    import cv2
    from scipy import ndimage
    import time
    
    
    def projection(image, theta):
        """
        计算投影值
        :param image: 原始图像
        :param theta: 射束旋转角度
        :return: 投影值矩阵
        projectionNum: 射束个数
        thetaNum: 角度个数
        randontansform: 拉登变换结果
        """
        print('正在计算投影...')
    
        projectionNum = len(image[0])
        thetaNum = len(theta)
        radontansform = np.zeros((projectionNum, thetaNum), dtype='float64')
        for i in range(len(theta)):
            # 进行离散拉登变换
            rotation = ndimage.rotate(image, -theta[i], reshape=False).astype('float64')
            radontansform[:, i] = sum(rotation)
        return radontansform
    
    
    def systeMartrix(theta, size, num, delta):
        """
        计算系统矩阵
        :param theta: 射束旋转角度
        :param size: 图片尺寸
        :param num: 射束条数
        :param delat: 网格边长
        :return: gridNum:穿过网格编号,gridLen:穿过网格长度
        """
        print('正在计算系统矩阵...')
    
        start = time.perf_counter()
        functionNum = len(theta) * num  # 方程个数,即为系统矩阵的行数
        # 射线最多穿过2*size个方格
        gridNum = np.zeros((functionNum, 2 * size))  # 系统矩阵:编号
        gridLen = np.zeros((functionNum, 2 * size))  # 系统矩阵:长度
        N = np.arange(-(size - 1) / 2, (size - 1) / 2 + 1)  # 射束
        for loop1 in range(len(theta)):
            th = theta[loop1]  # 射束角度
            for loop2 in range(size):
                u = np.zeros((2 * size))  # 编号
                v = np.zeros((2 * size))  # 长度
    
                # 垂直入射
                if th == 0:
                    # 射束未穿过图像时
                    if (N[loop2] >= size / 2 * delta) or (N[loop2] <= -size / 2 * delta):
                        continue
                    # 入射网格编号
                    kin = np.ceil(size / 2 + N[loop2] / delta)
                    # 穿过网格编号
                    kk = np.arange(kin, (kin + size * size), step=size)
                    u[0:size] = kk
                    v[0:size] = np.ones(size) * delta
    
                # 平行入射
                elif th == 90:
                    if (N[loop2] >= size / 2 * delta) or (N[loop2] <= -size / 2 * delta):
                        continue
                    # 出射网格编号
                    kout = size * np.ceil(size / 2 - N[loop2] / delta)
                    kk = np.arange(kout - size + 1, kout + 1)
                    u[0:size] = kk
                    v[0:size] = np.ones(size) * delta
    
                else:
                    # phi为射束与x轴所夹锐角
                    if th > 90:
                        phi = th - 90
                    elif th < 90:
                        phi = 90 - th
                    # 角度值换算为弧度制
                    phi = phi * np.pi / 180
                    # 截距
                    b = N / np.cos(phi)
                    # 斜率
                    m = np.tan(phi)
                    # 入射点纵坐标
                    y1 = -(size / 2) * delta * m + b[loop2]
                    # 出射点纵坐标
                    y2 = (size / 2) * delta * m + b[loop2]
    
                    # 射束未穿过图像
                    if (y1 < -size / 2 * delta and y2 < -size / 2 * delta) or (
                            y1 > size / 2 * delta and y2 > size / 2 * delta):
                        continue
    
                    # 穿过a、b边(左侧和上侧)
                    if (y1 <= size / 2 * delta and y1 >= -size / 2 * delta and y2 > size / 2 * delta):
                        """
                        (xin,yin): 入射点坐标
                        (xout,yout): 出射点坐标
                        kin,kout: 入射格子标号,出射格子编号
                        d1: 入射格子左下角与入射射束距离
                        """
                        yin = y1
                        yout = size / 2 * delta
                        # xin = -size / 2 * delta
                        xout = (yout - b[loop2]) / m
                        kin = size * np.floor(size / 2 - yin / delta) + 1
                        kout = np.ceil(xout / delta) + size / 2
                        d1 = yin - np.floor(yin / delta) * delta
    
                    # 穿过a、c边(左侧和右侧)
                    elif (
                            y1 <= size / 2 * delta and y1 >= -size / 2 * delta and y2 >= -size / 2 * delta and y2 < size / 2 * delta):
                        # xin = -size / 2 * delta
                        # xout = size / 2 * delta
                        yin = y1
                        yout = y2
                        kin = size * np.floor(size / 2 - yin / delta) + 1
                        kout = size * np.floor(size / 2 - yout / delta) + size
                        d1 = yin - np.floor(yin / delta) * delta
    
                    # 穿过d、b边(下侧和上侧)
                    elif (y1 < - size / 2 * delta and y2 > size / 2 * delta):
                        yin = - size / 2 * delta
                        yout = size / 2 * delta
                        xin = (yin - b[loop2]) / m
                        xout = (yout - b[loop2]) / m
                        kin = size * (size - 1) + size / 2 + np.ceil(xin / delta)
                        kout = np.ceil(xout / delta) + size / 2
                        d1 = size / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]
    
                    # 穿过d、c边(下侧和右侧)
                    elif (y1 < - size / 2 * delta and y2 >= -size / 2 * delta and y2 < size / 2 * delta):
                        yin = -size / 2 * delta
                        yout = y2
                        xin = (yin - b[loop2]) / m
                        # xout = size / 2 * delta
                        kin = size * (size - 1) + size / 2 + np.ceil(xin / delta)
                        kout = size * np.floor(size / 2 - yout / delta) + size
                        d1 = size / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]
    
                    else:
                        continue
    
                    # 计算穿过的格子编号和长度
                    """
                    k: 射线穿过的格子编号
                    c: 穿过格子的序号
                    d2: 穿过的格子的右侧与该方格右下角顶点的距离
                    """
                    k = kin  # 入射的格子即为穿过的第一个格子
                    c = 0  # c为穿过格子的序号
                    d2 = d1 + m * delta  # 与方格右侧交点
    
                    # 当方格数在1~n^2内迭代计算
                    while k >= 1 and k <= np.power(size, 2):
                        """
                        根据射线与方格的左右两侧的交点关系,来确定穿过方格的六种情况。
                        在每种情况中,存入穿过的方格编号,穿过方格的射线长度。
                        若该方格是最后一个穿过的方格,则停止迭代;若不是最后一个方格,则计算下一个穿过的方格的编号、左右边与射线的交点。
                        """
                        if d1 >= 0 and d2 > delta:
                            u[c] = k  # 穿过方格的编号
                            v[c] = (delta - d1) * np.sqrt(np.power(m, 2) + 1) / m  # 穿过方格的射线长度
                            if k > size and k != kout:  # 若不是最后一个方格
                                k -= size  # 下一个方格编号
                                d1 -= delta  # 下一个方格左侧交点
                                d2 = delta * m + d1  # 下一个方格右侧交点
                            else:  # 若是最后一个方格则直接跳出循环
                                break
    
                        elif d1 >= 0 and d2 == delta:
                            u[c] = k
                            v[c] = delta * np.sqrt(np.power(m, 2) + 1)
                            if k > size and k != kout:
                                k -= size + 1
                                d1 = 0
                                d2 = d1 + m * delta
                            else:
                                break
    
                        elif d1 >= 0 and d2 < delta:
                            u[c] = k
                            v[c] = delta * np.sqrt(np.power(m, 2) + 1)
                            if k != kout:
                                k += 1
                                d1 = d2
                                d2 = d1 + m * delta
                            else:
                                break
    
                        elif d1 <= 0 and d2 >= 0 and d2 <= delta:
                            u[c] = k
                            v[c] = d2 * np.sqrt(np.power(m, 2) + 1) / m
                            if k != kout:
                                k += 1
                                d1 = d2
                                d2 = d1 + m * delta
                            else:
                                break
    
                        elif d1 <= 0 and d2 > delta:
                            u[c] = k
                            v[c] = delta * np.sqrt(np.power(m, 2) + 1) / m
                            if k > size and k != kout:
                                k -= size
                                d1 -= delta
                                d2 = d1 + m * delta
                            else:
                                break
    
                        elif d1 <= 0 and d2 == delta:
                            u[c] = k
                            v[c] = d2 * np.sqrt(np.power(m, 2) + 1) / m
                            if k > size and k != kout:
                                k -= size + 1
                                d1 = 0
                                d2 = delta * m
                            else:
                                break
    
                        else:
                            print(d1, d2, "数据错误!")
    
                        c += 1
    
                    # 当射线斜率为负数时,利用对称性进行计算
                    if th < 90:
                        u_temp = np.zeros(2 * size)
                        # 排除掉未穿过图像的射束
                        if u.any() == 0:
                            continue
                        # 射束穿过的格子的编号
                        indexMTZero = np.where(u > 0)
    
                        # 利用对称性求得穿过网格编号
                        for loop in range(len(u[indexMTZero])):
                            """计算穿过的方格编号,利用方格编号与边长的取余的关系,得到对称的射线穿过的方格编号。"""
                            r = np.mod(u[loop], size)
                            if r == 0:
                                u_temp[loop] = u[loop] - size + 1
                            else:
                                u_temp[loop] = u[loop] - 2 * r + size + 1
                        u = u_temp
    
                # 方格编号
                gridNum[loop1 * num + loop2, :] = u
                # 穿过方格的射线长度
                gridLen[loop1 * num + loop2, :] = v
    
        end = time.perf_counter()
        print('系统矩阵计算耗时:%s 秒。'%(end - start))
    
        return gridNum, gridLen
    
    
    def iteration(theta, size, gridNum, gridLen, F, ite_num):
        """
        按照公式迭代重建
        :param theta: 旋转角度
        :param size: 图像边长
        :param gridNum: 射线穿过方格编号
        :param gridLen: 射线穿过方格长度
        :param F: 重建后图像
        :return: 重建后图像
        """
        print('正在进行迭代...')
    
        start = time.perf_counter()
        c = 0  # 迭代计数
        while (c < ite_num):
            for loop1 in range(len(theta)):  # 在角度theta下
                for loop2 in range(size):  # 第loop2条射线
                    u = gridNum[loop1 * size + loop2, :]
                    v = gridLen[loop1 * size + loop2, :]
                    if u.any() == 0:  # 若射线未穿过图像,则直接计算下一条射线
                        continue
                    # 本条射线对应的行向量
                    w = np.zeros(sizeSquare, dtype=np.float64)
                    # 本条射线穿过的网格编号
                    uLargerThanZero = np.where(u > 0)
                    # 本条射线穿过的网格长度
                    w[u[uLargerThanZero].astype(np.int64) - 1] = v[uLargerThanZero]
                    # 计算估计投影值
                    PP = w.dot(F)
                    # 计算实际投影与估计投影的误差
                    error = projectionValue[loop2, loop1] - PP
                    # 求修正值
                    C = error / sum(np.power(w, 2)) * w.conj()
                    # 进行修正
                    F = F + lam * C
            F[np.where(F < 0)] = 0
            c = c + 1
    
        F = F.reshape(size, size).conj()
    
        end = time.perf_counter()
        print('迭代耗时:%s 秒。'%(end - start))
    
        return F
    
    
    if __name__ == '__main__':
        image = cv2.imread("shepp-logan(256).png", cv2.IMREAD_GRAYSCALE)
        """
        theta: 射束旋转角度
        num: 射束条数
        size: 图片尺寸
        delta: 网格边长
        ite_num: 迭代次数
        lam: 松弛因子
        """
        theta = np.linspace(0, 180, 60, dtype=np.float64)
        num = np.int64(256)
        size = np.int64(256)
        delta = np.int64(1)
        ite_num = np.int64(10)
        lam = np.float64(.25)
        sizeSquare = size * size
    
        # 计算投影值
        projectionValue = projection(image, theta)
    
        # 计算系统矩阵
        gridNum, gridLen = systeMartrix(theta, size, num, delta)
    
        # 重建后图像矩阵
        F = np.zeros((size * size,))
    
        # 迭代法重建
        reconImage = iteration(theta, size, gridNum, gridLen, F, ite_num)
        print('共迭代%d次。'%ite_num)
    
    
        # 绘制原始图像、重建图像、误差图像
        plt.subplot(1, 3, 1)
        plt.imshow(image, cmap='gray')
        plt.title('Original Image')
    
        plt.subplot(1, 3, 2)
        plt.imshow(reconImage, cmap='gray')
        plt.title('Reconstruction Image')
    
        plt.subplot(1, 3, 3)
        plt.imshow(reconImage - image, cmap='gray')
        plt.title('Error Image')
    
        plt.savefig("ART.png", cmap='gray')
    
        plt.show()
    
    

    2.耗时:

    在这里插入图片描述

    3.重建结果:

    在这里插入图片描述

    2.SART方法:

    1.代码实现:

    对迭代过程进行修改:

    def iteration(theta, size, gridNum, gridLen, F, ite_num, num, sizeSquare, projectionValue, lam):
        """
        按照公式迭代重建
        :param theta: 旋转角度
        :param size: 图像边长
        :param gridNum: 射线穿过方格编号
        :param gridLen: 射线穿过方格长度
        :param F: 重建后图像
        :param sizeSquare: 图像边长平方
        :param projectionValue: 投影值
        :param lam: 松弛因子
        :return: 重建后图像
        """
        print('正在进行迭代...')
    
        start = time.perf_counter()
        c = 0  # 迭代计数
    
        while (c < ite_num):
            for loop1 in range(len(theta)):  # 在角度theta下
                u = gridNum[loop1 * num : (loop1 + 1) * num, :]  # 射线穿过网格编号
                v = gridLen[loop1 * num : (loop1 + 1) * num, :]  # 射线穿过网格长度
                w = np.zeros((num, sizeSquare))
                for loop2 in range(num):  # 对第loop2条射线
                    if u[loop2, :].any() == 0:  # 若射线未穿过图像,则直接计算下一条射线
                        continue
                    for loop3 in range(2 * size):
                        m = u[loop2, loop3]
                        if ((m > 0) and (m <= sizeSquare)):
                            w[loop2, np.int64(m)-1] = v[loop2, loop3]
                sumcol = np.sum(w, axis=0).reshape(-1, 1)
                sumrow = np.sum(w, axis=1).reshape(-1, 1)
                ind1 = np.where(sumrow > 0)
                corr = np.zeros((num, 1))
                err = projectionValue[:, loop1].reshape((256, 1)) - w.dot(F)
                corr[ind1] = err[ind1] * (1 / sumrow[ind1])
                backPro = (w.T).dot(corr)
                ind2 = np.where(sumcol > 0)
                delta = np.zeros((sizeSquare, 1))
                delta[ind2] = backPro[ind2] / sumcol[ind2]
                F += lam * delta
                F[np.where(F < 0)] = 0
            c += 1
                
        F = F.reshape(size, size)
        end = time.perf_counter()
        print('迭代耗时:%s 秒。' % (end - start))
    
        return F
    

    注意此时在main中创建F时应改为:

        # 重建后图像矩阵
        F = np.zeros((sizeSquare, 1))
    

    2.耗时:

    在这里插入图片描述

    3.重建结果:

    在这里插入图片描述

    展开全文
  • 针对运用最小二乘法求解DV-Hop定位算法带来的节点定位误差较大的问题,提出基于代数重建法的DV-Hop定位算法,运用一种由图象重建问题而引入的逐次迭代算法——代数重建法。仿真结果表明,改进算法能降低无线...
  • 该算法首先对CT稀疏投影数据采用联合代数重建算法(SART )进行重建,以获得满足数据一致性的重建图像,然后计算SART重建图像的离散梯度变换,并对其进行软阂值滤波,最后利用离散梯度变换的伪逆更新重建图像.由于在迭代...
  • 代数重建技术 参考 代数重建技术(ART) 代数重建技术(ART) 是使用的一类iterative algorithms 。 ART相对于其他重建方法(例如确定性方法 )的优势在于,将先验知识整合到重建过程中相对容易。 x^(k+1) = x^(k...
  • 平滑约束的OSEM代数重建算法激光与光电子学进展, ( ) « »54 0210062017 C2017 中国激光 杂志社Laser&OtoelectronicsPro...
  • 详细分析了投影噪声、投影方向数、场分布性质对重建精度的影响,并与代数迭代重建算法结果进行对比.结果表明,该算法以两个正交方向投影数据重建单峰余弦模拟场平均误差仅为0.03%,而代数迭代重建算法为3.81%;该算法以...
  • Python实现ART重建算法

    千次阅读 多人点赞 2020-06-26 21:28:30
    投影矩阵代数重建算法的基础,它将投影数据和断层图像联系了起来,投影矩阵的计算方法也将影响重建图像的质量,投影矩阵的模型可以分为以下几种: 把射束看为是宽度为0,间距为δ\deltaδ的一系列直线,将尺寸为N∗...
  • 1 简介 ART算法是一个不断迭代的图像重建方法,提高该算法的重建速度一直是研究的重要方面.针对ART算法简化权因子重建...英文名称为Algebraic reconstruction technique:即代数重建算法 代数重建技术(ART)是一种..
  • 在分析了相位衬度CT的特点之后,将压缩感知理论引入相位衬度CT重建中,并在该理论框架下将L1约束融入代数迭代重建(ART)算法中,提出了一种微分相位衬度CT重建算法。数值模拟和实际实验表明,该方法可以根据少量投影...
  • 采用代数迭代(ART)算法和最大似然期望最大化(MLEM)算法,利用开路傅里叶变换红外(OP-FTIR)光谱仪的测量结果,通过仿真模拟了高斯空间分布模型下的气体二维浓度场重建,并利用重建评价指标——逼近度和相关系数,分析了这...
  • 音视频-编解码-迭代法图像重建中系统矩阵构造算法的研究.pdf

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,848
精华内容 739
热门标签
关键字:

代数迭代重建算法