CN1096658C - 在不透明介质内部成象的方法和装置 - Google Patents

在不透明介质内部成象的方法和装置 Download PDF

Info

Publication number
CN1096658C
CN1096658C CN96191327A CN96191327A CN1096658C CN 1096658 C CN1096658 C CN 1096658C CN 96191327 A CN96191327 A CN 96191327A CN 96191327 A CN96191327 A CN 96191327A CN 1096658 C CN1096658 C CN 1096658C
Authority
CN
China
Prior art keywords
light source
images
light
detecting device
opaque medium
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 - Lifetime
Application number
CN96191327A
Other languages
English (en)
Other versions
CN1167537A (zh
Inventor
S·B·科拉克
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.)
Koninklijke Philips NV
Original Assignee
Koninklijke Philips Electronics NV
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 Koninklijke Philips Electronics NV filed Critical Koninklijke Philips Electronics NV
Publication of CN1167537A publication Critical patent/CN1167537A/zh
Application granted granted Critical
Publication of CN1096658C publication Critical patent/CN1096658C/zh
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/008Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47Scattering, i.e. diffuse reflection
    • G01N21/4795Scattering, i.e. diffuse reflection spatially resolved investigating of object in scattering medium
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/421Filtered back projection [FBP]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Biochemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Optics & Photonics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

为了生物医学图象的重建,一种不透明介质内部成象的方法和装置被公开。这种方法从得到第1幅图象的一种简单的反透影算法开始,然后,下一步该图象用空间变化卷积函数解卷积而被校正。该卷积函数是基于不透明介质的扩散和或吸收系数及光在不透明介质中的传播。第二次校正处理用坐标变换得到。根据本发明的方法的优点是这种方法相当容易实现,这种方法只需要少量的迭代步骤就能得到校正的图象,且这个方法不会遇到不稳定性问题。

Description

在不透明介质内部成象的方法和装置
发明领域
本发明涉及不透明介质内部成象的方法,本发明也涉及到实现这种方法的装置。
背景技术
这种方法可从国际专利申请WO 91/07655中了解到。在这个专利申请中,所知道的方法用于生物组织(biological tissue)的成象。在该专利申请中所谓不透明介质指的是大量高度扩散物质,如内酯溶液(intralipid solution)或生物组织。该熟知的方法,可在医学诊断中用来成象胸部组织中的肿瘤。在这个熟知的方法中,一种迭代处理用于根据测量强度确定不透明介质内部图象。在所说的迭代处理中,测量的强度同包括光源,不透明介质内部和检测器的仿真模型的计算强度比较。不透明介质在该仿真模型中表示成具有相关联的扩散和/或吸收系数的体元素集合。而且,通过迭代处理使误差函数最小,这个函数由计算的强度和测量的强度确定。体元素的扩散和/或吸收系数的最后值可相继以内部图象的形式表示。所熟知的方法的缺点是这种方法需要许多迭代步骤和相应的许多计算工作量。因而不可能实时成象。
发明概述
本发明的目标之一是提供减少迭代步骤的次数在不透明介质内部成象的一种图象重建算法。本发明的另一个目标是为医学成象处理提出一种可能的实时成象技术,如乳房X照相技术。
本发明的方法执行下列步骤:
a)用来自多个位置上光源的光连续照射不透明介质,并检测通过不透明介质沿着多个光路传播的光,
b)根据所测量的强度确定在不透明介质内部图象。
根据本发明,用于该目的方法其特征在于下面的步骤,
c)从测量的强度中确定光参数,
d)通过反投影(backprojection)确定所定光学参数的内部的第一幅图象,
e)从光学参数确定第一个空间变化卷积函数,
f)通过把第1个空间变化卷积函数同第1幅图象相组合,确定不透明介质内部的第2幅图象。
通过测量的光学参数所定值的反投影,得到第一幅图象。第一幅图象是第二幅图象的模糊后的版本,因为此图象是用正常发生的1/r作用和一个空间变化卷积函数卷积了的图象。正常1/r的作用出现在通过测量强度的反投影的目标重建中。因为在不透明介质内部的扩散和吸收,相对在第一例子中的卷积函数正常1/r的作用可忽略。当在强吸收的情况下,1/r的作用相对卷积函数不可忽略时,1/r的作用通过用一限带滤波器可予以补偿。
本发明基于识别光路分配函数可用于说明通过不透明介质光的传播。此光路分配函数接着用于确定卷积第二幅图象的卷积函数。从光路分配函数可以推出,在均匀的不透明介质中光的最可能的光路主要表现在光源和检测器之间的一个香蕉形空间或筒形空间内。光路分配函数特别从文章”对多散射媒质中的光子迁移路径分布的Monte Carlo模拟(Monte Carlosimulations of photon migration path distributions in multiplescattering media)”S.Feng等人,SPIE vol.1988,PP.78-88.可以了解。光路分配函数的空间变化参数是从所定的光学参数中确定。因为反投影及其组合是线性运算的,反投影及组合的次序可以颠倒。然后,通过所定的强度值形成第一幅图象。而且,在本专利申请中用的光是波长大约在200到1000nm之间的电磁辐射波。
根据本发明方法的实施例其特征在于通过用第一个空间变化卷积函数的第一个图象的解卷积确定第二个图象。
得到第二个图象的一种可能性是在伪时间域或空间域内,用卷积函数去解卷积第一个图象。
根据本发明方法的另一个实施例其特征在于通过对富里叶变换的第一个图象和富里叶变换的第一个空间变化卷积函数之商的富里叶逆变换,确定第二幅图象。在伪时间域内对解卷积的一种方法是在频率或富里叶域中的运算。富里叶域的优点是实现操作相对简单。
根据本发明方法的实施例的特征还在于用迭代处理,通过执行如下步骤得到第(n+1)幅图象:
a)从第n幅图象和所定的光学参数确定第n个空间变化卷积函数,及
b)通过把第n幅图象同第n个空间变化卷积函数的组合,确定目标的第n+1幅图象。
迭代处理执行,其中通过把第n幅图象同第n个空间变化卷积函数的组合,得到目标的第n+1幅改进图象。为了得到高质量的图象,二次迭代通常就足够。而在所引用的专利申请WO 91/07655中说明的那样,所熟知的方法需要多次迭代。需要很多迭代次数是因为每次迭代吸收和扩散系数只有很少的改变,否则出现稳定性问题,。本发明方法的优点是在迭代处理中没有这种不稳定出现。第n+1幅改进图象可用富里叶变换的第n幅图象和被富里叶变换的第n个卷积之商的富里叶逆变换确定。另外,第n+1幅改进图象也能够在伪时间域或空间域中,用第n个卷积函数去解卷积第n幅图象确定。
根据本发明方法的实施例其特征在于在这个方法中包括非线性坐标变换步骤。当不透明介质的衰减系数和浸入目标的衰减系数的差很大时,可用这种校正步骤。这种校正步骤补偿随着目标距离变化引起的光线弯曲。显然,这些校正从A.Yariv,Quantum Electronics,p107,1975及H.Schomberg,再现性超声层析X射线摄影的改进方法(An improvedapproach to reconstructive ultrasound tomography),Journal ofPhysics D.,Vol 11,page,1978中可以知道。这个校正是基于几何光学的镜象方程。
根据本发明方法的实施例其特征在于光学参数是一些衰减系数。使用这些衰减系数的优点是这些衰减系数是容易从测量强度和所知道的光路最短长度中得到。也可不用衰减系数,而用其他光学参数如吸收系数,扩散系数或散射系数。
根据本发明方法的实施例其特征在于这个方法包括降低高频噪声分量的步骤。在测量处理和反投影处理中,引入了高频噪声分量。这些高频分量可以通过低通滤波器的滤波而降低。滤波包括如通过低通滤波器频率响应和第一幅图象的多重富里叶频谱的乘积或根据低通滤波器二维脉冲响应和该图象的卷积。
根据本发明方法的实施例其特征在于这些光路在多个光源位置上的一个光源和在多个检测器位置上的一个检测器之间形成。在这种方法中,不透明介质被一个光源的光照射,用很多检测器位置上的一个检测器执行一定范围测量后,光源依次放置在很多光源位置上。
根据本发明方法的实施例其特征在于光路在位于多个位置上的多个光源中选择的一个光源和在多个位置上的多个检测器中选择的一个检测器之间形成。该方法的一个优点是所有的光源和检测器能够相对不透明介质处在同样的固定机械位置上,所以不需要机械定位装置。
根据本发明方法的实施例其特征在于这些光路是在一个位置的一个光源和通过多个位置上的多个检测器检测的光之间形成的。这个方法的优点是能够减少光源的数目。用于使不透明介质相对于光源和检测器定位的机械装置是需要的。
根据本发明方法的实施例其特征在于这些光路在多个位置上的多个光源中的一个选择的光源和在一个位置上的一个检测器之间形成。
本发明也涉及到一种装置,这种装置提供:
a)用于照射不透明介质的光源,
b)测量通过不透明介质传播的光的光强度的检测器,
c)从多个不同的位置照射不透明介质的装置,
d)在多个不同的位置上定位检测器的装置,
e)为装置产生控制信号的一个控制单元,以在多个位置上定位光源和检测器,及
f)根据测量的强度重建图象的重建单元,
特征在于准备好执行下列步骤:
a)从测量的光强度中确定衰减系数,
b)根据反投影,从子步a)中所定的衰减系数确定目标的第一个图象,
c)从衰减系数确定第1个空间变化卷积函数,
d)把第1幅图象同第1个空间变化卷积函数相组合,确定第二幅图象。
本发明的上述及其它的更详细的内容将通过例子参考附图在下面说明。
附图说明
在附图中:
图1显示了包括一个测量光源和一个检测器的第一个装置,
图2显示了在不透明介质中的测量光源和检测器之间的一个香蕉形光路分布函数,
图3显示了包括一个测量光源和一个检测器的第二个装置,
图4显示了在二个测量光源和检测器对中的可能光路,
图5显示了包括一个测量光源和一个检测器的第三个装置,
图6显示了包括一个测量光源和很多检测器的第四个装置,
图7显示了包括很多测量光源和很多检测器的第五个装置,
图8显示了在非线性坐标转换中所用的近似函数。
发明详述
图1所示的测量装置100用于在测量空间108内,对不透明介质实现测量。测量装置100包括一个由第一个光波导105的输出窗口104形成的测量光源103,光波导105光耦合到光源120,例如一产生670nm的波长及30mw的输出功率的激光二极管。测量光源103被固定在可移动座102上。可移动座被装在平台101上。检测器106也安排在可移动座102上的测量光源103的对面。检测器106由光耦合到光检测器121的第二光波导113的输入窗口107形成。光检测器121例如是光敏PIN二极管或光电倍增管。放置不透明介质或生物组织的测量空间108是在测量光源103和检测器106之间。测量空间108的直径是例如100mm。光波导可以是如光纤维,而且一种选择方案可能通过配有如透镜及镜面的几何光学系统形成。可移动座102根据耦合到控制单元122上的马达部件112的装置,对很多测量位置设定。测量位置的数目,如在16到256之间。为了实现强度测量,控制单元122对光源120将产生一个控制信号,以使在光源120中产生光。光源的光通过第一光波导105和输出窗口传到测量空间108中的不透明介质中。不透明介质包括例如1%的酯内溶液(intralipid solution),在此溶液中浸有5mm半径的Delrin圆柱。检测器106接收在测量空间108中通过不透明介质传送的部分光,且把接收的光通过第二光波导113传到光检测器量121。光检测器121把接收的光转成检测信号。检测信号经电导体127传给模数转换器123。在模数转换器中检测器信号转换成数字值I。数字强度值I接着通过数据总线128传到重建单元125。数据处理发生在重建单元125中。数据处理的结果以图象的形式显示在监示器126上,或存储在数据存储单元124中。
根据本发明方法的实施例,在M位置上的强度通过在第一步中检测器装置在第一个方向上,例如在正交坐标系统的X方向被测量,且存储在第一个数据表X(i)中。正交坐标系统的原点位于测量空间108的中心。在第二步中,强度在第二个方向的N位置上,例如在正交坐标系统的Y方向上被测量。在每个方向上的位置数如在16到256之间。实际上经常选择32。为了在第二个方向进行测量,移动座相对于出现在测量光源103和检测器106之间的测量空间108中的不透明介质旋转大约90度。第二个方向上的测量强度值通过重建单元125被接着存储在第二个表Y(j)中。而衰减系数Kx(i)的第一个表和衰减系数Ky(j)的第二个表根据下面的关系从测量第一个强度表X(i)和测量的第二个强度表Y(j)得到 ΔK = 1 L 1 n ( I I I 0 ) K , ΔK = K 0 - K i - - - - ( 1 )
这里的L表示测量光源和检测器之间的距离,Ii为第i检测器位置上的强度值,Io为参考位置上的测量的强度值,Ko为参考位置上的背景衰减系数,K1*是第i个检测位置上的衰减系数。
也可以用其他光学参数如吸收系数,扩散系数或散射距离代替衰减系数。此外,用不同波长执行二次连续的测量,也可以确定得到所述的各种光学参数的差。
根据本发明方法的下一步,不透明介质内部的第一幅图象用反投影从第一个表Kx(i)和第二个表Ky(j)计算。该技术例如从计算机层析X射线摄影术中了解,且“特别在图像处理基础”A.K.Jain,Prentice Hall,1989,pp 439-441中已说明。而且,第一幅图象用M*N矩阵描述。矩阵Ko(i,j)的元素(i,j)的值由第一个表Kx(i)中的第i个元素及第二个表Ky(j)中的第j个元素的和确定,因此Ko(i,j)=Kx(i)+Ky(j)。在归一化以后,第一幅图象被显示在监示器126上。下一步确定作为第一幅图象的起点的表Kx(i)和Ky(j)的衰减系数。
从文章“对多散射媒体中的光子迁移路径分布的Monte Carlo模拟”S.Feng等人,SPIE vol.1988,PP.78-88.可以了解到在均匀的不透明介质中光路分配函数相应于测量光源103和检测器106之间的一个香蕉形空间或筒形空间。在不透明介质中光源和检测器之间由于在位置r上点状吸收颗粒,联系着位置d上检测器的光路分配函数φ(r)如下给出 φ ( r ) = - a S 0 4 πD e - K | rd - r | | rd - r | e - K | r - rs | | r - rs | - - - - ( 2 )
这里的r是从一个目标的随机位置到极坐标系统的原点的距离,rs是从光源到坐标系统原点的距离,rd是从检测器到坐标系统的原点的距离,K是衰减系数。对类组织(tissue-like)介质的衰减系数大约为 K = 1 L a = μ a D - 3 μ a μ ′ s
这里μa是吸收系数,μs′是扩散系数,La是散射光穿透深度,D是不透明介质的散射系数。有关测量装置参考图1所述,三维空间降为二维空间。此外,当大致上认为是均匀的不透明介质时,在测量装置100中使用的测量光源103和检测器106之间的距离对所有的位置是相同的。光路的分配函数(2)适用于原点位于测量空间的中心的正交坐标系统X,Y,Z,简化为表达式 p ( y ) = exp [ - 2 K ( ( l 2 ) 2 + y 2 ) ] ( l 2 ) 2 + y 2 - - - - ( 3 )
这里的l是光源和检测器之间的距离,y是从垂直于线l并通过点l/2的线上的点y的距离,K是衰减系数。该表达式表示在所给的光源和检测器之间的光路分配函数的扩展,光源和检测器之间给定一个固定距离,同时该表达式给出卷积函数的形状,用该卷积函数,目标即第二幅图象被卷积。一个香蕉形光路分布函数的例子示于图2。
图2显示了衰减系数取值0.05mm-1,源检测器距离取值100mm的均匀不透明介质内,值为P(0.5)的光路分配函数的外形。值P(0.5)表示通过源S在不透明介质中注入总光子数的一半将通过不透明介质沿香蕉形曲线的封闭区的通路传播。对较高衰减系数的不透明介质,这个香蕉形状较窄,而对较小衰减系数的不透明介质,这个香蕉形状较宽。
所用的香蕉形强度函数(3)的衰减系数接着通过所有的元素K(i)和K(j)的平均值的计算,从所确定的衰减系数Kx(i)和Ky(j)估计。因此 K = 1 m + n ( Σ i = 1 m K x ( i ) + Σ j = 1 n K y ( j ) ) - - - - ( 4 )
对于高度非线性目标,用衰减系数Kx(i)和Ky(j)的最大值。高度非线性目标是这样一些目标,它具有的衰减系数接近50倍的不透明介质的衰减系数,其中,该目标被浸入。
在本发明的方法中,第一个空间变化卷积函数在第一幅图象通过衰减值反投影后被估计。在此可认识目标的重心从第一幅图象估算或从既定的知识中可得到。使重心例如处在位置(0,-Sobj),然后第一个卷积函数用G1(xi,yj)=F1(xi)·F2(yI)给出,这里F1(xi)对于
Figure C9619132700112
由等式(2)表示和F2(yj)对于
由等式(2)表示
这里xi是沿X轴的坐标,yi是沿Y轴的坐标,K确定衰减系数(4),l是光源检测器的距离,-Sobj是沿负y轴在可认识的目标和原点O之间的距离。这第一个卷积函数在原理上取决于目标到光源的距离和目标到检测器的距离。因此,它是空间变化卷积函数。现在,第1幅图象的富里叶变换等于第2幅图象的富里叶变换和第一个卷积函数的富里叶变换乘积,所以F{K1(i,j)}=F{K2(i,j)}.F{G1(i,j)}。因而富里叶变换的第2幅图象等于:F{K2(i,j)}=F{K1(i,j)}.F{G1(i,j)}.第2幅图象的重建也可以在另外的状态下执行,例如用卷积函数图象在伪时间域和空间域内解卷积。
确定第一幅图象的快速方法可以例如在感兴趣的区域不靠近光源或检测器的情况下得到。此时下面各步在投影域内执行。
1.从第一幅图象确定第一表K′x(i)和第二投影表K′y(j),确定第一个和第二个投影表的富里叶变换如F{K′x(i)}和F{K′y(j)},确定一维卷积函数的富里叶变换F{G1(i)},接着
2.富里叶变换的数据表的乘法如下:
{Kx2(i)}=F{Kx1(i)}F/F{F1(i)}                    (4)
{Ky2(i)}=F{Ky1(i)}F/F{F2(i)}                    (5)
接着
3.通过表Kx2(i)和表Ky2(i)的反投影确定第二幅图象K2(i,j)通过在步2和步3之间插入的噪声减少步骤,在第一图象中出现的任何噪声可以被减少。在噪声减少步中,二个表Kx2(i)和Ky2(j)乘以降噪滤波器,如低通滤波器的频率响应H(i)。用反投影和归一化的重建后,图象能显示在监视器126上。
按照本发明的方法的图象的进一步改进可通过非线性坐标变换和利用先验知识的下一个解卷积步骤得到,例如当浸入目标的衰减系数和周围不透明介质的衰减系数之差如超过50倍时。非线性坐标变换将参考图8予以解释。
在图8中,线800显示了在不透明介质中沿一个浸入目标图象K2(i,j)的中心线表示的衰减值。中心线800沿第一个坐标轴,如X轴的方向。沿中心线的衰减值用函数S(x)表示。然后得到函数S(x)对x变量的二阶导数S2(x)。二阶导数S2(x)在图8中用线801表示。然后而在第一个坐标轴上实现非线性坐标变换。这里,在位置x和x+Δa二个相邻的测量位置之间的空间间隔Δa可以从总的测量长度L除以测量的投影数M得到,所以
Figure C9619132700131
新的间隔Δa’,一个非均匀的间隔现在用
Δ a ′ = 1 1 + S ′ 2 ( x ) , x ∈ [ L 2 - W 2 , L 2 + W 2 ] , 确定,这里的S2′(x)表示S2(x)的近似函数。在此近似函数S2′用余弦函数的一个半周期给出。 S 2 ′ ( x ) = A cos ( 3 π W ( x - L 2 + W 2 ) + π 2 ) ,这里W表示曲线800的最大值的一半处的全宽度。在图8中虚线802表示近似函数S2′(x)。在非线性坐标变换后,一个解卷积步骤用已经从方程(3)得到的卷积函数实现以确定第三幅改进的图象。
根据本发明方法的实施例可用的另一安排参考测量装置予以解释,测量装置中测量光源和检测器之间距离的变化。这种安排显示在图3中。图3显示了测量光源303被固定在第一环309上,而检测器306置于第二个环310上。第一环309和第二个环310彼此相对旋转,且被固定到平台301上。而且,在第一环309上的测量光源303和在第二个环310上的检测器彼此相对这样固定,使得测量光源的输出窗口304,和检测器306的输入窗口307基本上处于同一平面。测量光源303和检测器306的位置通过马达部件312独立调整。马达部件312机械上连结到第一环309和第二环310上。第一环309和第二环310通过连结到控制单元322上的马达部件312,对很多的测量位置进行设定。测量位置数例如可以在16到256之间。测量空间308在第一环309和第二环310的中心。在测量空间308中可以放置如包含目标或生物组织的不透明介质。测量空间的直径例如是100mm。为了实现强度测量,控制单元322为光源320产生控制信号,以使在光源320中产生光。光源的光通过第一个光波导305和输出窗口304传到测量空间308的不透明介质中。检测器306接收一部分通过测量空间108中的不透明介质传送的光,且接收的光通过第二光波导313传向光检测器321。光检测器321把接收的光转成检测器信号。检测器信号经过电导器327传到模数转换器323。在模数转换器中检测器信号被转换成数值I。数字强度检测值I接着通过数据总线328送到重建单元325。在重建单元325发生数据处理。数据处理的结果以图象的形式显示在监示器326上或存储在数据存储单元324中。接着的步骤是根据本发明方法的第二实施例在不透明介质中实现成象目标。在第一步,通过不透明介质或生物组织传输的光的强度对测量光源303的第一M号位置同在第二N号位置上检测器306连续测量。这以后,测量光源303借助第一环309和马达部件312在选择位置上首先设定。然后,通过控制单元322和第二环310,对每个测量连续地以360/N角度旋转检测器位置。在N次测量后,第一环309作360/M度旋转,再执行M次测量。所测量的强度值通过重建单元325存储在表Rm(i)中。在执行了M*N次测量后,可得到每个N元素的M个数据表Rm。数据表接着通过下面的关系式转成衰减系数表Rm(Ki)。 ΔK = 1 l 1 n ( I I I 0 ) - - - - ( 6 )
下一步,通过扇形光束反投影从衰减表Rm(Ki)确定第一幅图象K1(i,j)。扇形光束反投影特别从“图像处理基础”A.K.KX,PrenticeHall,1989,pp 464-465中可了解。下一步,从所确定的衰减系数表Rm(Ki)中确定第一个卷积函数,不过由于测量光源303和检测器306之间的距离因检测器306相对于测量光源303的各种位置而不同,香蕉形光路分配函数的宽度可以因测量光源303和检测器306的各个位置而不同。根据本发明的方法中第一个卷积函数G1(i,j)现在被计算,其中的第一个卷积函数G1(i,j)的一个元素从起点计算,在M测量光源位置和N检测器位置之间对每个元素(i,j)属于光路Ls,d的衰减系数乘以光将要通过的包括源和检测器位置的正交二维空间中点i,j的概率,其关系为 K ′ ( i , j ) = Σ q , r = 1 M , N k ( q , r ) P q , r ( i , j ) Σ q , r = 1 M , N P q , r ( i , j ) - - - - ( 7 )
其中Pq,r(i,j)是测量光源位置Sq和检测器位置Dr之间的光路分配函数,这里的光路分配函数受位置i,j处微小的吸收颗粒影响。这里的Pq,r(i,j)可从方程(2)得到。关于元素i,j衰减系数的计算参见图4予以说明。
图4显示了在位置Sq1和Sq2的二个测量光源及在位置Dr1,Dr2上的二个检测器,检测器和测量光源的位置均置于圆周400中。而且点Pi,j的阵列,例如5*5阵列401是位于圆周400内。阵列401的每个点i,j对应于卷积阵列I的元素i,j。对应于第一卷积函数的元素i,j的衰减系数,现在用所有可能的测量光源和检测器光路之间光路分配衰减系数乘积乘以光路将要通过点Pi,j的概率给出。图4显示了在二个测量光源Sq1,Qq2和二个检测器Dr1,Dr2之间四个最短可能的光路402,403,404,405。对于二个在位置Sq1,Sq2的测量光源和位于Dr1,Dr2之间的二个检测器,衰减系数的计算结果为 K ( i , j ) = K q 1 , r 1 P q 1 , r 1 ( i , j ) + K q 2 , r 2 P q 2 , r 2 ( i , j ) + K q 1 , r 2 P q 1 , r 2 ( i , j ) + K q 2 , r 1 P q 2 , r 1 ( i , j ) P q 1 , r 1 ( i , j ) + P q 2 , r 2 ( i , j ) + P q 1 , r 2 ( i , j ) + P q 2 , r 1 ( i , j )
其中Pq,r(i,j)由方程(2)给出。在下一步,在确定了第一个卷积函数后,在频率域中执行的下面子步中,重建第一个图象:
1.确定第一幅图象的富里叶变换F{K1(i,j)}和第一个卷积函数的富里叶变换F{G1(i,j)},接着
2.通过F{K2(i,j)}=F{K1(i,j)}/F{G1(i,j)}确定第二幅图象的富里叶变换F{K2(i,j)},接着
3.通过富里叶逆变换确定第二幅图象K2(i,j)。
如果需要,在下一步可以减少出现的噪声,即在图象F{K2(i,j)}的富里叶逆变换之前执行一附加的滤波步骤,例如用低通滤波器把F{K2(i,j)}乘以低通滤波器的频率响应H(i,j)。这方法的另一个实施例可以用在测量光源和检测器之间的距离也变化的其他测量装置中。这种方法的实施例参见图5予以说明。
图5所示的测量装置500包括用于测量光源503的第一个移动座509,用于检测器506的第二个移动座510,这二个移动座相互围成一个角度α例如90度。第一个移动座509和第二个移动座510能彼此相对移动,且装在平台501上。另外,在第一个移动座509和第二个移动座510上的测量光源503和检测器506彼此相对这样固定:测量光源503的输出窗口504和检测器506的输人窗口507基本上处于同一平面上。测量光源503和检测器506的位置通过马达部件512独立调节。在第一个移动座509和第二个移动座510通过连结到控制单元522的马达部件512,能对很多的测量位置进行设定。测量位置数处于例如16到256之间。测量空间508在第一个移动座509和第二个移动座510之间。在测量空间508中,可以放置例如不透明介质或生物介质。测量空间的直径例如100mm。为了实现强度测量,控制单元522对光源520产生一个控制信号,以使在光源520中产生光。光源的光通过第一光波导505和输出窗口504传到测量空间508的不透明介质中。检测器506接收一部分在测量空间508中通过不透明介质传送的光,且把接收的光通过第二光波导513传到光检测器521。光检测器521把接收的光转成检测器信号。检测器信号经电导体527传给模数转换器523。在模数转换器中检测器信号转数值1。数字强度值1接着通过数据总线528传到重建单元525。数据处理发生在重建单元525中。数据处理的结果以图象的形式出现在监示器526上,或存储在数据存储单元524中。强度Io在测量光源的M位置和检测器的N位置中再被测量且列在表Pm(i)中,及通过下列方程转换成衰减系数表Pm(Ki) Δk = 1 l 1 n ( I i I 0 ) 。根据本发明的方法,G1(i,j)用M测量光源位置和N检测器位置之间属于光路Ls,d的衰减系数,在正交二维空间中对每个元素(i,j),用求和公式(7)再计算。这些衰减系数直接确定第一个卷积函数G1(i,j),用这个函数第二幅图象被卷积。下一步,在频率域内执行下列子步实现重建第二幅图象:
1.确定所定的衰减系数的富里叶变换F{K1(i,j)},和第一个卷积函数的富里叶变换F{G1(i,j)},确定第一个卷积函数的富里叶变换F{G1(i,j)},接着
2.通过F{K2(i,j)}=F{K1(i,j)}/F{G1(i,j)}确定第二幅图象的富里叶变换F{K2(i,j)},接着
3.通过富里叶逆变换F{K2(i,j)}确定第二幅图象K2(i,j)。
4.通过K2(i,j)的扇形光束反投影,重建第二幅图象K2(i,j)。
在此方法的实施例中反投影在数据解卷积以后执行。如果需要,在下一步可以减少出现的噪声,即,在图象F{K2(i,j)}的富里叶逆变换被执行之前附加一滤波步骤,例如用低通滤波器把F{K2(i,j)}乘以低通滤波器的频率响应Hi(i,j)。
应该注意,不是测量装置仅用一个测量光源和一个检测器,测量装置可以改用包含一个测量光源和多个检测器的测量装置。这种方法的实施例参见图6的说明。图6所示的测量装置600包括等距离置于平台601的圆周上的测量光源603和一些检测器606。这些检测器606例如是PIN光电二极管。在测量空间608中不透明介质被置于平台601内圆周中心的台面602上。台面602能通过马达部件612以相对于光源603的很多角度位置旋转。角度位置数可以在16到256位置之间变化,例如32个位置。马达控制部件612通过控制单元622控制。另外,控制单元622通过选择信号611控制选择装置610选择N个检测器606中的一个,其中检测器N的数量在值16到64之间,例如32。作为该选择的结果,强度信号613,所表示的在选择的检测器606上的强度值提供给AD转换器623。为了实现强度测量控制单元622为光源620产生一个控制信号,以使在光源620中产生光。测量光源620的光被引到第一个光波导605和输出窗口604进入测量空间608中不透明介质内。检测器606接收在测量空间608中通过不透明介质传播的一部分光。检测器606把接收的光转成检测器信号。检测器信号被传到选择装置610。控制单元622根据产生的选择信号611用选择装置610交替地选择全部检测器信号中的一个。所选择的检测器信号经电导体627传给模数转换器623。此模数转换器623把检测器信号转成数字强度值I。数字强度值I接着通过数据总线628传到重建单元625。然后,测量的强度存储在重建单元625中的强度表Pm(i)中。在一个位置上中测量了检测器的所有强度后台面旋转到下一个角度位置,这里向前面的角位置旋转角度360/M,其中M是为完成不透明介质测量所测量的表的角度位置数。在测量步骤以后,M*N强度的表可以进一步在重建单元625中处理和重建。进一步的处理类似于在图2中说明的处理步骤。重建后的图象可以在监示器626上显示或存储在存储单元624中。同该实施例密切相关的另一个实施例是包括一个检测器606及很多测量光源503的一个类似的装置。在这个装置中控制单元选择很多测量光源中的一个测量光源。在本发明的下一个实施例中用了很多测量光源和很多检测器。根据本发明方法用的测量装置同图7相关地被解释。图7显示了另一个测量装置700,该测量装置包括置于测量空间708的围着不透明介质在平台701上园周等距离的众多的M测量光源703和众多的N检测器706。在该实施例中,光源703包括例如光发射二极管及检测器包括例如PIN光电二极管。为了实现测量,控制单元722用形成的光源选择信号712和检测器选择信号711由光源选择装置713选择M个测量光源703中的一个测量光源703及根据检测选择装置710选择N个检测器中的一个检测器。所选择的检测器信号经电导体727传给模数转换器723。此模数转换器723把检测器信号转成数字强度值I。数字强度值I接着通过数据总线728传到重建单元725。通过连续选择测量器706所选择的测量光源703的强度测量结果存储在重建单元725中的强度表Pm(i)中。从所有可能的测量光源完成测量得到强度值矩阵P(i,j)。强度矩阵P(i,j)是可进一步在重建单元725中处理。重建后的图象可在监示器726上显示或存储到存储单元724中。处理方法类似于图2所述的方法。
而且,所有其它的源和检测器的位置的几何布局都可采用,只要所有源和检测器对之间的距离已知。因此,测量空间108的形状能够适应人类胸腔的形状对人类胸腔组织成图。

Claims (12)

1.在不透明介质内部成象的方法,这种方法执行下列步骤:
a)用来自多个位置上光源的光连续照射不透明介质及沿着通过不透明介质的多个光路检测传播的光,
b)根据所测量的强度确定在不透明介质内部图象,
其特征在于为确定目标的图象执行下列子步:
c)从测量的强度中确定光参数,
d)通过反投影确定所定光学参数的内部的第一幅图象,
e)从光学参数确定第一个空间变化卷积函数,
f)通过把第1个空间变化卷积函数同第1幅图象相组合,确定不透明介质内部的第2幅图象。
2.根据权利要求1的方法,其特征在于其中的第二幅图象是用第一个空间变化卷积函数对第一幅图象解卷积而确定的。
3.根据权利要求1的方法,其特征在于其中通过富里叶变换的第一个图象和富里叶变换的第一个空间变化卷积函数的商的富里叶逆变换,确定第二幅图象。
4.根据权利要求1,2或3的方法,其特征在于下列各步被执行以便在迭代处理中得到第n+1幅改进图象,其中的处理执行下列各步:
a)从第n幅图象和光学参数中确定第n+1个空间变化卷积函数,及
b)通过把第n幅图象同第n个空间变化卷积函数的组合,确定目标的第n+1幅图象。
5.根据权利要求4的方法,其特征在于该方法包括一个非线性坐标变换的步骤。
6.根据权利要求1的方法,其特征在于其中的光学参数是衰减系数。
7.根据权利要求1的方法,其特征在于该方法包括减少高频噪声元素的步骤。
8.根据权利要求1的方法,其特征在于在多个光源位置上的一个光源和在多个检测器位置上的一个检测器之间形成一些光路。
9.根据权利要求1的方法,其特征在于从位于多个位置上的多个光源中选择的一个光源和在很多位置上的多个检测器中选择的一个检测器之间形成一些光路。
10.根据权利要求1的方法,其特征在于在一个位置上的一个光源和在很多位置上的很多检测器中选择的一个检测器之间形成一些光路。
11.根据权利要求1的方法,其特征在于从多个位置上的很多光源选择的一个光源和在一个位置上的一个检测器之间形成一些光路。
12.一个不透明介质内部的成象装置,该装置提供:
a)用于照射不透明介质的光源,
b)测量穿过不透明介质传播来的光强度的检测器,
c)从多个不同的位置照射不透明介质的装置,
d)在多个不同的位置上定位检测量的装置,
e)为装置产生控制信号的控制单元,以在很多位置上定位光源和检测器,及
f)根据测量的强度重建图象的重建单元,
其特征在于准备好执行下列步骤:
a)从测量的光强度中确定衰减系数,
b)根据反投影,从子步a)中所定的衰减系数确定目标的第一个图象,
c)从衰减系数确定第1个空间变化卷积函数,
d)把第1幅图象同第1个空间变化卷积函数相组合,确定第二幅图象。
CN96191327A 1995-09-11 1996-08-29 在不透明介质内部成象的方法和装置 Expired - Lifetime CN1096658C (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
EP95202456 1995-09-11
EP95202456.0 1995-09-11

Publications (2)

Publication Number Publication Date
CN1167537A CN1167537A (zh) 1997-12-10
CN1096658C true CN1096658C (zh) 2002-12-18

Family

ID=8220625

Family Applications (1)

Application Number Title Priority Date Filing Date
CN96191327A Expired - Lifetime CN1096658C (zh) 1995-09-11 1996-08-29 在不透明介质内部成象的方法和装置

Country Status (6)

Country Link
US (1) US5719398A (zh)
EP (1) EP0791207B1 (zh)
JP (1) JP4376317B2 (zh)
CN (1) CN1096658C (zh)
DE (1) DE69616743T2 (zh)
WO (1) WO1997010568A1 (zh)

Families Citing this family (23)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11514096A (ja) * 1996-08-14 1999-11-30 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁状媒体の画像を形成する装置及び方法
US6100520A (en) * 1996-11-29 2000-08-08 Imaging Diagnostic Systems, Inc. Detector array for use in a laser imaging apparatus
US6130958A (en) * 1996-11-29 2000-10-10 Imaging Diagnostic Systems, Inc. Method for reconstructing the image of an object scanned with a laser imaging apparatus
EP1006866A1 (en) 1996-11-29 2000-06-14 Imaging Diagnostic Systems, Inc. Detector array for use in a laser imaging apparatus
US6150649A (en) * 1996-11-29 2000-11-21 Imaging Diagnostic Systems, Inc. Detector array with variable gain amplifiers for use in a laser imaging apparatus
EP0897282A2 (en) * 1996-12-03 1999-02-24 Koninklijke Philips Electronics N.V. Method and apparatus for imaging an interior of a turbid medium
WO1998051209A1 (en) 1997-05-09 1998-11-19 Koninklijke Philips Electronics N.V. Device for localizing an object in a turbid medium
JP2001524664A (ja) 1997-11-26 2001-12-04 イメージング・ダイアグノスティック・システムズ、インコーポレイテッド 時間分解胸部撮像装置
WO1999066832A1 (en) * 1998-06-25 1999-12-29 Koninklijke Philips Electronics N.V. Method of localizing an object in a turbid medium
US6332093B1 (en) 1998-08-06 2001-12-18 Art Recherches Et Technologies Avancees Inc./Art Advanced Research Technologies, Inc. Scanning module for imaging through scattering media
EP1207385B1 (en) * 1999-06-03 2011-02-09 Hamamatsu Photonics K.K. Optical ct device and method of image reformation
WO2001074238A1 (en) * 2000-03-31 2001-10-11 Koninklijke Philips Electronics N.V. Method and device for localizing a deviant region in a turbid medium
US6529770B1 (en) 2000-11-17 2003-03-04 Valentin Grimblatov Method and apparatus for imaging cardiovascular surfaces through blood
US7806827B2 (en) * 2003-03-11 2010-10-05 General Electric Company Ultrasound breast screening device
JP4595330B2 (ja) * 2004-01-19 2010-12-08 ソニー株式会社 画像処理装置および方法、記録媒体、並びにプログラム
EP1951108A2 (en) * 2005-11-18 2008-08-06 Koninklijke Philips Electronics N.V. Device for imaging an interior of a turbid medium
JP2009529948A (ja) * 2006-03-17 2009-08-27 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 混濁媒体画像形成装置
US8406846B2 (en) * 2006-03-31 2013-03-26 Shimadzu Corporation Mammographic apparatus
JP4808584B2 (ja) * 2006-10-18 2011-11-02 浜松ホトニクス株式会社 乳房撮像装置
EP1943941A1 (en) * 2006-12-21 2008-07-16 Koninklijke Philips Electronics N.V. Method for optically imaging an interior of a turbid medium, method for reconstructing an image of an interior of a turbid medium, device for imaging an interior of a turbid medium, medical image acquisition device and computer program products for use in said methods and devices
EP2224845A2 (en) * 2007-12-17 2010-09-08 Koninklijke Philips Electronics N.V. Method for detecting the presence of inhomogeneities in an interior of a turbid medium and device for imaging the interior of turbid media
CN103057261A (zh) * 2011-10-24 2013-04-24 致伸科技股份有限公司 用于边缘侦测装置的控制方法与控制装置
TWI605845B (zh) * 2016-12-21 2017-11-21 財團法人工業技術研究院 靜脈輸液用藥監測裝置及其方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5070455A (en) * 1989-11-22 1991-12-03 Singer Imaging, Inc. Imaging system and method using scattered and diffused radiation
US5396285A (en) * 1993-05-07 1995-03-07 Acuson Corporation Ultrasound imaging method and apparatus with dynamic non-linear filtering

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5079698A (en) * 1989-05-03 1992-01-07 Advanced Light Imaging Technologies Ltd. Transillumination method apparatus for the diagnosis of breast tumors and other breast lesions by normalization of an electronic image of the breast

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5070455A (en) * 1989-11-22 1991-12-03 Singer Imaging, Inc. Imaging system and method using scattered and diffused radiation
US5396285A (en) * 1993-05-07 1995-03-07 Acuson Corporation Ultrasound imaging method and apparatus with dynamic non-linear filtering

Also Published As

Publication number Publication date
WO1997010568A1 (en) 1997-03-20
EP0791207A1 (en) 1997-08-27
JP4376317B2 (ja) 2009-12-02
DE69616743T2 (de) 2002-07-25
EP0791207B1 (en) 2001-11-07
JPH10509243A (ja) 1998-09-08
CN1167537A (zh) 1997-12-10
US5719398A (en) 1998-02-17
DE69616743D1 (de) 2001-12-13

Similar Documents

Publication Publication Date Title
CN1096658C (zh) 在不透明介质内部成象的方法和装置
O'Leary Imaging with diffuse photon density waves
Mueller et al. Reconstructive tomography and applications to ultrasonics
Pogue et al. Initial assessment of a simple system for frequency domain diffuse optical tomography
CN108577876B (zh) 一种多边形静止ct及其工作方法
KR101252010B1 (ko) 라돈데이터로부터 (n+1)차원 영상 함수를 재구성하는방법과 장치
JP6139092B2 (ja) X線ct装置およびシステム
CN1200174A (zh) 基于荧光寿命的人体组织及其它无规则介质成象技术和光谱技术
JP6053770B2 (ja) 一体型マイクロトモグラフィおよび光学撮像システム
CN104013387B (zh) 一种太赫兹快速断层成像系统及方法
Kumavor et al. Target detection and quantification using a hybrid hand-held diffuse optical tomography and photoacoustic tomography system
Na et al. Transcranial photoacoustic computed tomography based on a layered back-projection method
CN105894562A (zh) 一种光学投影断层成像中的吸收和散射系数重建方法
CN102034266B (zh) 激发荧光断层成像的快速稀疏重建方法和设备
CN111915733B (zh) 基于LeNet网络的三维锥束X射线发光断层成像方法
Guggenheim et al. Multi-modal molecular diffuse optical tomography system for small animal imaging
Cao et al. Bayesian reconstruction strategy of fluorescence-mediated tomography using an integrated SPECT-CT-OT system
Feng et al. Bioluminescence tomography imaging in vivo: recent advances
Wang et al. A novel finite-element-based algorithm for fluorescence molecular tomography of heterogeneous media
JP2001346894A (ja) 線量分布測定装置
Ramm et al. Optical CT scanner for in-air readout of gels for external radiation beam 3D dosimetry
Mudeng et al. Computational image reconstruction for multi-frequency diffuse optical tomography
Zhang et al. Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy
Cai et al. Three-dimensional radiative transfer tomography for turbid media
Lo et al. Small-animal 360-deg fluorescence diffuse optical tomography using structural prior information from ultrasound imaging

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CX01 Expiry of patent term

Granted publication date: 20021218

EXPY Termination of patent right or utility model