专利名称：图像重构装置、图像重构方法、图像重构程序、ct装置的制作方法
技术领域：
本发明涉及从对象物的放射线投影像，重构截面像的技术。
背景技术：
计算机层析成像(CT)是一种使用X线等光线，对物体从各个 角度投影的投影像进行摄影，其后通过计算机获得截面像的技术。图 1表示典型的X线CT装置的模式图。
在典型的CT装置中，通过X线光源的旋转移动，取得测定对象 的从各个角度投影的投影像。从由此得到的投影像，通过计算操作， 获得截面像的技术即为CT。通常，根据Filtered Back-projection法 (FBP)，通过变换操作从投影像变换为截面像而获得。在FBP中， 在投影像上实施过滤(基本上为微分过滤)后，在原来的投影方向实 施"反投影"等操作，由此获得截面像。此时，微分过滤具有放大噪 音、误差，容易产生伪像(原来不存在的误差、虛像)等特点。并且， 反投影操作具有使其伪像向截面像整体传播的效果。因此，CT中的 伪像并不仅仅止于有问题的部分，会破坏截面像的整体，这点构成致
伪像的多半是由包含在FBP中的过滤操作、反投影操作引起， 因此，若不使用FBP，就应该能够得到没有伪像的截面像。Algebraic Reconstruction Technique ( ART )作为除FBP以外的计算截面像的方 法，在历史上非常重要。ART在提出FBP之前就被利用，是一种有 历史的方法。在ART中将截面像的计算作为拟合(fitting)问题加以 对待，将截面像作为参数，将投影像作为拟合对象，逐次修正截面像， 使从截面像计算出的投影像(p)与通过实验求得的投影像(pQ) —
12致。在ART中，对截面像渐进地进行改进，使(p-po)趋于0，这点 即为特征。ART完全不使用FBP之外，但在求解截面像的计算上耗 费时间，因此如今只用在特殊的用途(地震波的解析等)。ART不会 像FBP那样出现极端的伪像，但一般来说FBP所得出的截面像更自 然一些。
另外，除了过滤操作、反操作之外，作为引起伪像的原因，还可 以列举投影像中的数据的脱落、过少。若数据少，则作为结果而得出 的截面像的画质也跟着降低，这是显而易见的。人们发现在FBP中， 数据的脱落、减少同样会产生致命的伪像，这被认为是重大的问题。 基于拟合的CT (包括ART)，关于数据的脱落、过少，比FBP强。 但是，CT中的数据的脱落，被认为是极端地"条件差"的问题，就 算使用ART等，也很难得到改善。在ART中将(p-p。)作为误差加 以考虑，但条件差时拟合容易失败。所谓最小化二乘误差(p-pj的方 法，作为拟合而言稳定一些。作为最小化二乘误差的方法，最小二乘 法最为标准。在最小二乘法中，求解一边具有相当参数个数的方阵的 逆矩阵。CT中的参数是构成截面像的各个像素值，因此，参数的个 数十分庞大。比如1000 x 1000像素的截面像，参数的个数达100万 个，矩阵的因子数达l兆个，非常庞大。因此，在正攻法的最小二乘 法中矩阵变得过于庞大，从而导致失败。因此，避开正攻法的最小二
法,提出了 Simultaneous iterative reconstruction technique (SIRT)、
Iterative least squaretechnique (ILST )。 与ART相同，将截面像的计
算作为拟合问题加以对待，但如上所述，由于不得不避开正攻法的最 小二乘法，因此SIRT与ILST在计算中也同样作为投影操作的逆演 算利用FBP。因此，不能从根本上解决由FBP引起的过滤操作和有 关反投影的问题。或许，出于这种原因，只报导过说，在SIPT、 ILST 中仅具有"降低"拟合的效果。 [非专利文献1]Yazdi M, Gingras L, BeaulieuL. An adaptive approach to metal artifact reduction in helical computedtomography for radiation therapy treatment planning: experimental and clinicalstudies. Int J Radiat Oncol Biol Phys，62: 1224-1231,2005。
发明内容
从以上论述可知，为了得到没有伪像的截面像，不使用FBP，抗 数据脱落强等性质非常重要。因此，像ART那样，不使用FBP,而作 为拟合的评价函数利用二乘误差是谁都能想到的。尽管如此，CT中 的伪像问题不能得到根本上的解决，这是因为这种单纯的想法不能够 简单地实现。第一问题在于，截面像的计算作为拟合问题，规模过于 庞大，计算耗费时间。有必要适当进行高速化。第二问题在于，最小 二乘法作为拟合算法，不太抗"[条件差"问题。因此，可以认为， 根据最小二乘法的已知的算法(ILST等)的延伸并不能解决伪像问题。 第三问题在于，就像SIRT、 ILST—样，已知的算法不能完全排除FBP。
本发明的第一新颖性在于，作为拟合方法釆用了 Simulated Annealing法(SA )。 SA为公知的技术，拟合耗费时间，本质上与CT 不符。但是，即使通过硬将SA适用于CT，从而完全不使用FBP,也 实现了基于二乘误差最小化的截面像的计算。并且，SA对于数据脱 落等条件差的拟合，具有较为稳定的性质，即使在这点上，也有利于 伪像的根本上的解决。鉴于以上的情况，通过在CT上适用SA，从而 发现了在根本性地解决伪像上的效果，这点即为本发明重要的新颖 性。单纯考虑在CT上适用SA时，需无数次的反复进行从截面像求解 投影像的计算。真正这样计算，从截面像求解投影像的计算，与FBP 大致相同，当在CT上适用SA时，则相比FBP需耗费数百万倍的时间。 就算用现在的高速计算机，也需要以年为单位的计算时间。因此，在 本发明中，通过式子的变形，大幅降低了计算量，解决了该问题。该 办法在将SA适用于CT时为必须的技术，是本发明的重要因素。第二新颖性在于，作为拟合的评价函数，除了二乘误差，还导入
了积极破坏伪像的平滑(smothing)项和熵(entropy)项。在以往的 CT中，仅仅将p与po的差或二乘误差作为对象。此时，将伪像模糊化， 则可能得到较好的拟合结果(误差小)，这种可能性时常存在，而现 实中的大多情况为这种情况。即，伪像成为被其他伪像抵消的状态。 因此，很难消除伪像。在本发明中，仅仅将二乘误差作为评价函数的 情况下，同样也没有完全消除伪像。这种事实意味着，为了获取没有 伪像的截面像，除了二乘误差以外，还需要其他条件。在本发明中， 从统计力学出发，提出了平滑项、熵项。这些项是将截面像必须为光 滑自然的浓淡图像等理所当然的条件用数式形式表现出来的。熵项对 拟合进行诱导，使之破坏伪像，从而使截面像整体画质均勻化。平滑 项具有对用熵项诱导出的截面像的粗糙度进行控制的效果。通过导入 这些项，能够消除伪像，并且能够得到自然的截面像。熵项和平滑项 虽然分别单独具有降低伪像的效果，但将两个进行组合的话更有效。
并且，在本发明中，作为比一般定义更为广泛的概念，将包含X 线、可视光线、电波的电磁波，由电子、荷电粒子构成的粒子线，作 为介质振动的音波等所有的情况，概括称为"放射线"。
作为本发明的第一效果，可以列举起因于数据脱落的伪像的减少 带来的效果。作为数据脱落构成问题的情况，可以列举在观察对象上 存在不透明部位的情况、投影角度存在限制的情况、投影角度的刻紋 不均匀的情况、锥束(conebeam )、螺旋状扫描等三维CT的情况。
尤其能够看出效果的是，在观察对象上存在不透明部位的情况下 所表现出来的金属伪像的减少效果。金属伪像是指，在观察对象中， 相对X线存在不透明部位(多数情况下为金属部分)时，用CT得到的 截面像整体(不仅仅是不透明部位)遭到破坏性错乱的情况。金属伪 像产生的原因在于，在不透明部位投影像的亮度的不连续变化、和不 透明部位的信息脱落。当FBP中适用微分过滤器时，亮度的不连续变
15化成为异常值，由反投影操作，生成以不透明部位为中心的筋状的伪像。并且，由于信息的脱落，在与不透明部位没有直接关系的部位上，产生预料之外的明暗。本发明不使用FBP而利用对于信息的脱落比较稳定的SA,对金属伪像的去除有效，这是不难想象的。
事实上，重要的是在于，向锥束、螺旋状扫描的应用。两者称为
3维CT,是一种目前急速普及的技术。但是公知的是3维CT上会出现特有的伪像，虽然知道其产生的原因，但还没有解决办法。其原因，在锥東的情况下，为数据的脱落。在锥束中无法满足获取完整截面像的条件，从而出现起因于数据脱落的伪像。螺旋状扫描中的伪像产生的根本原因在于反投影操作上。在螺旋状扫描中系的几何学异向性(螺旋)影响过滤操作、反投影操作，出现风车伪像。本发明无需进行过滤操作、反投影操作，抗数据脱落强，因此能够解决3维CT中的诸多问题。
作为投影角度不均匀的情况可以列举的是，根据地震波的地球内部的CT的解析、根据使用来自人工卫星的电波的大气状态的CT的解析。这些被认为是无法利用FBP的典型。通过利用本发明，能够期望提高解析精度。
本发明的第二效果在于，投影像摄影的迅速化、X线曝光量的减少上。SA具有不仅对于数据的脱落比较稳定，对于过少也比较稳定的特点，本发明也一样。在CT的情况中，相当于数据过少的情况是投影角度的数目少。即，利用本发明，能够比以往的CT减少投影角度的数目。投影角度的数目是指根据放射线的投影像的摄影张数。摄影张数与摄影时间和曝光率成比例，因此，若能减少摄影张数，就能缩短摄影时间，降低曝光率。
并且，作为数据过少的其他形态，可以列举投影像的画质不好(S/N比低)的情况。可以看出，本发明在这种情况下也较稳定。若能允许投影像画质的下降，则这也能缩短摄影时间、降低曝光量。或者，在S/N比极其低的SPECT、 PET中，本发明能够提高画质。
本发明的第三效果在于，截面像亮度值的决定精度高。本发明中能够得到相当如实地再现了所测定的投影像的截面像。再现精度比FBP高两个数量级左右。这是由所谓二乘误差最小化的拟合算法带来的。担高亮度值的决定精度，保证亮度值的定量性，从而能够进行使用亮度值的密度测定。这能提高骨密度测定等的测定精度。并且，也能提高异常部位(有肿瘤的脏器等)的检测精度。
本发明的第四效果在于，根据本发明所得到的截面像的对比度高于FBP。如第三效果所述，虽然本发明的亮度值的决定精度高，但作为其次要效果，截面像的对比度变高。若对比度高的话，表面上看来，空间分解能也变高。其结果，由本发明所得到的截面像其画质比以往高。要强调一点，即本发明不仅在特殊的条件、用途的CT中有效，而且在普通的CT上也有效。
图1为X线CT装置的模式图2为表示截面像f(xj)、投影像p(/^)的关系的模式图；图3为典型p。(r")的例子(横轴取检测器的信道位置、纵轴取投影角度的图像)；
图4为表示本发明基本步骤的流程图5 (a)为步骤(I )至(VI)的拟合过程的模式图5 (b)为步骤(1)至(9)的拟合过程的模式图6为表示根据"非专利文献l"刊载的方法所得出的结果的图7为从图6(a)通过模拟求解的sinogram(白色部分为不可视领
域)；
图8 (a)为将图6 (c)进一步放大的图；图8 (b)为表示本发明所得出的结果的图；图9为表示假想能量E的各项的效果的图；图10为表示具有角度限制的情况下的伪像和本发明的伪像减少效果的图。
具体实施例方式
(第l实施方式)
图1表示根据本发明实施方式的X线CT装置的概略结构。该X线CT装置具有摄影部和计算机。摄影部包括用于将X线照射在测定对象的光源，和检测穿透测定对象的X线的检测器，将测定对象的周围用X线从多个方向进行投影，从而取得投影像。计算机包括控制X线CT装置整体的控制部，和根据由摄影部得到的X线投影像，生成测定对象的所关心部位的截面像的图像重构处理部。另外，图l所示的结构在以下的各实施方式中是共通的，而在各实施方式中图像重构处理部所进行的处理分别不同。
本实施方式的图像重构处理部釆用Simulated Annealing法(SA )，作为用于从投影像求解截面像的拟合方法。因此，首先对SimulatedAnnealing法(SA )的架构进行说明。SA为从Monte Carlo法派生的拟合算法，其特征在于，根据随机数推进拟合，并模仿热力学处理假想能量、假想温度等。SA本身为公知的技术，按照以下步骤(i) (vi)进行。
(i) 根据随机数选择一个参数，并且根据随机数改变其参数(随
(ii) 进行变更后的评价。作为评价函数考虑假想能量E。在通常的SA中，E为二乘误差的总和。变更前后E的变化记为AE (评价)。
(iii) 若评价得到改善(A￡<0),则接受该变更(变更)。
(iv) 以exp(-AE/r)的概率接受改坏。
(v) —点点减少T。Ui)从(i)反复。
在SA中，如(iv)，按照波尔茨曼统计接受改坏，由此确保能够
18脱离局部解的可能性。因此，不会拘于局部解，能够达到大域解，对于条件不好的伪像也能稳定操作。并且，通过逐步减少T，软着陆于大域解。自(i)到(V)执行一次称为一个蒙特卡洛步骤。在SA中，无限次反复进行蒙特卡洛步骤，以此进行拟合。计算所需的时间为一个蒙特卡洛步骤所需时间x所需蒙特卡洛步骤数。所需的蒙特卡洛步骤数与参数的个数、自由度成比例。
接着，按照权利要求的形式，对本实施方式进行说明。在将CT考虑成拟合问题的情况下，拟合参数即为截面像(f(x，")。构成拟合
对象的数据即为所测定的投影像(P。(r") )。 r为摄影投影像的l维检测器的信道位置，大为投影角度。坐标的定义表示在图3。 p。(r")为通过改变角度0测定投影像而得出的数据组，当横轴取r、纵轴取^时，可以看作2维图像。这种数据称为sinogram。典型的sinogram在图3中表示。
极论而言，CT为从sinogram求解截面像的技术。在本实施方式中，考虑了临时截面像f (x,与所测定的投影像p。 (r, e)的二乘误差，但为了计算二乘误差，使用数101从f(xj)进行计算，求解投影像P。(r")。
<formula>formula see original document page 19</formula>
在图2中可知，数101为在s方向求的和、对f(x,少)的x方向的投影像进行的计算。使用由此得出的p(r，"， 二乘误差H能够按照以下进
行计算。[数102]
在通常的拟合中，作为假想能量E,直接利用数102，但在本实施方式中，其特征在于除H以外导入平滑项和熵项。E的定义如下。
1[数103]
T为表示假想温度(温度参数)的系数，S为表示熵项的系数，CT
为表示像素值标准偏差的系数，c为表示平滑项强度的系数。TS即为
熵项，CCT即为平滑项。S、C7的定义、计算法将在后续中论述。在本
实施方式中，使用这些算式，按照如图4所示的步骤计算截面像。 步骤(a):对以某种形式求解的暂定截面像f(x,"，求解评价函
数(以下称为"能量")(Eo) (ST10及ST20)。
步骤(b):用随机数选择(x。,少。)及A//,改变截面像f(;c,y)的一部
分(ST30 )。
步骤(c):对改变后的截面像f(x，y),求解评价函数E, (ST40及 ST50 )。
步骤(d):求解上述能量(E。)与上述能量(E,)的差(AE )( ST60)。 步骤(e):根据使用上述差(AE )及温度参数(T)的受理函数 (典型的为波尔茨曼统计)，判断是否釆用上述改变(ST70 ST110)。 步骤(f):返回上述步骤(a) (ST120)。
步骤(g):上述步骤(a) ~ (f)的重复每达到规定次数，就变 更一次假想温度(T)的值(ST130)。
步骤(h):判断上述步骤(e)的判断结果是否满足规定的停止 条件，若满足则结束处理。其中，执行了ST80、 ST100的情况设为"成 功"，执行了ST110的情况设为"失败"，当成功概率低于10% (该值 可适当调整)时结束。然后将结東时所釆用的推定截面像作为测定对 象的截面像，将其表示在计算机的显示器或存储到存储介质。
从步骤(a) 步骤(h)的一系列操作对应于权利要求12，在本 实施方式中，这些 一系列操作由图1的图像重构处理部来执行。
若考虑与SA的基本步骤(i) (vi)的对应，步骤(b)对应(i)， 步骤(a) (c) (d)对应(ii)，步骤(e)包含(iii)和(iv)，步骤(g)对应(V),步骤(f)对应(Vi)。从这点可以看出，本实施方式在CT
上如实地适用了SA。
并且，对于上述步骤(h)中的处理结東的判断，还可以不是由 图l的图像重构处理部进行，而是将推定截面像逐次表示在计算机的 显示器上，用户看见后指示计算机结東处理。
(第2实施方式)
为了计算数103所示的E，会经过数101和数102,执行关于s、 r、 0的级数和。由于计算三重积分(级数和)，会耗费大量时间。即， 意味着一个蒙特卡洛步骤所需时间很长。并且，在CT中参数的个数 非常庞大，其结果，实施SA所需的总计算时间即使利用现在的计算 机，也需要"年单位"。
因此，在本实施方式中，不对E进行计算，而仅仅是主要计算进 行变更时的变化量AE。 [数104]
现将有关暂定截面像的变更设为Af(x,y) 。 Af (x,少)为仅仅在坐标 (x。,少。)上取A//的值，除此之外为取O的截面像。该Af(x,少)的投影像 Ap(r")可以用与数101相同的方法计算。使用该Ap(r"), AH可以如下
计算。 [数105]
<formula>formula see original document page 21</formula>
将其写下去成为，
<formula>formula see original document page 21</formula>
现在，Af(x，y)只在(x。,少。)取值，因此Ap(r")只在一卜x。coW +少。sine上取值A^，除此之外为0。因此，数106的U中的 内容也仅仅在r的二x。cose + y。sin^上取值。因此，级数和无需对r和^的 双方进行，可如下表示。 [数107]<formula>formula see original document page 22</formula>
在数107中，仅仅构成关于e的级数和，这点非常重要。由此， 能够大幅减少计算量。p、 p。为数字图像，因此计算数107时，通常 有必要进行关于r(W的内插。因此，不能将数107的U的内容进一步 展开。但是，当允许误差而展开数107时可得出下式。 [数108]<formula>formula see original document page 22</formula>
但是，M为投影角度的个数。当使用数108来替代数107时，AH的 值会产生1%左右的误差。但是，数108的速度高于数107,在本发明 中的利用价值高。
以下对Acr的计算进行说明。a为坐标(x。,jO附近的亮度值的标准 偏差。作为(x。，凡)附近，考虑(x。,y。)周围dxd像素(本实施例中使用 d=5)。这些像素的标准偏差可通过下式求解。 [数109]<formula>formula see original document page 22</formula>
但是，
<formula>formula see original document page 22</formula>[数lll]
1 c /2
/( Wo》=7 I]/(x。 + z'， 7。 +力
-<i/2
若使用数i09 in,则A(j可如下计;
2 , 2
/ 2\ , /.
Aa= 〈/(x0，_y0)2〉-~'」2 ;」〈/(xo,}o)〉-~
但是，f,， f,为变更前及后的/(x。,y。)的值。
以下是对于AS的计算，在此之前，对平滑项S进行定义。 一般计 算机所处理的图像为数字图像，不仅是像素的坐标(x， y)，就连像素 的值也是数字值。因此， 一个像素看作一个量子，将像素值视为量子 状态。那么，图像可以看成由多个量子构成的系。按照统计力学，对 这种系的熵进行如下定义。
TV!
S二yHn
N为图像中像素的总数，Ni是像素值以数字值为i的像素的总数。 K在通常的物理中为波尔茨曼常数，但本发明与物理没有直接关系，
因此任取。在本发明中，与"一样，S也被定义为(x。,y。)周围的dxd
对于被定义成数113的S考虑AS。通过变更，像素值从数字值i变 更为j。那么AS变成如下式子。 [数114]
… W！
AV^2!…(M — +1)-"A^! AVAy…AV…^'!…气！
即，通过变更，N,减l，而Nj加l。若展开数114，则可以得到非常简单的式子。
M = Hn TV, - Hn(7V; +1)
综上，本实施方式的顺序如下
(I )用随机数选择(x。,y。)及A^。
(II )使用数104、 108(或107)、 112、 115计算AE。
(III )若评价得到改善(AE<0 )，则在/(x。,y。)上加Ap。
(IV)AE〉0时，也同样以概率exp(-AE/T)，在/(x。，y。)上加A//。
(V ) —点点减小T。
(VI) 从(I )重复。
(VII) 判断是否满足规定的停止条件，若满足则结束处理。
在这里，将在(m)、 (iv)中附加了A^的情况设为"成功"，除 此之外的情况设为[失败]，当成功概率低于规定值(例如10%，该值 可适当调整)时结束。然后，在结束时所采用的推定截面像作为测定 对象的截面像，将其显示在计算机的显示器上，或存储到存储介质。
可以看出，步骤(I ) (VI)完全遵循SA的基本步骤(i) (vi)。 然而，步骤(I ) (VII)的处理由图l中的图像重构处理部来执行。 另外，对于上述步骤(vn)的处理结束的判断，还可以不是由图
l中的图像重构处理部来执行，而是逐次表示在计算机的显示器上， 用户肉眼观察后指示计算机结束处理。
(第3实施方式)
继而，对上述方法的重要变形进行论述。在上述的方法中，利用 数1O7或数1O8求解了有关0的级数和。因此，研究一个像素的变更需 进行M次的反复演算。本实施方式进一步减少了该部分的计算量。
首先，考虑P(W)的反投影像g(xj)。 [数116]<formula>formula see original document page 25</formula>在数116中，没有进行过滤操作，因此g(x,y)虽与截面像/(x,3;)相 似，但会成为模糊的图像。使用g(x,"时，数108变成如下。<formula>formula see original document page 25</formula>但是，g。(x,y)为p。(r,e)的反投影像，可以与数116相同的方法进 行计算。数117在没有反复进行运算的这一点上优于数108。 g。(u)在 计算过程中完全没有改变，因此可以事先计算。另外，g(^)在每更 改/(x,少)的l点时，就会产生微小的变化，因此严格来说，每经过一
次(i ) (vn)步骤，有必要进行重新计算。但是，当假定/(w)的
改变带来的g(x,y)的变化微小时，可以适用其他步骤，这就是本实施 方式。以下表示步骤。然而，以下的步骤(1 ) (9)的处理由图1 中的图像重构处理部来执行。
(1 )从p。(r刀)求解g。(x,y)。
(2) 从f(x,少)求解p(/^),并且从p(r,e)求解g(x，少)。
(3) 用随机数生成将f^^的变更值作为像素值的图像A/4,少;)。
(4) 对各图像值适用数117，计算图像AH(xj)。
(5) 同样对各像素值适用数112、数115，计算图像A"x,少)， △d
(6) 根据数104计算AE(xj)。
(7) 对于AE为正的坐标(x,少)，将A/,(x，;;)设为0。
(8 )向f(x，y)加上A/'(x,y)。
(9) 将丁放大《倍("<1 ),从(2)重复。
(10) 判断是否满足规定的停止条件，若满足则结束处理。其中， 在(7)中的A/z(x,y)成为0的情况设为"失败"，不是0的情况设为"成 功"，在成功概率低于规定值(例如10%，该值可适当调整)时结束
25处理。然后，将在结束时所釆用的推定截面像作为测定对象的截面像， 将其显示在计算机的显示器上，或存储到存储介质。
然而，对于上述步骤(10)的处理结束的判断，可以不是由图l 中的图像重构处理部来执行，而是将推定截面像逐次显示在计算机的 显示器上，用户肉眼观察后指示计算机结東处理。
图5为表示步骤(I ) (VII)和步骤(1) ~ (10)的模式图。
在图5中将参数值(本件中为像素值)作为轴，以浓淡表示假想能E (E越小越浓)。按照图示的情况，参数的个数(像素的个数)取两 种情况。在图5 (a)所示的步骤(I ) (VII)中，弯弯曲曲地曲折 前进，时而完全偏离，最终探索E的最小值。另外，在图5(b)所 示的步骤(1) (10)中，考虑E的倾斜的同时探索E的最小值。 可以直观地理解图5 (b)的效率高。数117带来的计算量的降低及 图5 (b)所示的探索效率的提高，有利于本实施方式的计算速度。 步骤(1) (10)虽然没有完全遵循SA的基本算法(i) (vi),但 (I ) (VII)步骤为自然的发展形式。 (其他实施方式)
执行上述各实施方式中所说明的处理的图像重构处理部，可以通 过以下方式实现用于执行这些处理的计算机程序，安装有该程序的 计算机，及执行这些处理的专用LSI等。 (实施例)
作为实施例，给出模拟的金属伪像的去除效果。首先，为了比较， 在图6 (非专利文献的fig.5)示出刊载在"非专利文献l"的方法所 带来的结果。
图6中，左栏(a) (d) (g)为没有Metal Artifact的原版模型图 像。中央(b) (e) (h)为在原版模型的规定位置设定不可视领域， 通过通常的FBP法重构的图像。右栏(c) (f) (i)为与上述相同， 设定不可视领域，通过"非专利文献l"的方法重构的图像，在以往技术中去除Metal Artifact最好的 一个例子。
因此，从"非专利文献l"中取图6(a)的图像，在与文献相同 的位置设定不可视领域，模拟计算投影像。结果表示在图7。图7对 应本实施例的P。(。W。对于图7,将根据本发明重构的结构及用于比 较的图6(c)的扩大图表示在图8。效果一目了然。在本实施方式(b) 中即使想找出Metal Artifact的特征都困难。
但是在本实施方式(b)中边缘附近有些模糊。这是因为"使截 面像光滑的因子"(熵项c(T )有些过强。可以看出，在Metal Artifact 的减少上，需稍强地设置该因子。由此按照本实施方式，在计算假想 能E时的各因子的平衡上，还有改进的余地，但Metal Artifact去除 效果非常高。
另外，为参考起见，将没有对平滑项w、熵项TS进行设定而执 行本实施方式算法的情况的结果表示在图9中。图9 (a)表示没有 平滑项^ (仅有熵项TS)时的结果，图9(b)表示没有熵项TS(仅 有平滑项ca)时的结果，图9(c)表示即没有熵项TS也没有平滑项 ^时的结果，当比较图9和图8 (b)时，可以看出，在本实施方式 中，作为能量E的因子设定平滑项ccr和熵项TS的两方的情况，比仅 仅设定某一个，和两方都不设定的情况，Metal Artifact去除效果更好。 只有平滑项时，具有在伪像中以不透明为中心残留有呈放射状扩展的 条紋的特点。另外，只有熵项时，具有截面像粗糙，S/N低的特点。 并且，最终的画质，慢冷过程中的平滑项的系数c和假想温度T之比 非常重要。
接着，将具有旋转角度限制的情况的实施例表示在图10。图10 U)表示仅根据FBP (通常的CT)的重构结果。图10 (b)表示适 用ILST (根据最小二乘法的图像复原法)的结果。图10 (c)表示根 据本实施方式的重构结果。将这些进行比较时，可以发现，相对限制 角度时所出现的伪像，本实施方式更为有效。根据角度限制的伪像具有以下特点，即圆状领域成为两端尖的杏仁领域，杏仁状领域的不尖 锐的躯干周围的亮度进行反转。若按照本实施方式，则能够降低这方 面的特征。
产业上的可利用性
本发明相对Metal Artifact具有显著的效果，因此在Metal Artifact 严重的领域，例如牙齿的CT、包含金属制植入管的CT中尤其有用。
并且， 一般地，本发明对具有信息脱落的投影像的装置有效。作 为信息脱落显著的情况中的一个，可以列举具有投影角度限制的情 况。投影角度限制构成问题的领域为三维电子显微镜法、乳腺摄影的 CT版、平行移动CT法(日本特开2006-71472号公报)等。作为信 息脱落构成问题的其他情况，可以列举锥束CT、螺旋状扫描CT等3 维CT。本发明对于去除3维CT上所出现的伪像也有效。
并且，也可利用在信息极少的系。例如，可用在荧光X线CT、 根据地震波的地球内部的CT等上。
并且，作为本发明的次要的效果，可以列举能够将截面像的亮度 值(在X线中与线吸收系数对应)以比以前高的精度进行决定。该 效果可以应用在骨密度测定等的精度改善上。
若使用本发明，则能够得到对比度比以前的方法更高的重构像。 因此，在伪像等不成问题的通常的X线CT中，本发明的利用价值也 高。并且，本发明对于数据过少，也处于稳定，因此在缩短投影像测 定的时间，进而降低X线曝光量上有效。综上，可以看出，本发明 隐藏了能够替换所有已有的X线CT技术的可能性。
权利要求
1、一种图像重构装置，从向对象物照射放射线而得到的投影像(以下称为“放射线投影像”)求解所述对象物的截面像，其特征在于，包括单元(a)，求解评价函数(以下称为“能量”)(E0)，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述放射线投影像的差；单元(b)，改变所述目前的推定截面像的一部分；单元(c)，求解能量(E1)，该能量(E1)包含从所述改变后的推定截面像，通过演算而得到的投影像和所述放射线投影像的差；单元(d)，求解所述能量(E0)与所述能量(E1)的差(ΔE)；单元(e)，根据使用所述差(ΔE)与控制受理概率的温度参数(T)的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；单元(f)，所述由单元(a)至(e)的一系列的反复处理，每达到规定次数，就改变一次所述温度参数(T)的值。
2、如权利要求1所述的图像重构装置，其特征在于，所述单元(a)求解能量(E0)，该能量(E0)包含从所述目前的推定截面像，通过演算得到的投影像与所述放射线与所述改变后的推定截面像的局部领域的标准偏差的和。
3、如权利要求1所述的图像重构装置，其特征在于，所述单元(a)求解能量(E0),该能量(E0)包从所述目前的推定截面像，通过演算得到的投影像与所述放射线投影像的差，与所述目前的推定截面像的局部领域的熵项的和；所述单元(c)求解能量(E。，该能量(Ei)包含从所述改变后的推定截面像，通过演算得到的投影像与所述放射线投影像的差，与所述改变后的推定截面像的局部领域的熵项的和。
4、 如权利要求1所述的图像重构装置，其特征在于，所述单元(a)求解能量(Eo)，该能量(Eq)包含从所述目前的推定截面像，通过演算得到的投影像与所述放射线与所述目前的推定截面像的局部领域的标准偏差，与所述目前的推定截面像的局部领域的熵项的和；所述单元(c)求解能量(E,)，该能量(E!)包含从所述改变后的推定截面像，通过演算得到的投影像与所述放射线投影像的差，与所述改变后的推定截面像的局部领域的标准偏差，与所述改变后的推定截面像的局部领域的熵项的和。
5、 如权利要求1所述的图像重构装置，其特征在于，用单元(h)替代所述单元(a)、 (c)、 (d),所述单元(h)根据[数l]算出AH，并求解AE (AE-AH +…)，该△E (AE = AH + ...)作为构成要素包含通过[数l]算出的AH,i:数i]其中，当将所述对象物的目前的推定截面像设为f(:c,y)、将根据所述单元(b)的改变部分设为Af(xj)时，Af(x，y)仅在坐标(x。,少。)取A//的值，除此之外为0的截面像，p(r,e)为从所述对象物的目前的推定截面像，通过演算得到的投影像，p。(r")为所述对象物的放射线投影像，r为摄影投影像的1维检测器的信道位置，^为投影角度，<formula>formula see original document page 4</formula>
6、 如权利要求l所述的图像重构装置，其特征在于，用单元(h)替代所述单元(a)、 (c)、 (d),所述单元(h)根据[数2]算出AH，并求解<formula>formula see original document page 4</formula>该<formula>formula see original document page 4</formula>作为构成要素包含通过[数2]算出的AH，<formula>formula see original document page 4</formula>其中，当将所述对象物的目前的推定截面像设为f(x,y)、将根据所述单元(b )的改变部分设为Af(u)时，Af(x,少)仅在坐标(x。,少。)取A^的值，除此之外为0的截面像，p(r，S)为从所述对象物的目前的推定截面像通过演算得到的投影像，p。(r")为所述对象物的放射线投影像，r为摄影投影像的1维检测器的信道位置，^为投影角度，<formula>formula see original document page 4</formula> M为投影角度的数目。
7、 如权利要求5或6所述的图像重构装置，其特征在于，所述单元(h)根据[数3]算出Ac7,并求解<formula>formula see original document page 4</formula>该<formula>formula see original document page 4</formula>作为构成要素包含cA(j与所述AH ，所述cAcx为根据[数3]算出的Aa与系数C的相乘，<formula>formula see original document page 4</formula>并，a为坐标Oc。,y。)的周围dxd像素的亮度值的标准偏差，由[数4J表示，f,、 f,分别为根据所述单元(b)的变更前后的f(x。,y。)的值,[数4]<formula>formula see original document page 5</formula>
8、如权利要求5或6所述的图像重构装置，其特征在于，所述单元(h)根据[数7]算出AS，并求解AE ( AE = AH-TAS +…)，所述AE(AE = AH-TAS +…)作为构成要素包含-TAS与AH的和，所述-TAS为通过[数7]算出的AS与温度参数(T)的相乘，[数7]AS41nA^ - Hn(7Vy. +1)并，S为坐标(x。,y。)的周围"d像素的局部领域图像的熵项，由1数[数8]7V!并且，N:所述局部领域图像中的像素的总数，N,:像素值以数字值为i的像素的总数，N|:像素值以数字值为j的像素的总数，K:常数，通过所述单元(b)的改变，像素值以数字值从i改变为j。
9、如权利要求l所述的图像重构装置，其特征在于，用单元(el)和单元(fl)代替所述单元(e)、 (f)，单元(el)进行如下预约，即根据使用所述差(AE )和控制受理概率的温度参数(T)的受理函数，判断是否受理所述改变，将判断结果反映在所述目前的推定截面像上；单元(fl)，由所述单元(a)至(d)及(el)的一系列的反复处理，每达到规定次数，就使所述单元(el)中的预约反映在所述目前的推定截面像中，并变更一次所述温度参数(T)的值。
10、 一种图像重构装置，从向对象物照射放射线而得到的投影像求解所述对象物的截面像，其特征在于，包括单元(ml),通过没有过滤器的反投影操作，计算所述对象物的放射线投影像P。(^)的反投影像g。(x,y);单元(m2)，从所述对象物的目前的推定截面像f(xj)计算投影像P(r"),并且通过没有过滤器的反投影操作，计算该投影像pO^)的反投影像g(x，少)；单元(m3)，生成将所述对象物的目前的推定截面像f(x,少)的变更值作为像素值的图像AMx，y);单元(m4)，对于各像素值，适用[数9]，生成图像AH(x，y)[数9]A/f =她//2 + 2 A// fe(x。 ， _y0)- go (x0, y0)}其中，M:投影角度的数目；单元(m5)，使用所述AH(x，少)，计算AE(x，y),其中，AE(x，;/)为评价函数E。 (x j)与E, (x, y)的差，E。(x,少)为包含从所述推定截面像f(xj)，通过演算得到的投影像P(r")与所述放射线投影像p。 (r")的差，E,(x,jO为包含从(f(x,少)+A/z(x,;^通过演算得到的投影像Ip(r")+Ap(r,外与所述放射线投影像p。(r^)的差，所述(f(xj)+A^(x，少》将所述推定截面像f(x,y)与由所述单元(m3)得到的图像A/4，j;)相加而得；单元(m6)，对所述AE取正的坐标(x,y),将所述A^(x,y)设为0;单元(m7)，将所述推定截面像f(x,y)与由所述单元(m6)得到的图像A/4,y)相加而得到的结果作为新的推定截面像f(x,少)，反复进行根据所述单元(m2)至(m6)的处理。
11、 一种图像重构装置，从向对象物照射放射线而得到的投影像,求解所述对象物的截面像，其特征在于，包括单元(ml)，通过[数10]，求解所述对象物的放射线投影像p。(r,力的反投影像g。(x,少)，[数10]<formula>formula see original document page 7</formula>其中，r:摄影投影像的1维检测器的信道位置，&投影角度；单元(m2)，从所述对象物的目前的推定截面像f(;cj)，通过演算求解投影像p(r刀)，并且通过[数11j求解该投影像pO^)的反投影像[数ii]<formula>formula see original document page 7</formula>单元(m3 )，生成将相对所述对象物的目前的推定截面像f(x,少)的变更之作为像素值的图像A/z(x,jO;单元(m4),对各像素值适用[数12],生成图像AH(x,少)，[数12]<formula>formula see original document page 7</formula>其中，M:投影角度的数目；单元(m5),对各像素值适用[数13],生成图像AcT(;c,y)，[数13]<formula>formula see original document page 8</formula>然而，cr为坐标(x。,y。)的周围dxd像素的亮度值的标准偏差，由[数14]表示，f,.、 f,分别为变更前后的f(x。,少。)的值，[数14]<formula>formula see original document page 8</formula>但是，[数l习为' 、2\ 1:数16]为<formula>formula see original document page 8</formula>(m6)，对各像素值适用[数17]，生成图像AS(x，少)，<formula>formula see original document page 8</formula>并，S为坐标(x。j。)的周围dxd像素的局部领域图像的熵项，由数18]表示，[数18]<formula>formula see original document page 8</formula>并且:N:所述局部领域图像中的像素的总数，N」像素值以数字值为j的像素的总数，K:常数；单元(m7),根据[数19]计算AE(x，力，[数19]<formula>formula see original document page 9</formula>其中，c:系数，T:假想温度(温度参数)；单元(m8)，对所述AE取正的坐标(x,y),将所述A/z(;c,;;)设为0;单元(m9)，将所述推定截面像fOc，少)与由所述单元(m8)得到的图像A/^(x，少)相加得到的结果作为新的推定截面像f (x, _>');单元(m10)，将所述T放大a倍(《<1 ),反复进行由所述单元(m2)至(m9)的处理。
12、 一种图像重构方法，从向对象物照射放射线而得到的投影像(以下称为"放射线投影像")求解所述对象物的截面像，其特征在于，包括步骤(a),求解评价函数(以下称为"能量")(E。)，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述放射线投影像的差；步骤(b),改变所述目前的推定截面像的一部分；步骤(c)，求解能量(E,)，该能量(Ei)包含从所述改变后的推定截面像，通过演算而得到的投影像和所述放射线投影像的差；步骤(d)，求解所述能量(Eq)与所述能量(E,)的差(ae );步骤(e)，根据使用所述差(ae )与控制受理概率的温度参数(T)的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；步骤(f),所述由步骤(a)至(e)的一系列的反复处理，每达到规定次数，就改变一次所述温度参数(T)的值。
13、 一种图像重构程序，从向对象物照射放射线而得到的投影像(以下称为"放射线投影像")求解所述对象物的截面像的计算机程序，其特征在于，包括步骤(a)，求解评价函数(以下称为"能量")(E。)，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述放射线投影像的差；步骤(b)，改变所述目前的推定截面像的一部分；步骤(c)，求解能量(EJ,该能量(Ei)包含从所述改变后的推定截面像，通过演算而得到的投影像和所述放射线投影像的差；步骤(d)，求解所述能量(Eo)与所述能量(Ei)的差(ae);步骤(e)，根据使用所述差(ae )与控制受理概率的温度参数(T)的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；步骤(f)，所述由步骤(a)至(e)的一系列的反复处理，每达到规定次数，就改变一次所述温度参数(T)的值。
14、 一种CT装置，包括单元(A)，向对象物照射放射线，获取投影像；单元(B)，从所述投影像求解所述对象物的截面像，其特征在于，所述单元(B)包括单元(bl)，求解评价函数(以下称为"能量")(E。)，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述放射线投影像的差；单元(b2)，改变所述目前的推定截面像的一部分；单元(b3)，求解能量(E,)，该能量(E,)包含从所述改变后的推定截面像，通过演算而得到的投影像和所述放射线投影像的差；单元(b4),求解所述能量(Eq)与所述能量(EJ的差(ae);单元(b5)，根据使用所述差(AE )与控制受理概率的温度参数(T)的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；单元(b6)，所述由单元(bl)至(b5)的一系列的反复处理，每达到规定次数，就改变一次所述温度参数(T)的值。
全文摘要
一种层析成像装置及程序，相对通过某种形式求解的暂定截面像f(x，y)，定义评价函数E，该评价函数E包含从f(x，y)计算的投影像与测定的投影像的差，变更f(x，y)，使E大致减少，由此求解对应投影像的截面像。特征在于，本质上无需在通常的计算机层析成像中所需的反投影操作，在金属伪像、人为现象等的去除、降低上尤其效果好。
文档编号A61B6/03GK101553171SQ20078003849
公开日2009年10月7日 申请日期2007年11月13日 优先权日2006年11月13日
发明者西川幸宏 申请人:国立大学法人京都工艺纤维大学