CN113536557B - 成像系统中探测器布局的优化方法 - Google Patents
成像系统中探测器布局的优化方法 Download PDFInfo
- Publication number
- CN113536557B CN113536557B CN202110751354.XA CN202110751354A CN113536557B CN 113536557 B CN113536557 B CN 113536557B CN 202110751354 A CN202110751354 A CN 202110751354A CN 113536557 B CN113536557 B CN 113536557B
- Authority
- CN
- China
- Prior art keywords
- image
- time
- detection
- detector
- interest
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000003384 imaging method Methods 0.000 title claims abstract description 42
- 238000001514 detection method Methods 0.000 claims abstract description 81
- 238000005457 optimization Methods 0.000 claims abstract description 43
- 238000011835 investigation Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 26
- 238000007476 Maximum Likelihood Methods 0.000 claims description 11
- 238000001914 filtration Methods 0.000 claims description 11
- 238000004422 calculation algorithm Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 5
- 230000014759 maintenance of location Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 5
- 230000006870 function Effects 0.000 description 32
- 230000000875 corresponding effect Effects 0.000 description 10
- 238000013461 design Methods 0.000 description 6
- 238000012986 modification Methods 0.000 description 4
- 230000004048 modification Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000002591 computed tomography Methods 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000002596 correlated effect Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000004069 differentiation Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 210000001015 abdomen Anatomy 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 210000004185 liver Anatomy 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000002603 single-photon emission computed tomography Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- 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/20—Special algorithmic details
- G06T2207/20024—Filtering details
-
- 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/20—Special algorithmic details
- G06T2207/20092—Interactive image processing based on input by user
- G06T2207/20104—Interactive definition of region of interest [ROI]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computer Hardware Design (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computing Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
Abstract
本公开实施例涉及一种成像系统中探测器布局的优化方法，其方法包括：计算设备获取扫描对象的侦查图像，并接收操作者基于侦查对象确定的感兴趣区域；计算设备基于预先建立的探测设备中各探测模块的扫描时间和服从泊松分布的探测数据的模型，获取侦查图像中感兴趣区域和感兴趣区域关联的区域中图像值的协方差信息，将感兴趣区域及与感兴趣区域关联的区域的图像值的协方差作为目标函数；计算设备对目标函数进行优化处理，获取各探测模块的最优采集时间；根据各探测模块的最优采集时间对探测设备的探测器进行布局优化。通过本申请的方法，降低探测器数目和数据冗余，可以有效的降低成像系统的硬件成本并提高数据处理速度，满足实时重建的需求。
Description
技术领域
本申请属于医学成像技术领域，具体涉及一种成像系统中探测器布局的优化方法。
背景技术
随着高分辨率成像需求的提高，传统的高端医学影像诊断设备，比如电子计算机断层扫描(Computed Tomography，CT)、正电子发射正电子发射断层显像(PositronEmission Tomography，PET)、单光子发射计算机断层成像术(Single-Photon EmissionComputed Tomography，SPECT)等需要使用更小尺寸的探测器单元以提升本征空间分辨率，再结合高精度的系统建模和重建算法，才可以实现对扫描物体的高分辨率成像。同时，高端医学影像诊断设备也在向着更大扫描范围的方向发展，以适应越来越广泛的扫描应用。
探测器单元的小型化和探测器扫描覆盖范围的扩大化不可避免地导致了所需探测器数量的大大增加，然而探测器数目的增加对扫描物体信息的获取和图像质量的提升是存在边际效应的，通常采用的全空间均匀布满探测器的系统设计往往不是最佳选择。成像系统采集数据时存在信息冗余，传统的全空间均匀覆盖探测器的布局方案存在优化空间，也使得节省探测器硬件成本具有可能性。但是目前成像系统很少有成熟的探测器布局的优化方案，即使有些成像系统凭借经验尝试着去掉一些探测器，一方面很难准确评估这种非均匀探测器分布的效果，另一方面这种系统设计也很难达到图像质量和硬件成本的最佳平衡。
发明内容
(一)要解决的技术问题
鉴于现有技术的上述缺点、不足，本申请提供一种成像系统中探测器布局的优化方法。
(二)技术方案
为达到上述目的，本申请采用如下技术方案：
本申请提供一种成像系统中探测器布局的优化方法，所述探测设备包括多个探测模块，多个探测模块成环状分布，所述探测模块包括一个或多个探测器，该方法包括：
S10、计算设备获取扫描对象的侦查图像，所述侦查图像为预先通过数字建模方式建立的能够确定感兴趣区域的图像；
S20、所述计算设备接收操作者基于所述侦查对象确定的感兴趣区域；
S30、所述计算设备基于预先建立的探测设备中各探测模块的扫描时间和服从泊松分布的探测数据的模型，获取侦查图像中感兴趣区域和感兴趣区域关联的区域中图像值的协方差信息，将感兴趣区域及与感兴趣区域关联的区域的图像值的协方差作为目标函数；
S40、所述计算设备对所述目标函数进行优化处理，获取各探测模块的最优采集时间；优化目标为在所有探测模块的采集时间总和为固定值时，所述感兴趣区域的图像值的协方差为最小值；
S50、根据各探测模块的最优采集时间对所述探测设备的探测器进行布局优化。
可选地，当所述探测模块包括一个探测器时，步骤S20包括：
通过以下公式定义所述感兴趣区域：
其中，
为所述感兴趣区域的图像值，zT为所述感兴趣区域的指示向量，/>
为侦查图像。
可选地，当所述探测模块包括一个探测器时，步骤S30包括：
S31、对所述成像系统中的图像采集过程建模：
其中，y表示探测数据的平均值，A＝[Anj]为系统矩阵，xj表示为未知图像的像素，M为数字化图像空间像素的数目，N表示为探测器的数目；
S32、通过后滤波带惩罚项的最大似然法，定义一定时间下、各探测模块的扫描时间分布为非均匀时间分布t＝[t1,t2,…,tN]T时所述成像系统中待优化的探测器扫描得到的扫描对象图像估计：
其中，
为带惩罚项的最大似然法表示的扫描对象图像估计，/>
为经过后滤波的扫描对象图像估计，F表示后滤波函数，L(x,y)表示对数似然函数，R(x)表示标量惩罚函数，β表示权重因子，用来权衡对数似然函数和惩罚函数的重要性，x＝[x1，x2,…,xM]T表示为未知图像；y＝[y1，y2,…,yN]T表示探测到的数据；
S33、以费歇尔信息矩阵来表征最大似然估计的准确性，获取以下公式表示的所述感兴趣区域的图像协方差，将所述图像协方差的最小化作为目标函数：
其中，J(n)(0)为在第n个位置的探测器单位时间测量到的数据对应的费歇尔信息矩阵。
可选地，所述标量惩罚函数为粗糙惩罚方程，如下公式所示：
其中，i，j表示成像系统扫描视野中任意两个离散空间位置，R为惩罚方程的黑塞矩阵，R＝▽2R(x)；对于相邻的像素wij＝1，对于平面对角像素
对于空间对角像素/>
其他位置的像素wij＝0。
可选地，步骤S40中，对所述目标函数进行优化处理的方法为：通过基于梯度的优化算法迭代计算出在所有探测模块的采集时间总和为固定值时成像系统中各探测模块的最优采集时间。
可选地，当所述探测模块包括一个探测器时，步骤S40包括：
S41、通过拉格朗日乘子法对所述目标函数进行转换，得到优化目标函数；
其中，T为各扫描模块的扫描总时间；
S42、迭代更新时间分布，直至所述图像协方差相对于不同探测器时间分布的梯度达到均匀梯度，并将最后一次迭代时的时间分布作为最优采集时间。
可选地，通过以下公式迭代更新时间分布：
其中，h表示迭代次数，Δt是迭代步长。
可选地，根据各探测模块的最优采集时间对所述探测设备的探测器进行布局优化的方法为：
基于预设时间阈值，得到优化后的探测器分布向量：
I＝[I1,I2,…,IN]T
其中，Tthreshold为预设时间阈值，tn为第n个探测器的最优采集时间，In＝1表示探测器保留，In＝0表示探测模器取消；
将所述探测器分布向量作为探测器布局优化的结果。
可选地，步骤S40之后，步骤S50之前还包括：
对所述最优采集时间通过预设的滤波方法进行滤波，所述预设的滤波方法为样条插值、线性插值、双线性插值、多项式插值中的一种或多种。
(三)有益效果
本申请的有益效果是：本发明提出了一种成像系统中探测器布局的优化方法。针对扫描物体的感兴趣区域最小化此区域内图像值估计的方差，得到最优采集时间，并根据最优采集时间中的扫描时间进行探测模块的取舍，达到优化探测器分布的目的。通过优化探测器分布，可降低探测器数目和成像系统采集数据的冗余，从而有效地降低了成像系统的硬件成本；通过降低数据冗余还可提高数据处理速度，满足图像实时重建的需求。
附图说明
本申请借助于以下附图进行描述：
图1为本申请一个实施例中的成像系统中探测器布局的优化方法流程示意图；
图2为本申请另一个实施例中的侦查图像示例图；
图3为本申请另一个实施例中的探测器布局迭代优化流程示意图。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。可以理解的是，以下所描述的具体的实施例仅仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合；为了便于描述，附图中仅示出了与发明相关的部分。
为了在保证图像质量不变的情况下，尽可能的减少探测器数目和降低成本，本发明提出了一种用于成像系统中探测器布局的优化方法，以便对成像系统中探测器布局进行优化。以下通过实施例对本发明作详细描述。
本申请提出了一种成像系统中探测器布局的优化方法，该方法中涉及的探测设备包括多个探测模块，多个探测模块成环状分布，探测模块包括一个或多个探测器。图1为本申请一个实施例中的成像系统中探测器布局的优化方法流程示意图，如图1所示，该方法包括：
S10、计算设备获取扫描对象的侦查图像，侦查图像为预先通过数字建模方式建立的能够确定感兴趣区域的图像；
S20、计算设备接收操作者基于侦查对象确定的感兴趣区域；
S30、计算设备基于预先建立的探测设备中各探测模块的扫描时间和服从泊松分布的探测数据的模型，获取侦查图像中感兴趣区域和感兴趣区域关联的区域中图像值的协方差信息，将感兴趣区域及与感兴趣区域关联的区域的图像值的协方差作为目标函数；
S40、计算设备对目标函数进行优化处理，获取各探测模块的最优采集时间；优化目标为在探测模块的采集时间总和为固定值时，感兴趣区域的图像值的协方差为最小值；
S50、根据各探测模块的最优采集时间对探测设备的探测器进行布局优化。
通过本实施例的成像系统中探测器布局的优化方法，可降低探测器数目和成像系统采集数据的冗余，从而有效地降低了成像系统的硬件成本；通过降低数据冗余还可提高数据处理速度，满足图像实时重建的需求。
实施例一
针对探测模块中仅包括一个探测器的情况，本实施例提出了一种成像系统中探测器布局的优化方法，以下对本实施例的方法进行说明。
本实施例方法包括以下步骤：
S100、计算设备获取扫描对象的侦查图像，侦查图像为预先通过数字建模方式建立的能够确定感兴趣区域的图像。
图2为本申请另一个实施例中的侦查图像示例图，如图2所示，侦查图像(scoutimage)是利用蒙卡或数字建模生成扫描患者的代表图像，作为扫描患者的一个良好近似。
S200、计算设备接收操作者基于侦查对象确定的感兴趣区域。
利用侦查图像，操作者决定并勾划感兴趣区域，这个区域可以通过一个感兴趣区域指示向量来表示：
其中，z表示感兴趣区域指示向量，ROI表示感兴趣区域目。
感兴趣区域既可以是连续区域，比如患者的胸部和下腹部，也可以是特定不相邻部位的集合，比如心脏和肝部。
这样在感兴趣区域内的图像值的估计为：
其中，
为感兴趣区域的图像值，/>
为侦查图像。
S300、计算设备基于预先建立的探测设备中各探测器的扫描时间和服从泊松分布的探测数据的模型，获取侦查图像中感兴趣区域和感兴趣区域关联的区域中图像值的协方差信息，将感兴趣区域及与感兴趣区域关联的区域的图像值的协方差作为目标函数，包括：
S301、通过后滤波带惩罚项的最大似然法表示未知图像。
令x＝[x1,x2,…,xM]T表示为未知的扫描对象图像，y＝[y1,y2，…,yN]T表示探测到的数据，M表示为数字化图像空间像素的数目，N表示为物理探测器的数目。实际探测数据可能还会有其他的测量维度，比如时间分辨率信息等等，因为与探测器扫描时间优化无关，因此不予考虑。未知图像通过后滤波带惩罚项的最大似然法来进行估计：
其中，
为带惩罚项的最大似然法表示的扫描对象图像的估计，/>
为经过后滤波的扫描对象图像估计，F表示后滤波函数，L(x,y)表示对数似然函数(log-likelihood function)，R(x)为标量惩罚函数，β为权重因子，用来权衡对数似然函数和惩罚函数的重要性。
成像系统探测数据服从泊松分布，忽略与未知数无关的项，探测数据的log-likelihood函数表示为：
其中，
表示探测数据的平均值。
PET采集过程可以被建模为如下公式：
其中，A＝[Anj]为系统矩阵。
系统矩阵用数学的形式表达了成像系统扫描视野中离散空间位置j的点源被探测器n探测到的概率，反映了系统的物理特性。
为了简化推导，不考虑散射噪声、随机噪声、衰减校正、归一化校正等因素，公式(5)写成矩阵形式为：
其中，A为为系统矩阵，x为未知图像。
将公式(5)带入到公式(4)，则对数似然函数可以写作：
本实施例标量惩罚函数R(x)使用但不局限于粗糙惩罚方程(roughness penaltyfunction)，其定义为：
其中，i，j表示成像系统扫描视野中任意两个离散空间位置，R为惩罚方程的黑塞矩阵，R＝▽2R(x)；对于相邻的像素wij＝1，对于平面对角像素
对于空间对角像素/>
其他位置的像素wij＝0。
需要说明的是，标量惩罚函数也有其他的定义方法，都可以用在优化算法中，本实施例不对标量惩罚函数做限定。
本实施例中每个探测模块中仅包含一个探测器，由于优化后的成像系统不同探测器扫描时间不均匀，因此探测器对应的扫描时间分布表示为t＝[t1,t2,…,tN]T，其中N为系统中的探测器数目。在非均匀时间分布情况下，测量数据的平均值改写为：
其中，n表示为对应的探测器，
表示为第n个探测器单位时间测量数据的平均值，
表示为第n个探测器在tn时间内测量数据的平均值，An为系统矩阵的子矩阵模块，表示系统中所有空间位置点源被探测器n探测到的概率。
如果扫描时间都相等，即t1＝t2＝…＝tN，则公式(9)回归到公式(5)所示的传统的探测器均匀分布的成像系统模型。
在非均匀时间分布情况下，log-likelihood函数可以改写为：
S302、建立时间分布优化的目标函数。
费歇尔信息矩阵(Fisher Information matrix)用来表征最大似然估计的准确性，其定义为：
公式(11)可以写成矩阵形式：
其中，J(n)(0)为在第n个位置的探测器单位时间测量到的数据对应的费歇尔信息矩阵：
其中，系统矩阵的子矩阵模块An表示为系统中所有空间位置点源被探测器n探测到的概率。
根据费歇尔信息矩阵，后滤波带惩罚项的最大似然法评估的图像协方差矩阵可以近似为：
利用图像协方差矩阵可以确认目标函数的优化指标。将步骤S10获取的侦查图像进行投影，可以作为公式(11)中的单位时间测量数据平均值
的一个良好近似。
此感兴趣区域的方差可以近似为：
S400、计算设备对目标函数进行优化处理，获取各探测器的最优采集时间；优化目标为在探测器的采集时间总和为固定值时，感兴趣区域的图像值的协方差为最小值。
公式(15)给出了针对时间分布t优化的目标函数，其中逆矩阵操作的计算量近似于一次重建的计算量。由于公式的对称性，经过一次重建的时间就可以得到特定时间分布下感兴趣区域内图像值估计的方差。但是最终的目的是寻找最佳的时间组合toptimum＝[toptimum1,toptimum2,…,toptimumN]T，最小化感兴趣区域内的图像值估计方差，而时间分布是连续可变的，由于无法遍历所有可行的时间组合，所需计算量也会增长为不可接受。因此即便有了如公式(15)所示的感兴趣区域内图像值方差的良好估计，仍然需要一个快速优化的算法来有效解决这个最小化问题。
由于扫描时间的总时间通常是固定值，即
因此这个优化问题可以转换为变量约束优化问题。利用拉格朗日乘子法可以将时间优化问题定义为：
其中，λ为拉格朗日乘子。
为了求解公式(16)，有：
公式(17)表明如果感兴趣区域内的图像值估计方差相对于不同探测器的时间分布的梯度均匀，则对应的为最佳扫描时间分布。因此可以利用感兴趣区内的图像值对于扫描时间的梯度：
来指导最佳扫描时间分布的优化过程。
感兴趣区域内的图像值估计针对扫描时间的梯度成为优化任务的关键，因此需要一种有效的方法来对其进行计算。对第n个时间tn引入一个小的扰动δ，则对应的感兴趣区域的放射性活度方差对于第n个时间tn的偏微分为：
其中：
对于小的扰动δ，逆矩阵可以做如下近似：
[J+β·R+δ·J(n)(0)]-1≈[J+β·R]-1-δ·[J+β·R]-1·J(n)(0)·[J+β·R]-1 (21)
将公式(21)带入公式(19)和(20)，可以得到感兴趣区域的方差对于第n个时间tn的偏微分的近似解：
公式(22)中计算量最大的部分是两个逆矩阵操作，因此其计算量近似等于两个图像重建。而矩阵J(n)(0)相乘操作的计算量相比很小，因此最终只需要近似两个重建的时间就可以计算出感兴趣区域内图像值估计的方差对于所有探测器的时间分布的梯度。
根据梯度，可以沿着有效的方向来迭代更新时间分布，最终的目标是达到均匀的梯度。时间更新的原则是，以梯度的平均值为基准，梯度振荡大的探测器表示其扫描时间需要明显的变化，梯度接近平均值的探测器表示其扫描时间仅需要微小的变化，时间更新公式如下：
其中，h表示迭代次数，Δt表示迭代步长。
Δt的选择需要进行权衡，如果Δt选的太大，会造成时间分布随着迭代振荡而无法收敛至最佳解；如果Δt选的太小，会降低收敛速度。时间迭代优化算法的最终目的是使得感兴趣区域内的图像值估计方差针对所有时间分布的梯度均匀，对应的即为最优采集时间。
图3为本申请另一个实施例中的探测器布局迭代优化流程示意图，如图3所示，以均匀时间分布为起始值，通过公式(15)和(22)分别计算感兴趣区域的方差和梯度，若梯度不均匀则通过公式(23)进行时间分布更新，最终得到梯度均匀时的时间分布作为最优采集时间。在试验中通常迭代十次左右可以得到最优采集时间，对应的计算量近似为二十次重建的计算量。结合GPU的应用和其他速度优化算法，优化算法变得可行。
对最优采集时间进行滤波，降低数据噪声的影响。滤波方法通常选择样条插值，但是线性插值，双线性插值、多项式插值等也都可以完成同样的操作。
S500、根据各探测器的最优采集时间对探测设备的探测器进行布局优化。
对时间分布设置阈值，得到优化后的探测器分布向量I＝[I1,I2,…,IN]T，In＝1表示探测器保留，In＝0表示探测器可以取消。采集时间少于阈值Tthreshold认为对应的探测器对于感兴趣区域成像重要性很低，可以取消，高于阈值的即为有用探测器，可以保留。
其中，Tthreshold为时间阈值，tn为第n个探测器的最优采集时间。
需要说明的是，探测器布局的优化效果取决于侦查图像的经典图像，患者的身高和体重会影响采集时间的优化，一般是选择传统身体质量指数(Body Mass Index，BMI)的患者图像。实际分析，患者的BMI对探测器布局优化影响不大，可以通过改变扫描时间阈值来适应。
本实施例先依据经验预先确认扫描物体的感兴趣区域，然后通过基于梯度的优化算法分析感兴趣区域与病人其他扫描区域的相关性，最终目标是在总时间固定的情况下，尽量降低与感兴趣区域不相关区域的扫描时间，增加与感兴趣区域相关性大的区域的扫描时间，最终最小化感兴趣区域内的图像值估计的统计噪声，这意味着感兴趣区域最大化了采集信息，迭代计算出最优采集时间，使得感兴趣区域内的图像值估计的方差最小化，最后以优化过的扫描时间分布为重要性权重来对探测器进行取舍。优化探测器分布，降低探测器数目和数据冗余，可以有效的降低成像系统的硬件成本。此外，系统的可靠性是至关重要的，探测器的减少可以简化电子学和传输线路的设计，也减少了探测器损坏的概率，增加了系统的设计冗余。最后，探测器分布优化降低了数据冗余，可以提高数据处理速度，满足实时重建的需求。
实施例二
实际探测器系统通常是模块化设计，将探测器晶体按一定阵列方式进行排布，共享后端电子信号处理电路，这样既便于探测器的维护，又降低了复杂电路的设计难度。因此，针对探测模块中包括多个探测器的情况，本实施例提出了一种成像系统中探测器布局的优化方法，以下对本实施例的方法进行说明。
当探测模块中包括多个探测器时，同一个探测器模块内所有探测器具有相同的扫描时间，因此本实施例仅仅给出求解思路，对求解过程不再再开具体描述。
针对探测器模块化设计，可以在实施例一的公式(9)中将探测器数目N改为探测器模块数目K，则探测器模块对应的扫描时间扫描时间表示为t＝[t1,t2,…,tK]T。后续推导相同，只是更改了一个变量。另一种求解方法是仍保证所有探测器单元时间可变，但是在迭代过程中，如公式(23)，每次迭代后对相同模块内的探测器时间均匀化，用平均值赋值。
这种基于优化算法的探测模块布局优化方案，在探测模块非均匀分布的情况下保证了感兴趣区域的重建图像质量不受影响，有效的降低成像系统的硬件成本。
应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。
尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。
Claims (9)
1.一种成像系统中探测器布局的优化方法，探测设备包括多个探测模块，多个探测模块成环状分布，所述探测模块包括一个或多个探测器，其特征在于，该方法包括：
S10、计算设备获取扫描对象的侦查图像，所述侦查图像为预先通过数字建模方式建立的能够确定感兴趣区域的图像；
S20、所述计算设备接收操作者基于侦查对象确定的感兴趣区域；
S30、所述计算设备基于预先建立的探测设备中各探测模块的扫描时间和服从泊松分布的探测数据的模型，获取侦查图像中感兴趣区域和感兴趣区域关联的区域中图像值的协方差信息，将感兴趣区域及与感兴趣区域关联的区域的图像值的协方差作为目标函数；
S40、所述计算设备对所述目标函数进行优化处理，获取各探测模块的最优采集时间；优化目标为在所有探测模块的采集时间总和为固定值时，所述感兴趣区域及与感兴趣区域关联的区域的图像值的协方差为最小值；
S50、根据各探测模块的最优采集时间对所述探测设备的探测器进行布局优化。
2.根据权利要求1所述的方法，其特征在于，当所述探测模块包括一个探测器时，步骤S20包括：
通过以下公式定义所述感兴趣区域：
其中，
为所述感兴趣区域的图像值，zT为Z的向量转置，Z为所述感兴趣区域的指示向量，/>
为侦查图像。
3.根据权利要求2所述的方法，其特征在于，当所述探测模块包括一个探测器时，步骤S30包括：
S31、对所述成像系统中的图像采集过程建模：
其中，
表示探测数据的平均值，A＝[Anj]为系统矩阵，xj表示为侦查图像的像素，M为数字化图像空间像素的数目，N表示为探测器的数目；/>
表示为第n个探测器在tn时间内测量数据的平均值，
S32、通过后滤波带惩罚项的最大似然法，定义一定时间下、各探测模块的扫描时间分布为非均匀时间分布t＝[t1,t2,…,tN]T时所述成像系统中待优化的探测器扫描得到的扫描对象图像估计：
其中，
为带惩罚项的最大似然法表示的扫描对象图像估计，/>
为经过后滤波的扫描对象图像估计，F表示后滤波函数，L(x,y)表示对数似然函数，R(x)表示标量惩罚函数，β表示权重因子，用来权衡对数似然函数和惩罚函数的重要性，x＝[x1,x2,…,xM]T表示为未知图像；y＝[y1,y2,…,yN]T表示探测到的数据；
S33、以费歇尔信息矩阵来表征最大似然估计的准确性，获取以下公式表示的所述感兴趣区域的图像协方差，将所述图像协方差
的最小化作为目标函数：
其中，J(n)(0)为在第n个位置的探测器单位时间测量到的数据对应的费歇尔信息矩阵，R为惩罚方程的黑塞矩阵，
4.根据权利要求3所述的方法，其特征在于，所述标量惩罚函数为粗糙惩罚方程，如下公式所示：
其中，i，j表示成像系统扫描视野中任意两个离散空间位置，对于相邻的像素wij＝1，对于平面对角像素
对于空间对角像素/>
其他位置的像素wij＝0。
5.根据权利要求3所述的方法，其特征在于，步骤S40中，对所述目标函数进行优化处理的方法为：通过基于梯度的优化算法迭代计算出在所有探测模块的采集时间总和为固定值时成像系统中各探测模块的最优采集时间。
6.根据权利要求5所述的方法，其特征在于，当所述探测模块包括一个探测器时，步骤S40包括：
S41、通过拉格朗日乘子法对所述目标函数进行转换，得到优化目标函数；
其中，T为各扫描模块的扫描总时间，
λ为拉格朗日乘子，toptimum为最佳的时间组合；
S42、迭代更新时间分布，直至图像协方差相对于不同探测器时间分布的梯度达到均匀梯度，并将最后一次迭代时的时间分布作为最优采集时间。
7.根据权利要求6所述的方法，其特征在于，通过以下公式迭代更新时间分布：
其中，h表示迭代次数，Δt是迭代步长，tn为第n个时间。
8.根据权利要求6所述的方法，其特征在于，根据各探测模块的最优采集时间对所述探测设备的探测器进行布局优化的方法为：
基于预设时间阈值，得到优化后的探测器分布向量：
I＝[I1,I2,…,IN]T
其中，Tthreshold为预设时间阈值，tn为第n个探测器的最优采集时间，In＝1表示探测器保留，In＝0表示探测模器取消；
将所述探测器分布向量作为探测器布局优化的结果。
9.根据权利要求8所述的方法，其特征在于，步骤S40之后，步骤S50之前还包括：
对所述最优采集时间通过预设的滤波方法进行滤波，所述预设的滤波方法为样条插值、线性插值、双线性插值、多项式插值中的一种或多种。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110751354.XA CN113536557B (zh) | 2021-07-02 | 2021-07-02 | 成像系统中探测器布局的优化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110751354.XA CN113536557B (zh) | 2021-07-02 | 2021-07-02 | 成像系统中探测器布局的优化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113536557A CN113536557A (zh) | 2021-10-22 |
| CN113536557B true CN113536557B (zh) | 2023-06-09 |
Family
ID=78126602
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110751354.XA Active CN113536557B (zh) | 2021-07-02 | 2021-07-02 | 成像系统中探测器布局的优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113536557B (zh) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8305575B1 (en) * | 2008-06-23 | 2012-11-06 | Spectral Sciences, Inc. | Adaptive spectral sensor and methods using same |
| CN102803893A (zh) * | 2009-06-04 | 2012-11-28 | 瑞尼斯豪公司 | 视觉测量探针和操作方法 |
| CN109507717A (zh) * | 2018-12-28 | 2019-03-22 | 江苏赛诺格兰医疗科技有限公司 | 一种SiPM探测器架构 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7899253B2 (en) * | 2006-09-08 | 2011-03-01 | Mitsubishi Electric Research Laboratories, Inc. | Detecting moving objects in video by classifying on riemannian manifolds |
| WO2010151809A1 (en) * | 2009-06-26 | 2010-12-29 | University Of Virginia Patent Fondation | Time-domain estimator for image reconstruction |
| US8705831B2 (en) * | 2009-09-24 | 2014-04-22 | Koninklijke Philips N.V. | System and method for generating an image of a region of interest |
| JP5759161B2 (ja) * | 2010-12-16 | 2015-08-05 | キヤノン株式会社 | 物体認識装置、物体認識方法、学習装置、学習方法、プログラム、および情報処理システム |
| US9911208B2 (en) * | 2016-04-11 | 2018-03-06 | Toshiba Medical Systems Corporation | Apparatus and method of iterative image reconstruction using regularization-parameter control |
| US10803350B2 (en) * | 2017-11-30 | 2020-10-13 | Kofax, Inc. | Object detection and image cropping using a multi-detector approach |
-
2021
- 2021-07-02 CN CN202110751354.XA patent/CN113536557B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8305575B1 (en) * | 2008-06-23 | 2012-11-06 | Spectral Sciences, Inc. | Adaptive spectral sensor and methods using same |
| CN102803893A (zh) * | 2009-06-04 | 2012-11-28 | 瑞尼斯豪公司 | 视觉测量探针和操作方法 |
| CN109507717A (zh) * | 2018-12-28 | 2019-03-22 | 江苏赛诺格兰医疗科技有限公司 | 一种SiPM探测器架构 |
Non-Patent Citations (1)
| Title |
|---|
| 利用块间耦合稀疏贝叶斯学习的建筑物布局成像方法;晋良念;冯飞;刘庆华;欧阳缮;;电子与信息学报(第04期);853-859 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113536557A (zh) | 2021-10-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9123098B2 (en) | Medical image processing device and medical image processing method, applying weighted penalty term in iterative approximation | |
| US8058601B2 (en) | Determining a multimodal pixon map for tomographic-image reconstruction | |
| US5440647A (en) | X-ray procedure for removing scattered radiation and enhancing signal-to-noise ratio (SNR) | |
| US9031300B1 (en) | System and method reconstructing a nuclear medicine image using deformed attenuation image | |
| JP2939749B2 (ja) | 画像化方法および装置 | |
| CN101681520B (zh) | Pet局部断层摄影 | |
| JP4792463B2 (ja) | 断層撮影画像における時間的アーチファクトの訂正のためのシステム及び方法 | |
| US20090190807A1 (en) | Reconstruction Stabilizer and Active Vision | |
| US10475215B2 (en) | CBCT image processing method | |
| WO2022000192A1 (zh) | 一种ct图像的构建方法、ct设备以及存储介质 | |
| CN103026379A (zh) | 推算图像噪音水平的方法 | |
| Chun et al. | Noise properties of motion-compensated tomographic image reconstruction methods | |
| US11605162B2 (en) | Systems and methods for determining a fluid and tissue volume estimations using electrical property tomography | |
| Watson | An evaluation of image noise variance for time-of-flight PET | |
| CN107348969A (zh) | 一种pet数据处理方法、系统及pet成像设备 | |
| Zhou et al. | Theoretical analysis and simulation study of a high-resolution zoom-in PET system | |
| CN111161182B (zh) | Mr结构信息约束的非局部均值引导的pet图像部分容积校正方法 | |
| CN112798654A (zh) | 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法 | |
| EP3951434A1 (en) | Estimating background radiation from unknown sources | |
| CN112529977B (zh) | 一种pet图像重建的方法和系统 | |
| CN113536557B (zh) | 成像系统中探测器布局的优化方法 | |
| CN111968192A (zh) | 一种ct图像的构建方法、ct设备以及存储介质 | |
| CN113177991A (zh) | 一种基于计划ct校正cbct中散射伪影的方法 | |
| Marin et al. | Deformable left‐ventricle mesh model for motion‐compensated filtering in cardiac gated SPECT | |
| Mesquita et al. | Choosing the ART relaxation parameter for Clear-PEM 2D image reconstruction |
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 |