CN103026379A - 推算图像噪音水平的方法 - Google Patents

推算图像噪音水平的方法 Download PDF

Info

Publication number
CN103026379A
CN103026379A CN2012800007264A CN201280000726A CN103026379A CN 103026379 A CN103026379 A CN 103026379A CN 2012800007264 A CN2012800007264 A CN 2012800007264A CN 201280000726 A CN201280000726 A CN 201280000726A CN 103026379 A CN103026379 A CN 103026379A
Authority
CN
China
Prior art keywords
noise
image
psd
technique according
projectional technique
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.)
Granted
Application number
CN2012800007264A
Other languages
English (en)
Other versions
CN103026379B (zh
Inventor
杨智
A·扎米亚京
大石悟
A·姆修-纳塔拉詹
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.)
Canon Medical Systems Corp
Original Assignee
Toshiba Corp
Toshiba Medical Systems Corp
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 Toshiba Corp, Toshiba Medical Systems Corp filed Critical Toshiba Corp
Publication of CN103026379A publication Critical patent/CN103026379A/zh
Application granted granted Critical
Publication of CN103026379B publication Critical patent/CN103026379B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Facsimile Image Signal Circuits (AREA)

Abstract

通过简易方法实质上针对所有图像诊断医用诊断摄影装置都无需人工介入而推算高可靠性的噪音指数。计算与图像内的多个局部区域分别对应的多个局部噪音指数(S10-S50)。根据计算出的噪音指数作成表示噪音水平与次数的对应关系的直方图(S60)。根据作成的直方图推算与图像有关的噪音指数(S70)。

Description

推算图像噪音水平的方法
技术领域
本发明涉及一种推算图像噪音水平的方法。
背景技术
噪音水平的可靠性对图像质量评价至关重要。在噪音水平的评价中需要主观性视觉性评价。在测量主观性评价所需的关心区域内的噪音水平中需要复杂的手动操作。主观性评价耗时且价格昂高,缺乏统一的判断基准。换而言之,主观性评价依赖于操作者。并且,噪音评价的手动方法一般不与用于主观性评价的自动处理互换。
某种噪音不足以表示图像整体。例如,标准偏差(SD)等在某种图像上可靠性低。对此,一个例示性方法是对图像应用拉普拉斯算子消除边缘像素,并根据除去了边缘的图像,利用平均值或方差中的任一个来计算。拉普拉斯法是经由图像处理计算噪音而不是根据原图像直接计算噪音。
在其他以往技术的自动噪音评定方法中,假设噪音依赖于图像诊断医用图像摄影装置,但该噪音医用图像摄影装置复杂。双图像相减法(twin image subtraction method)等其他更简单的现有技术的自动噪音评定方法需要来自CT收集序列的奇数视图和偶数视图的两个重建,但这对于在每次迭代中图像噪音发生变化的迭代处理将使得该方法不具有可行性。
发明内容
目的在于提供一种通过简易方法实质上对于所有图像诊断医用图像摄影装置无需人工介入而推算高可靠性的噪音指数的方法。
本实施方式涉及的图像噪音推算方法首先计算与图像内的多个局部区域分别对应的多个局部噪音指数。根据计算出的噪音指数作成表示噪音水平与次数之间的对应的直方图。根据作成的直方图推算与图像有关的噪音指数。
附图说明
图1为表示根据本发明的一个实施方式的多切片X射线CT装置的结构的图。
图2为表示由根据本实施方式的噪音推算装置执行的处理步骤的整体流程的图。
图3A为表示由根据本实施方式的噪音推算装置执行的处理步骤的详情的流程图。
图3B为表示在本实施方式中,在计算标准偏差时取原图像的平均值从而使原信号平滑化的一些步骤的流程图。
图4A为表示在本实施方式中,具有心肌区域内的2个ROI(关心区域)的肺部CT图像的图。
图4B为在本实施方式中,用于根据利用图4A的肺部CT图像的PSD计算单个噪音指数的直方图。
图5A为在本实施方式中,具有3个ROI(关心区域)的头部CT图像。
图5B为在本实施方式中,用于根据利用图5A的头部CT图像的PSD计算单个噪音指数的直方图。
图6A为在本实施方式中,具有1个ROI(关心区域)的例示性的腹部CT图像。
图6B为在本实施方式中,用于根据利用图6A的腹部CT图像的PSD计算单个噪音指数的直方图。
图7为表示在本实施方式中,比较横跨头部区域内的CT切片的PSD值与SD值的图表。
图8为表示在本实施方式中,比较横跨腹部区域内的CT切片的PSD值与SD值的图表。
图9为表示在本实施方式中,比较横跨心脏区域内的CT切片的PSD值与SD值的图表。
图10A为表示在本实施方式中,与图像整体的尺寸相等的区域有关的多个中间PSD值的图。
图10B为表示在本实施方式中,与尺寸不相等的区域有关的多个中间PSD值的图。
图11为表示在本实施方式中,利用被计算的基于PSD的噪音推算值的这一整体概念的图。
图12为表示在本实施方式中,图像(原图像)的一个例子的图。
图13为表示在本实施方式中,通过对图12的原图像的局部滤波而生成的局部平均值图像的图。
图14为表示在本实施方式中,图13的局部平均值图像与原图像的差图像的图。
图15为表示在本实施方式中,与图14的差图像有关的平方图像的一个例子的图。
图16为表示在本实施方式中,通过对图15的平方图像的卷积运算而生成的卷积图像的图。
图17为表示在本实施方式中,与图16的卷积图像有关的平方根图像的图。
具体实施方式
本实施方式为了生成与图像有关的高可靠性的噪音指数,将抑制使噪音指数的可靠性降低的区域内的影响或提高噪音指数可靠性的区域(特征区域)作为对象来决定噪音指数。特征区域在图像内未必清晰地存在。通过局部的平均值或其时间变化来搜索特征区域。例如,平均值为图像内或局部的像素值中央。因此,本实施方式为了计算在大多类图像中实质上能够信赖的噪音指数,将在最能感知噪音的像素中评价噪音水平的方法与系统作为对象。
在本实施方式中,希望以系统噪音低的组织区域为对象来执行噪音指数的计算。组织区域一般在CT图像或超声波图像中比较容易确定为软组织。在CT图像中,软组织区域对应于HU范围(-1000HU~1000HU)的中央内(-100HU~100HU)。另一方面,肺部图像内的极小构造在噪音指数的计算上是困难的。肿瘤及外伤一般具有相对于周围的健全组织低的对比度,所以从诊断观点来看要求更狭窄的显示范围或显示窗口。CT图像适合于确定图像内的组织区域。在是超声波图像时,以明确地展开统计上稳定的斑点噪音的软组织区域内的像素为对象来计算噪音指数。
可靠的噪音指数的计算对噪音滤波器的有效性评价很重要。噪音滤波器一般适用于减少图像内的某种噪音或伪影。但是,噪音滤波器的有效性评价常常被主观所左右。在其他情况下,减少某种噪音的滤波器实质上减少目标噪音,但有时同一噪音滤波器在图像中产生其它种类的噪音或伪影。由这些理由可以看出,能够信赖的噪音指数的计算对评价减少噪音的滤波器的有效性很重要。
能够信赖的噪音指数也可以用于迭代降噪(iterative noisereduction)处理的结束基准。噪音指数使用具有以中央水平亮度缓慢变化的平均值的像素(有效像素)来计算。这些有效像素一般而言是图像内的像素的大多数。
本实施方式为了计算(推算)可靠性高且实用的噪音指数,引入伪标准偏差(PSD:pseudo-standard deviation)和直方图方法。
直方图方法识别上述有效像素。与边缘、线、或小斑点等快速变化的区域相比,缓慢变化的区域通常形成大的、更多的像素总体。这些因素可以使用于分离缓慢变化的区域的直方图方法正当化,并使分离的区域接受噪音水平分析。在直方图方法中,频度计数的划分(面元尺寸)影响与噪音图像有关的最终指数计算结果的稳健性及准确度。一般,大面元使直方图方法稳定,但使得更不准确。另一方面,小面元使计算准确,但使得更不稳定。在本实施方式中,最初利用比较小的面元尺寸,然后,通过在维持面元的精度的同时利用移动平均,从而放大面元尺寸。
一般,受到不规则尺寸的影响很大的数据在本实施方式中为重建图像数据。在本实施方式中,使用重建图像数据,计算PSD作为噪音指数。计算图像数据内的噪音指数(PSD)的方法实质上可以适用于例如计算测量数据内的噪音指数(PSD)的方法。
图像数据的使用根据本实施方式的几个形态具有对于测量出的数据的几个实用的益处。图像数据不管检测器等输入装置如何都具有像素值的集合体这一统一的性质。因此图像数据相比测量数据一般化,且多用途。至少从这些理由可以看出,图像区域内的图像数据的使用在实施全自动、功能简单的与噪音图像有关的最终指数计算方法时有利。
如图1所示,本实施方式的多切片X摄影CT装置具有机架100和其他装置或单元。机架100具备X射线管101、环状架102、多排型(多列型)或二维阵列型X射线检测器103。X射线管101及X射线检测器103隔着被检体S相对地安装在环状架102上。环状架102围绕轴RA旋转。在沿旋转轴RA移动被检体S期间,旋转单元107以0.4秒/旋转等高速连续地旋转架102。
多切片X射线CT装置具有电流调整器113、为了从X射线管101中产生X射线而对X射线管101施加灯电压的高电压发生器109。例如,高电压发生器109安装于架102。X射线朝被检体S照射。被检体S的断面近似圆形。X射线检测器103为了检测透过了被检体S的X射线,经由被检体S与X射线管101相对。
设置了用于处理来自X射线检测器103的检测信号的其他装置。数据收集电路或数据收集系统(DAS)104针对每一通道将来自X射线检测器103的输出信号转换为电压信号,并将其放大,再将其转换为数字信号。X射线检测器103及DAS104被构成为对每一旋转的庞大的投影数据集进行处理。
DAS104的输出数据经由非接触数据发送器105被发送至收纳在机架100外部的控制台内的预处理装置106。预处理装置106对于原始数据执行灵敏度校正等某种校正。然后,存储装置112保存处于紧接重建处理之前阶段的被称为投影数据的数据。存储装置112与重建装置114、显示装置116、输入装置115及扫描计划支持装置200一起经由数据/控制总线连接于系统控制器110。扫描计划支持装置200包含为了开发扫描计划而支持图像诊断技师的功能。
噪音推算装置117包括各种软件模块与硬件组件的组合。噪音推算装置117具备噪音指数和伪标准偏差(PSD)和直方图的关联功能。噪音推算装置117经由数据/控制总线连接于重建装置114与存储装置112。重建装置114根据存储装置112中保存的投影数据,重建图像。投影数据根据检测器103的输出由数据收集系统104和预处理装置106来生成。图像数据和测量数据分别与图像信号和测量信号同义。在单单称为数据、信号时,分别表示图像数据和测量数据。
噪音推算装置117从重建装置114及/或存储装置112中获取重建图像数据。由于被重建的图像的数据已处于图像区域,因此噪音推算装置117不需要使用了图像数据的与噪音图像有关的最终指数计算处理以外的追加处理。
在噪音推算装置117从重建装置114及/或112中获取测量数据时,因为测量数据为重建处理前的数据,所以噪音推算装置117可能需要与噪音图像有关的最终指数计算处理以外的追加处理。
噪音推算装置117实质上无需人工介入,在计算作为噪音指数的PSD时,对获取到的信号执行任务的规定的设置。自动处理或自动方法需要直方图分析的面元尺寸及用于除去不希望的背景数据的阈值等参数的规定的条件数据集。并且,参数在计算PSD作为噪音指数时还定义规定对信号或数据滤波的处理单位(局部)的内核的尺寸或特性。
图2示出了根据本实施方式的噪音推算装置117的处理步骤。在步骤S10中,原信号D1被卷积规定的均匀内核KⅠ。由此,原信号D1被平滑化。对以构成原图像D1的各个像素为中心的局部应用内核KⅠ。对局部内的像素的像素值乘以内核KⅠ的构成要素,并合计。该计算结果从该局部的中心像素的原像素值置换来的。对所有像素执行相同的局部处理。原信号D1为测量数据或图像数据中的任一个。内核KI典型地为N×N等正方形内核。内核KI的其他例子为N×M等长方形内核。N与M均为整数。通过步骤S10的卷积使图像平滑化。在步骤S20中,计算由步骤S10生成的平滑化信号(一次平滑化图像)与原信号(原图像)D1之间的差。分别对构成平滑化图像的多个像素在与位于原图像D1上的同一位置的像素之间,计算像素值的差。该计算结果根据需要称为“差图像”。
接着,噪音推算装置117执行步骤S30。在步骤S30中,使步骤S20计算出的差信号平方。如果是图像,则使构成差图像的像素的像素值分别平方。该计算结果根据需要称为“平方图像”。
步骤S40是第2次平滑化处理。在步骤S40中,对步骤S30的平方结果(平方图像),通过使用规定的均匀内核KⅡ卷积从而再次使其平滑化。对以构成平方图像的各个像素为中心的局部应用内核KⅡ。对局部内的像素的像素值乘以内核KⅡ的构成要素,并合计。该计算结果从该局部的中心像素的原像素值置换来的。对所有像素执行相同的局部处理。通过第2次平滑化处理,生成二次平滑化图像。内核KⅡ典型地为P×P等正方形内核。内核KⅡ的其他例子为P×Q等长方形内核。P与Q均为整数。内核KⅠ与内核KⅡ典型的情况是不同的,但并不否定相同的情况。
在步骤S50中,计算构成来自步骤S40的二次平滑化图像的各个像素的像素值的平方根。通过这些一系列的S10-S50的计算处理,计算作为噪音指数的伪标准偏差(PSD)。计算结果作为整体示出了PSD分布。
为了生成直方图,构成PSD分布的多个PSD例如按照升序重新排列。PSD数据被记录为一维(1D)向量数据(D3)。接着,对在均等地分割整个PSD范围的多个直方图面元(也称为多个分组、多个划分)各个中所包含的PSD的个数进行计数。垂直轴为PSD的个数或出现次数,水平轴为直方图面元列、即各直方图面元的PSD范围或中心值的排列轴。
在步骤S60中,从直方图的对象中消除表示图像的背景区域的像素。背景像素与其他有用的组织图像相比,对于噪音指数变动的影响大。换而言之,这些不希望的背景像素的选择性消除的失败可能带来直方图内的几个错误的峰值,也就是说由于反映临床上无关联的背景的这些错误的波峰可能降低与图像有关的最终噪音指数的可靠性。最后,在步骤S70中,将作为与图像有关的最终噪音指数的基于PSD的噪音水平决定为基于直方图的分布模式。
噪音指数推算方法的上述方法假设缓慢变化的区域内的像素总体在提高噪音指数的可靠性上比边缘等快速变化的区域内的像素总体有用。通过选择分布模式,从而选择最普通的或频繁产生的伪标准偏差(PSD)的值作为图像整体或信号整体的实质上表示真正的噪音水平的单个噪音指数。
图3A示出了由根据本实施方式的噪音推算装置117执行的噪音推算处理步骤的详情。一般,噪音由于光子计数中的量子噪音、数据收集系统(DAS)104中的电子噪音、模数转换(A/D)中的量化噪音、重建处理的近似误差、其他产生源而产生。假设I(x,y)表示由某加法性噪音破坏的图像,则图像数据I(x,y)由下面的公式(1)来表示。
I(x,y)=Io(x,y)+n(x,y)(1)
这里,Io为无噪音的真值图像,n为由于各种噪音产生源而加在真值上的噪音分量。x及y为像素的二维坐标。图3A的流程图中所包含的步骤为了计算规定的噪音指数PSD,用公式(2)来概括。
PSD ( x , y ) = w ( x , y ) ⊗ [ I ( x , y ) - w ( x , y ) ⊗ I ( x , y ) ] 2 - - - ( 2 )
这里,
Figure BDA00002030487600082
表示卷积运算符,w(x,y)为标准化的均匀的移动平均内核。括号内部的移动平均计算各像素附近的平均值,但括号外部的移动平均计算均方误差(MSE)的附近平均值。希望注意到的是MSE并不是通过像通常进行标准偏差公式那样地减去采样的信号平均值来计算得到的,而是通过减去被滤波后的样本来计算得到的。PSD在(1)的假设有效时接近“真的”标准偏差。PSD用于改善标准偏差(SD)的直接计算的计算效率。
直方图方法用于分离区域并分析噪音水平。该方法基于在大多图像或信号中,缓慢变化的区域通常大,并且具有比边缘、线或小斑点等快速变化的区域大得多的像素总体这一一般发现。如果使用直方图方法,则面元尺寸对与噪音图像有关的最终指数计算结果的稳健性与准确度来说很重要。如上所述,关于面元尺寸,由于具有与噪音图像有关的最终指数计算的稳定性与准确度的权衡,因此比较小的面元尺寸一开始被使用于取得直方图,其后,移动平均被应用于一边保持面元位置的精度一边放大面元尺寸。对此,如果假设h(u)为直方图,m(u)为1维标准化后的移动平均内核,则被滤波后的直方图H(u)通过卷积被定义为下述公式(3)。
H ( u ) = m ( u ) ⊗ h ( u ) - - - ( 3 )
这里,u为PDS值内的面元位置。移动平均窗口的长度(L)为直方图内的最大值的位置(Umax)的单调增加函数。
长度L使用下述公式(4)来计算。
L=Larger odd Nearest Integer of{Umax/k}(4)
另外,L为超过{Umax/K}的最接近的奇数整数。这里,系数k还是数据噪音水平的函数。在Umax小于10时,k=5保证稳定的解,在Umax超过90的期间内,k=3保证稳定的解。另一方面,Umax在10与90之间时,k必须使用接下来的单调减少的公式(5)来选择。
k = 11 60 + 1 600 U max - - - ( 5 )
k的选择可以通过轻的内核或无需平滑化内核而处理比较无变化的直方图,同时通过更重的平滑化内核处理噪音多的直方图。
并且,为了保证准确的与噪音图像有关的最终指数计算,在本实施方式的某种例示性方法与系统中可选地执行以上步骤的多个路径。一般,2个路径就足够了。在第1路径中,第1直方图峰值中的计数用于判断平滑化掩模w(u,v)的最佳尺寸。对于大的计数选择更大的掩模,相反亦然。与掩模尺寸有关的经验性导出的集合的一个例子如下所示。
计数<5000、掩模尺寸=7×7
计数>10000、掩模尺寸=11×11
并且,参照图3A,在步骤S80中,原图像I通过取该图像的移动平均而被平滑化。即,应用规定尺寸的被标准化了的均匀移动平均内核。虚线表示附近以外的图像的展度,实线表示包括中心像素及其附近的像素的局部像素组。对此,该例子中的D1的局部是含有9个像素的3×3附近像素。通过对整个图像I取移动平均,从而作成新的图像IMA。构成该图像IMA的各个像素的值为该3×3附近的平均值。对于特定的附近区域D1,对应的移动平均结果为D4。如图2中已述的那样,步骤S80中的使用标准化的均匀内核的卷积处理的结果在这里是各像素中的图像的附近平均值。在步骤S20中,如减号所示,计算来自步骤S80的D4内的值中的对应值与原信号D1之间的差。并且,在步骤S30中,使差进行平方。在步骤S80A中,被平方后的差的附近平均通过应用标准化的均匀移动平均内核来计算。最后,在步骤S50中,计算来自步骤S80A的结果的平方根。因此,计算出图像内的各像素的伪标准偏差(PSD)。步骤S80及S80A中的标准化的均匀移动平均内核可选地在一个方法中相同,在另一个方法中不同。
参照图3B,图像I被分割成多个块。块尺寸等价于D1。在步骤S90中,计算各块内的像素平均值。平均值被分配给D5内的所有3×3附近像素。在步骤S100中,计算来自步骤S90的平均值与原信号D1之间的差(误差)。并且,在步骤S110中,使误差进行平方。在步骤S90A中,计算被平方的误差的块平均。最后,在步骤S120中,计算来自步骤S90A的结果的平方根。也就是说对局部区域每一个计算通常的标准偏差(SD)。
相对于通常的SD,PSD有优点。PSD通过从标准偏差的计算处理对象中排除边缘变动因素,从而可以抑制SD本质上具备的由边缘移动引起的可靠性的降低。PSD在图像具有固定平均值或缓慢变化的平均值的区域内,接近SD。
对包含模型数据及患者数据的多个数据集应用上述噪音推算方法而推算出的噪音指数表示与介入谨慎的手动操作所计算出的标准偏差等价的结果。在手动评价中,通过手动方式选择出的最佳关心区域(ROI)作为最终的判断基准来被使用。如以上所说明的那样,下一例示性图像为重建的CT图像,但根据本实施方式的噪音推算方法及实施方式不限于重建的CT图像,也可以应用于来自CT的其他测量出的数据及其他医用图像摄影装置。
在此,参照图4A及图4B,对肺部CT图像应用在计算根据本实施方式的基于PSD的通用噪音推算值时所说明的例示性方法。图4A是具有由心肌区域内的白线所示的2个谨慎选择出的ROI(关心区域)的例示性缩放的肺部CT图像。在包含ROI11和ROI112的2个ROI内进行与手动噪音图像有关的最终指数计算。测量出的噪音的标准偏差(SD)在ROI11内为33.25HU,在ROI12内为33.79HU。将这些SD值与基于本实施方式上述所说明的例示性方法的同一肺部图像的PSD的噪音指数进行比较。
这里,参照图4B,生成基于在本实施方式所说明的例示性方法的、用于根据利用图4A的肺部CT图像的PSD计算单个噪音指数的直方图。该直方图根据针对肺部CT图像内的每个像素所计算出的PSD来生成。由于33HU的弱峰值表示最频繁生成的噪音值,因此该峰值自动选择为单个推算出的噪音指数。自动计算的噪音水平是实质上与图4A的ROI11和ROI12的手动测量出的结果一致的33HU。
这里,参照图5A及图5B,对肺部CT图像应用在计算根据本实施方式的基于PSD的通用噪音推算值时所说明例示性方法。图5A为具有由白线所示的3个谨慎选择出的ROI(关心区域)的例示性头部CT图像。在包含ROI11和ROI12和ROI13的3个ROI区域内进行与手动噪音图像有关的最终噪音指数计算。测量出的噪音的标准偏差(SD)在ROI11内为21.49HU,在ROI12内为21.76HU,在ROI13内为23.70HU。将这些SD值与基于在本实施方式所说明的例示性方法的同一头部图像的PSD的噪音指数进行比较。
这里,参照图5B,生成基于在本实施方式中所说明的例示性方法的、用于根据利用图5A的头部CT图像的PSD计算单个噪音指数的直方图。该直方图根据针对头部CT图像内的每个像素计算出的PSD来生成的。由于22HU的弱峰值表示最频繁生成的噪音值,因此该峰值自动选择为单个推算出的噪音指数。自动计算的噪音水平为实质上与图5A的ROI11和ROI12和ROI13的手动测量出的结果一致的22HU。
这里,参照图6A及图6B,对腹部CT图像应用在计算根据本实施方式的基于PSD的通用噪音推算值时所说明的例示性方法。图6A为具有由白线所示的1个谨慎选择出的ROI(关心区域)的例示性腹部CT图像。在ROI内进行与手动噪音图像有关的最终指数计算。测量出的噪音的标准偏差(SD)为11.77HU。将该SD值与基于在根据本实施方式所说明的例示性方法的同一腹部图像的PSD的噪音指数进行比较。
这里,参照图6B,生成基于在本实施方式所说明的例示性方法的、用于根据利用图6A的腹部CT图像的PSD计算单个噪音指数的直方图。该直方图根据针对腹部CT图像内的每个像素计算出的PSD来生成的。由于13HU的弱峰值表示最频繁生成的噪音值,因此该峰值自动选择为单个推算出的噪音指数。自动计算的噪音水平是实质上与图6A的ROI的手动测量出的结果一致的13HU。
除了以上所说明的ROI方法,还在CT图像内的不同的身体区域内的整个体积内将PSD值与切片上的对应的标准偏差(SD)进行比较。图7为表示横跨头部区域内的CT切片的PSD值与SD值的比较的图表。x轴表示切片编号,y轴表示PSD值或SD值中的任一个。实线表示对应的切片上的自动计算出的PSD值,伴有正方形的实线表示手动计算出的SD值。在整个体积的头部内的整个切片内,自动计算出的PSD值实质上与手动计算出的SD值一致。
图8为表示横跨腹部区域内的CT切片的PSD值与SD值的比较的图表。x轴表示切片编号,y轴表示PSD值或SD值中的任一个。实线表示对应的切片上的自动计算出的PSD值,伴有正方形的实线表示手动计算出的SD值。在整个体积的腹部内的整个切片内,自动计算出的PSD值实质上与手动计算出的SD值一致。
图9为表示横跨心脏区域内的CT切片的PSD值与SD值的比较的图表。x轴表示切片编号,y轴表示PSD值或SD值中的任一个。实线表示对应的切片上的自动计算出的PSD值,伴有正方形的实线表示手动计算出的SD值。在整个体积的心脏内的整个切片内,自动计算出的PSD值实质上与手动计算出的SD值一致。
在计算单个噪音推算指数时所说明的方法中,通过根据本实施方式的例示性方法或系统针对整个图像、整个体积或数据的完整的集合自动选择最具有代表性的PSD值。为了保护单个选择方法,在根据本实施方式的另一例示性方法或系统中,在更特定的关心区域或关心范围内根据PSD值计算多个噪音推算指数值,而不是在整个图像、整个体积或数据的完整的集合内根据PSD值计算多个噪音推算指数值。
这里,参照图10A及10B,图中示出了通过根据本实施方式的另一例示性方法或系统在更特定的关心区域或关心范围内根据PSD选择的多个噪音推算指数值。图10A示出了整个图像的相等尺寸的区域的多个中间PSD值。整个图像像素区域I1被分割成9个相等的尺寸的部分A1至A9。一开始按照在本实施方式所说明的方法与系统对图像I1内的每个像素来计算PSD值。其后,在规定的相等尺寸的部分A1至A9中的各个中,而不是在整个图像I1中自动地选择单个最具有单表性的PSD值。由于多个PSD值保持部分A1至A9各个的实质上局部的噪音特性,因此可选地在某种情况下计算选择出的PSD噪音推算值的各个。多个中间噪音推算值可以有利地去除以手动方式消除背景噪音的需要。有利的情况是,多个中间噪音推算值可能处于最可能包含边缘等快速变化的像素总体的图像内。
这里,参照图10B,图中示出了尺寸不相等的区域的多个中间PSD值。整个图像像素区域I2被分割成规定个数的尺寸不等的部分B1和B2。例如,所谓尺寸不等的部分B1和B2分别是一个关心区域(ROI),并且分别包含在医疗上关心的特定构造。一开始按照在本实施方式说明的方法与系统对图像I2内的每个像素来计算PSD值。其后,在规定的尺寸不相等的部分B1和B2各个中,而不是在整个图像I2内自动地选择单个最具有代表性的PSD值。由于多个PSD值保持医疗上重要的区域B1、B2的各自的实质上局部的噪音特性,因此可选地在某种情况下计算选择出的PSD噪音推算值的各个。多个中间噪音推算值可以有利地去除以手动方式消除背景噪音的需要。有利的情况是,多个中间噪音推算值可能处于最可能包含边缘等快速变化的像素总体的图像内。
这里,参照图11,图中示出了由在本实施方式所说明的例示性方法和系统而计算的基于PSD的噪音推算值的利用这一整体概念。该概念按照本实施方式针对整个图像的像素各个计算出PSD值后,以处理这些PSD值的流程来表示。在步骤S100中,按照规定的像素附近尺寸针对各个像素计算PSD值。由于这些PSD值中的每一个实质上都对局部数据敏感,因此如步骤S100所示那样将这些PSD值称为局部PSD值。
在后续的步骤S110中,可选地根据步骤S100的局部PSD值计算中间PSD值。如图10A与图10B所说明的那样,中间PSD值按照规则的规定的集合来计算。如例示的那样,可选地针对图像的相等尺寸或不等尺寸的部分来选择中间PSD值。不论在哪种情况下,中间PSD值都不是整个图像、整个体积、整个切片或整个数据的单个代表性噪音推算值。在实施中间PSD值的计算时,按照本实施方式,针对整个图像、整个体积、整个切片或整个数据的各部分或各区域生成多个直方图。多个基于PSD的值不限定于以上判定基准,可选地根据测量出的信号的角度依赖性等其他数据特性来计算。
并且,参照图11,单个代表性PSD值可选地根据步骤S100的局部PSD值及/或步骤S110的中间PSD值来计算。如图10A和10B所说明的那样,中间PSD值按照规则的规定集合来计算。如在例子中所示的那样,多个中间PSD值可选地针对图像的相等尺寸或不等尺寸的部分来选择。对比而言,单个代表性PSD值按照本实施方式的另一实施方式,作为广域的基于PSD的噪音指数被分配给整个图像、整个体积、整个切片或整个数据。广域的基于PSD的噪音指数可选地与中间PSD值组合来进行计算。在代替方案中,广域的基于PSD的噪音指数可选地作为单个噪音推算指数来计算。
对上述公式(2)具体进行说明。这里,例示了图12所示的图像(原图像)。对原图像I(x,y)的局部卷积滤波器w(x,y)。滤波器w(x,y)的尺寸为3×3个像素,各要素为1/9。被滤波的图像的各像素表示3×3个像素的局部内的平均值。图12中虚线所示的局部的平均值为7.7,如图13所示,该平均值是从该局部的中心像素(被画阴影线的像素)的原值置换来的。如图14所示,计算局部平均值与该中心像素的原值的差。差在公式(2)中表示[]内的值。
图15示出了由该差的平方所得到的平方图像。通过滤波器w(x,y)对平方图像进行卷积处理,取得图16所示的图像。图17示出了作为该平方根的PSD分布。
在图12至图17的例子中,组织边界横穿PSD计算单位(虚线所示的局部)。假设通过通常的步骤来计算SD则示出了2.35。这因为边界也与噪音一样起误差作用。但是,由于PSD对每个像素求取周围的平均值,因此存在所谓更难受到边界影响的特征。
也可以求取PSD以外的噪音拟合。例如,也可以使用SD来代替PSD。通过使用SD,即使在组织边界SD值也变大,但如果计算直方图,那么由于组织边界只是整个图像中的一部分,因此通过PSD得到的结果与通过SD得到的结果大致相同。或者,也可以对图像使用降噪滤波器、以高阶函数进行拟合等取得噪音分布。降噪滤波器从扩散滤波器(Diffusion Filter)、干涉滤波器(Coherent Filter)、中央值滤波器(Median Filter)、高斯滤波器(Gauss Filter)、平均值滤波器等中任意选择。作为拟合方法,采用例如通过B样条(B-spline)拟合、通过三维函数拟合等任意种类的拟合。也可以通过这样基于噪音分布计算SD来求取局部的噪音指数。
但是,在上述说明中与本发明的构造和功能的详情一起示出了本发明的多个特性和特征,但是本公开只是例示,在详情中,特别是在部分形状、尺寸及配置问题与软件、硬件或两者的组合中的实施中,可以进行变更,但希望理解的是其变更包含在由其中表示所附权利要求范围的用语的广义的一般意思所示的完整的范围内的本发明的原理内。

Claims (16)

1.一种图像噪音推算方法,其特征在于,包括:
计算与图像内的多个局部区域分别对应的多个局部噪音指数的步骤;
根据计算出的所述噪音指数,作成表示噪音水平与次数的对应关系的直方图的步骤;以及
根据作成的所述直方图,推算与所述图像有关的噪音指数的步骤。
2.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述局部噪音指数是所述局部区域的标准偏差。
3.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述局部噪音指数是所述局部区域的伪标准偏差。
4.根据权利要求3所述的图像噪音推算方法,其特征在于:
在计算所述局部噪音指数的步骤中包括根据所述图像生成平滑化图像的步骤。
5.根据权利要求4所述的图像噪音推算方法,其特征在于:
通过对所述图像与所述平滑化图像之间的差图像针对每个像素依次实施平方处理、平滑化处理、平方根处理而产生所述伪标准偏差。
6.根据权利要求1所述的图像噪音推算方法,其特征在于,还包括:
根据作成的所述直方图决定所述局部区域的尺寸的好坏的步骤。
7.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述多个局部区域具有相同尺寸。
8.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述多个局部区域具有不同的尺寸。
9.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述多个局部区域包含临床上重要的构造。
10.根据权利要求1所述的图像噪音推算方法,其特征在于:
在计算所述局部噪音指数的步骤中包括:将所述图像转换为降噪图像的步骤;将所述降噪图像和所述图像相减,生成噪音图像的步骤;以及计算所述噪音图像的标准偏差的步骤。
11.根据权利要求1所述的图像噪音推算方法,其特征在于:
通过对所述图像施用降噪滤波器而生成所述降噪图像。
12.根据权利要求1所述的图像噪音推算方法,其特征在于:
通过对所述图像拟合为特定的函数而生成所述降噪图像。
13.根据权利要求1所述的图像噪音推算方法,其特征在于:
在推算与所述图像有关的噪音指数的步骤中,从作成的所述直方图中排除与所述图像的背景部分对应的次数。
14.根据权利要求13所述的图像噪音推算方法,其特征在于:
所述背景部分包含视野以外的部分和与骨骼对应的部分。
15.根据权利要求1所述的图像噪音推算方法,其特征在于:
在作成直方图的步骤中,根据所述直方图的最高次数调整所述局部区域的尺寸。
16.根据权利要求1所述的图像噪音推算方法,其特征在于:
所述图像为三维图像,
对构成所述三维图像的多个二维图像计算所述多个局部噪音指数,
以对所述多个二维图像所计算出的全部所述多个局部噪音指数为对象而作成所述直方图;
根据作成的所述直方图来推算与所述三维图像有关的噪音指数。
CN201280000726.4A 2011-06-14 2012-06-14 推算图像噪音水平的方法 Active CN103026379B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US13/160,111 US8761540B2 (en) 2011-06-14 2011-06-14 Method and system for estimating noise level
US13/160,111 2011-06-14
PCT/JP2012/065279 WO2012173205A1 (ja) 2011-06-14 2012-06-14 画像のノイズレベルを推定する方法

Publications (2)

Publication Number Publication Date
CN103026379A true CN103026379A (zh) 2013-04-03
CN103026379B CN103026379B (zh) 2016-08-31

Family

ID=47353699

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201280000726.4A Active CN103026379B (zh) 2011-06-14 2012-06-14 推算图像噪音水平的方法

Country Status (4)

Country Link
US (1) US8761540B2 (zh)
JP (2) JPWO2012173205A1 (zh)
CN (1) CN103026379B (zh)
WO (1) WO2012173205A1 (zh)

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140056536A1 (en) * 2012-08-27 2014-02-27 Toshiba Medical Systems Corporation Method and system for substantially removing dot noise
US9031297B2 (en) * 2013-02-01 2015-05-12 Kabushiki Kaisha Toshiba Alternative noise map estimation methods for CT images
US9208591B2 (en) 2013-04-18 2015-12-08 International Business Machines Corporation Providing user controlled ability to determine data level of detail in a graph
CN103442266B (zh) * 2013-08-16 2017-02-08 北京视博数字电视科技有限公司 一种获得图像源数据的方法与装置
JP6252750B2 (ja) * 2013-12-03 2017-12-27 株式会社島津製作所 骨梁解析装置
US9406107B2 (en) * 2013-12-18 2016-08-02 General Electric Company System and method of computed tomography signal restoration via noise reduction
JP6333551B2 (ja) * 2013-12-27 2018-05-30 キヤノンメディカルシステムズ株式会社 医用画像診断装置及び画像処理装置
JP6341059B2 (ja) 2014-10-31 2018-06-13 オムロン株式会社 文字認識装置、文字認識方法、およびプログラム
JP6713690B2 (ja) * 2015-04-02 2020-06-24 イービーエム株式会社 血管形状構築装置、その方法及びコンピュータソフトウエアプログラム
WO2018104349A1 (en) 2016-12-06 2018-06-14 Koninklijke Philips N.V. Image noise estimation using alternating negation
US10475214B2 (en) * 2017-04-05 2019-11-12 General Electric Company Tomographic reconstruction based on deep learning
US10674045B2 (en) 2017-05-31 2020-06-02 Google Llc Mutual noise estimation for videos
JP7233850B2 (ja) * 2018-04-25 2023-03-07 キヤノンメディカルシステムズ株式会社 医用画像処理装置、x線診断装置及び医用画像処理プログラム
US11449979B2 (en) * 2020-10-09 2022-09-20 Applied Materials Israel Ltd. Measuring a pattern
JP7395523B2 (ja) * 2021-02-02 2023-12-11 富士フイルムヘルスケア株式会社 医用画像処理装置および医用画像処理方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020071600A1 (en) * 2000-10-17 2002-06-13 Masahiko Yamada Apparatus for suppressing noise by adapting filter characteristics to input image signal based on characteristics of input image signal
US20030161548A1 (en) * 2002-02-22 2003-08-28 Pieter Vuylsteke Noise reduction method
CN1698540A (zh) * 2004-05-17 2005-11-23 Ge医疗系统环球技术有限公司 图像处理方法、图像处理系统以及x-射线ct系统
CN101479762A (zh) * 2006-06-23 2009-07-08 皇家飞利浦电子股份有限公司 用于在包含图像值的图像中确定阈值的方法、系统和计算机程序

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE69331719T2 (de) * 1992-06-19 2002-10-24 Agfa-Gevaert, Mortsel Verfahren und Vorrichtung zur Geräuschunterdrückung
JP4679710B2 (ja) 2000-10-25 2011-04-27 富士フイルム株式会社 ノイズ抑制処理装置並びに記録媒体
CN1305011C (zh) 2001-04-19 2007-03-14 株式会社东芝 图像处理方法和图像处理设备
EP1345170B1 (en) 2002-02-22 2005-01-12 Agfa-Gevaert N.V. Noise reduction method.
EP1566935B1 (fr) * 2004-02-19 2012-07-25 St Microelectronics S.A. Dispositif et procédé de suppression d'interférences impulsionnelles dans un signal
JP4869653B2 (ja) 2005-08-03 2012-02-08 オリンパス株式会社 画像処理装置
JP2007260292A (ja) 2006-03-29 2007-10-11 Toshiba Corp 画像処理装置及びプログラム
JP5224768B2 (ja) * 2007-10-12 2013-07-03 インフォコム株式会社 胸部正面x線画像の検像方法及び検像装置
US7706497B2 (en) * 2008-03-14 2010-04-27 General Electric Company Methods and apparatus for noise estimation for multi-resolution anisotropic diffusion filtering
BR112012012231A2 (pt) * 2009-11-25 2017-12-19 Koninl Philips Electronics Nv método e sistema

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020071600A1 (en) * 2000-10-17 2002-06-13 Masahiko Yamada Apparatus for suppressing noise by adapting filter characteristics to input image signal based on characteristics of input image signal
US20030161548A1 (en) * 2002-02-22 2003-08-28 Pieter Vuylsteke Noise reduction method
CN1698540A (zh) * 2004-05-17 2005-11-23 Ge医疗系统环球技术有限公司 图像处理方法、图像处理系统以及x-射线ct系统
CN101479762A (zh) * 2006-06-23 2009-07-08 皇家飞利浦电子股份有限公司 用于在包含图像值的图像中确定阈值的方法、系统和计算机程序

Also Published As

Publication number Publication date
JP2015024155A (ja) 2015-02-05
JPWO2012173205A1 (ja) 2015-02-23
CN103026379B (zh) 2016-08-31
WO2012173205A1 (ja) 2012-12-20
US20120321157A1 (en) 2012-12-20
US8761540B2 (en) 2014-06-24
JP6038850B2 (ja) 2016-12-07

Similar Documents

Publication Publication Date Title
CN103026379A (zh) 推算图像噪音水平的方法
US12053319B2 (en) Image quality compliance tool
US11399779B2 (en) System-independent quantitative perfusion imaging
CN109389655B (zh) 时变数据的重建
JP6026214B2 (ja) 連続マルチスケール再構成において詳細画像を補うx線コンピュータ断層撮像装置(x線ct装置)、医用画像処理装置及び医用画像処理方法
CN112651885A (zh) 用于减少图像记录噪声的方法和装置
JP2005197792A (ja) 画像処理方法、画像処理装置、プログラム、記憶媒体及び画像処理システム
CN105321155A (zh) 一种cbct图像环形伪影消除方法
CN111540025A (zh) 预测用于图像处理的图像
CN110533738B (zh) 重建数据处理方法、装置、医学成像系统及存储介质
Gao et al. A novel software tool for semi-automatic quantification of thoracic aorta dilatation on baseline and follow-up computed tomography angiography
KR20170117324A (ko) 화상 처리장치, 화상 처리방법, 및 기억매체
US8077954B2 (en) System and method for image processing
US9867586B2 (en) Stereo X-ray tube based suppression of outside body high contrast objects
KR20200021156A (ko) 사전 영상을 이용하는 고품질의 4차원 콘빔 전산화 단층 촬영 시스템
US10217248B2 (en) Method for removing streak from detector cell with performance difference
US10600182B2 (en) Image processing apparatus and image processing method, that determine a conformable image
CN105488824B (zh) 一种重建pet图像的方法和装置
CN117203671A (zh) 基于机器学习的迭代图像重建改进
KR101870890B1 (ko) 콘빔 ct 영상의 금속으로 인한 아티팩트 보정방법 및 장치
CN112446931A (zh) 一种重建数据处理方法、装置、医学成像系统及存储介质
JP2020080913A (ja) 非造影ct画像からの3次元メディアルアクシスモデルに基づく関心臓器画像自動セグメンテーション装置、及び自動セグメンテーション方法
Giacalone et al. An unsupervised spatio-temporal regularization for perfusion MRI deconvolution in acute stroke
Passand Quality assessment of clinical thorax CT images
Fok et al. Cross-modality Training Approach for CT Super-resolution Network

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C41 Transfer of patent application or patent right or utility model
TA01 Transfer of patent application right

Effective date of registration: 20160705

Address after: Japan Tochigi

Applicant after: Toshiba Medical System Co., Ltd.

Address before: Tokyo, Japan, Japan

Applicant before: Toshiba Corp

Applicant before: Toshiba Medical System Co., Ltd.

C14 Grant of patent or utility model
GR01 Patent grant