CN101553171B - 图像重构装置、图像重构方法、图像重构程序、ct装置 - Google Patents

Abstract

一种层析成像装置及程序，相对通过某种形式求解的暂定截面像f(x，y)，定义评价函数E，该评价函数E包含从f(x，y)计算的投影像与测定的投影像的差，变更f(x，y)，使E大致减少，由此求解对应投影像的截面像。特征在于，本质上无需在通常的计算机层析成像中所需的反投影操作，在金属伪像、人为现象等的去除、降低上尤其效果好。

Description

图像重构装置、图像重构方法、图像重构程序、CT装置

技术领域

本发明涉及从对象物的放射线投影像，重构截面像的技术。

背景技术

计算机层析成像(CT)是一种使用X线等光线，对物体从各个角度投影的投影像进行摄影，其后通过计算机获得截面像的技术。图1表示典型的X线CT装置的模式图。

在典型的CT装置中，通过X线光源的旋转移动，取得测定对象的从各个角度投影的投影像。从由此得到的投影像，通过计算操作，获得截面像的技术即为CT。通常，根据Filtered Back-projection法(FBP)，通过变换操作从投影像变换为截面像而获得。在FBP中，在投影像上实施过滤(基本上为微分过滤)后，在原来的投影方向实施“反投影”等操作，由此获得截面像。此时，微分过滤具有放大噪音、误差，容易产生伪像(原来不存在的误差、虚像)等特点。并且，反投影操作具有使其伪像向截面像整体传播的效果。因此，CT中的伪像并不仅仅止于有问题的部分，会破坏截面像的整体，这点构成致命的缺陷。

伪像的多半是由包含在FBP中的过滤操作、反投影操作引起，因此，若不使用FBP，就应该能够得到没有伪像的截面像。AlgebraicReconstruction Technique(ART)作为除FBP以外的计算截面像的方法，在历史上非常重要。ART在提出FBP之前就被利用，是一种有历史的方法。在ART中将截面像的计算作为拟合(fitting)问题加以对待，将截面像作为参数，将投影像作为拟合对象，逐次修正截面像，使从截面像计算出的投影像(p)与通过实验求得的投影像(p₀)一致。在ART中，对截面像渐进地进行改进，使(p-p₀)趋于0，这点即为特征。ART完全不使用FBP之外，但在求解截面像的计算上耗费时间，因此如今只用在特殊的用途(地震波的解析等)。ART不会像FBP那样出现极端的伪像，但一般来说FBP所得出的截面像更自然一些。

另外，除了过滤操作、反操作之外，作为引起伪像的原因，还可以列举投影像中的数据的脱落、过少。若数据少，则作为结果而得出的截面像的画质也跟着降低，这是显而易见的。人们发现在FBP中，数据的脱落、减少同样会产生致命的伪像，这被认为是重大的问题。基于拟合的CT(包括ART)，关于数据的脱落、过少，比FBP强。但是，CT中的数据的脱落，被认为是极端地“条件差”的问题，就算使用ART等，也很难得到改善。在ART中将(p-p₀)作为误差加以考虑，但条件差时拟合容易失败。所谓最小化二乘误差(p-p₀)²的方法，作为拟合而言稳定一些。作为最小化二乘误差的方法，最小二乘法最为标准。在最小二乘法中，求解一边具有相当参数个数的方阵的逆矩阵。CT中的参数是构成截面像的各个像素值，因此，参数的个数十分庞大。比如1000×1000像素的截面像，参数的个数达100万个，矩阵的因子数达1兆个，非常庞大。因此，在正攻法的最小二乘法中矩阵变得过于庞大，从而导致失败。因此，避开正攻法的最小二乘法，提出了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 artifactreduction in helical computedtomography for radiation therapy treatmentplanning：experimental and clinicalstudies.Int J Radiat Oncol BiolPhys，62：1224-1231，2005。

发明内容

从以上论述可知，为了得到没有伪像的截面像，不使用FBP，抗数据脱落强等性质非常重要。因此，像ART那样，不使用FBP，而作为拟合的评价函数利用二乘误差是谁都能想到的。尽管如此，CT中的伪像问题不能得到根本上的解决，这是因为这种单纯的想法不能够简单地实现。第一问题在于，截面像的计算作为拟合问题，规模过于庞大，计算耗费时间。有必要适当进行高速化。第二问题在于，最小二乘法作为拟合算法，不太抗“[条件差”问题。因此，可以认为，根据最小二乘法的已知的算法(ILST等)的延伸并不能解决伪像问题。第三问题在于，就像SIRT、ILST一样，已知的算法不能完全排除FBP。

本发明的第一新颖性在于，作为拟合方法采用了SimulatedAnnealing法(SA)。SA为公知的技术，拟合耗费时间，本质上与CT不符。但是，即使通过硬将SA适用于CT，从而完全不使用FBP，也实现了基于二乘误差最小化的截面像的计算。并且，SA对于数据脱落等条件差的拟合，具有较为稳定的性质，即使在这点上，也有利于伪像的根本上的解决。鉴于以上的情况，通过在CT上适用SA，从而发现了在根本性地解决伪像上的效果，这点即为本发明重要的新颖性。单纯考虑在CT上适用SA时，需无数次的反复进行从截面像求解投影像的计算。真正这样计算，从截面像求解投影像的计算，与FBP大致相同，当在CT上适用SA时，则相比FBP需耗费数百万倍的时间。就算用现在的高速计算机，也需要以年为单位的计算时间。因此，在本发明中，通过式子的变形，大幅降低了计算量，解决了该问题。该办法在将SA适用于CT时为必须的技术，是本发明的重要因素。

第二新颖性在于，作为拟合的评价函数，除了二乘误差，还导入了积极破坏伪像的平滑(smothing)项和熵(entropy)项。在以往的CT中，仅仅将p与p₀的差或二乘误差作为对象。此时，将伪像模糊化，则可能得到较好的拟合结果(误差小)，这种可能性时常存在，而现实中的大多情况为这种情况。即，伪像成为被其他伪像抵消的状态。因此，很难消除伪像。在本发明中，仅仅将二乘误差作为评价函数的情况下，同样也没有完全消除伪像。这种事实意味着，为了获取没有伪像的截面像，除了二乘误差以外，还需要其他条件。在本发明中，从统计力学出发，提出了平滑项、熵项。这些项是将截面像必须为光滑自然的浓淡图像等理所当然的条件用数式形式表现出来的。熵项对拟合进行诱导，使之破坏伪像，从而使截面像整体画质均匀化。平滑项具有对用熵项诱导出的截面像的粗糙度进行控制的效果。通过导入这些项，能够消除伪像，并且能够得到自然的截面像。熵项和平滑项虽然分别单独具有降低伪像的效果，但将两个进行组合的话更有效。

并且，在本发明中，作为比一般定义更为广泛的概念，将包含X线、可视光线、电波的电磁波，由电子、荷电粒子构成的粒子线，作为介质振动的音波等所有的情况，概括称为“放射线”。

作为本发明的第一效果，可以列举起因于数据脱落的伪像的减少带来的效果。作为数据脱落构成问题的情况，可以列举在观察对象上存在不透明部位的情况、投影角度存在限制的情况、投影角度的刻纹不均匀的情况、锥束(conebeam)、螺旋状扫描等三维CT的情况。

尤其能够看出效果的是，在观察对象上存在不透明部位的情况下所表现出来的金属伪像的减少效果。金属伪像是指，在观察对象中，相对X线存在不透明部位(多数情况下为金属部分)时，用CT得到的截面像整体(不仅仅是不透明部位)遭到破坏性错乱的情况。金属伪像产生的原因在于，在不透明部位投影像的亮度的不连续变化、和不透明部位的信息脱落。当FBP中适用微分过滤器时，亮度的不连续变化成为异常值，由反投影操作，生成以不透明部位为中心的筋状的伪像。并且，由于信息的脱落，在与不透明部位没有直接关系的部位上，产生预料之外的明暗。本发明不使用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(x，y)、投影像p(r，θ)的关系的模式图；

图3为典型p₀(r，θ)的例子(横轴取检测器的信道位置、纵轴取投影角度的图像)；

图4为表示本发明基本步骤的流程图；

图5(a)为步骤(I)至(VI)的拟合过程的模式图；

图5(b)为步骤(1)至(9)的拟合过程的模式图；

图6为表示根据“非专利文献1”刊载的方法所得出的结果的图；

图7为从图6(a)通过模拟求解的sinogram(白色部分为不可视领域)；

图8(a)为将图6(c)进一步放大的图；

图8(b)为表示本发明所得出的结果的图；

图9为表示假想能量E的各项的效果的图；

图10为表示具有角度限制的情况下的伪像和本发明的伪像减少效果的图。

具体实施方式

(第1实施方式)

图1表示根据本发明实施方式的X线CT装置的概略结构。该X线CT装置具有摄影部和计算机。摄影部包括用于将X线照射在测定对象的光源，和检测穿透测定对象的X线的检测器，将测定对象的周围用X线从多个方向进行投影，从而取得投影像。计算机包括控制X线CT装置整体的控制部，和根据由摄影部得到的X线投影像，生成测定对象的所关心部位的截面像的图像重构处理部。另外，图1所示的结构在以下的各实施方式中是共通的，而在各实施方式中图像重构处理部所进行的处理分别不同。

本实施方式的图像重构处理部采用Simulated Annealing法(SA)，作为用于从投影像求解截面像的拟合方法。因此，首先对SimulatedAnnealing法(SA)的架构进行说明。SA为从Monte Carlo法派生的拟合算法，其特征在于，根据随机数推进拟合，并模仿热力学处理假想能量、假想温度等。SA本身为公知的技术，按照以下步骤(i)～(vi)进行。

(i)根据随机数选择一个参数，并且根据随机数改变其参数(随机数)。

(ii)进行变更后的评价。作为评价函数考虑假想能量E。在通常的SA中，E为二乘误差的总和。变更前后E的变化记为ΔE(评价)。

(iii)若评价得到改善(ΔE＜0)，则接受该变更(变更)。

(iv)以exp(-ΔE/T)的概率接受改坏。

(v)一点点减少T。

(vi)从(i)反复。

在SA中，如(iv)，按照波尔茨曼统计接受改坏，由此确保能够脱离局部解的可能性。因此，不会拘于局部解，能够达到大域解，对于条件不好的伪像也能稳定操作。并且，通过逐步减少T，软着陆于大域解。自(i)到(v)执行一次称为一个蒙特卡洛步骤。在SA中，无限次反复进行蒙特卡洛步骤，以此进行拟合。计算所需的时间为一个蒙特卡洛步骤所需时间×所需蒙特卡洛步骤数。所需的蒙特卡洛步骤数与参数的个数、自由度成比例。

接着，按照权利要求的形式，对本实施方式进行说明。在将CT考虑成拟合问题的情况下，拟合参数即为截面像(f(x，y))。构成拟合对象的数据即为所测定的投影像(p₀(r，θ))。r为摄影投影像的1维检测器的信道位置，θ为投影角度。坐标的定义表示在图3。p₀(r，θ)为通过改变角度θ测定投影像而得出的数据组，当横轴取r、纵轴取θ时，可以看作2维图像。这种数据称为sinogram。典型的sinogram在图3中表示。

极论而言，CT为从sinogram求解截面像的技术。在本实施方式中，考虑了临时截面像f(x，y)与所测定的投影像p₀(r，θ)的二乘误差，但为了计算二乘误差，使用数101从f(x，y)进行计算，求解投影像p₀(r，θ)。

[数101]

$p(r,\mathrm{\θ})=\underset{s}{\mathrm{\Σ}}f(r\cos \mathrm{\θ}-s\sin \mathrm{\θ},r\sin \mathrm{\θ}+s\cos \mathrm{\θ})$

在图2中可知，数101为在s方向求的和、对f(x，y)的s方向的投影像进行的计算。使用由此得出的p(r，θ)，二乘误差H能够按照以下进行计算。

[数102]

$H=\underset{r,\mathrm{\θ}}{\mathrm{\Σ}}{\{p(r,\mathrm{\θ})-p_0(r,\mathrm{\θ})\}}^2$

在通常的拟合中，作为假想能量E，直接利用数102，但在本实施方式中，其特征在于除H以外导入平滑项和熵项。E的定义如下。

[数103]

E＝H-TS+cσ

T为表示假想温度(温度参数)的系数，S为表示熵项的系数，σ为表示像素值标准偏差的系数，c为表示平滑项强度的系数。TS即为熵项，cσ即为平滑项。S、σ的定义、计算法将在后续中论述。在本实施方式中，使用这些算式，按照如图4所示的步骤计算截面像。

步骤(a)：对以某种形式求解的暂定截面像f(x，y)，求解评价函数(以下称为“能量”)(E₀)(ST10及ST20)。

步骤(b)：用随机数选择(x₀，y₀)及Δμ，改变截面像f(x，y)的一部分(ST30)。

步骤(c)：对改变后的截面像f(x，y)，求解评价函数E₁(ST40及ST50)。

步骤(d)：求解上述能量(E₀)与上述能量(E₁)的差(ΔE)(ST60)。

步骤(e)：根据使用上述差(ΔE)及温度参数(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)中的处理结束的判断，还可以不是由图1的图像重构处理部进行，而是将推定截面像逐次表示在计算机的显示器上，用户看见后指示计算机结束处理。

(第2实施方式)

为了计算数103所示的E，会经过数101和数102，执行关于s、r、θ的级数和。由于计算三重积分(级数和)，会耗费大量时间。即，意味着一个蒙特卡洛步骤所需时间很长。并且，在CT中参数的个数非常庞大，其结果，实施SA所需的总计算时间即使利用现在的计算机，也需要“年单位”。

因此，在本实施方式中，不对E进行计算，而仅仅是主要计算进行变更时的变化量ΔE。

[数104]

ΔE＝ΔH+cΔσ-TΔS

现将有关暂定截面像的变更设为Δf(x，y)。Δf(x，y)为仅仅在坐标(x₀，y₀)上取Δμ的值，除此之外为取0的截面像。该Δf(x，y)的投影像Δp(r，θ)可以用与数101相同的方法计算。使用该Δp(r，θ)，ΔH可以如下计算。

[数105]

$\mathrm{\ΔH}=\underset{r,\mathrm{\θ}}{\mathrm{\Σ}}{\{p(r,\mathrm{\θ})+\mathrm{\Δp}(r,\mathrm{\θ})-p_0(r,\mathrm{\θ})\}}^2-\underset{r,\mathrm{\θ}}{\mathrm{\Σ}}{\{p(r,\mathrm{\θ})-p_0(r,\mathrm{\θ})\}}^2$

将其写下去成为，

[数106]

$\mathrm{\ΔH}=\underset{r,\mathrm{\θ}}{\mathrm{\Σ}}\{\mathrm{\Δp}{(r,\mathrm{\θ})}^2+2\mathrm{\Δp}(r,\mathrm{\θ})[p(r,\mathrm{\θ})-p_0(r,\mathrm{\θ})\left]\right\}$

现在，Δf(x，y)只在(x₀，y₀)取值，因此Δp(r，θ)只在r(θ)＝x₀ cosθ+y₀ sinθ上取值Δμ，除此之外为0。因此，数106的{}中的内容也仅仅在r(θ)＝x₀ cosθ+y₀ sinθ上取值。因此，级数和无需对r和θ的双方进行，可如下表示。

[数107]

$\mathrm{\ΔH}=\underset{\mathrm{\θ}}{\mathrm{\Σ}}\{{\mathrm{\Δ\μ}}^2+2\mathrm{\Δ\μ}[p(r\left(\mathrm{\θ}\right),\mathrm{\θ})-p_0(r\left(\mathrm{\θ}\right),\mathrm{\θ})\left]\right\}$

在数107中，仅仅构成关于θ的级数和，这点非常重要。由此，能够大幅减少计算量。p、p₀为数字图像，因此计算数107时，通常有必要进行关于r(θ)的内插。因此，不能将数107的{}的内容进一步展开。但是，当允许误差而展开数107时可得出下式。

[数108]

$\mathrm{\ΔH}=M{\mathrm{\Δ\μ}}^2+2\mathrm{\Δ\μ}\underset{\mathrm{\θ}}{\mathrm{\Σ}}\{p(r\left(\mathrm{\θ}\right),\mathrm{\θ})-p_0(r\left(\mathrm{\θ}\right),\mathrm{\θ})\}$

但是，M为投影角度的个数。当使用数108来替代数107时，ΔH的值会产生1％左右的误差。但是，数108的速度高于数107，在本发明中的利用价值高。

以下对Δσ的计算进行说明。σ为坐标(x₀，y₀)附近的亮度值的标准偏差。作为(x₀，y₀)附近，考虑(x₀，y₀)周围d×d像素(本实施例中使用d＝5)。这些像素的标准偏差可通过下式求解。

[数109]

$\mathrm{\σ}=\sqrt{\⟨f{(x_0,y_0)}^2\⟩-{\⟨f(x_0,y_0)\⟩}^2}$

但是，

[数110]

$\⟨f{(x_0,y_0)}^2\⟩=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\Σ}}}f{(x_0+i,y_0+j)}^2$

[数111]

$\⟨f(x_0,y_0)\⟩=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\Σ}}}f(x_0+i,y_0+j)$

若使用数109～111，则Δσ可如下计算。

[数112]

$\mathrm{\Δ\σ}=\sqrt{\⟨f{(x_0,y_0)}^2\⟩-\frac{{f_i}^2-{f_j}^2}{d^2}-{\{\⟨f(x_0,y_0)\⟩-\frac{f_i-f_j}{d^2}\}}^2}$

但是，f_i，f_j为变更前及后的f(x₀，y₀)的值。

以下是对于ΔS的计算，在此之前，对平滑项S进行定义。一般计算机所处理的图像为数字图像，不仅是像素的坐标(x，y)，就连像素的值也是数字值。因此，一个像素看作一个量子，将像素值视为量子状态。那么，图像可以看成由多个量子构成的系。按照统计力学，对这种系的熵进行如下定义。

[数113]

$S=k\ln \frac{N!}{N_1!N_2!\·\·\·N_i!\·\·\·N_n!}$

N为图像中像素的总数，N_i是像素值以数字值为i的像素的总数。K在通常的物理中为波尔茨曼常数，但本发明与物理没有直接关系，因此任取。在本发明中，与σ一样，S也被定义为(x₀，y₀)周围的d×d像素。

对于被定义成数113的S考虑ΔS。通过变更，像素值从数字值i变更为j。那么ΔS变成如下式子。

[数114]

$\mathrm{\ΔS}=k\ln \frac{N!}{N_1!N_2!\·\·\·(N_i-1)!\·\·\·(N_j+1)!\·\·\·N_n!}-k\ln \frac{N!}{N_1!N_2!\·\·\·N_i!\·\·\·N_j!\·\·\·N_n!}$

即，通过变更，N_i减1，而N_j加1。若展开数114，则可以得到非常简单的式子。

[数115]

ΔS＝klnN_i-kln(N_j+1)

综上，本实施方式的顺序如下：

(I)用随机数选择(x₀，y₀)及Δμ。

(II)使用数104、108(或107)、112、115计算ΔE。

(III)若评价得到改善(ΔE＜0)，则在f(x₀，y₀)上加Δμ。

(IV)ΔE＞0时，也同样以概率exp(-ΔE/T)，在f(x₀，y₀)上加Δμ。

(V)一点点减小T。

(VI)从(I)重复。

(VII)判断是否满足规定的停止条件，若满足则结束处理。

在这里，将在(III)、(IV)中附加了Δμ的情况设为“成功”，除此之外的情况设为[失败]，当成功概率低于规定值(例如10％，该值可适当调整)时结束。然后，在结束时所采用的推定截面像作为测定对象的截面像，将其显示在计算机的显示器上，或存储到存储介质。

可以看出，步骤(I)～(VI)完全遵循SA的基本步骤(i)～(vi)。然而，步骤(I)～(VII)的处理由图1中的图像重构处理部来执行。

另外，对于上述步骤(VII)的处理结束的判断，还可以不是由图1中的图像重构处理部来执行，而是逐次表示在计算机的显示器上，用户肉眼观察后指示计算机结束处理。

(第3实施方式)

继而，对上述方法的重要变形进行论述。在上述的方法中，利用数107或数108求解了有关θ的级数和。因此，研究一个像素的变更需进行M次的反复演算。本实施方式进一步减少了该部分的计算量。

首先，考虑p(r，θ)的反投影像g(x，y)。

[数116]

$g(x,y)=\underset{\mathrm{\θ}}{\mathrm{\Σ}}p(x\cos \mathrm{\θ}+y\sin \mathrm{\θ},\mathrm{\θ})$

在数116中，没有进行过滤操作，因此g(x，y)虽与截面像f(x，y)相似，但会成为模糊的图像。使用g(x，y)时，数108变成如下。

[数117]

ΔH＝MΔμ²+2Δμ{g(x₀，y₀)-g₀(x₀，y₀)}

但是，g₀(x，y)为p₀(r，θ)的反投影像，可以与数116相同的方法进行计算。数117在没有反复进行运算的这一点上优于数108。g₀(x，y)在计算过程中完全没有改变，因此可以事先计算。另外，g(x，y)在每更改f(x，y)的1点时，就会产生微小的变化，因此严格来说，每经过一次(I)～(VII)步骤，有必要进行重新计算。但是，当假定f(x，y)的改变带来的g(x，y)的变化微小时，可以适用其他步骤，这就是本实施方式。以下表示步骤。然而，以下的步骤(1)～(9)的处理由图1中的图像重构处理部来执行。

(1)从p₀(r，θ)求解g₀(x，y)。

(2)从f(x，y)求解p(r，θ)，并且从p(r，θ)求解g(x，y)。

(3)用随机数生成将f(x，y)的变更值作为像素值的图像Δμ(x，y)。

(4)对各图像值适用数117，计算图像ΔH(x，y)。

(5)同样对各像素值适用数112、数115，计算图像Δσ(x，y)，ΔS(x，y)。

(6)根据数104计算ΔE(x，y)。

(7)对于ΔE为正的坐标(x，y)，将Δμ(x，y)设为0。

(8)向f(x，y)加上Δμ(x，y)。

(9)将T放大α倍(α＜1)，从(2)重复。

(10)判断是否满足规定的停止条件，若满足则结束处理。其中，在(7)中的Δμ(x，y)成为0的情况设为“失败”，不是0的情况设为“成功”，在成功概率低于规定值(例如10％，该值可适当调整)时结束处理。然后，将在结束时所采用的推定截面像作为测定对象的截面像，将其显示在计算机的显示器上，或存储到存储介质。

然而，对于上述步骤(10)的处理结束的判断，可以不是由图1中的图像重构处理部来执行，而是将推定截面像逐次显示在计算机的显示器上，用户肉眼观察后指示计算机结束处理。

图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)示出刊载在“非专利文献1”的方法所带来的结果。

图6中，左栏(a)(d)(g)为没有Metal Artifact的原版模型图像。中央(b)(e)(h)为在原版模型的规定位置设定不可视领域，通过通常的FBP法重构的图像。右栏(c)(f)(i)为与上述相同，设定不可视领域，通过“非专利文献1”的方法重构的图像，在以往技术中去除Metal Artifact最好的一个例子。

因此，从“非专利文献1”中取图6(a)的图像，在与文献相同的位置设定不可视领域，模拟计算投影像。结果表示在图7。图7对应本实施例的p₀(r，θ)。对于图7，将根据本发明重构的结构及用于比较的图6(c)的扩大图表示在图8。效果一目了然。在本实施方式(b)中即使想找出Metal Artifact的特征都困难。

但是在本实施方式(b)中边缘附近有些模糊。这是因为“使截面像光滑的因子”(熵项cσ)有些过强。可以看出，在Metal Artifact的减少上，需稍强地设置该因子。由此按照本实施方式，在计算假想能E时的各因子的平衡上，还有改进的余地，但Metal Artifact去除效果非常高。

另外，为参考起见，将没有对平滑项cσ、熵项TS进行设定而执行本实施方式算法的情况的结果表示在图9中。图9(a)表示没有平滑项cσ(仅有熵项TS)时的结果，图9(b)表示没有熵项TS(仅有平滑项cσ)时的结果，图9(c)表示即没有熵项TS也没有平滑项cσ时的结果，当比较图9和图8(b)时，可以看出，在本实施方式中，作为能量E的因子设定平滑项cσ和熵项TS的两方的情况，比仅仅设定某一个，和两方都不设定的情况，Metal Artifact去除效果更好。只有平滑项时，具有在伪像中以不透明为中心残留有呈放射状扩展的条纹的特点。另外，只有熵项时，具有截面像粗糙，S/N低的特点。并且，最终的画质，慢冷过程中的平滑项的系数c和假想温度T之比非常重要。

接着，将具有旋转角度限制的情况的实施例表示在图10。图10(a)表示仅根据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技术的可能性。

Claims (13)

1.一种图像重构装置，从向对象物照射放射线而得到的投影像求解所述对象物的截面像，其特征在于，包括：

单元a，求解评价函数E₀，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述从向对象物照射放射线而得到的投影像的差；

单元b，改变所述目前的推定截面像的一部分；

单元c，求解评价函数E₁，该评价函数E₁包含从改变后的推定截面像，通过演算而得到的投影像和所述从向对象物照射放射线而得到的投影像的差；

单元d，求解所述评价函数E₀与所述评价函数E₁的差ΔE；

单元e，根据使用所述差ΔE与控制受理概率的温度参数T的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；

单元f，由所述单元a至e的一系列的反复处理，每达到规定次数，就改变一次所述温度参数T的值。

2.如权利要求1所述的图像重构装置，其特征在于，所述单元a求解评价函数E₀，该评价函数E₀包含：

从所述目前的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述目前的推定截面像的局部领域的标准偏差的和；

所述单元c求解评价函数E₁，该评价函数E₁包含：

从所述改变后的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述改变后的推定截面像的局部领域的标准偏差的和。

3.如权利要求1所述的图像重构装置，其特征在于，所述单元a求解评价函数E₀，该评价函数E₀包含：

从所述目前的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述目前的推定截面像的局部领域的熵项的和；

所述单元c求解评价函数E₁，该评价函数E₁包含：

从所述改变后的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述改变后的推定截面像的局部领域的熵项的和。

4.如权利要求1所述的图像重构装置，其特征在于，所述单元a求解评价函数E₀，该评价函数E₀包含：

从所述目前的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述目前的推定截面像的局部领域的标准偏差，

与所述目前的推定截面像的局部领域的熵项的和；

所述单元c求解评价函数E₁，该评价函数E₁包含：

从所述改变后的推定截面像，通过演算得到的投影像与所述从向对象物照射放射线而得到的投影像的差，

与所述改变后的推定截面像的局部领域的标准偏差，

与所述改变后的推定截面像的局部领域的熵项的和。

5.如权利要求1所述的图像重构装置，其特征在于，用单元h替代所述单元a、c、d，

所述单元h根据数1算出ΔH，并求解ΔE，该ΔE作为构成要素包含通过数1算出的ΔH，

数1

$\mathrm{\ΔH}=\underset{\mathrm{\θ}}{\mathrm{\Σ}}\{\mathrm{\Δ}{\mathrm{\μ}}^2+2\mathrm{\Δ\μ}[p(r\left(\mathrm{\θ}\right),\mathrm{\θ})-p_0(r\left(\mathrm{\θ}\right),\mathrm{\θ})\left]\right\}$

其中，当将所述对象物的目前的推定截面像设为f(x，y)、将根据所述单元b的改变部分设为Δf(x，y)时，Δf(x，y)仅在坐标(x₀，y₀)取Δμ的值，除此之外为0的截面像，p(r，θ)为从所述对象物的目前的推定截面像，通过演算得到的投影像，p₀(r，θ)为从向对象物照射放射线而得到的投影像，r为摄影投影像的1维检测器的信道位置，θ为投影角度，r(θ)＝x₀cosθ+y₀sinθ。

6.如权利要求1所述的图像重构装置，其特征在于，用单元h替代所述单元a、c、d，

所述单元h根据数2算出ΔH，并求解ΔE，该ΔE作为构成要素包含通过数2算出的ΔH，

数2

$\mathrm{\ΔH}=\mathrm{M\Δ}{\mathrm{\μ}}^2+2\mathrm{\Δ\μ}\underset{\mathrm{\θ}}{\mathrm{\Σ}}\{p(r\left(\mathrm{\θ}\right),\mathrm{\θ})-p_0(r\left(\mathrm{\θ}\right),\mathrm{\θ})\}$

其中，当将所述对象物的目前的推定截面像设为f(x，y)、将根据所述单元b的改变部分设为Δf(x，y)时，Δf(x，y)仅在坐标(x₀，y₀)取Δμ的值，除此之外为0的截面像，p(r，θ)为从所述对象物的目前的推定截面像通过演算得到的投影像，p₀(r，θ)为从向对象物照射放射线而得到的投影像，r为摄影投影像的1维检测器的信道位置，θ为投影角度，r(θ)＝x₀cosθ+y₀sinθ，M为投影角度的数目。

7.如权利要求5或6所述的图像重构装置，其特征在于，所述单元h根据数3算出Δσ，并求解ΔE，该ΔE作为构成要素包含cΔσ与所述ΔH的和，所述cΔσ为根据数3算出的Δσ与系数c的相乘，

数3

$\mathrm{\&Delta;\&sigma;}=\sqrt{<f{(x_0,y_0)}^2>-\frac{f_i^2-f_j^2}{d^2}-{\{<f(x_0,y_0)>-\frac{f_i-f_j}{d^2}\}}^2}$

并，σ为坐标(x₀，y₀)的周围d×d像素的亮度值的标准偏差，由数4表示，f_i、f_j分别为根据所述单元b的变更前后的f(x₀，y₀)的值，

数4

$\mathrm{\&sigma;}=\sqrt{<f{(x_0,y_0)}^2>-{<f(x_0,y_0)>}^2}$

但是，数5为，

$<f{(x_0,y_0)}^2>=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\&Sigma;}}}f{(x_0+i,y_0+j)}^2$

数6为，

$<f(x_0,y_0)>=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\&Sigma;}}}f(x_0+i,y_0+j).$

8.如权利要求5或6所述的图像重构装置，其特征在于，所述单元h根据数7算出ΔS，并求解ΔE，所述ΔE作为构成要素包含-TΔS与ΔH的和，所述-TΔS为通过数7算出的ΔS与温度参数T的相乘，

数7

ΔS＝klnN_i-kln(N_j+1)

并，S为坐标(x₀，y₀)的周围d×d像素的局部领域图像的熵项，由数8表示，

数8

$S=k\ln \frac{N!}{N_1!N_2!...N_i!...N_n!}$

并且，

N：所述局部领域图像中的像素的总数，

N_i：像素值以数字值为i的像素的总数，

N_j：像素值以数字值为j的像素的总数，

K：常数，

通过所述单元b的改变，像素值以数字值从i改变为j。

9.如权利要求1所述的图像重构装置，其特征在于，用单元e1和单元f1代替所述单元e、f，

单元e1进行如下预约，即根据使用所述差ΔE和控制受理概率的温度参数T的受理函数，判断是否受理所述改变，将判断结果反映在所述目前的推定截面像上；

单元f1，由所述单元a至d及e1的一系列的反复处理，每达到规定次数，就使所述单元e1中的预约反映在所述目前的推定截面像中，并变更一次所述温度参数T的值。

10.一种图像重构装置，从向对象物照射放射线而得到的投影像求解所述对象物的截面像，其特征在于，包括：

单元m1，通过没有过滤器的反投影操作，计算所述从向对象物照射放射线而得到的投影像p₀(r，θ)的反投影像g₀(x，y)；

单元m2，从所述对象物的目前的推定截面像f(x，y)计算投影像p(r，θ)，并且通过没有过滤器的反投影操作，计算该投影像p(r，θ)的反投影像g(x，y)；

单元m3，生成将所述对象物的目前的推定截面像f(x，y)的变更值作为像素值的图像Δμ(x，y)；

单元m4，对于各像素值，适用数9，生成图像ΔH(x，y)

数9

ΔH＝MΔμ²+2Δμ{g(x₀，y₀)-g₀(x₀，y₀)}

其中，M：投影角度的数目；

单元m5，使用所述ΔH(x，y)，计算ΔE(x，y)，

其中，

ΔE(x，y)为评价函数E₀(x，y)与E₁(x，y)的差，

E₀(x，y)为包含从所述推定截面像f(x，y)，通过演算得到的投影像p(r，θ)与所述从向对象物照射放射线而得到的投影像p₀(r，θ)的差，

E₁(x，y)为包含从{f(x，y)+Δμ(x，y)}通过演算得到的投影像{p(x，θ)+Δp(r，θ)}与所述从向对象物照射放射线而得到的投影像p₀(r，θ)的差，所述{f(x，y)+Δμ(x，y)}将所述推定截面像f(x，y)与由所述单元m3得到的图像Δμ(x，y)相加而得；

单元m6，对所述ΔE取正的坐标(x，y)，将所述Δμ(x，y)设为0；

单元m7，将所述推定截面像f(x，y)与由所述单元m6得到的图像Δμ(x，y)相加而得到的结果作为新的推定截面像f(x，y)，反复进行根据所述单元m2至m6的处理。

11.一种图像重构装置，从向对象物照射放射线而得到的投影像，求解所述对象物的截面像，其特征在于，包括：

单元m1，通过数10，求解所述从向对象物照射放射线而得到的投影像p₀(r，θ)的反投影像g₀(x，y)，

数10

$g_0(x,y)=\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}{p}_0(x\cos \mathrm{\&theta;}+y\sin \mathrm{\&theta;},\mathrm{\&theta;})$

其中，

r：摄影投影像的1维检测器的信道位置，

θ：投影角度；

单元m2，从所述对象物的目前的推定截面像f(x，y)，通过演算求解投影像p(r，θ)，并且通过数11求解该投影像p(r，θ)的反投影像g(x，y)；

数11

$g(x,y)=\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}p(x\cos \mathrm{\&theta;}+y\sin \mathrm{\&theta;},\mathrm{\&theta;})$

单元m3，生成将相对所述对象物的目前的推定截面像f(x，y)的变更值作为像素值的图像Δμ(x，y)；

单元m4，对各像素值适用数12，生成图像ΔH(x，y)，

数12

ΔH＝MΔμ²+2Δμ{g(x₀，y₀)-g₀(x₀，y₀)}

其中，

M：投影角度的数目；

单元m5，对各像素值适用数13，生成图像Δσ(x，y)，

数13

$\mathrm{\&Delta;\&sigma;}=\sqrt{<f{(x_0,y_0)}^2>-\frac{f_i^2-f_j^2}{d^2}-{\{<f(x_0,y_0)>-\frac{f_i-f_j}{d^2}\}}^2}$

然而，σ为坐标(x₀，y₀)的周围d×d像素的亮度值的标准偏差，由数14表示，f_i、f_j分别为变更前后的f(x₀，y₀)的值，

数14

$\mathrm{\&sigma;}=\sqrt{<f{(x_0,y_0)}^2>-{<f(x_0,y_0)>}^2}$

但是，数15为

$<f{(x_0,y_0)}^2>=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\&Sigma;}}}f{(x_0+i,y_0+j)}^2$

数16为

$<f(x_0,y_0)>=\frac{1}{d^2}\underset{i,j=-d/2}{\overset{d/2}{\mathrm{\&Sigma;}}}f(x_0+i,y_0+j);$

单元m6，对各像素值适用数17，生成图像ΔS(x，y)，

数17

ΔS＝klnN_i-kln(N_j+1)

并，S为坐标(x₀，y₀)的周围d×d像素的局部领域图像的熵项，由数18表示，

数18

$S=k\ln \frac{N!}{N_1!N_2!...N_i!...N_n!}$

并且，

N：所述局部领域图像中的像素的总数，

N_i：像素值以数字值为i的像素的总数，

N_j：像素值以数字值为j的像素的总数，

K：常数；

单元m7，根据数19计算ΔE(x，y)，

数19

ΔE＝ΔH+cΔσ-TΔS

其中，

c：系数，

T：假想温度，所述假想温度为温度参数；

单元m8，对所述ΔE取正的坐标(x，y)，将所述Δμ(x，y)设为0；

单元m9，将所述推定截面像f(x，y)与由所述单元m8得到的图像Δμ(x，y)相加得到的结果作为新的推定截面像f(x，y)；

单元m10，将所述T放大α倍，其中，α＜1，反复进行由所述单元m2至m9的处理。

12.一种图像重构方法，从向对象物照射放射线而得到的投影像求解所述对象物的截面像，其特征在于，包括：

步骤a，求解评价函数E₀，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述从向对象物照射放射线而得到的投影像的差；

步骤b，改变所述目前的推定截面像的一部分；

步骤c，求解评价函数E₁，该评价函数E₁包含从改变后的推定截面像，通过演算而得到的投影像和所述从向对象物照射放射线而得到的投影像的差；

步骤d，求解所述评价函数E₀与所述评价函数E₁的差ΔE；

步骤e，根据使用所述差ΔE与控制受理概率的温度参数T的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；

步骤f，由所述步骤a至e的一系列的反复处理，每达到规定次数，就改变一次所述温度参数T的值。

13.一种CT装置，包括：

单元A，向对象物照射放射线，获取投影像；

单元B，从所述投影像求解所述对象物的截面像，

其特征在于，所述单元B包括：

单元b1，求解评价函数E₀，该评价函数包含从所述对象物的目前的推定截面像，通过演算而得到的投影像与所述从向对象物照射放射线而得到的投影像的差；

单元b2，改变所述目前的推定截面像的一部分；

单元b3，求解评价函数E₁，该评价函数E₁包含从所述改变后的推定截面像，通过演算而得到的投影像和所述从向对象物照射放射线而得到的投影像的差；

单元b4，求解所述评价函数E₀与所述评价函数E₁的差ΔE；

单元b5，根据使用所述差ΔE与控制受理概率的温度参数T的受理函数，判断是否受理所述改变，并将判断结果反映在所述目前的推定截面像上；

单元b6，由所述单元b1至b5的一系列的反复处理，每达到规定次数，就改变一次所述温度参数T的值。

Images