CBP(卷积反投影)实现
CBP公式如下:
f(x1,x2)=∫0πf^(r,ϕ)∗H(r)∣r=(x1,x2)⋅ϕdϕf(x_1,x_2)=\int_0^{\pi}\hat{f}(r,\phi)*H(r)|_{r=(x_1,x_2)\cdot{\phi}}d\phif(x1,x2)=∫0πf^(r,ϕ)∗H(r)∣r=(x1,x2)⋅ϕdϕ
求f(x1,x2)f(x_1,x_2)f(x1,x2)就简单地分成以下几步:
- 求滤波器h(w)对应的卷积核。需要注意的是h(w)在频域有限,在时域就是无穷的,类似一系列冲击函数的叠加。
函数表示为:
显示h图像为 - 对于每一个ϕ\phiϕ,将f^(r,ϕ)\hat{f}(r,\phi)f^(r,ϕ)与H(r)H(r)H(r)做卷积,取中间(N)个卷积结果
- 此时得到图像结果是离散的,也就是对于r=nd取不到的点,所以这时需要对r方向进行线性插值,把r变成连续的,可取值的;
这时我们随便给一个r值,就都有对应的灰度值了 - 接下来我们就需要把原图像(x1,x2)(x_1,x_2)(x1,x2)的信息提取出来。已知r=(x1,x2)⋅ϕr=(x_1,x_2)\cdot\phir=(x1,x2)⋅ϕ那么x1=rcosϕx_1=rcos\phix1=rcosϕx2=rsinϕx_2=rsin\phix2=rsinϕ
经过上述操作,我们知道了,所谓f(x1,x2)f(x_1,x_2)f(x1,x2)的值,就是经过这一点的所有射线radon值的叠加。
射线由标准位置ϕ\phiϕ从0度到180度,r与直线夹角就是π/2−ϕ\pi/2-\phiπ/2−ϕ度,那r是很好求的r=x1cos(π/2−ϕ)+x2sin(π/2−ϕ)r=x_1cos(\pi/2-\phi)+x_2sin(\pi/2-\phi)r=x1cos(π/2−ϕ)+x2sin(π/2−ϕ)
具体代码
import shelve
import numpy as np
import math
import matplotlib.pyplot as plt #plt用于显示图片
tuoyuanlist=shelve.open('test')
m=tuoyuanlist['2Dlist']#此时m由二维数组转化为矩阵,好操作,存储了r,phi对应位置的值。#以下为一些常量定义,d长度为探测器元长度
d=0.02
phi = np.linspace(0.0314,math.pi,500)#将(0,pi)平均分为500份#首先应该有一个长度为n的h对f做卷积操作
#定义R_L卷积核为:
def R_L(n):h=np.zeros(n)for i in range(n):if(i== (n/2-1)):#n=0时取值h[i]=1/(4*pow(d,2))print(h[i])elif((i-(n/2-1))%2==0):#n为偶数取值h[i]=0elif((i-(n/2-1))%2==1):#n为奇数取值h[i]=-1/(pow(math.pi*(i-(n/2-1))*d,2))return hdef shiyu(Array,filter,n):#进行时域滤波,对r做操作,n是phi分的角度g=np.zeros((n,2*n-1))#角度for i in range(n):g[i,:]=np.convolve(Array[i,:],filter,'full')#每一行(每一个phi)进行时域滤波(卷积)return g#线性插值连续化
def Continue_Line(xn, yn):def line(x):index = -1#找出x所在的区间for i in range(1, len(xn)):if x <= xn[i]:index = i-1breakelse:i += 1if index == -1:return -100#插值result = (x-xn[index+1])*yn[index]/float((xn[index]-xn[index+1])) + (x-xn[index])*yn[index+1]/float((xn[index+1]-xn[index]))return resultreturn linedef Continue_Figure(figure,n,phi,xr):#xr是r的取值,phi是取滤波后图像的第几行xn = [i for i in range(n)]yn = figure[phi,:]lin = Continue_Line(xn, yn)y=lin(xr+5)#挪一下相对应return ydef transport_xy(figure_filtited,a,b,xn,yn,n,phi_n):theta=math.pi-phi[phi_n]xr=np.linspace(-a,a,xn)yr=np.linspace(-b,b,yn)value=np.zeros((xn,yn))for i in range(xn):for j in range(yn):r=xr[i]*math.cos(theta)+yr[j]*math.sin(theta)value[i][j]=Continue_Figure(figure_filtited,n,phi_n,r)return valuedef transport_phi(figure_filtited,a,b,xn,yn,n):Result_Figure=np.zeros((xn,yn))for i in range(n):Result_Figure=Result_Figure+transport_xy(figure_filtited,a,b,xn,yn,n,i)return Result_Figurey=R_L(500)#卷积核长度
# # 进行离散化求值
c=shiyu(m,y,500)[:,250:750]
# #求x,y坐标下的灰度值Resultxy=transport_phi(c,6,4,60,40,500)# print(Result_Figure[i])
plt.imshow(Resultxy,cmap='gray')
plt.axis('off')
plt.show()
# print(Continue_Figure(c,500,0,-1))
CBP(卷积反投影)实现相关推荐
- OpenCV2:开头篇 介绍
一.OpenCV简介 OpenCV所有的类和函数都在cv命名空间里面,可以用 using namespace cv; #include "opencv2/opencv.hpp" 1 ...
- 如何使用计算机模拟函数图像,模拟图像
本词条缺少概述图,补充相关内容使词条更完整,还能快速升级,赶紧来编辑吧! 又称连续图像,是指在二维坐标系中连续变化的图像,即图像的像点是无限稠密的,同时具有灰度值(即图像从暗到亮的变化值).连续图像的 ...
- 【遥感数字图像处理】基础知识:第一章 绪论
第一章 绪 论 ◆ 课程学习要求 主要教学内容:遥感数字图像处理的概念和基础知识,遥感数字图像的几何处理,遥感图像的辐射校正,遥感数字图像的增强处理,遥感图像的计算机分类,遥感数字图像的分析方法, ...
- 计算机系统结构相关的论文,计算机系统结构毕业论文题目.doc
计算机系统结构毕业论文题目计算机系统结构毕业论文题目 毕业论文(设计) 题 目 学 院 学 院 专 业 学生姓名 学 号 年级 级 指导教师 毕业教务处制表毕业 二〇一五年 九月毕业二十 日 一.论文 ...
- matlab 画图 断层显示,MATLAB编程实现连续断层工业CT图像的三维重建_张爱东
第26卷 第4期核电子学与探测技术 V ol.26 N o.4 2006年 7月Nuclear Electr onics &Detection T echnolo gy July 2006 M ...
- 计算机图像处理怎么学,计算机图像处理在全息学中的应用
[摘 要]全息技术的发现是物理学的一项重大突破,尤其是计算机技术的普遍应用为全息技术的发展提供了一定的平台.在全息技术中包含着两个重要的部分,一是CCD,一是计算机的图像处理技术.因此,本文主要通过对 ...
- 数字图像处理(Digital Image Processing)
数字图像处理(Digital Image Processing) 是指用计算机处理图像,主要包括: (1)点运算:针对图像的像素进行基本数学运算.点运算可以有效的改变图像的直方图分布,可以有效提高图像 ...
- 传统基本图像处理方法:图像增强(灰度变换、直方图增强、空间域滤波、频率域滤波)、图像分割、图像配准等
图像处理设计主要有以下几种处理:图像增强(灰度变换.直方图增强.空间域滤波.频率域滤波).图像分割.图像配准等等. 图像增强: 图像增强作为基本的图像处理技术,目的在于通过对图像进行加工使其比原始图像 ...
- 计算机图像处理数据 流行病学,漫谈计算机图像处理在全息学中的应用.docx
漫谈计算机图像处理在全息学中的应用 摘要全息技术是物理学中的重大发现,近年来在各个行业得到广泛的应用. 作为全息技术中的两个重要部分,和计算机图像处理技术,在推动数字全息新一轮发展中起到至关重要的作用 ...
最新文章
- Java从SFTP服务器下载文件一
- 网络推广策略之网站文章收录少的时候都是哪些因素导致的?
- 引用可以是void类型吗?
- jquery 手指滑动多半屏_JS拖拽专题(五)——「玩出花儿来」移动端滑动事件的封装...
- 将SQL Server中所有表的列信息显示出来
- java main函数_都知道Java程序的入口方法是main,那你知道为什么是main方法吗?
- Retrofit 注解参数详解
- MongoDB 里面日期查询的问题
- 搜索+回溯问题(DFS\BFS详解)
- 新华智云基于MaxCompute建设媒体大数据开放平台
- oracle安装 插件的执行方法失败_解决 VS Code 中 golang.org 被墙导致的 Go 插件安装失败问题...
- 面试官:ca证书存储在哪的
- 抽奖活动软件 html,APP怎么制作抽奖活动,制作APP抽奖活动有何亮点
- MySQL 数据库命名规范.PDF
- java表格标题栏_java使用poi自定义excel标题头并导出(springmvc+poi)
- 算法---木头切割问题(二分法)
- Unity学习笔记——TimeLine的简单使用方法(一)
- 力扣刷题 DAY_73 回溯
- ArGIS Engine专题(6)之利用GP水文分析工具实现基于DEM的山脊线提取
- css写阴影颜色渐变,css3——阴影(立体感,层次效果),渐变色按钮