CN106408616B - 一种ct成像中投影图背景不一致的校正方法 - Google Patents
一种ct成像中投影图背景不一致的校正方法 Download PDFInfo
- Publication number
- CN106408616B CN106408616B CN201611046838.XA CN201611046838A CN106408616B CN 106408616 B CN106408616 B CN 106408616B CN 201611046838 A CN201611046838 A CN 201611046838A CN 106408616 B CN106408616 B CN 106408616B
- Authority
- CN
- China
- Prior art keywords
- image
- bright
- coordinate
- value
- perspective
- 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.)
- Expired - Fee Related
Links
- 238000013170 computed tomography imaging Methods 0.000 title claims abstract description 12
- 238000000034 method Methods 0.000 claims abstract description 40
- 238000013519 translation Methods 0.000 claims abstract description 19
- 239000011159 matrix material Substances 0.000 claims description 35
- 230000009466 transformation Effects 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 238000005314 correlation function Methods 0.000 claims description 4
- 238000011426 transformation method Methods 0.000 claims description 3
- 102000002274 Matrix Metalloproteinases Human genes 0.000 claims description 2
- 108010000684 Matrix Metalloproteinases Proteins 0.000 claims description 2
- 238000010219 correlation analysis Methods 0.000 abstract description 2
- 238000001514 detection method Methods 0.000 abstract 1
- 238000002247 constant time method Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 description 1
- 229910052782 aluminium Inorganic materials 0.000 description 1
- 239000004411 aluminium Substances 0.000 description 1
- 238000003705 background correction Methods 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 239000000428 dust Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000010358 mechanical oscillation Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
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
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
Landscapes
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Veterinary Medicine (AREA)
- Heart & Thoracic Surgery (AREA)
- High Energy & Nuclear Physics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Pulmonology (AREA)
- General Physics & Mathematics (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种CT成像中投影图背景不一致的校正方法,具体是一种基于图像灰度信息的CT投影图背景互相关配准法,包括如下步骤:1)在采集样品CT投影图像之前,采集亮场图像一张。确定亮场图像的光斑中心;2)计算亮场图像补偿函数;3)亮场图像照度补偿;4)投影图像照度补偿;5)投影图背景与亮场图背景的实际偏移量计算;6)图像平移配准。本发明基于图像灰度信息,利用互相关原理实现对CT投影图的配准,具有不需要预先知道探测器具体偏转参数,操作简单,计算效率高的优点。
Description
技术领域
本发明涉及一种CT成像中投影图背景不一致的校正方法,具体是一种基于图像灰度信息的CT投影图背景互相关配准法。
背景技术
CT成像技术能够快速、无损地探测物体的三维结构,目前已经成为一种重要的材料三维结构表征手段。CT重构过程中,为了避免探测器上的灰尘、缺陷或电子噪声对CT投影图像的影响,经常在CT实验过程中采集亮场图像和暗场图像,通过对投影数据图像扣除背景以提高CT图像精度。高精度的CT图像重构要求在实验过程中保证X射线源、样品台轴心、探测器三者保持相对静止,从而保证探测器在样品转过不同角度时所采集到的投影图具有相同的背景图像。然而在实际成像过程中,由于射线源漂移或机械振动等原因导致的射线源、样品台轴心及探测器三者之间的相对位移很难完全避免,这使得不同角度采集的CT投影图像背景不一致,从而影响CT投影图像的背景扣除以及后续的CT图像重建。
目前,针对CT实验扫描过程中样品台的移动造成CT重建产生的伪影,常用旋转中心测量校正方法校正。常见的旋转中心测量方法有投影边界法、几何关系法、针模法、中心射线法、对称关系法和迭代收敛法等。针对CT实验扫描过程中探测器偏转造成的CT重建图像伪影,常见校正方法有利用原始投影数据正弦图计算校正算法所必需的几何参数、利用含探测器偏转参数的重建公式直接由原始投影数据重建CT图像。但现有方法需要预先知道探测器偏转参数且存在操作复杂和计算量大的技术问题。
发明内容
本发明针对CT扫描重建过程中,由于投影图背景不一致导致的CT重建图像伪影的校正方法存在操作复杂和计算量大的技术问题,提供了一种CT成像中投影图背景不一致的校正方法,这一方法基于互相关原理,不需要预先知道探测器偏转的具体参数,具有操作简便,计算效率高的优点。
为解决上述技术问题,本发明采用的技术方案为:
一种CT成像中投影图背景不一致的校正方法,包括以下步骤:
1)在采集样品CT投影图像之前,采集一张亮场图像,确定亮场图像的光斑中心O(i0,j0);
2)计算亮场图像补偿函数f(r);
3)亮场图像照度补偿;
4)投影图像照度补偿;
5)投影图背景与亮场图背景的实际偏移量计算;
6)图像平移配准。
所述步骤1)中确定亮场图像的光斑中心O(i0,j0)的具体操作步骤为:将亮场图像每点的灰度数值读出,表示为一个灰度数值矩阵{pi,j}(i=0,1,2...,imax;j=0,1,2,...,jmax),其中i、j分别表示矩阵{pi,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应;imax是图像列方向像素点的最大坐标值,jmax是图像行方向像素点的最大坐标;pi,j表示坐标为(i,j)处的灰度数值;
从亮场图像中光斑边界清晰的图像区域选择第j1,j2,...,jq行,共q行的灰度数据,对每一行灰度数据求其对列坐标i的差分f'j(i):
其中Δi取自然数,其数值要在保证光斑边界处的f'j(i)能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值f'j(i)min和最大值f'j(i)max时,对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标
其中表示取数值的整数部分;计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i0,即:
其中,表示取数值的整数部分;
采用同样的办法得到亮场图像列方向光斑中心坐标j0,则亮场图像光斑中心O点坐标确定为O(i0,j0)。
所述步骤2)中亮场图像补偿函数的计算方法为:以光斑中心O(i0,j0)为圆心,半径为(r+0.5)作圆,其中r=[1,2,...,min(i0,j0)],单位为一个像素,min(i0,j0)表示取i0,j0这两个数值中较小的一个数值。r取不同值时,分别求出圆周上所有像素点灰度值的算术平均值,表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值;根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r),即为亮场图像补偿函数。
所述步骤3)中亮场图像照度补偿的具体操作步骤为:根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿:亮场图像灰度数值矩阵{pi,j}中,坐标为(i,j)的像素点与光源中心O(i0,j0)的距离表示为R,则:
将r=R代入亮场图像补偿函数f(r)中,计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值ai,j:
对亮场图像中所有像素点执行公式(5)的计算,得到补偿后的亮场图像灰度数值矩阵{ai,j},根据{ai,j}利用计算机编程可输出补偿后的亮场图像。
所述步骤4)中投影图像照度补偿的具体操作步骤为:
将每张投影图像每点的灰度数值读出,表示为一个灰度数值矩阵(m=0,1,2...,mmax;n=0,1,2,...,nmax,k=1,2,...,kmax),其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,mmax是图像列方向像素点的最大坐标值,nmax是图像行方向像素点的最大坐标值,k表示投影图的编号,kmax是投影图总的张数;表示第k张投影图上坐标为(m,n)处的灰度数值;确定每张投影图的光斑中心坐标分别表示第k张投影图光斑中心的行坐标和列坐标;利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿,表示补偿后的第k张投影图的灰度数值矩阵,表示照度补偿后的第k张投影图在列坐标为m,行坐标为n处的灰度数值;根据利用计算机编程可输出照度补偿后的投影图像。
所述确定每张投影图的光斑中心坐标的具体操作步骤为:
从第k张投影图像中光斑边界清晰的图像区域选择第n1,n2,...,nz行,共z行的灰度数据,对每一行灰度数据求其对列坐标m的差分
其中Δm取自然数,其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值和最大值时,对应自变量分别记为和根据公式(7)计算得到光斑在第k张投影图像第n行的中心坐标
其中表示取数值的整数部分;计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即:
其中,表示取数值的整数部分;
采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心Ok点坐标确定为
所述对每张投影图的灰度数值进行补偿的具体操作步骤为:
根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿:第k张投影图像灰度数值矩阵中,坐标为(m,n)的像素点与光斑中心的距离表示为D,
则:
将r=D代入亮场图像补偿函数f(r)中,计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值
对第k张图像中所有像素点执行公式(10)的计算,得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出补偿后的第k张投影图像。
所述步骤5)中投影图背景与亮场图背景的实际偏移量计算方法为:以照度补偿后的亮场图像作为模板图像,照度补偿后的第k张投影图像作为待配准图像,在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2,这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χk(e,f),其中e,f分别表示所选图像区域的列坐标和行坐标,e=0,1,2,...,emax,f=0,1,2,...,fmax,emax和fmax分别表示所选图像区域中列方向和行方向的最大坐标值,其数值根据所选图像区域大小确定;根据公式(11)计算互相关函数ck(u,v)的数值,
公式(11)中,u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量,u=0,1,2,...,umax,v=0,1,2,...,vmax,其中umax与vmax分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量,其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定,代表序列ξ(e,f)的算术平均值,代表序列χk(e,f)的算术平均值;当ck(u,v)取到最大值时,对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量;采用相同的方法将所有投影图像背景配准,从而得到每张投影图像背景与亮场图像背景的实际偏移量。
所述步骤6)中图像平移配准的具体操作步骤为:根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换,实现CT投影图像背景与亮场图像背景的配准;根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。
所述投影图像平移变换方法如公式(12)所示,
其中,k表示照度补偿后的第k张投影图,表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。
所述步骤1)中,所选择的总行数q可以根据图像情况选择不同数目,选择的行数越多,计算得到亮场图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。
步骤1)中,在采集样品CT投影图像之前,采集亮场图像一张。从亮场图像中光斑边界清晰的图像区域选择第j1,j2,...,jq行,共q行的灰度数据,对每一行灰度数据求其对列坐标i的差分f'j(i)。亮场图像采集方法以及确定亮场图像中光斑边界清晰的图像区域为本领域常规方法。
步骤2)中,在根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r)。多项式拟合法为本领域常规使用方法。
步骤4)中,从第k张投影图像中光斑边界清晰的图像区域选择第n1,n2,...,nz行,共z行的灰度数据,对每一行灰度数据求其对列坐标m的差分所选择的总行数z可以根据图像情况选择不同数目。选择的行数越多,计算得到投影图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。确定投影图像中光斑边界清晰的图像区域是本领域常规方法。
步骤5)中,在其中umax与vmax分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量,其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定。确定umax与vmax的数值是本领域常规使用方法。
本发明中提到的将图像灰度数值读出表示为一个灰度数值矩阵、将灰度数值矩阵输出为图像均是利用计算机编程实现,是本领域常规使用手段。本发明中提到的对灰度数值矩阵中全部或部分坐标处的灰度值根据不同公式进行计算都是利用计算机编程实现,属于本领域常规使用手段。
本发明采用以上技术方案,与背景技术相比,具有以下优点:
1、不需要预先知道探测器偏转参数;
2、操作简单,计算量小。
附图说明
图1为亮场图像;
图2为图1中横线位置处灰度值;
图3为照度补偿后亮场图像;
图4为图3中横线位置处灰度值分布;
图5为第一张投影图;
图6为图5中横线位置处的灰度值分布;
图7为照度补偿后的第一张投影图;
图8为图7中横线位置处灰度值分布;
图9为背景校准后的第一张投影图;
图10为利用未进行背景配准的CT投影图重建得到的某张CT切片;
图11为利用背景配准后的CT投影图重建得到的某张CT切片。
具体实施方式
本发明所述的CT投影图像背景修正方法可修正不同CT设备获取的不同样品的CT投影图像背景。本实施例所用投影数据是在三英精密仪器有限公司制造的实验室CT机上采集到的。为了验证互相关配准法的有效性,数据采集过程中人为引入探测器抖动。样品为直径2.5mm的圆柱状人造砂岩与直径2.85mm的圆柱状纯铝样品粘合在一起。亮场图像、暗场图像、样品CT投影图的像素均为2048×2048,本次实验共采集投影图像901张。该实施例只是为了进一步阐述本发明,但并不限制本发明所保护的范围。
一种CT成像中投影图背景不一致的校正方法,包括以下步骤:
1)在采集样品CT投影图像之前,采集一张亮场图像,确定亮场图像的光斑中心O(i0,j0):将亮场图像每点的灰度数值读出,表示为一个灰度数值矩阵{pi,j}(i=0,1,2...,imax;j=0,1,2,...,jmax),其中i、j分别表示矩阵{pi,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应;imax是图像列方向像素点的最大坐标值,jmax是图像行方向像素点的最大坐标;pi,j表示坐标为(i,j)处的灰度数值;
从亮场图像中光斑边界清晰的图像区域选择第j1,j2,...,jq行,共q行的灰度数据,对每一行灰度数据求其对列坐标i的差分fj'(i):
其中Δi取自然数,其数值要在保证光斑边界处的fj'(i)能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值fj'(i)min和最大值fj'(i)max时,对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标
其中表示取数值的整数部分;计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i0,即:
其中,表示取数值的整数部分;
采用同样的办法得到亮场图像列方向光斑中心坐标j0,则亮场图像光斑中心O点坐标确定为O(i0,j0)。
这一步中,所选择的总行数q可以根据图像情况选择不同数目,选择的行数越多,计算得到亮场图像在行方向的光斑中心位置越精确。在确定列方向光斑中心时同样如此。
图1所示为CT实验中采集的亮场图像,图像像素为2048×2048,将图像中所有点的灰度值读出可以表示为一个灰度数值矩阵{pi,j}(i=0,1,2...,2047;j=0,1,2,...,2047)。图2为图1中横线位置处灰度值,由图2可看到亮场图像灰度值沿图1中横线方向至左向右明显减小,灰度值分布不均匀。本文选取图像中光斑边界清晰的第100-150行(图中a点与b点之间)共51行的灰度数值数据,对每一行灰度数据求其对列坐标i的差分f'j(i):
以j=100(即第100行)为例,若Δi=1,差分涨落参差不齐,很难找到f'100(i)的极值。若Δi=2,除光斑边界处发生突变外,其他区域f'100(i)的涨落均匀,即选Δi=2较为理想。根据公式(1)计算得到,当差分函数取到最小值f'100(i)min=-200.5时,对应自变量为光斑在第100行的左边界;当差分函数取到最大值f'100(i)max=236.75时,对应自变量为光斑在第100行的右边界。根据公式(2)计算得到光斑在亮场图像第100行的中心坐标利用同样的方法计算出图像中第100行到150行中所有行的中心坐标,根据公式(3)计算得到亮场图像在行方向的光斑中心坐标为1018。采用同样的方法,计算得到亮场图像在列方向的光斑中心坐标为1017,这样确定亮场图像光斑中心O点坐标(1018,1007)。
2)计算亮场图像补偿函数f(r):以光斑中心O(i0,j0)为圆心,半径为(r+0.5)作圆,其中r=[1,2,...,min(i0,j0)],单位为一个像素,min(i0,j0)表示取i0,j0这两个数值中较小的一个数值。r取不同值时,分别求出圆周上所有像素点灰度值的算术平均值,表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值;根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r),即为亮场图像补偿函数。
以图1所示亮场图像中坐标为(1018,1007)处的像素点为圆心,半径分别取r=1.5,2.5,...,1007.5,在图像上做圆,并分别计算出每个半径对应圆周上像素点的平均灰度值利用多项式拟合法得到关于r的函数:f(r)=4.90051r3-0.03556r2+0.0000118r
3)亮场图像照度补偿。根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿:亮场图像灰度数值矩阵{pi,j}中,坐标为(i,j)的像素点与光源中心O(i0,j0)的距离表示为R,则:
将r=R代入亮场图像补偿函数f(r)中,计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值ai,j:
对亮场图像中所有像素点执行公式(5)的计算,得到补偿后的亮场图像灰度数值矩阵{ai,j},根据{ai,j}利用计算机编程可输出补偿后的亮场图像。
根据公式(4)计算出图1中除圆心外,其他所有点到圆心O的距离R,将R的数值代入函数f(r)=4.90051r3-0.03556r2+0.0000118r中,计算出f(R)后,代入公式(5)可计算出照度补偿后(i,j)点处的灰度值ai,j。对图1中所有点进行上述操作后输出补偿后的亮场图像,如图3所示。图4为图3中横线位置处灰度值分布。通过对比图1图2与图3图4可看出照度补偿后亮场图像灰度分布明显变得较为均匀。
4)投影图像照度补偿:将每张投影图像每点的灰度数值读出,表示为一个灰度数值矩阵(m=0,1,2...,mmax;n=0,1,2,...,nmax,k=1,2,...,kmax),其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,mmax是图像列方向像素点的最大坐标值,nmax是图像行方向像素点的最大坐标值,k表示投影图的编号,kmax是投影图总的张数;表示第k张投影图上坐标为(m,n)处的灰度数值;确定每张投影图的光斑中心坐标分别表示第k张投影图光斑中心的行坐标和列坐标;利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿,表示补偿后的第k张投影图的灰度数值矩阵,表示照度补偿后的第k张投影图在列坐标为m,行坐标为n处的灰度数值;根据利用计算机编程可输出照度补偿后的投影图像。
确定每张投影图的光斑中心坐标的具体操作步骤为:从第k张投影图像中光斑边界清晰的图像区域选择第n1,n2,...,nz行,共z行的灰度数据,对每一行灰度数据求其对列坐标m的差分
其中Δm取自然数,其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值和最大值时,对应自变量分别记为和根据公式(7)计算得到光斑在投影图像第n行的中心坐标
其中表示取数值的整数部分;计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即:
其中,表示取数值的整数部分;
采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心Ok点坐标确定为
对每张投影图的照度进行补偿的具体操作步骤为:根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿:第k张投影图像灰度数值矩阵中,坐标为(m,n)的像素点与光斑中心的距离表示为D,则:
将r=D代入亮场图像补偿函数f(r)中,计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值
对第k张投影图像中所有像素点执行公式(10)的计算,得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出照度补偿后的第k张投影图像。
以CT实验中的第一张投影图为例,如图5所示,图6为图5中横线位置处灰度数值分布,由图6可看出图5中沿横线方向灰度数值至左向右明显减小,分布不均匀。将图5中所有点灰度数值读出可表示为一个灰度数值矩阵(m=0,1,2...,2047;n=0,1,2,...,2047),从第1张投影图像中光斑边界清晰的图像区域选择第100,101,...,150行(图中c点和d点之间),共51行的灰度数据,对每一行灰度数据求其对列坐标m的差分
以第100行的数据为例,当Δm=1时,计算得到第100行的差分涨落参差不齐,很难找到的极值。当Δm=2时,除光斑边界处发生突变外,其他区域的涨落均匀,即选Δm=2较为理想。当差分函数取到最小值时,对应自变量为光斑左边界;当差分函数取到最大值时,对应自变量为光斑右边界。根据公式(7)计算得到光斑在投影图像第100行的中心坐标利用同样的方法计算出第一张投影图像中第100行到150行中所有行的中心坐标,根据公式(8)计算得到第一张投影图像在行方向的光斑中心坐标为992。采用同样的方法,计算得到第一张投影图像在列方向的光斑中心坐标为1006,这样确定第一张投影图像光斑中心点坐标(992,1006)。
以第一张投影图为例说明投影图像的照度补偿,根据公式(9)计算出图5中除圆心外,其他所有点到圆心O的距离D,将D的数值代入函数f(r)=4.90051r3-0.03556r2+0.0000118r中,计算出f(D)后,代入公式(10)可计算出照度补偿后(m,n)点处的灰度值对图5中所有点进行上述操作后输出补偿后的第一张投影图像,如图7所示。图8为图7中横线位置处灰度值分布。通过对比图4图5与图7图8可看出经过照度补偿后的第一张投影图像灰度数值分布变得较为均匀。
5)投影图背景与亮场图背景的实际偏移量计算:以照度补偿后的亮场图像作为模板图像,照度补偿后的第k张投影图像作为待配准图像,在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2,这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χk(e,f),其中e,f分别表示所选图像区域的列坐标和行坐标,e=0,1,2,...,emax,f=0,1,2,...,fmax,emax和fmax分别表示所选图像区域中列方向和行方向的最大坐标值,其数值根据所选图像区域大小确定;根据公式(11)计算互相关函数ck(u,v)的数值,
公式(11)中,u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量,u=0,1,2,...,umax,v=0,1,2,...,vmax,其中umax与vmax分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量,其数值可以根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定,代表序列ξ(e,f)的算术平均值,代表序列χk(e,f)的算术平均值;当ck(u,v)取到最大值时,对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量;采用相同的方法将所有投影图像背景配准,从而得到每张投影图像背景与亮场图像背景的实际偏移量。
以照度补偿后的亮场图(图3)为模板,照度补偿后的投影图(图7)为待配准图像,在两张图中相同位置处选取大小相同的区域(图3中的区域1和图7中的区域2)。把这两个区域中的灰度数值读出形成两组离散序列ξ(e,f)和χk(e,f),其中e,f分别表示所选图像区域的列坐标和行坐标,e=0,1,2,...,34,f=0,1,2,...,691,k=1,根据公式(11)计算互相关函数c1(u,v)的数值。
公式(11)中,u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量,其中u=0,1,2,...,umax,v=0,1,2,...,vmax,根据本实施例中X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定投影图与亮场图像在列方向和行方向可能的最大背景偏移量umax=100,vmax=100。代表序列ξ(e,f)的算术平均值,代表序列χ1(e,f)的算术平均值。
通过计算发现,当c1(u,v)取到最大值0.988时,对应的即为亮场图像和第1张投影图像的实际背景偏移量。
6)图像平移配准:根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换,实现CT投影图像背景与亮场图像背景的配准;根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。投影图像平移变换方法如公式(12)所示,
其中,k表示照度补偿后的第k张投影图,表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。
以照度补偿后的第一张投影图像为例,根据亮场图像和第1张投影图像的实际背景偏移量按照公式(12)对照度补偿后的第一张投影图进行平移变换,因为在第一张投影图中总是有所以平移变换公式(12)的具体形式变为:
其中,表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。对投影图进行平移变换后,灰度数值矩阵即为背景校准后的第1张投影图灰度数值矩阵,利用计算机编程输出该图像如图9所示。图10和图11分别为利用未进行背景配准的CT投影图和背景配准后的CT投影图重建得到样品相同位置处CT切片。对比图10和图11可看出,对投影图像背景配准后明显提高了CT切片重构质量。
Claims (4)
1.一种CT成像中投影图背景不一致的校正方法,其特征在于:包括以下步骤:
1)在采集样品CT投影图像之前,采集一张亮场图像,确定亮场图像的光斑中心O(i0,j0);
所述确定亮场图像的光斑中心O(i0,j0)的具体操作步骤为:将亮场图像每点的灰度数值读出,表示为一个灰度数值矩阵{pi,j}(i=0,1,2...,imax;j=0,1,2,...,jmax),其中i、j分别表示矩阵{pi,j}的列坐标和行坐标,与亮场图像上各像素点坐标一一对应;imax是图像列方向像素点的最大坐标值,jmax是图像行方向像素点的最大坐标;pi,j表示坐标为(i,j)处的灰度数值;
从亮场图像中光斑边界清晰的图像区域选择第j1,j2,...,jq行,共q行的灰度数据,对每一行灰度数据求其对列坐标i的差分f′j(i):
其中Δi取自然数,其数值要在保证光斑边界处的f′j(i)能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值f′j(i)min和最大值f′j(i)max时,对应自变量分别记为和根据公式(2)计算得到光斑在亮场图像第j行的中心坐标
其中表示取数值的整数部分;计算出不同j值时的算术平均值作为亮场图像在行方向的光斑中心坐标i0,即:
其中,表示取数值的整数部分;
采用同样的办法得到亮场图像列方向光斑中心坐标j0,则亮场图像光斑中心O点坐标确定为O(i0,j0);2)计算亮场图像补偿函数f(r);
所述亮场图像补偿函数的计算方法为:以光斑中心O(i0,j0)为圆心,半径为(r+0.5)作圆,其中r=[1,2,...,min(i0,j0)],单位为一个像素,min(i0,j0)表示取i0,j0这两个数值中较小的一个数值,r取不同值时,分别求出圆周上所有像素点灰度值的算术平均值,表示半径为(r+0.5)的圆周上所有像素点灰度值的算术平均值;根据不同r值得到的用多项式拟合法得到关于r的函数,表示为f(r),即为亮场图像补偿函数;
3)亮场图像照度补偿;
所述亮场图像照度补偿的具体操作步骤为:根据亮场图像补偿函数f(r)对亮场图像中所有像素点进行照度补偿:亮场图像灰度数值矩阵{pi,j}中,坐标为(i,j)的像素点与光源中心O(i0,j0)的距离表示为R,则:
将r=R代入亮场图像补偿函数f(r)中,计算得到亮场图像补偿后在坐标(i,j)处像素点的灰度值ai,j:
对亮场图像中所有像素点执行公式(5)的计算,得到补偿后的亮场图像灰度数值矩阵{ai,j},根据{ai,j}利用计算机编程可输出补偿后的亮场图像;
4)投影图像照度补偿;
所述投影图像照度补偿的具体操作步骤为:将每张投影图像每点的灰度数值读出,表示为一个灰度数值矩阵其中m、n分别表示矩阵的列坐标和行坐标,与投影图像上各像素点坐标一一对应,mmax是图像列方向像素点的最大坐标值,nmax是图像行方向像素点的最大坐标值,k表示投影图的编号,kmax是投影图总的张数;表示第k张投影图上坐标为(m,n)处的灰度数值;确定每张投影图的光斑中心坐标 分别表示第k张投影图光斑中心的行坐标和列坐标;利用亮场图像补偿函数f(r)对每张投影图的灰度数值进行补偿,表示补偿后的第k张投影图的灰度数值矩阵,表示照度补偿后的第k张投影图在列坐标为m,行坐标为n处的灰度数值;根据利用计算机编程可输出照度补偿后的投影图像;
5)投影图背景与亮场图背景的实际偏移量计算;
所述投影图背景与亮场图背景的实际偏移量计算方法为:以照度补偿后的亮场图像作为模板图像,照度补偿后的第k张投影图像作为待配准图像,在模板图像和待配准图像中的相同位置处选取大小相同的区域1和区域2,这两个区域所包含像素点的灰度值分别作为两组离散序列ξ(e,f)和χk(e,f),其中e,f分别表示所选图像区域的列坐标和行坐标,e=0,1,2,...,emax,f=0,1,2,...,fmax,emax和fmax分别表示所选图像区域中列方向和行方向的最大坐标值,其数值根据所选图像区域大小确定;根据公式(11)计算互相关函数ck(u,v)的数值,
公式(11)中,u、v分别表示亮场图像和投影图像在列方向和行方向可能的背景偏移量,u=0,1,2,...,umax,v=0,1,2,...,vmax,其中umax与vmax分别表示投影图与亮场图像在列方向和行方向可能的最大背景偏移量,其数值根据X射线CT实验中X射线源、样品台、探测器三者之间可能发生的最大偏移量确定,代表序列ξ(e,f)的算术平均值,代表序列χk(e,f)的算术平均值;当ck(u,v)取到最大值时,对应的背景偏移量(分别记为)即为亮场图像和第k张投影图像的实际背景偏移量;采用相同的方法将所有投影图像背景配准,从而得到每张投影图像背景与亮场图像背景的实际偏移量;
6)图像平移配准;
所述图像平移配准的具体操作步骤为:根据照度补偿后每张投影图背景与亮场图像背景的实际偏移量可以将所有投影图像进行平移变换,实现CT投影图像背景与亮场图像背景的配准;根据背景配准后的第k张投影图灰度数值矩阵利用计算机编程可输出背景配准后的投影图。
2.根据权利要求1所述的一种CT成像中投影图背景不一致的校正方法,其特征在于:所述确定每张投影图的光斑中心坐标的具体操作步骤为:
从第k张投影图像中光斑边界清晰的图像区域选择第n1,n2,...,nz行,共z行的灰度数据,对每一行灰度数据求其对列坐标m的差分
其中Δm取自然数,其数值要在保证光斑边界处的能明显区别于图像其他位置的前提下取最小值;当差分函数取到最小值和最大值时,对应自变量分别记为和根据公式(7)计算得到光斑在第k张投影图像第n行的中心坐标
其中表示取数值的整数部分;计算出不同n值时的算术平均值作为第k张投影图像在行方向的光斑中心坐标即:
其中,表示取数值的整数部分;
采用同样的办法得到第k张投影图像列方向光斑中心坐标则第k张投影图像光斑中心Ok点坐标确定为
3.根据权利要求2所述的一种CT成像中投影图背景不一致的校正方法,其特征在于:所述对每张投影图的灰度数值进行补偿的具体操作步骤为:
根据亮场图像补偿函数f(r)对第k张投影图像中所有像素点进行照度补偿:第k张投影图像灰度数值矩阵中,坐标为(m,n)的像素点与光斑中心的距离表示为D,则:
将r=D代入亮场图像补偿函数f(r)中,计算得到第k张投影图像照度补偿后在坐标(m,n)处像素点的灰度值
对第k张图像中所有像素点执行公式(10)的计算,得到补偿后的第k张投影图像灰度数值矩阵根据利用计算机编程可输出补偿后的第k张投影图像。
4.根据权利要求3所述的一种CT成像中投影图背景不一致的校正方法,其特征在于:所述投影图像平移变换方法如公式(12)所示,
其中,k表示照度补偿后的第k张投影图,表示平移变换后图像在列坐标为m,行坐标为n处的灰度数值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201611046838.XA CN106408616B (zh) | 2016-11-23 | 2016-11-23 | 一种ct成像中投影图背景不一致的校正方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201611046838.XA CN106408616B (zh) | 2016-11-23 | 2016-11-23 | 一种ct成像中投影图背景不一致的校正方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106408616A CN106408616A (zh) | 2017-02-15 |
CN106408616B true CN106408616B (zh) | 2019-02-26 |
Family
ID=58082181
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201611046838.XA Expired - Fee Related CN106408616B (zh) | 2016-11-23 | 2016-11-23 | 一种ct成像中投影图背景不一致的校正方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106408616B (zh) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108956658B (zh) * | 2017-05-27 | 2021-10-22 | 上海西门子医疗器械有限公司 | X射线系统和校准x射线管的偏转电流的方法 |
CN108364325B (zh) * | 2018-01-30 | 2021-03-30 | 山西大学 | 规则样品x射线ct投影图位置平移偏差检测及校正方法 |
CN109493293A (zh) * | 2018-10-30 | 2019-03-19 | 京东方科技集团股份有限公司 | 一种图像处理方法及装置、显示设备 |
CN109636874B (zh) * | 2018-12-17 | 2023-05-26 | 浙江科澜信息技术有限公司 | 一种三维模型透视投影方法、系统及相关装置 |
CN109870730B (zh) * | 2018-12-28 | 2020-11-20 | 中国科学院重庆绿色智能技术研究院 | 一种用于x光机图像解析度测试体定检的方法及系统 |
CN110310313B (zh) * | 2019-07-09 | 2021-10-01 | 中国电子科技集团公司第十三研究所 | 一种图像配准方法、图像配准装置及终端 |
CN114757853B (zh) * | 2022-06-13 | 2022-09-09 | 武汉精立电子技术有限公司 | 平场校正函数的获取方法及系统、平场校正方法及系统 |
CN117953022A (zh) * | 2024-03-27 | 2024-04-30 | 杭州邦杰星医疗科技有限公司 | 基于cuda卷积算法的医学影像配准处理系统及方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102542569A (zh) * | 2011-12-21 | 2012-07-04 | 武汉市兑尔科技有限公司 | 一种快速图像配准及其标定方法及实现其的系统 |
CN103116879A (zh) * | 2013-03-15 | 2013-05-22 | 重庆大学 | 一种基于邻域加窗的非局部均值ct成像去噪方法 |
CN103390272A (zh) * | 2013-07-16 | 2013-11-13 | 西安应用光学研究所 | 实现多光谱伪彩色图像的配准融合方法 |
DE102012216272A1 (de) * | 2012-09-13 | 2014-03-13 | Siemens Aktiengesellschaft | Röntgenfokusjustage |
CN105869108A (zh) * | 2016-03-23 | 2016-08-17 | 北京环境特性研究所 | 一种用于动平台移动目标侦测中的图像配准方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB201019694D0 (en) * | 2010-11-19 | 2011-01-05 | Setred As | Camera and visualisation of images and video for autostereoscopic display |
-
2016
- 2016-11-23 CN CN201611046838.XA patent/CN106408616B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102542569A (zh) * | 2011-12-21 | 2012-07-04 | 武汉市兑尔科技有限公司 | 一种快速图像配准及其标定方法及实现其的系统 |
DE102012216272A1 (de) * | 2012-09-13 | 2014-03-13 | Siemens Aktiengesellschaft | Röntgenfokusjustage |
CN103116879A (zh) * | 2013-03-15 | 2013-05-22 | 重庆大学 | 一种基于邻域加窗的非局部均值ct成像去噪方法 |
CN103390272A (zh) * | 2013-07-16 | 2013-11-13 | 西安应用光学研究所 | 实现多光谱伪彩色图像的配准融合方法 |
CN105869108A (zh) * | 2016-03-23 | 2016-08-17 | 北京环境特性研究所 | 一种用于动平台移动目标侦测中的图像配准方法 |
Non-Patent Citations (3)
Title |
---|
A Flat-Field Correction Method for Photon-Counting-Detector-Based Micro-CT;So E. Park 等;《Medical Imaging: Physics of Medical Imaging》;20140331;第9033卷(第90335N期);第1-7页 |
Micro CT投影图像噪声的去除;董歌 等;《医疗卫生装备》;20090228;第30卷(第2期);第8页第3节 |
转台中心偏移与探测器偏转造成的图像伪影分析与校正;陈明 等;《第十一届中国体视学于图像分析学术年会论文集》;20061001;摘要 |
Also Published As
Publication number | Publication date |
---|---|
CN106408616A (zh) | 2017-02-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106408616B (zh) | 一种ct成像中投影图背景不一致的校正方法 | |
CN105849537B (zh) | 校准设备和计算机断层扫描方法 | |
US10481110B2 (en) | Radiographic image generating device | |
US8353628B1 (en) | Method and system for tomographic projection correction | |
KR20120099471A (ko) | 엑스선 이미지들 내의 마커 식별 및 처리 | |
WO2018218611A1 (zh) | 锥形束计算机断层成像系统的几何参数确定方法 | |
KR101850320B1 (ko) | 단층 촬영 장치 및 오차 제거 방법 | |
US7912174B2 (en) | Computed tomography system and method | |
US9253449B2 (en) | Mosaic picture generation | |
KR20120087995A (ko) | 모션이 존재하는 삼차원 공간에서의 엑스선 마커 위치측정 방법 | |
KR20120098807A (ko) | 연속 ct 투영 이미지들 내의 엑스선 마커들을 추적하기 위한 방법 | |
CN106488744A (zh) | X射线拍摄装置以及图像重建方法 | |
CN109358068B (zh) | 一种基于线扫描和环带拼接的大口径平面镜的瑕疵检测装置和方法 | |
CN110621985A (zh) | X线计算机断层装置 | |
CN103218834A (zh) | 一种基于点扩展函数的工业ct图像重建中心定位方法 | |
Kingston et al. | An auto-focus method for generating sharp 3D tomographic images | |
CN112955735B (zh) | X射线相位摄像系统 | |
Illemann et al. | An efficient procedure for traceable dimensional measurements and the characterization of industrial CT systems | |
CN105319225B (zh) | 一种实现板状样品高分辨率大视野cl成像的扫描方法 | |
Zhao et al. | Method for correction of rotation errors in Micro-CT System | |
CN110308165B (zh) | 使用Pi线优化的计算机断层扫描投影的几何对齐、样本移动校正以及强度归一化 | |
CN110292390A (zh) | Ct层敏感度曲线测试体模摆位误差校正方法 | |
Blumensath et al. | Calibration of robotic manipulator systems for cone-beam tomography imaging | |
Franco et al. | Error sources analysis of computed tomography for dimensional metrology: an experimental approach | |
Zemek et al. | Voxel size calibration for high-resolution CT |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190226 Termination date: 20211123 |