CN110400253B - 一种基于双线性插值原理确定发射层析权重矩阵的方法 - Google Patents

一种基于双线性插值原理确定发射层析权重矩阵的方法 Download PDF

Info

Publication number
CN110400253B
CN110400253B CN201910591999.4A CN201910591999A CN110400253B CN 110400253 B CN110400253 B CN 110400253B CN 201910591999 A CN201910591999 A CN 201910591999A CN 110400253 B CN110400253 B CN 110400253B
Authority
CN
China
Prior art keywords
projection
dimensional
determining
weight matrix
coordinate system
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.)
Active
Application number
CN201910591999.4A
Other languages
English (en)
Other versions
CN110400253A (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.)
Xian Technological University
Original Assignee
Xian Technological University
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 Xian Technological University filed Critical Xian Technological University
Priority to CN201910591999.4A priority Critical patent/CN110400253B/zh
Publication of CN110400253A publication Critical patent/CN110400253A/zh
Application granted granted Critical
Publication of CN110400253B publication Critical patent/CN110400253B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/06Topological mapping of higher dimensional structures onto lower dimensional surfaces

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)
  • Complex Calculations (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明涉及一种基于双线性插值原理确定发射层析权重矩阵的方法,首先获得发射层析的三维投影模型;然后将投影积分离散化表示;最后利用双线性插值原理确定层析系统的权重矩阵。本发明方法得到了用于准确描述发射层析投影的三维空间二维Radon变换模型,为三维层析重建提供了理论基础;该方法基于双线性插值原理确定权重矩阵的原理简单,计算效率高,可以快速获得权重矩阵数值;其计算结果精确度高,可以减少计算过程中的离散引起的干扰,大大增加计算结果的准确性;本发明方法不仅能够用于三维空间中的二维Radon变换,对于传统的二维层析技术中的权重矩阵计算同样适用。

Description

一种基于双线性插值原理确定发射层析权重矩阵的方法
技术领域
本发明属于光学三维成像与燃烧场三维测量技术领域,更具体的是燃烧场发射光谱层析重建技术,涉及一种基于双线性插值原理确定发射层析权重矩阵的方法。
背景技术
计算层析技术(Computed Tomography,CT)是一种重要的三维成像方法,其理论基础为Radon变换,已广泛应用于医学诊断、工业静态和动态无损检测、地震波和地质构造检测等。在计算层析的各种应用中,燃烧场发射光谱层析技术(Emission SpectrumTomography,EST),简称发射层析技术,是将发射光谱测量和层析原理相结合的一种非接触的新型三维燃烧场诊断测量技术。该技术直接利用CCD相机或光纤探测器接收被测燃烧场自身发射光强度的积分投影信息,并结合计算层析重建算法,可以对燃烧场的三维结构、燃烧组分的三维分布等重要物理参数进行测量。
现有的发射层析技术主要基于二维Radon变换原理,实现的是二维切片层析重建。被测三维物体被分割成一系列相互平行的二维切片,然后对每个二维切片进行单独重建,再将所有重建得到的二维数据进行堆积形成三维重建结果。由于二维层析技术需要单独重建多个二维平面,三维成像的计算效率低、速度慢。此外,实现二维“切片式”层析的前提是不同切片之间的投影相对独立。因此,发射层析系统中各探测器要垂直放置在被测燃烧场周围的同一水平面内,探测器中每一行像素记录不同切片的投影信息。然而,在实际的应用中,系统装置和投影测量都是在三维空间中实现的。由于系统中各相机的安装中不可避免地存在误差,其探测到的投影不再是二维切片面内的Radon变换,而是发射光强度在三维空间中沿投影方向的积分。
以三维Radon变换为基础的三维层析将被测场作为一个整体来考虑,跨越二维切片过程,可实现真正意义上的三维成像,已广泛应用在医学CT中。三维Radon变换对投影矢量确定的二维平面进行积分,其正变换模型和用于重建的傅里叶切片定理、滤波反投影、代数迭代算法等均比较成熟。然而,层析系统中探测器测量的投影不是被测场在某一平面的积分,其实质是三维空间中的二维Radon变换。三维空间中的二维Radon变换与传统的三维Radon变换不同,到目前为止,并没有相应的数学模型对三维空间中的二维Radon变换过程进行描述。
在层析重建算法中,代数迭代算法是层析技术中最常用的重建方法。它将层析投影过程进行离散化表示,用权重矩阵表示被测物体中不同离散网格的物函数对投影信息的贡献量,将投影过程转化为一系列线性方程组,重建过程为利用不同形式的代数迭代方法求解该线性方程组。其中,权重矩阵的计算精度严重影响了层析重建精度和质量。在已有的层析系统中,通过计算投影射线与被测区域离散网格的相交长度确定权重矩阵,其算法复杂,且精度和效率都比较低。
发明内容
针对以上问题,本发明涉及一种基于双线性插值原理确定发射层析权重矩阵的方法,解决现有技术存在的算法复杂,且精度和效率都比较低的问题。
为了达到上述目的,本发明的技术方案如下:
一种基于双线性插值原理确定发射层析权重矩阵的方法,其特征在于,首先获得发射层析的三维投影模型;然后将投影积分离散化表示;最后利用双线性插值原理确定层析系统的权重矩阵。
进一步的,具体包括以下步骤:
步骤一、获得发射层析的三维投影模型:
发射层析的三维投影模型为三维空间中的二维Radon变换,表示为:
Figure BDA0002115351220000021
其中,(x,y,z)为被测物体所处的世界坐标系,(ξ,η,ζ)为投影坐标系,r4、r5、r6、r7、r8、r9为世界坐标系变换至投影坐标系的旋转矩阵中的参数。
步骤二:投影积分离散化:
S1:将投影表示为投影方向的积分;
S2:将积分区间沿投影射线方向均匀划分为N个微元区间;
S3:根据世界坐标系(x,y,z)和投影坐标系(ξ,η,ζ)之间的转换关系,确定微元在世界坐标系中的位置;
S4:将重建区域划分为X×Y×Z个离散网格,确定微元所在的最小相邻离散网格序号以及该微元与最小相邻离散网格中心的距离;
S5:根据双线性插值原理,由相邻的8个离散网格中心点物函数确定该微元的物函数值,得到每个离散网格对微元的权重因子。
步骤三:对所有的微元进行求和运算,并根据投影模型确定的微元与投影的关系,确定离散网格对投影的权重矩阵。
本发明的有益效果如下:
1.本发明方法得到了用于准确描述发射层析投影的三维空间二维Radon变换模型,为三维层析重建提供了理论基础;该方法基于双线性插值原理确定权重矩阵的原理简单,计算效率高,可以快速获得权重矩阵数值;其计算结果精确度高,可以减少计算过程中的离散引起的干扰,大大增加计算结果的准确性;
2.本发明方法不仅能够用于三维空间中的二维Radon变换,对于传统的二维层析技术中的权重矩阵计算同样适用。
附图说明
图1是三维空间中的二维Radon变换投影及坐标系定义;
图2是以二维Radon变换为例的双线性插值计算权重矩阵示意,其中,图2(a)是投影的离散化表示,图2(b)是双线性插值原理。
图3是本发明提供的双线性插值权重矩阵计算方法与传统相交长度方法的计算投影结果对比,其中,图3(a)是采样率为1时的投影对比,图3(b)是采样率为2时的投影对比,图3(c)是采样率为4时的投影对比。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。其中不同实施方式中类似元件采用了相关联的类似的元件标号。在以下的实施方式中,很多细节描述是为了使得本申请能被更好的理解。然而,本领域技术人员可以毫不费力的认识到,其中部分特征在不同情况下是可以省略的,或者可以由其他元件、材料、方法所替代。在某些情况下,本申请相关的一些操作并没有在说明书中显示或者描述,这是为了避免本申请的核心部分被过多的描述所淹没,而对于本领域技术人员而言,详细描述这些相关操作并不是必要的,他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。
本发明方法建立了三维空间中的二维Radon变换数学模型,并以该模型为基础提出了一种基于双线性插值原理确定权重矩阵的方法,包括下述步骤:
一、获得发射层析三维投影模型,具体包括:
参见图1,以原点O为转动定点,通过欧拉角
Figure BDA0002115351220000047
可以将世界坐标系(x,y,z)旋转至投影坐标系(ξ,η,ζ),其中ψ为章动角,θ为旋进角,φ自转角,旋转角度的正方向为从正半轴向原点看时角度为逆时针方向转动。投影方向平行于ξ轴,(η,ζ)平面为投影探测平面。投影过程为待测物体沿投影方向的直线积分,结果为探测平面上的二维数组。三维空间中的二维Radon变换投影及坐标系定义如图1所示。
世界坐标系(x,y,z)和投影坐标系(ξ,η,ζ)的转换关系可以表示为
Figure BDA0002115351220000041
其中旋转矩阵
Figure BDA0002115351220000042
简单起见,将公式(2)表示为
Figure BDA0002115351220000043
三维空间中的二维Radon变换与传统的三维Radon变换不同。传统的三维Radon变换是由投影矢量确定一个平面,在整个平面内进行积分,其结果为一维的数组。而三维空间中的二维Radon变换是指沿平行于投影方向的直线积分,其结果为探测平面上的二维数组。
空间中与投影方向ξ平行,且与投影平面相交于(η,ζ)位置的直线可以表示为:
Figure BDA0002115351220000044
因此,三维空间中的二维Radon变换可以表示为:
Figure BDA0002115351220000045
二、将投影积分离散化,具体为以下几个步骤:
(1)从投影方向可以将三维空间中的二维Radon变换表示为:
Figure BDA0002115351220000046
其中ξmin和ξmax分别为投影射线与被测区域相交的上下边界,其可以由被测区域的顶点位置确定。
(2)将积分区间[ξminmax]均匀划分为N个微元区间,即积分的步数为N,步长为
Figure BDA0002115351220000051
假设每个微元区间的物函数为一个常数,因此
Figure BDA0002115351220000052
(3)根据世界坐标系(x,y,z)和投影坐标系(ξ,η,ζ)之间的转换关系,该积分微元在世界坐标系中的位置可以确定:
Figure BDA0002115351220000053
因此,
Figure BDA0002115351220000054
(4)将重建区域划分为X×Y×Z个离散网格,每个网格的大小为Δg×Δg×Δg。根据双线性插值原理,三维空间中任意位置的物函数值可以由其相邻的离散网格中心点的物函数值确定。具体过程为:
a.确定该微元所在的最小相邻离散网格序号
Figure BDA0002115351220000055
b.确定该微元与最小相邻离散网格中心的距离分别为Δxm、Δym和Δzm
(5)由其相邻的8个离散网格中心点物函数确定该微元的物函数值
Figure BDA0002115351220000061
根据相邻离散网格的距离所确定的系数,可以得到每个离散网格对微元的权重因子。
三、利用双线性插值原理计算层析系统的权重矩阵:
三维空间中的二维Radon变换可以写成以下形式
Figure BDA0002115351220000062
根据公式(13),对所有的微元进行积分运算,并根据投影模型确定微元与投影间的关系,可以得到离散网格对投影的权重矩阵。
实施例:
传统的二维平面Radon变换是三维空间二维Radon变换的一个特例。以二维平面Radon变换为例对本发明的方法和效果进行说明。
二维平面Radon变换投影的坐标系定义如图2所示。当章动角ψ=0,旋进角θ=0时,探测平面垂直于水平面位于被测物体周围,探测平面上每行像素的投影值为相同二维水平面上的投影结果。根据公式(2)可得此时的旋转矩阵为
Figure BDA0002115351220000063
并且投影表达式为
Figure BDA0002115351220000064
根据狄拉克函数的性质,上式可写为
Figure BDA0002115351220000071
即探测器每行像素探测到的投影为同一水平面上物体平面的二维Radon变换,这与传统二维Radon变换的定义相同。验证了本发明中投影模型的正确性。
利用双线性插值方法计算二维Radon变换权重矩阵的过程如图2所示。从投影方向看,以上积分公式可以表示为
Figure BDA0002115351220000072
其中ξmin和ξmax分别为投影射线与被测区域相交的上下边界。将积分区间均匀划分为500个微元区间,每个微元在世界坐标系中的位置为
xm=cosφ(ξmin+mΔξ)-sinφη
ym=sinφ(ξmin+mΔξ)+cosφη
该微元所处的离散网格序号为
Figure BDA0002115351220000073
Figure BDA0002115351220000074
并且该微元与网格中心的距离为Δxm和Δym,如图2(b)所示。利用双线性插值方法确定的该微元的函数值为
Figure BDA0002115351220000075
最终的投影为
Figure BDA0002115351220000076
通过计算正方形在不同角度上的投影,验证了基于双线性插值原理计算权重矩阵的准确性。这个正方形被分成64×64的离散网格,网格大小为1mm,正方形区域中的物函数值都等于1。积分区间被划分为500个微元。图3给出了72°时不同采样率条件下利用双线性插值方法和传统相交长度方法计算的投影值与理论投影值的对比。理论投影值为投影射线与正方形区域的相交长度。图3(a)表明,当采样率R=1,两种方法的计算结果与理论值吻合度较高。然而,随着采样率的增加,传统相交长度方法的误差增大,而双线性插值的权重矩阵计算方法保持了较高的精度,如图3(b)和3(c)所示。该实例说明了基于双线性插值的投影矩阵计算方法具有较高的精度。
以上应用了具体个例对本发明进行阐述,只是用于帮助理解本发明,并不用以限制本发明。对于本发明所属技术领域的技术人员,依据本发明的思想,还可以做出若干简单推演、变形或替换。

Claims (1)

1.一种基于双线性插值原理确定发射层析权重矩阵的方法,其特征在于,首先获得发射层析的三维投影模型;然后将投影积分离散化表示;最后利用双线性插值原理确定层析系统的权重矩阵;
步骤一、获得发射层析的三维投影模型:
发射层析的三维投影模型为三维空间中的二维Radon变换,表示为:
Figure FDA0003998648850000011
其中,(x,y,z)为被测物体所处的世界坐标系,(ξ,η,ζ)为投影坐标系,r4、r5、r6、r7、r8、r9为世界坐标系变换至投影坐标系的旋转矩阵中的参数;
步骤二:投影积分离散化:
S1:将投影表示为投影方向的积分;
S2:将积分区间沿投影射线方向均匀划分为N个微元区间;
S3:根据世界坐标系(x,y,z)和投影坐标系(ξ,η,ζ)之间的转换关系,确定微元在世界坐标系中的位置;
S4:将重建区域划分为X×Y×Z个离散网格,确定微元所在的最小相邻离散网格序号以及该微元与最小相邻离散网格中心的距离;
S5:根据双线性插值原理,由相邻的8个离散网格中心点物函数确定该微元的物函数值,得到每个离散网格对微元的权重因子;
步骤三:对所有的微元进行求和运算,并根据投影模型确定的微元与投影的关系,确定离散网格对投影的权重矩阵。
CN201910591999.4A 2019-07-02 2019-07-02 一种基于双线性插值原理确定发射层析权重矩阵的方法 Active CN110400253B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910591999.4A CN110400253B (zh) 2019-07-02 2019-07-02 一种基于双线性插值原理确定发射层析权重矩阵的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910591999.4A CN110400253B (zh) 2019-07-02 2019-07-02 一种基于双线性插值原理确定发射层析权重矩阵的方法

Publications (2)

Publication Number Publication Date
CN110400253A CN110400253A (zh) 2019-11-01
CN110400253B true CN110400253B (zh) 2023-03-17

Family

ID=68323878

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910591999.4A Active CN110400253B (zh) 2019-07-02 2019-07-02 一种基于双线性插值原理确定发射层析权重矩阵的方法

Country Status (1)

Country Link
CN (1) CN110400253B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111208567B (zh) * 2020-01-07 2020-10-20 中国科学院地理科学与资源研究所 一种矿层成像方法、设备及计算机可读存储介质
CN111950639B (zh) * 2020-08-14 2024-03-19 四川维思模医疗科技有限公司 一种同步实时显示超声与断层解剖图像的成像方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104899827A (zh) * 2015-05-26 2015-09-09 大连理工大学 基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
CN106600687A (zh) * 2016-12-08 2017-04-26 南京理工大学 一种多方向火焰发射层析系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1677253A1 (en) * 2004-12-30 2006-07-05 GSF-Forschungszentrum für Umwelt und Gesundheit GmbH Method and device of reconstructing an (n+1)-dimensional image function from radon data

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104899827A (zh) * 2015-05-26 2015-09-09 大连理工大学 基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
CN106600687A (zh) * 2016-12-08 2017-04-26 南京理工大学 一种多方向火焰发射层析系统

Also Published As

Publication number Publication date
CN110400253A (zh) 2019-11-01

Similar Documents

Publication Publication Date Title
US8000513B2 (en) System and method for 3D time of flight PET forward projection based on an exact axial inverse rebinning relation in fourier space
US8803943B2 (en) Formation apparatus using digital image correlation
Hesselink Optical tomography
CN100587391C (zh) 一种适用于2d-ct扫描系统的投影旋转中心测量方法
EP1031943B1 (en) Efficient cone-beam reconstruction system using circle-and-line orbit data
CN110400253B (zh) 一种基于双线性插值原理确定发射层析权重矩阵的方法
US20110091007A1 (en) System and Method for Image Reconstruction by Using Multi-Sheet Surface Rebinning
WO2003027954B1 (en) Versatile cone-beam imaging apparatus and method
WO1997015841A1 (en) Multi-slice limited projection pet
CN104809750A (zh) 一种直线扫描ct系统及图像重建方法
CN102809494A (zh) 数字x射线成像系统调制传递函数的刀口法测量方法
CN111899344B (zh) 基于相机阵列的火焰发射层析重建装置及方法
CN106228601A (zh) 基于小波变换的多尺度锥束ct图像快速三维重建方法
CN104268938A (zh) 一种三维温度场的重建方法
CN101825433B (zh) 扇束2d-ct扫描系统转台旋转中心偏移量的测量方法
CN102865814A (zh) 植物群体三维重建误差测量方法
CN108364326B (zh) 一种ct成像方法
CN107870159B (zh) 用于可调谐半导体激光吸收光谱的气体浓度二维重建方法
Lin et al. Calibration method of center of rotation under the displaced detector scanning for industrial CT
Yu et al. Interior SPECT—exact and stable ROI reconstruction from uniformly attenuated local projections
EP2453798B1 (en) System and method for image reconstruction by using multi-sheet surface rebinning
CN115797493A (zh) 基于一维系统矩阵稀疏采样的磁场自由线磁粒子成像方法
Cai et al. Three-dimensional combustion diagnostics based on computed tomography of chemiluminescence
Petrou et al. Full tomographic reconstruction of 2D vector fields using discrete integral data
Zhang et al. Design and Research of Low‐Cost and Self‐Adaptive Terrestrial Laser Scanning for Indoor Measurement Based on Adaptive Indoor Measurement Scanning Strategy and Structural Characteristics Point Cloud Segmentation

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