对投影值进行线性插值之后再进行滤波反投影的Python实现

前面一篇文章中我介绍了滤波反投影，实际中我们的扫描都是分立而非连续的，因此我们通常需要对投影值进行插值之后再进行滤波反投影，这样能够获得更好的效果。我现在先把代码贴上来，具体的数学过程过几天再详细讲。

#作者：CSDN用户：宋体的微软雅黑
#时间：2020年7月12日
#脚本任务：线性插值之后进行滤波反投影
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from cv2 import cv2
import numba as nb
from numba import jitdef DiscreteRadonTransform(image, steps):res = np.zeros((len(image[0]), steps), dtype='float64')for s in range(steps):rotation = ndimage.rotate(image, -s*180/steps, reshape=False).astype('float64')res[:,s] = np.sum(rotation, axis=0)return res@jit("double[:](int64, int64)", nopython=True, fastmath=False, cache=True)
def SLFilter(N, d):filterSL = np.zeros(N)for i in range(N):filterSL[i] = - 2.0 / (np.pi**2.0 * d**2.0 * (4.0 * (i - N / 2)**2.0 - 1.0))return filterSL@jit("double[:](int64, int64)", nopython=True, fastmath=False, cache=True)
def RLFilter(N, d):filterRL = np.zeros(N)for i in range(N):#判断除数是否为0，如果为0，直接开始下一个循环if np.mod(i - N / 2, 2) == 0:filterRL[i] = 0else:filterRL[i] = - 1.0 / np.power((i - N / 2) * np.pi * d, 2.0)filterRL[np.int64(N/2)] = 1 / (4 * np.power(d, 2.0))return filterRL@jit("double[:, :](int64, int64, double[:, :], double)", nopython=True, fastmath=False, cache=True)
def RLIRadonTransform(steps, pictureSize, projectionValue, delta):res = np.zeros((pictureSize, pictureSize))filter = SLFilter(pictureSize, 1)for step in range(steps):pm = projectionValue[:, step]filterPmWithoutCut = np.convolve(filter, pm)#, "same")#numba仅支持convolve前两个参数，这也就意味着无法将mode设置为"same"，#所以我们需要对卷积的结果进行一个裁剪lenFirst = len(filter)lenSecond = len(pm)lenCore = max(lenFirst, lenSecond)lenRes = len(filterPmWithoutCut)sideLen = np.int64((lenRes - lenCore) / 2.0)firstIndex = np.int64(sideLen)lastIndex = np.int64(lenRes - 1 - sideLen)#print(sideLen, firstIndex, lastIndex)filterPm = filterPmWithoutCut[firstIndex:lastIndex+1]deg = (step - 1) * deltabias = (pictureSize / 2.0) * (1 - np.cos(deg) - np.sin(deg))for row in range(pictureSize):for col in range(pictureSize):pos = bias + (col - 1) * np.cos(deg) + (row - 1) * np.sin(deg)n = np.int64(np.floor(pos))t = pos - nn = np.int64(max(0, n))n = np.int64(min(n, pictureSize - 2))p = (1 - t) * filterPm[n] + t * filterPm[n+1]res[np.int64(pictureSize - 1 - row), np.int64(col)] = res[np.int64(pictureSize - 1 - row), np.int64(col)] + preturn resN = 256
step = 256
theta = np.linspace(0, 180, step, endpoint=False)
print(theta)
delta = np.pi / step
d=1image = cv2.imread("matlabSheppLogan.png", cv2.IMREAD_GRAYSCALE)
P = DiscreteRadonTransform(image, step)
#P = np.array(image).astype(np.float64)
print(P.shape)df = pd.DataFrame(P)
df.to_csv("ProjectioValue.csv", header=False, index=False, mode="w")
res = RLIRadonTransform(step, N, P, delta)plt.figure()
plt.imshow(P, cmap="gray")plt.figure()
plt.title("SL Filter")
plt.imshow(res, cmap="gray")
plt.savefig("iradonWithoutCorrection.png")plt.show()

重建结果如下：
可以看到成像的结果略有改善，但是边缘还是有些“伪影”，我推测是因为我自己写的Radon函数的问题，因为我是直接按列求和产生的投影值，相当于探测器宽度和像素边长相等，而MATLAB中的iradon函数得到的投影值数组的列长度大于Phantom的体素的列的数目，这就相当于使用了更小尺寸的探测器。

希望大家批评指正。

对投影值进行线性插值之后再进行滤波反投影的Python实现相关推荐

【转】由投影重建图像：滤波反投影、FDK、TFDK三维重建算法理论基础
转自:由投影重建图像:滤波反投影.FDK.TFDK三维重建算法理论基础_m0_37357063的博客-CSDN博客_fdk算法 1. 基础理论从: [1] RafaelC.Gonzalez, Rich ...
PDF转图片再转长图、python、pil
PDF转图片再转长图.python.pil 直接贴代码运行环境直接贴代码 # -*- coding: utf-8 -*- """ 1.安装库pip install p ...
使用权值衰减算法解决神经网络过拟合问题、python实现
使用权值衰减算法解决神经网络过拟合问题.python实现一.what is 过拟合二.过拟合原因三.权值衰减四.实验验证 4.1制造过拟合现象 4.2使用权值衰减抑制过拟合一.what is ...
C语言中枚举变量与枚举值，枚举类型变量再赋值问题
1.枚举定义及其使用 1.1 定义枚举是一种特殊的整型,关键词为enum,将变量的值一一列举出来,变量的值只限于列举出来的值的范围内 1.2 使用枚举的定义使用 enum msgtype { eo ...
使用MP进行等值eq查询前，先判断这个值不为null再操作。然后添加多个排序条件
用MybatisPlus进行等值查询和排序操作的代码如下: 这是根据分类type查询对应的分类数据并返回给前端显示.查询到的数据先按照Sort进行升序排序,如果sort相同的情况下再按照updateT ...
mysql 字段a减字段b_SQL 数据库如何实现第一行字段A减字段B得到值C，然后再用C减去第二行字段B，以此类推，求高手解答...
select * ,0 as 缺货 into tmp_r from table_1 --创建结果表 declare @i int declare @j int declare @q int selec ...
indexbar 将数组键值转化为字母再通过字母分组，排序。
这段时间需要写一个手机通讯录的功能,后端只给了数据,没有分组,也没有字母分类,下面将自己转化字母,分类,排序,实现这个功能. 效果图首先新建俩个js文件. pinying.js export con ...
字典删除多个键值对方法_Life is short，you need Python——Python序列（元组、字典、集合）...
一.元组 tuple 列表属于可变序列,可以任意修改列表中的元素. 元组属于不可变序列,不能修改元组中的元素.因此,元组没有增加元素.修改元素.删除元素相关的方法. 下面只介绍元组的创建和删除,元组中 ...
goupby 两个值结果变了_一道问题引出的python中可变数据类型与不可变数据类型...
一. 问题的提出我们先来看两个对比第一道题,当对象为整数时,最终结果:b = 2, a = 1,b的变化没有引起a的变化第二道题,当对象为字典时,最终结果:a = {"name&quo ...

对投影值进行线性插值之后再进行滤波反投影的Python实现

对投影值进行线性插值之后再进行滤波反投影的Python实现相关推荐

最新文章

热门文章