CN103069432A - 医学图像的非线性分辨率降低方法 - Google Patents

医学图像的非线性分辨率降低方法 Download PDF

Info

Publication number
CN103069432A
CN103069432A CN2011800322137A CN201180032213A CN103069432A CN 103069432 A CN103069432 A CN 103069432A CN 2011800322137 A CN2011800322137 A CN 2011800322137A CN 201180032213 A CN201180032213 A CN 201180032213A CN 103069432 A CN103069432 A CN 103069432A
Authority
CN
China
Prior art keywords
image
pixel
resolution
noise
value
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.)
Pending
Application number
CN2011800322137A
Other languages
English (en)
Inventor
埃泽尔·巴拉维夫
伊利兰·达昂
兹维·德维尔
塔尔·克尼希
盖伊·罗斯曼
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.)
MEDIC VISION BRAIN TECHNOLOGIE
Medic Vision Imaging Solutions Ltd
Original Assignee
MEDIC VISION BRAIN TECHNOLOGIE
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 MEDIC VISION BRAIN TECHNOLOGIE filed Critical MEDIC VISION BRAIN TECHNOLOGIE
Publication of CN103069432A publication Critical patent/CN103069432A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00Manipulating 3D models or images for computer graphics
    • G06T19/20Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
    • 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
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20016Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30016Brain

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Architecture (AREA)
  • Computer Graphics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

一种用于从三维(3D)高分辨率医学图像产生3D低分辨率图像的方法。所述3D高分辨率医学图像是多个2D图像组成,这些2D图像在轴向方向上具有轴向分辨率,所述多个2D图像是沿着所述轴向方向而获得的,所述方法包括以下步骤:分解多个2D图像中的每个图像,产生多个转换的数据组,每个转换的数据组对应于所述多个2D图像中的一个图像;加权在转换的数据组内的每个像素,以产生多个加权的转换的数据组;合并多个加权的转换的数据组成为单独的新的转换的数据组;以及产生3D低分辨率图像,该图像在所述3D低分辨率图像的轴向方向上具有第一分辨率,所述第一分辨率是小于或等于所述轴向分辨率。

Description

医学图像的非线性分辨率降低方法
相关申请(仅用于美国)
本专利申请要求享有申请号为61/359,845、申请日为2010年6月30日美国临时专利申请的权益。本专利申请也是申请号为13/142,282、申请日为2011年6月27日的美国专利申请的部分延续申请,该美国专利申请是PCT申请号为PCT/IL2009/001216、申请日为2009年12月24的PCT专利申请的国家阶段申请,该PCT申请要求了以下优先权:PCT申请号为PCT/IL2008/001679、申请日为2008年12月25的PCT专利申请,以及申请号为61/219,857、申请日为2009年6月24日的美国临时专利申请。
本发明的技术领域与背景技术
本发明,在它的一些具体实施方式中,涉及一种用于处理医学图像的方法与系统,以产生具有降低的噪声与其他期望的特征的图像,更特别地,但不排除,涉及一种处理CT图像的方法,该方法考虑到噪声在图像中的非均一分布和/或采用非线性滤波器来保留精细的图像细节。
文献E.H.Adelson,C.H.Anderson,J.R.Bergen,P.J.Burt and J.M.Ogden,“Pyramidmethods in image processing”,RCA Engineer,29-6,Nov.1984描述了一种将一个场景以不同照相机焦距设置拍摄的两个图像融合的方法,采用了拉普拉斯金字塔图像分解技术。
文献Hui Li,B.S.Manjunath,and S.K.Mitra,″Multi-sensor image fusion using thewavelet transform″in Proceedings of IEEE International Conference on Image Processing,1994描述了相同区域的不同图像的融合,采用不同类型的传感器来获得,该方法采用小波变换代替拉普拉斯金字塔技术。
由北京大学重离子物理研究所(中国北京,邮编100871)于2001年所发表的报告Yu Lifeng,Zu Donglin,Wan Weidong与Bao Shanglian,″Multi-Modality MedicalImage Fusion Based on Wavelet Pyramid and Evaluation″描述了一种融合两个医学图像的方法,这些图像是采用不同的图像形态特征来获得的,例如,CT和MRI,或PET和MRI,采用小波金字塔技术。
文献Hassam El-Din Moustafa and Sameh Rehan,″Applying Image FusionTechniques for Detection of Hepatic Lesions″Proceedings of the 6th WSEASInternational Conference on Wavelet Analysis & Multirate Systems,Bucharest,Romania,October 16-18,2006,pages 41-44比较了采用不同方法融合CT和MRI图像的结果,这些方法包括拉普拉斯金字塔、小波变换、计算有效像素级图像融合方法、以及基于空间频率的多焦点技术。
文献Richard Alan Peters II,″A New Algorithm for Image Noise Reduction usingMathematical Morphology″,IEEE Trans.Image Processing 4,554-568(1995)描述了一种形态图像清洗算法,当去噪时保留薄特征。该方法通过形态大小分布而计算了在一定数量的不同场景上的残留图像,并丢弃了在不同残留图像内判断为包含噪声的区域,所提供的噪声具有比薄特征更小的动态范围。
Garnier等人的美国专利申请US2008/0118128描述了一种产生模拟图像的方法,该图像具有设定大小的人工产生的噪声,该噪声加入到图像中。
Toth的美国专利申请US2008/0118128描述了一种产生终端图像的方法,该图像具有预定量的人工产生的噪声。
下列出版物和专利涉及图像处理噪声降低、图像采集和/或计算机视觉:
US 2007/053477——锥形束CT成像与扇束CT成像的全局去噪方法与装置;
KR 2005-0031210——图像去噪的方法与装置;
JP 2000-050109——用于去除噪声的非线性图像滤波器;
US 6,459,755——用于低剂量给药CT扫描的方法与装置;
US 2003/099405——带有计算有效实施的CT剂量减薄滤波器;
EP 1 774 837——活性剂量减薄装置与方法;
JP 2001-39874——用于MRI的磁场发生器;
WO 2007/047599——用于高增益磁共振的方法与装置;
Steven Haker,Lei Zhu,Allen Tannenbaum,and Sigurd Angenent,“Optimal MassTransport for Registration and Warping”,International Journal of Computer Vision,Volume 60,Issue 3(December 2004),Pages 225–240;
Yossi Rubner,Carlo Tomasi,and J.Leonidas Guibas,“A Metric for Distributionswith Applications to Image Databases”,ICIP 1998,Pages 59–66;
Belongie Serge,Jitendra Malik,and Puzicha Jan,“Shape Matching and ObjectRecognition Using Shape Contexts”,IEEE T-PAMI,Volume 24,No.4,(April 2002);
Robert Osada,Thomas Funkhouser,Bernard Chazelle,and David Dobkin,“Matching3D Models with Shape Distributions”,Proceedings of the International Conference onShape Modeling & Applications 2001,Pages 154–166;
P.J.Burt and E.H.Adelson,“The Laplacian Pyramid as a Compact Image Code”,IEEE Trans.on Communications,pp.532–540,April 1983;
Iddo Drori,Daniel Cohen-Or,and Hezy Yeshurun,“Fragment based imagecompletion”,ACM Transactions on Graphics 22(3),(Proc.of SIGGRAPH 2003),303–312;
John Goutsias and Henk J.A.M.Heijmans,“Nonlinear Multiresolution SignalDecomposition Schemes—Part I:Morphological Pyramids”,IEEE Trans.on ImageProcessing,Vol.9,No.11,November 2000;
John Goutsias and Henk J.A.M.Heijmans,“Nonlinear Multiresolution SignalDecomposition Schemes—Part II:Morphological Wavelets”,IEEE Trans.on ImageProcessing,Vol.9,No.11,November 2000;
Jean Serra,“Image Analysis and Mathematical Morphology”,1982;
A.J.Britten,M.Crotty,H.Kiremidjian,A.Grundy,and E.J.Adam,“The addition ofcomputer simulated noise to investigate radiation dose and image quality in images withspatial correlation of statistical noise:an example application to X-ray CT of the brain”,The British Journal of Radiology,77(2004),323–328;
C.Tomasi and R.Manduchi,“Bilateral filtering for gray and color images”,inProceedings of the 6th International Conference in Computer Vision(ICCV),1998,pp.839–846;
J.Weickert,“Coherence-Enhancing Diffusion Filtering”,International Journal ofComputer Vision,31(2-3),pp.111–127,1999;
A.Buades,B.Coll,and J.-M.Morel,“On image denoising methods”,Centre deMathématiques et de Leurs Applications(CMLA)publication No.2004-15,2004;
P.Coupéet al,“Fast Non Local Means Denoising for 3D MR Images,”9thInternational Conference on MICCAI 2006,R.Larsen,M.Nielsen,J.Sporring(eds.),Lecture Notes in Computer Science,Vol.4191,pp.33-40,Copenhagen,Denmark,Oct.2006;
M.Mahmoudi and G.Sapiro,“Fast Image and Video Denoising via Nonlocal MeansOf Similar Neighborhoods,”IEEE Signal Processing letters 12,839-842(2005);
A.Heiderzadeh and A.N.Avanaki,“An Enhanced Nonlocal-Means Algorithm forImage Denoising,”Proc.IEEE 9th International Symposium on Signal Processing and itsApplications(ISSPA’07),pp.1-4,Sharjah,UAE,2007;
N.Azzabou et al,“Image Denoising Based on Adaptive Dictionary Computation,”Proceedings of IEEE International Conference on Image Processing,2007。
发明内容
本发明的一些具体实施方式的一些方面涉及一种用于处理医学图像的方法与系统,以致当噪声和/或分辨率被降低时,产生带有特定的期望特征的输出图像,包括一个或多个降低的噪声、期望水平和空间分布的噪声、期望的分辨率、和/或保持精细细节和结构。
根据本发明的教导,这里提供了一种用于从三维(3D)高分辨率医学图像产生3D低分辨率图像的方法,所述3D高分辨率医学图像是多个2D图像组成,这些2D图像在轴向方向上具有轴向分辨率,所述多个2D图像是沿着所述轴向方向而获得的,所述方法包括以下步骤:采用可逆带通分解技术分解多个2D图像中的每个图像,以产生多个转换的数据组,每个转换的数据组对应于所述多个2D图像中的一个图像;采用非均一权重向量对转换的数据组内的每一部分中的每个像素进行加权,以产生多个加权的转换的数据组;以非线性方式将对于每个部分的所述多个加权的转换的数据组从每个部分合并成为单独的新的转换的数据组;以及采用可逆带通分解技术的逆向技术从所述单独的新的转换的数据组的每个数据组中产生3D低分辨率图像,该图像在所述轴向方向上具有第一分辨率,所述第一分辨率是小于或等于所述轴向分辨率。
在一个可选的具体实施方式中,所述3D高分辨率医学图像是CT图像。
在另一个可选的具体实施方式中,所述3D高分辨率医学图像是MRI图像。
在另一个可选的具体实施方式中,所述可逆带通分解技术是拉普拉斯金字塔分解技术。
在另一个可选具体实施方式中,所述可逆带通分解技术是小波变换分解技术。
在另一个可选的具体实施方式中,所述多个加权的转换的数据组的合并包括:基于加权的像素的期望值,在所述部分的转换的数据组内从对应的位置选择像素。
在另一个可选的具体实施方式中,所述选择是选择具有最高的加权的绝对值的像素。
在另一个可选的具体实施方式中,如果所选择的像素具有比特定极限更大的加权值,在所述部分的转换的数据组内留在对应位置的所有像素的加权平均值是被用于代替所选择的像素值。
在另一个可选的具体实施方式中,所述特定的极限是硬极限。
在另一个可选的具体实施方式中,所述特定的极限是软极限。
在另一个可选的具体实施方式中,所述3D高分辨率医学图像在第二方向上的至少一个第二分辨率是等于所述3D低分辨率图像在对应方向上的对应分辨率,且所述3D高分辨率医学图像的对比度是等于所述3D低分辨率医学图像的对应的对比度。
除非另作定义,这里所采用的所有技术术语和/或科学术语都具有本发明所涉及领域的技术人员所熟知的普通含义。与这里所描述的方法和材料的类似或等同的方法和材料都可被用于本发明的具体实施方式的实施或测试,下面给出了示例性的方法和/或材料。在发生冲突时,本专利说明书(包括这些定义)都会控制冲突。另外,所有材料、方法和例子都仅是示例性的,而不试图作为对本发明的限制。
本发明的实施方式所述的方法和/或系统的执行,可涉及手动、自动或两者结合地执行或完成所选择的任务。而且,根据本发明所述的方法和/或系统的实际设备和装置的实施例,可通过硬件、软件或固件或它们的结合,采用操作系统来执行几种选择的任务。
例如,用于执行根据本发明所述的实施例的选择任务的硬件可以是以芯片或电路的方式实施。对于软件,本发明所述的实施例的选择任务可采用任意合适的操作系统通过计算机来执行多个软件指令来实施。在本发明的一个示例性实施例中,根据这里所描述的方法和/或系统的示例的一个或多个任务是通过数据处理器来执行,例如,用于执行多个指令的计算平台。可选地,该数据处理器包括用于存储指令和/或数据的易失性存储器,和/或用于存储指令和/或数据的非易失性存储器,例如,磁性硬盘和/或可移动介质。可选地,还提供了网络连接。显示器和/或用户输入设备例如键盘或鼠标也可选地提供。
附图的简要说明
本发明的一些具体实施方式在这里仅通过实施例并结合所附的附图来描述。现在特别详细参考有关附图,需要强调的是,通过实施例的方式来显示特别内容,是用于说明本发明的实施例的示例性讨论的目的。在这方面,本说明书采用附图使本领域技术人员能清楚地了解本发明的实施方式。
在这些附图中:
图1A显示了根据本发明的一个示例性实施方式所述的医学图像的去噪方法的流程图。
图1B图解地显示了根据图1A所示的方法用于获取医学图像和对该图像去噪的一种医学成像系统。
图2显示了范围压缩函数的图,该函数可选地用于采用如图1A所示的方法对医学图像预处理。
图3A图解地显示了一个没有噪声的二维图像,而图3B图解地显示了噪声增加的同一图像,该图像带有选择的像素和邻域,以显示用于降低图像噪声的算法,该算法可被用于图1A所示的方法。
图4A图解地显示了图3B所示的图像和选择的像素和邻域,与类似于根据用于在图像中去噪的一个算法所选择的其它像素在一起,这些像素可被用于图1A所示的方法中。
图4B图解地显示了图3B所示的图像和选择的像素和邻域,与类似于根据用于在图像中去噪的另一个算法所选择的其它像素在一起,这些像素可被用于图1A所示的方法中。
图5是用于在图像中去噪的方法的一个流程图,该流程图可被用于图1A所示的方法中。
图6A是采用相对低的X-射线剂量制作的噪声CT图像;
图6B显示了图6A所示的图像在采用图5所示方法降低噪声后的图像;以及
图6C是一个低噪声CT图像,类似于在图6A中所示的图像,但该图像是采用相对高的X-射线剂量来制作的。
图7A是通过整数因子的下采样切片的图。
图7B是采用分数权重的下采样切片的图。
图8A是一系列图像,它们是从3D高分辨率医学图像的输入切片图像。
图8B是以线性下采样方法产生的3D低分辨率图像。
图8C是以非线性下采样方法产生的3D低分辨率图像。
具体实施方式详述
本发明,在它的一些具体实施方式中,涉及一种用于处理医学图像的方法与系统,以产生降低的噪声以及其它期待特征的图像,尤其是,但不排除,涉及一种处理CT图像的方法,考虑到噪声在图像中的非均一分布和/或采用非线性滤波器来保持精细的图像细节。
本发明的具体实施方式的一个特征是从三维(3D)高分辨率医学图像产生3D低分辨率图像,该3D高分辨率医学图像包括一系列2D图像,其中分辨率是在该系列的方向上降低的,且该分辨率在其他方向上被保持。本发明的一些实施例的一个方面是涉及一种对医学图像去噪的方法,其中,该去噪是在高空间分辨率下进行的,例如,原始图像是在该高空间分辨率下获取的,例如,通过CT机获得的图像。然后,去噪图像被向下转换为较低分辨率,可选地采用非线性下采样算法,该算法保持精细的细节,比线性下采样算法更好。较低的分辨率是指,例如,放射科医生通常检查CT图像的分辨率。在现有技术中,在通常的去噪程序中,去噪是在图像已经被转换为低分辨率之后进行的。
本发明的一些实施例的一个方面是涉及一种复原结构为去噪图像的方法,其中,去噪算法在去除噪声的同时已经从该图像中移除了一些结构。残留的图像,在原始图像与去噪图像之间的差异,可选地被过滤以致去除噪声和增强边缘,被完全或部分地加回到已去噪的图像中。在本发明的一个示例性实施例中,在图像中的每个位置,用于残留图像的相对权重取决于在该残留图像中的结构程度,或者在原始或去噪图像中的该位置的结构程度。例如,在一个低的结构程度的位置,可给予残留图像小的权重或不给权重,可选地,该权重轻微或者根本不增加结构的程度;而在另一个高的结构程度的位置,可给予该残留图像相对大的权重,可选地,给予该残留图像最大权重。可选地,用于残留图像的相对权重也取决于在原始图像或残留图像中的每个位置的本地噪声水平,例如,对具有较高噪声水平的残留图像采用较低权重。
本发明的一些实施例的一个方面是涉及一种产生带有特定量值和分布的噪声的去噪图像的方法。例如,特定量值和分布的噪声可以是由放射科医生在检查CT图像时所期望的量值和分布的噪声,以便使该图像比完全去噪的图像看起来更自然。该特定量值和分布的噪声是通过加入噪声到该图像来获得的,带有空间包络,这样会获得想要的量值和分布的噪声。可选地,该噪声是通过采用空间变化加权参数对原始图像与去噪图像进行平均而获得的。
本发明的一些实施例的一个方面是涉及一种使图像去噪的方法,其中,去噪算法测量与在图像中的位置成函数关系的噪声水平。例如,噪声水平是通过围绕该体素的大窗口查看而在给出的体素处查明的。在该大窗口内,仅考虑在一定范围内的体素的灰度值。该范围可选地取决于正被成像的组织,例如,仅考虑对于该组织的在中间范围内的体素。对于这些体素中的每个体素,在小窗口内查明变化的本地测量的灰度值,例如,标准偏差。测量变化的一组子集,变化的测量结果低于特定的分位数,对在子集的变化的测量值进行平均,以便找到对于大窗口的本地噪声水平。作为与在图像中的位置成函数关系的噪声水平,它是通过采用对于分布在该图像中的多个体素的这个方法来查明的,并可选地内插以覆盖整个图像。一旦该噪声水平是已知的与在图像中的位置成函数关系,可采用这里所描述的其他程序,来根据本地噪声水平最优化或改善程序本身的性能。
本发明的一些实施例的一个方面是涉及一种使图像去噪的方法,其中,与在图像中的位置成函数关系的噪声水平是被用于去噪算法中。例如,该去噪算法是一个“非本地平均”(NLM)类型的算法,其中,围绕给出体素的补码是与围绕在该图像的另一处的对照体素的补码进行比较,并对于这连个补码进行类似测量。基于在给出体素的本地噪声水平、在对照体素的本地噪声水平,或者它们两者,计算由于噪声导致的在类似测量中所期待的变化,并用于使该类似测量标准化。在现有技术的NLM算法中,取决于体素的位置的恒定噪声水平已经被用于使类似测量标准化。对于对照体素,基于标准化的类似测量,查明一个权重,也查明对照体素的灰度值的加权平均。然后,基于对照体素的灰度值的加权平均,改变该给出体素的灰度值。可选地,通过上述的程序查明与在图像中位置成函数关系的噪声水平。
本发明的一些实施例的一个方面是涉及一种锐化图像的方法,其中,变化的本地程度的灰度值,例如本地标准偏差的灰度值,是作为在图像中位置的函数关系来查明的。然后,对该图像应用锐化滤波器,但锐化的程度是取决于变化的本地程度的灰度值。当变化的本地程度的灰度值越大时,采用越小的锐化程度。所得到的图像在锐化边缘可具有相对少的或者没有视觉假象,与线性锐化滤波器相比,后者对于给出的平均锐化程度会产生更多的视觉假象。
总的来说,在图像处理文献中,术语“像素(pixel)”是用于表示二维图像的一个元素,而“体素(voxel)”是用于表示二维图像的一个元素。因为这里所描述的方法总体上既可用于二维图像,也可用于三位图像,这里所采用的术语“像素”和“体素”不应理解为将本说明书限制为二维图像或三维图像的例子。除非另有特别说明,这里所采用的术语“像素”和“体素”应被理解为应用于任意例子的通用术语,它们通常可交替使用。
在医学成像模式中,例如CT或MRI,它是定制的,当显示一个图像时,绘制该图像密度,例如在一个CT图像中以豪森菲尔德单位(HU)表示的密度,以表明亮度水平或从黑到白的灰度值范围,对于感兴趣的特定的密度范围。术语“灰度值”也可用于指图像密度,用于成像模式的图像密度的范围之外,该图像密度是在图像处理的中间步骤中产生的,例如一个负数。此外,这里所使用的术语“灰度值”不仅是指黑白图像的亮度,还是指彩色图像中任意色彩变量的程度,例如,在彩色图像中的红、绿或蓝,或者彩色图像的亮度或饱和度。在诸如CT或MRI图像的医学图像中,通常只显示一种单一密度变量,例如HU密度、T1或T2加权密度,它被绘制为灰度显示的亮度水平,而在本例中,“灰度值”是特别倾向于采用,但这里所述的方法并不限制于它们对医学图像或对黑白图像的适用性。
这里所描述的降低噪声方法可被特别用于医学图像,因为医学图像通常具有相对高的噪声水平,由于在噪声水平和图像采集参数(例如,X-射线剂量或MRI采集时间)之间通常有折中,施加经济或安全惩罚以用于降低噪声。而且,因为医学图像总体上没有在“光照”差别,一个像素的邻域的特征通常是它的真实灰度值的好的指示。这是特别真实的,因为医学图像倾向于具有类似结构,在该图像的不同部分重复出现,有时带有改变的尺度或方向。
举例来说,这里所描述的方法可被应用在图像采集设备或它的工作站(例如,CT机、MRI机、超声成像机)上,应用在图像处理站和/或通过网络连接到远程位置或远程服务器。
在详细解释本发明的至少一个实施例之前,需要明确的是,本发明并不必然限制其申请在以下说明中的详细描述的内容。本发明能采用多种不同方式实施或执行的其他实施例。
去噪方法的概述
现在参考附图,图1A显示了一个流程图1000,显示了根据本发明的示例性实施例所述的用于产生降低噪声的图像的程序。该程序是对于采用CT图像而开发的,除了可用于处理CT图像外,也可被用于处理其他医学图像,或者其他类型的图像。不同类型的医学图像,以及非医学图像,具有不同的特征,这些特征可以使该程序的特定实施更适合于它们。诸如CT和MRI图像等的医学图像不依赖于普通非医学图像处理的“光照”方法,它们的噪声水平和噪声的空间分布可以是一致的,并可从一个图像预期到另一个图像。
在步骤1002,通过一个医学成像设备(例如CT机)获得高分辨率和有噪声的3D原始图像I。在图1B中原理地显示了用于获得这样的图像和对图像去噪的系统500。采用控制器506,例如计算机,或者与成像设备相关联的指定电路,从由成像设备502从一个病人504上获得的原始数据来重构该图像。该控制器也可以执行任意或者全部下面所述的去噪程序。在本发明的一些实施例中,物理上分离的计算机或者电路可执行不同部分的去噪和/或图像重构程序,但它们在这里都是全部用作控制器。可选地,诸如显示器等输出设备508被用于显示去噪图像,而诸如键盘或控制台等输入设备510被用于控制成像程序和/或成像设备。
再参见图1A,在步骤1004,原始图像I被可选地采用非线性滤波器来过滤,导致获得标示为C的预处理图像,以后将用于比较在去噪程序中的补码。当在常规CT图像中最感兴趣的密度范围是相对小时,例如–50HU至+100HU豪森菲尔德单位(HU),这里全范围是–1000HU至+4000HU,可在图像I的非线性过滤之前进行可选范围的压缩程序。根据预期的噪声模式对图像的“预变白”,定义在下面的标题为“预处理的图像”的章节A,也是可选地进行的。
在CT图像中,由于获取机制和光束硬化,噪声总体上在空间不是均匀分布的,通过一种方法使X-射线的非单色光束增加穿透身体的平均能量,因为更柔和的X-射线被优先吸收了。可选地,该去噪程序采用根据在该CT图像的每个点的本地噪声特征的空间依赖参数,这些参数在步骤1006中被评估。可选地,在预处理该图像之前,或者在预处理该图像的同时,评估本地噪声特征。然而,先进行预处理是有优势的,例如,如果本地噪声特征是基于在预处理的图像中的“变白噪声”,而不是基于在原始图像中的“有色噪声”时。这些术语在下面标题为“预处理的图像”的章节A中进行定义。
在步骤1008中,在图像上进行非本地平均(NLM)类型的去噪程序。总的来说,NLM程序是通过找到类似于在给出的体素周围的补码来对在图像中的给出体素进行去噪的。该体素的去噪值是计算为这些补码的中间体素的加权平均值,其中,加权的权重是与在这些补码之间的一些类似度量成比例关系的。除了NLM程序之外,也可采用本领域已知的其他去噪程序,取代NLM去噪程序或者作为NLM去噪程序的补充。
在基于NLM程序的特征中,类似度量是基于从这些补码中抽取的多种特征,其中,这些特征可基于图像C和图像I。该类似度量也可取决于在步骤1006中评估的参数。图像D是在执行了基于NLM程序之后的去噪图像。在步骤1008中采用的NLM程序可以是某个特征类型的NLM程序,例如基于下面标题为“在NLM去噪算法中采用的示例性类型的特征”中描述的特征的NLM程序,以及相关的公开于PCT申请WO2009/081410的内容,或者在例如上面所引用的早期文献中关于去噪图像所描述的类型的NLM程序。
在步骤1010中,已经在去噪程序中从图像中移除的结构被可选地复原到该图像。即使在去噪后,也可能从残留图像R=I–D中复原附加的特征,该图像含有大多数已经移除的噪声,但也包含一些空间上可区别的结构。相干增强扩散(CED)或类似滤波器是可选地用于找到那些结构,这些滤波器对残留图像进行非线性平滑处理。可选地,过滤的残留图像S不是直接加到图像D,但可以乘以参数α,该参数可取决于图像D的本地特征。图像D′=D+αS是在复原后的去噪图像。
在步骤1012中,噪声是可选地加入到图像中。对于习惯于在CT图像中特定水平的噪声的放射科医生来说,去噪图像有时看起来不自然。在图像D′与图像I之间的平均对于放射科医生看起来更自然。这样一个图像的例子是E=βD′+(1–β)I,其中,β是精选的以致所得到的图像E具有一些特定标准偏差的噪声模式,而β可在该图像上变化。可选地,图像D′也可被存储,用于后面的可能的应用,在该图像中,附加的噪声可以是不利的,例如,如果放射科医生想要在更高分辨率缩小一部分图像时。
在步骤1014中,可选地降低了该图像的分辨率。放射科医生通常在比由CT扫描仪所获得的原始图像的分辨率更低的分辨率检查常规的CT图像。分辨率的最终降低是可选地采用非线性降低程序来执行的。该分辨率降低程序是被设计为保持比图像的线性过滤和采样更多的细节。
去噪方法的细节
下面提供了流程图100的每个程序的更多细节。
图像获得
在日常实践中,由放射科医生检查的CT图像通常具有约2.5mm的轴分辨率(切片厚度)。然而,大多数现代CT机的天然轴分辨率是约0.6mm。因此,由CT机获得的原始数据是比由临床医生所看到的数据具有更高的分辨率。通常不在日常实践中利用该高分辨率数据,由于它包含高的噪声水平,且它需要更长的阅读时间。在临床常规中检查的每个低分辨率切片包含从在原始图像中一定数量的切片整合的数据。
这里所述的去噪程序和相关的图像处理方法可选地采用高分辨率数据,这些数据通常是在CT机中获得,并包含更多信息。在完成处理之后,可选地产生较低分辨率图像,供放射科医生检查。
在本发明的示例性实施例中,在临床轴分辨率与原始轴分辨率之间的差异是以以下两种途径之一或同时来充分利用:第一,较丰富的数据组被用作输入到图像处理算法。由于高信息内容,该算法是潜在地可产生更精确的结果。第二,可选地对处理的图像以非线性方式进行下采样,以产生较低分辨率的图像,潜在地保留在原始图像中可用的精细的细节。
预处理的图像
在预处理阶段,预处理的图像C是可选地从原始图像I中计算的。该预处理的图像是可选地用于在去噪程序中的补码的比较。一个基础的去噪算法可写成下式:
x ^ i = Σ j x j · w ( x i , x j ) Σ j w ( x i , x j ) ,
其中,
Figure BDA00002670283000122
是第i个体素的去噪值,而权重w(xi,xj)是基于在围绕图像I的体素xi和xj的补码之间的距离的特征的函数。更多常规去噪算法可选地比较在原始图像I的补码与预处理的图像C的补码:
x ^ i = Σ j x j · w ( { x i , c i } , { x j , c j } ) Σ j w ( { x i , c i } , { x j , c j } ) ,
预处理的图像C是可选地采用下面的步骤来从图像I进行重构:
A.根据噪声模式的有色噪声的预变白
CT噪声可被近似地模式化为附加的有色高斯噪声,也就是,由一些滤波器过滤的白色高斯噪声,该滤波器被表示为噪声—色彩滤波器。优选地,可通过应用维纳滤波器对图像进行预变白处理,该滤波器试图转化噪声着色的操作,因此产生含有白色噪声的图像,该噪声不是空间相关的。这种类型的噪声已经在文献中广泛研究,且可以是比有色噪声更容易去除的。
可选地,通过计算来自均一仿真的CT图像的协方差矩阵来评估噪声—色彩滤波器,类似于在文献A.J.Britten,M.Crotty,H.Kiremidjian,A.Grundy,and E.J.Adam,“The addition of computer simulated noise to investigate radiation dose and imagequality in images with spatial correlation of statistical noise:an example application to X-ray CT of the brain”,The British Journal of Radiology,77(2004),323–328中描述的方法。
B.范围压缩
CT图像的动态范围通常是在约–1000HU(空气)至+4000HU(金属)之间。然而,在0HU至100HU之间的范围是更为重要的考虑范围,相比于例如在1000HU至1500HU之间的范围,因为第一个范围区分了软组织,而第二个范围表示更致密的组织,例如骨。此外,诸如骨等致密组织通常是用更宽的密度窗口来检查的,这导致噪声不那么可视。放射科医生通常检查脑CT图像,其密度窗口约为0HU至+80HU,而肝CT图像的密度窗口为–15HU至+155HU。在图像C采用灰度转换是有优势的,这样扩展了软组织范围,而压缩了–1000HU至–200HU的范围,并在+300HU之上。例如,灰度水平转换可以是原始密度x的下列函数y:
x ′ ← 2 · x - a b - a - 1
 y′←sign(x′)·log(1+|x′|)
y ← y ′ · b - a 2 + a
其中,a和b是范围常数,而y是转换的密度。图2显示了范围压缩函数y(x)的图200。
需要注意的是,其他灰度水平转换函数也可被构建,并也能被使用。例如,合适的灰度水平转换可从伸展CT图像的柱状图或者平衡该CT图像的柱状图中获得,可选地仅对图像中感兴趣的部分,例如,没有周围空气的部分。
C.Robust非线性噪声移除
可选地,将噪声移除滤波器应用到图像C。本发明人已经发现,即使噪声移除滤波器使得图像C显得过度平滑可视,损失了一些分辨率,过滤后的图像C仍是优于未过滤的图像C,以用于比较补码的目的。可选地,图像C仅可用于比较补码的目的,以评价类似度量,并在原始图像I而不是图像C上执行平均操作。
可选地,图像C是采用双边滤波器来过滤,正如文献C.Tomasi and R.Manduchi,“Bilateral filtering for gray and color images”,in Proceedings of the 6thInternational Conference in Computer Vision(ICCV),1998,pp.839–846所描述的。可选地或附加地,图像C可采用相干增强扩散(CED)滤波器来过滤,正如文献J.Weickert,“Coherence-Enhancing Diffusion Filtering”,International Journal of ComputerVision,31(2-3),pp.111–127,1999所描述的,也是在保留边缘的同时使图像平滑。也可采用本领域已知的其他类型的噪声移除滤波器。
本地噪声水平的评估
根据本发明的一个示例性实施例所述的本地平均去噪方案的一个重要参数是σR,它可控制在不同邻域之间的加权,例如由下式给出:
w ( { x i , c i } , { x j , c j } ) ≡ W ij = exp ( - d P 2 ( i , j ) σ P 2 - | | C i - C j | | 2 σ R 2 ) ,
其中,Cx表示在图像C中围绕体素x的图像补码,而dP(i,j)是在体素i和j之间的空间距离。这里||Ci-Cj||表示在围绕图像C中体素i和j的补码之间的差异测量。
在文献中,以前曾有提议,参数σR的值可以在由于噪声的测量的变化中的标准偏差的量值等级上;例如参见上面所引用的A.Buades、B.Coll与J.-M.Morel的文献。加权的选择是根源于下面的启发式观察。类似于在图像中的统计噪声的量值的两种图像补码将接收大的权重,指示它们的类似性。与之相反,在所期待的来自图像中的噪声差异之外的不相同的两种补码将接收低的权重。
传统上,在图像中的噪声是被考虑为同等分布的,也就是,在该图像中任意空间位置具有类似统计学特征。因此,σR是通常每个图像一次确定的,而相同值被用于计算在图像内所有空间位置的权重。
在CT图像中,噪声通常在空间不是均匀的,由于获得机制和线断层重构算法。本发明人已经发现,对于σR穿过整个图像的恒定值不能很好地用于CT图像。这个选择可导致更低噪声图像区域的超平滑,以及具有更高噪声水平的不足去噪的图像区域。替代的,一种空间变化值可选地用于σR,其中,例如,在图像中的每个体素是分配对于σR的合适的值,根据由于噪声的测量的本地标准偏差。
下面总结了用于评估由于噪声的空间依赖性标准偏差的示例性方法。首先,对于在原始噪声图像I中的每个体素,计算本地标准偏差的灰度值,或者本地变化的另一测量的灰度值。这是通过计算小邻域例如3×3×3体素的标准偏差(或变化的其他测量值),与正被考量的体素相关联,例如在它的中心含有正被考量的体素。然后,围绕图像中的每个体素,检查更大的邻域,例如33×33×5体素。对于每个大邻域,在特定密度范围内的一些或所有体素,例如–50HU至+150HU被抽取。最后,所抽取的体素的子样本的本地标准偏差值被平均化以产生由于噪声的本地标准偏差的评估值,对应于例如本地分位数,低于在大窗口中抽取的体素的本地标准偏差的特定值(例如0.3)。
由于采用大的邻域,由于噪声的本地标准偏差的评估不是期望在相邻体素之间变动太大。可选地,仅在所有图像体素的子样本上进行计算,而可选地找到在其他位置的值,如果它们是需要的,采用标准内插技术来插入所计算的值。
可选地或附加地,采用其他方法来评估本地噪声水平,例如,仿真成像,直接从在仿真的均一部分的灰度水平的标准偏差中确定噪声水平。
σR的值被设为特定分数,例如由于噪声的本地标准偏差的评估为0.54。可选地,在去噪算法中可采用其他参数,或在基于由于噪声的本地标准偏差的评估的其他相关图像处理程序中也可采用其他参数。
在本发明的一些实施例中,用于去噪算法或相关程序的σR和/或任意其他参数,具有在图像中任意位置的相同值,但该值取决于由于噪声的本地标准偏差的整个图像的平均值,它是可选地按上述方法计算。例如,当查明本地标准偏差时,仅考虑在特定密度范围内的体素。可选地,当计算平均值时,仅包括低于特定分位数的本地标准偏差。
示例性NLM去噪算法
图3A显示了二维图像100,包括像素阵列,每个像所有数值,该数值映射到在黑与白之间的一个灰度值。在CT图像中灰度值表示所成像的对象的真实密度的方便映射,定制为豪森菲尔德单位(HU)。例如,在脑的CT图像中,该图像通常是可视的以致0HU,这表示水的密度,被映射到黑色,而70HU被映射到白色。
图像100包括亮区102和暗区104,在它们之间有相当尖锐的边界。在图3B中,图像108是将噪声加入图像100得到的图像。在现有技术中,通过将像素的灰度值与邻域像素的灰度值进行平均化,使噪声有些降低,对位置最接近的像素给予最大权重。这在没有精细细节的均匀区域做得很好,例如在图像108中的区域102和104,但在它们之间会导致边界的模糊。另一个现有技术中的去噪方法,采用双边滤波器,是一种非线性滤波器,试图避免由对像素i的灰度值Ii与首要的在灰度值上类似的其他像素j的灰度值Ij的平均化带来的问题。例如,当在位于(xi,yi)的像素i上进行操作时,对于位于(xj,yj)的另一像素j的灰度值的权重Wj是由下式给出:
W j = exp ( - ( x i - x j ) 2 σ P - ( y i - y j ) 2 σ P - ( I i - I j ) 2 σ R ) = exp ( - d P 2 σ P - ( I i - I j ) 2 σ R )
其中,dp是两个像素之间在空间上的欧几里得距离,而|Ii–Ij|可以被考虑为在这两个像素之间的抽象“距离”,测量它们互相类似的程度。对于像素i的新的灰度值是由下式定义:
Figure BDA00002670283000162
其中,N是围绕像素i的检索窗口,而总和是在该检索窗口的对于所有像素j的。
用于去噪的另一种类型的非线性滤波器描述在文献L.Rudin,S.Osher,and E.Fatemi,“Nonlinear total variation based noise removal algorithms,”Physica D60,259-268(1992)。
在非本地平均滤波器中,两个像素的类同取决于这两个像素的邻域的像素对像素的比较。例如,为降低像素i的噪声水平,在图3B中标记为110,邻域Mi,在图3B中标记为112,限定在像素110周围。然后检索其他像素j,带有相同大小和形状的邻域Mj,围绕每个这样的检索像素j,在这些像素之间的邻域112查明平均方差MSE(Mi,Mj),以及每个检索像素j的邻域的相应像素。对于检索像素,在它们的邻域与像素110的邻域之间的平均方差是小的,这些像素被给予最大权重,当对这些检索像素的灰度值进行平均以获得对于像素110的降低噪声灰度水平时。该权重Wj是由下式给出:
W j = exp ( - d P 2 σ P - MSE ( M i , M j ) 2 σ R ) ,
对于像素i的新值是由下式确定:
Figure BDA00002670283000164
图4A显示了图像200,类似于在图3B中的图像108,带有一组像素202,这些像素具有类似于像素110的邻域112的邻域。每个像素202具有类似的邻域,因为像素202都是与在亮区102和暗区104之间的边缘有相同距离,导向几乎相同的方向。
在其他采用非线性滤波器的去噪方法中,两个邻域的类同是基于在邻域的所有像素的平均灰度值,或者在邻域的像素的灰度值的梯度的方向上,正如由文献M.Mahmoudi,and G.Sapiro,“Fast image and video denoising via nonlocal means of similarneighborhoods,”IEEE,Signal Proc.,vol.12,no.12,pp.839-842,Dec.2005所描述的。在另一不同方法,描述在A.Heidarzadeh and A.N.Avanaki,“An Enhanced NonlocalMeans Algorithm for Image Denoising,”9th ISSPA,Feb.2007,这两个邻域的类同取决于这两个邻域的双边映射的平均方差,如采用Canny边缘探测器所确定的,也取决于在这两个邻域的原始图像的平均方差。
图4B显示了图像204,类似于图像108。根据本发明的一个示例性实施例,采用不同标准计算权重Wj,找到更好的一组检索像素206,这些检索像素具有类似于足够接近的像素110的邻域112的邻域。在图2B所示的特定实施例中,标准(将在下面详细描述)不取决于邻域的相对方向,因此所有像素与暗区104的距离相等,像素110将具有非常类同的邻域112的邻域,根据这些标准。具有高权重的放大组的检索像素206,与采用非本地平均方法具有高权重的检索像素202相比,可使得噪声进一步降低,因为有更多像素的灰度值平均化。在本发明的一些实施例中,在两个邻域之间的类同标准可取决于这两个邻域的相对方向,而在本发明的这些或其他实施例中,带有高权重的检索像素的量可以比现有技术的方法更多,而这些检索像素的质量会更好,显然它们提供了像素110的真实灰度值的更好的评估。
图5显示了根据本发明的一个示例性实施例所述的一种降低图像中噪声的方法的流程图300。流程图300所示的方法是图4B所示方法的总结,具有用于在邻域之间的相似性的不同标准。在步骤302中,获得有噪声的图像。噪声降低算法在一次检查一个像素,初始在步骤304中设定像素i等于1。在步骤306中,考量像素i,在步骤308中找到像素i的特征向量F1。特征向量是一个或多个特征的排序的一组值,每个值取决于所考量的像素的灰度值和/或在周围邻域内的其他像素的灰度值。该邻域不需要是邻接的,也不围绕像素,但可以在一侧。被考量的像素的坐标,例如xi和yi(在二维图像),或xi、yi和zi(在三维图像),也可被处理为特征。本领域已知的特征的例子包括:像素i的灰度值,采用上述的双边滤波器,以及在围绕像素i的特定尺寸的邻域内的每个像素的灰度值,采用非本地滤波器装置。如上所述,本领域已知的其他特征包括:在像素i的邻域内的所有像素的平均灰度值,在像素i的邻域内的灰度值梯度的方向,以及在像素i的邻域的二进制边缘图的每个像素的灰度值,正如采用Canny边缘检测器确定的。有宽的变化范围的其他特征可被使用,正如在公开的PCT申请WO2009/081410中的示例性描述,本申请要求该PCT申请的优先权,在第14-17页。
在步骤310起,检查一组检索像素,标记为检索像素j,以便找到具有与像素i类似特征值的像素。检索像素j的灰度值非常近似于像素I,这将最大程度地有助于评估像素i的没有噪声的真实灰度值。在步骤310,指数j被初始设为等于1。在步骤312,考量检索像素j。检索像素可选地包括在该图像中的所有像素,或除了像素i以外的所有像素。另外,检索像素仅包括在该图像内的像素的子集,例如,仅包括在检索窗口内围绕像素i的像素,或仅包括随机选择的一些像素,或者在该检索窗口内规则间距的像素,和/或具有足够接近于像素i的灰度值的像素。可选地,例如在医学图像中,该图像是分割为不同类型的组织,采用任意已知的分割技术,仅或优选地,从如像素i的同类组织的像素中选择检索像素。
在本发明的一些实施例中,取代或附加地采用选自正被去噪的图像的检索像素,检索像素可在数据库中从其他图像中选择。例如,从被预期为类似于这个图像的其他图像中,预先可完成一个可能的检索像素和它们的邻域的字典。例如,如果本图像时一个医学图像,字典包括来自从同一病人的身体的相同部分制备的早期图像的像素,或者来自其他病人的身体的相同部分的像素。
在步骤314,评估对于检索像素j的特征向量F2。特征向量F2是一个或多个特征的排序的一组值,每个值对应于特征向量F1的一个特征值。可选地,以同样方式定义在F1和F2中的对应特征,采用在像素i和像素j的邻域的对应像素的灰度值。在本发明的一些实施例中,在F1和F2中的对应特征值是不同定义的,例如,围绕其中一个像素的邻域可被指向不同的角度,或以不同尺寸缩放,相对于围绕其他像素的邻域,如果需要,在计算特征值时插入灰度值。在任何情况下,在F1和F2中的对应特征都是可选地以类似的足够方式来定义,以致能够比较它们,也能采用它们的值的差来计算在像素i和像素j之间的理论距离测度,可知它们之间的相似度如何,用于降低噪声的目的。
如果检索像素j是从之前存储的检索像素字典中取得,而不是从正被检查的图像中取得,则对于像素j的特征向量F2,或者它的一些分量,也被存储在该字典中,在使用时不需要每次计算。类似地,如果检索像素j之前是作为对于另一像素i的检索像素来使用的,则它的特征向量F2可选地是存储在内存中,不需要再次计算。可选地,对于在图像中的所有像素,特征向量F2是提前评估的,并存储在内存中,因此F2不需要在循环中对于检索像素j和像素i而被评估。以至于F2的特征值是以相同方式定义的,作为F1的对应特征值,特征向量F2或它的一些分量,也可从内存中取回,而不是再次计算,如果检索像素j之前采用为在步骤306中检查的像素i时。
在步骤316中,距离测度d(F1,F2)是可选地计算的,它是反映像素j与像素i的相似度的理论距离,由它们的灰度值和它们邻域的灰度值来定义,也可能由它们的位置来定义。距离测度d取决于构成特征向量F1和F2的每个对应特征值的差异。如果特征向量F1和F2的每个特征值具有由
Figure BDA00002670283000191
Figure BDA00002670283000192
给出的k分量(特征值),则该距离测度可由下式定义:
d ( F 1 , F 1 ) = ( α 1 | f 1 1 - f 1 2 | β + α 2 | f 2 1 - f 2 2 | β + . . . + α k | f k 1 - f k 2 | β ) 1 / β
其中,(α12,…αk)是加权向量,给予在计算距离测度中用于不同特征的权重。参数β通常是同一级的正数,通常设为等于2,这使d(F1,F2),正交分量的欧几里得距离,每个等于在两个像素i和j的特征值之间的加权绝对差。正如下面所描述的,加权向量(α12,…αk)是可选地采用遗传算法找到的,这试图找到一个最佳的加权向量以使降低噪声方法的效率最大化。
对于d(F1,F2)的另一种表达,考虑在不同特征值之间的相关性,例如在一个邻域中的不同像素的灰度值之间,是描述在公开的PCT申请WO2009/081410中,参见该申请的图8的说明。对于d(F1,F2)的表达式可包括交叉项,例如(f1 1–f1 2)(f 21–f2 2),可提供对于在不同邻域之间的相似度的更有用的测度,在不同特征值是相关的情形下。
在步骤318中,对于像素j的权重Wj是可选地从d(F1,F2)计算的,并存储在内存。当像素i和j的邻域彼此最大相似时,也就是,当d是小的时,权重Wj是最大的,而当d是大的时,Wj是小的。例如,Wj=exp(-d2N)。如果特征值仅取决于像素和它的邻域的灰度值,而不取决于该像素的位置,则Wj是是由Wj=exp(-d2N–dp 2p)定义的,这里dp是在像素i和j之间物理距离的测度,例如,欧几里得距离。其中,σN和σp是决定伴随在像素i和j之间理论距离d和空间距离dp的增加而Wj减弱的程度的参数。此外,Wj具有与d与dp的不同的相关性,但仍随着d与dp的值的增大而减弱。可选地,为节省计算时间,或者为增强性能,权重Wj是被设为零,当它小于某些极限时,或当d和/或dp大于某些极限时。
在步骤320,检索像素j被增加1,以检查下一个检索像素。在步骤322,确定所有检索像素是否已经被检查完。如果否,在步骤312考量下一个检索像素。如果是,即所有检索像素已经被检查完,将计算检索像素j的加权平均灰度值(由Wj加权)。
在步骤326,评估像素i的没有噪声的真实灰度值,基于检索像素的灰度值,可选地也基于像素i的原始灰度值,检索像素具有对于评估的真实灰度值的更大影响,如果它们被视为更近似于像素i,基于所具有的类似特征值。例如,特征值的相似度被用于计算理论上的距离测度d(F1,F2),如前所述,每个检索像素j被分配一个权重Wj,该权重基于它的从像素i的距离测度,并从检索像素j的加权平均灰度值(由Wj加权)中找到像素i的评估的真实灰度值。该平均可以是平均数、中值、模式、去除轮廓的平均数,或者任意其他类型的平均。
可选择地,像素i的真实灰度值的评估是以不同方式计算的:从检索像素的灰度值、检索像素的特征向量F2,以及像素i的特征向量F1。例如,检索像素被分为不同级别,表示不同组织类型,基于它们的特征向量F2的簇,只有在同类的检索像素(如像素i)被用于评估像素i的真实灰度值,或者具有比像素i的评估的真实灰度值更大的效果。另外,只有顶部的几个检索像素j具有接近于F1的特征向量F2(通过一些测量),这几个检索像素被用于评估像素i的真实灰度值。可选地,代替采用检索像素的平均灰度值,从基于几个检索像素的灰度的检查表找到像素i的评估的真实灰度值。
校正的灰度值可选地是像素i的原始灰度值与检索像素的加权平均的线性结合。可选地,不明确考虑像素i的原始灰度值,但像素i自身被处理类似另一个检索像素,并包括在加权平均。在本例中,如果F2的特征值是以与F1的对应特征值相同的方式来定义,对于像素i自身的权重Wj将是1,如果F2的特征值是以不同的方式(例如,邻域旋转或缩放)来定义,则对于像素i自身的权重Wj将会小于1。
需要明确的是,这里所指的像素的灰度值不必然是该图像的原始灰度值,还可以是变换的图像或过滤的图像的灰度值,例如,高斯过滤后的图像,其σ等于或小于仅几个像素的宽度。
在步骤328,像素i增加1,在步骤330中,将确定是否仍留有任何像素需要考量。如果有,将回到步骤306考量下一个像素i。如果没有,该流程在步骤332结束,带有降低噪声的图像,采用在步骤326中找到的校正的灰度值,作为输出。
图6A显示了噪声图像400,头部切片的CT图像,以说明在图5中描述的方法。该图像比普通图像有更多噪声,因为它是采用减少的X-射线剂量获得的。图6B显示了降低噪声的图像402,是采用图5所示的方法从图像400获得的,带有一组特征和加权向量,将在下面描述。为了比较,图6C显示了一个低噪声的图像404,采用正常的对这类图像的X-射线剂量获得的。减少噪声的图像402具有比原始图像400显著更低的噪声,可看到更多细节,尤其是在脑中,在不同组织之间有相对低的对比度。与图像400相比,图像402在质量上显得更接近与低噪声图像404。
从残留图像中复原细节
图像D是原始图像I的去噪版本。该去噪图像试图被平滑化,但有时小的结构和器官也被平滑掉。在应用去噪算法后,可选地采用复原程序来进行复原和增强由去噪算法部分地去除或平滑掉的结构。
残留图像R取决于在原始图像与去噪图像之间的差异,例如,R=I–D。残留图像是主要噪声的,所附加的噪声是从原始图像中移除的,但除非有空间连贯性的结构通常可从残留图像中采用非线性边缘保存滤波器来复原。
在本发明的一个示例性实施例中,去噪算法将类似的区域一起平均化,因此当以精细的尺度检查时,残留图像是主要噪声的。该残留图像仍然可保留一些隐藏信息,当以精细的尺度检查时,隐藏信息看起来类似噪声,但当以稍微大的尺度检查时,展示出空间连贯性。为了暴露这些隐藏的结构,可选地采用非线性边缘保留滤波器,试图将边缘定位在3D残留图像的2D表面,或者在2D图像的1D边缘,并平行于这些边缘使残留图像平滑化。合适的滤波器的例子包括:非线性各向异性扩散滤波器,例如,Beltrami流滤波器和相干增强扩散(CED)滤波器。过滤后的残留图像称为S,可暴露空间连贯的结构,同时平滑掉残留噪声。在本发明的一些实施例中,在图像S中增强了边缘或其他结构,但噪声未被平滑掉。可选地,残留图像是未过滤的,在下面所述的程序中,残留图像R被用于取代过滤的残留图像S。无论采用哪种方法,可选地获得复原的去噪图像,该图像是比残留图像R的灰度值更敏感,或者对于过滤的残留图像S,在残留图像具有更大程度的结构的位置的灰度值比在残留图像具有更小程度的结构的位置的灰度值更敏感。
过滤的残留图像S是可选地加回到去噪图像D,以获得带有平滑的细节恢复原处的复原的去噪图像D′。可选地,采用一个适应性参数α,该参数是取决于图像S的本地“结构信息”。有许多测量方法可被用于评估包含连贯结构或相干性的特定区域。例如,可选地采用结构张量的特征值或图像S的海赛矩阵的特征值,以便确定在每个位置的参数α。可选地,在残留图像R的信息被复原到去噪图像D,仅在结构能在残留图像中找到的区域,因此小的噪声也被带回到去噪图像中。
在实践中,放射科医生经常喜欢小量的噪声留在图像中,因此该图像看起来更“现实”。可选地,为获得这样的图像,参数α会保持大于0,即使在残留图像中看起来仅有噪声的区域。
采用D′=D+α·S能增加误导的结构到图像的相对平滑的区域。可选地,图像S的密度范围是包边的。一个替代方法将参数α看作一个“信封”,这样限制了图像S可在每个位置改变图像密度的量。采用这种方法,正切或压缩版本的S图像被加入图像D中,
D′=D+α·正切(S/α)或者
D′=D+正切(α)·S,
其中,正切函数是例如双曲线正切(tanh)函数,或者例如上面给出的范围压缩函数y(x),或者直到极限值并在该值上恒定的的线性函数。参数α的本地值确定了允许的密度范围,以致图像S能在给出的体素处从去噪图像D中移除或增加。
可选地,将锐化滤波器应用到图像D或图像D′。当应用时,这个滤波器是可选地受限于一个类似的信封,例如基于该图像的本地标准偏差。可选地,只应用该锐化滤波器,或者更强烈地应用在具有较低标准偏差的区域。没有这个信封,锐化滤波器可导致视觉假象,特别是近似尖锐的边缘,例如头骨的边缘。采用这个信封,附加的尖锐不会超过相对极限,而视觉假象也可能减少或者完全避免。在本发明的一些实施例中,仅应用滤波器,或者更强烈地应用在具有中间标准偏差的区域,而不应用在高标准偏差的区域,滤波器在该区域可导致视觉假象,而可应用于低标准偏差的区域,滤波器在该区域仅会放大噪声。
CED方案的新变化是从去噪图像D中提取结构信息,而应用滤波器到残留图像R。复原程序是可选地分为三个步骤。首先,计算图像的结构张量,然后计算每个像素的扩散张量,基于结构张量的特征分解,然后进行扩散步骤。可选地,采用迭代程序,对于给出的迭代从(1-λ)·D+λ·S计算结构张量,其中,S是跟着前面的迭代的过滤的残留图像,而当迭代持续进行时,系数λ逐渐接近于1。
增加噪声到去噪图像
去噪图像有时在放射科医生看来显得不自然,因为放射科医生习惯于检查带有特定量的噪声的CT图像,即使当这些图像是在好的成像条件下拍摄的。为了产生看起来更自然的图像,可采用图像D′和I的加权平均图像。例如,找到一个图像E=βD′+(1–β)I,其中β是空间变化加权参数,它是精选的以致所得到的图像E包含类似于在输入图像I的噪声分布的噪声,但比输入图像I的噪声量值低。需要明确的是,如果不采用复原程序,则这里可采用图像D代替D′。
在一些实施方式中,加权参数β是如下计算的。首先,对于图像D′和I,计算本地标准偏差图。这些标准偏差图分别标示为STDD’和STDI。接着,根据以下公式计算β的初步值:
β = ( STD I - σ N ) max ( STD I - STD D ′ , ϵ )
其中,σR是控制所增加的噪声的量的参数,例如σR=4,而ε是一个小的正的常数,例如ε=10-3。其次,β是可选地范围压缩的,例如采用上面所定义的范围压缩函数y(x),例如a=0及b=1。最后,β是可选地平滑化的,例如通过高斯滤波器进行过滤。
分辨率降低
放射科医生通常在较低分辨率检查CT图像,而由CT机获得的原始图像具有较高的分辨率。通常,所检查的图像的切片厚度是约四倍于由CT机所重建的原始图像的切片厚度。
可选地,原始图像I是从CT机获得的高分辨率原始图像,上述的去噪算法和相关程序都可应用于该图像。为了获得在低分辨率的由放射科医生检查的CT图像,可选地采用最终分辨率降低程序。在现有技术中,这样的分辨率降低程序通常是这样进行:首先,以低带通滤波器过滤该图像,例如抗锯齿滤波器,然后,子采样到想要的分辨率。
在本发明的一个示例性实施例中,采用了一种非线性分辨率降低方法,该方法是设计为比现有技术的的去噪方法保持更多细节。例如,在上面所引用的E.H.Adelson、C.H.Anderson、J.R.Bergen、P.J.Burt和J.M.Ogden的文献中所描述的方法通常也可被采用。
最初,Adelson等人的方法被发展为用于将在不同照相机焦距设置所拍摄的两个图像融合为一个单一图像,包含来自两个图像的重要数据特征。类似的方法可被用于融合多个CT切片,可选地可被应用于分离模式的其他图像维度,以获得检查的图像分辨率。在此之后,描述了这样的非线性下采样程序。可选地,当应用到包含安排在一个切片方向上的多个切片的三维图像时,该程序导致在该切片方向上具有降低的分辨率的图像,但该分辨率不会沿着每个切片的方向降低,例如在平面切片的平面上,垂直于切片方向。
首先,在高分辨率的原始CT图像的每个切片上进行带通分解。例如,切片被转换为拉普拉斯金字塔。在前述引用的P.J.Burt与E.H.Adelson的文献中描述了拉普拉斯金字塔。然后,对于每组ns连贯切片,例如每组4个切片,单个金字塔结构,对应于单个较厚的切片,是可选地以下述方式创设的。对于金字塔的至少一些水平面,采用非线性的程序来合并ns切片以形成对于较厚的切片的金字塔的水平面,虽然这不需要完成金字塔的所有水平面。例如,对于金字塔的顶层(包含低频率内容的图像),较好的切片是可选地通过一个线性程序对ns拉普拉斯金字塔的所有顶层图像进行平均化来创建的。但对于该金字塔的其他水平面,采用非线性程序来形成较厚的切片,至少在一些例子中,将更多权重给予较薄的切片,图像对于该切片局部具有更大振幅。例如,该金字塔对于较厚的切片的其他水平面是通过从ns金字塔的一层的每个定位分别取值来形成的。所选择的值是可选地在所有ns金字塔中具有最高绝对值的一个值,除非该值高于特定的极限,例如10HU,在这个例子中,穿过所有ns金字塔取平均。最后,较厚的CT切片是通过重构已形成的金字塔结构来形成的。所有较厚的切片在一起形成较低分辨率的图像,以非线性程序试图在沿着每个切片的方向上保持精细的细节,该图像垂直于切片方向。
在本发明的一些实施例中,采用另一类型的带通分解,取代或者附加地采用拉普拉斯金字塔作为带通分解。例如,可采用小波变换,正如在上面所引用的LiHui、B.S.Manjunath和S.K.Mitra的文献中所描述。附加地或可选地,也可使用任意融合图像的方法,如上面所引用的Moustafa和Rehan的文献中所描述的方法。
可选地,由这种非线性方案创设的低分辨率图像是被合并为加权平均的线性子采样的低分辨率图像。
用于NLM去噪算法中的示例性类型的特征
几种类型的特征可被用于特征向量F1和F2
在本发明的一些实施例中,计算一个或多个特征值的步骤包括:找到在邻域内像素的灰度值的分布特征。可选地,该特征值是灰度值的分布的矩,或者一个或多个矩的函数,其中,分布的第一矩是平均数,第二矩是标准背离,第三矩是偏斜度,等等。分布的第k矩,k>1,可被定义为
Figure BDA00002670283000251
其中In是在该邻域内第n像素的灰度值,总和是在该邻域内N个像素的和,而M1是第一矩,也就是灰度值的平均数。另外,特征值是(或者取决于)分布的次序统计,其灰度值对应于该分布的给定百分数。例如,特征值是中间灰度值,它是在50%百分数处的灰度值。另外,采用不同百分数的灰度值,例如25%、37.5%、62.5%或75%。可选地,采用中间百分数,例如在25%至75%中间,它具有潜在优势:特征值将是作为一个整体的邻域的特征,而不仅是在该邻域的几个轮廓像素。可选地,如果从字典中选择检索像素,该字典包括其他图像的检索像素,带有不同的标准化的灰度值,则这两个图像的灰度值是标准化的,因此它们可做有意义的比较,例如,在基于次序统计的特征的比较。
仅取决于在邻域内像素的灰度值的分布特征的一个特征,尤其是如果该邻域是正方形或者是相当各向同性的形状,具有潜在的优势:该特征值对于在图像中的结构的方向是相对不敏感的。例如,采用在图4B中图像204这样的一个特征,很可能产生一组类似像素206的像素,这些像素具有接近于像素110的特征值,因为该特征值将主要取决于该像素与暗区104的距离,而不取决于在暗区与亮区之间的边缘的局部方向。在另一方面,如果已知特定的一部分图像具有面向某一特定方向边缘或纹理,例如来自身体组织的分割图,则采用对结构方向敏感的特征是有优势的。
可选地,从在邻域内像素的原始灰度值的分布中未找到该特征值,但从该图像已被平滑化或其他方式处理之后的灰度值分布中可找到该特征值。在评估该特征值之前平滑处理该图像,具有潜在的优势:该特征值可更多地取决于在邻域的图像的结构特征,并对在邻域的噪声更低敏感。可选地,对于这里所描述的任意类型的特征,所述平滑或其他图像处理时在评估该特征值之前完成的,而不仅用于取决于灰度值分布的那些特征。例如,可以通过高斯滤波器、二进制滤波器或者全变差滤波器(如前面引用的Rudin等人的文献)来完成平滑处理。可选地,该平滑以某种方式完成,在最大尺寸规模的用于该特征的邻域内不会平滑掉大多数结构,或者甚至用在最小尺寸规模的邻域内。例如,如果以宽度参数σ使用高斯滤波器,则σ可以是比邻域的最大尺寸或者邻域的最小尺寸更小,或者至少不太大。另外,以某种方式完成平滑处理,该方式有效地平滑去除邻域的所有空间结构,而特征值是不在该邻域内的结构的测量,但特征值是围绕该邻域的更大尺度的结构的测量,或者是围绕该邻域的平均梯度的测量。
可选地,在找到灰度值的分布之前,该图像是以不同的方式来处理的。例如,应用导数算子到该图像,代替每个像素的灰度值,通过在特定方向与该图像的导数成正例的值,或者通过与该图像的梯度量值成比例的值。如果完成,则例如在邻域内像素的平均分布值将是在该邻域内平均梯度的测量。可选地,该图像是在找到该梯度之前被平滑的,足够平滑以致在该邻域的大多数像素具有几乎相同的梯度,使该特征值对噪声低敏感。
在本发明的一些实施例中,计算特征值的步骤包括:将变换或过滤应用到该图像,至少在一个邻域内,优先选择在该邻域的最大尺寸与几个像素之间的中间尺度范围内的结构。另外,该变换或过滤优先选择在对于其他方向的一些方向上的结构。以这种方式定义的特征可以是在该图像中有利于挑选期望具有在特定范围(例如,血管)内的尺寸和/或方向的结构,而忽视由于噪声导致的精细的密度变化。
本发明的这些实施例可采用任意大变化范围的滤波器或变换器(线性或非线性)所应用的特征,这些特征已经被应用于诸如计算机手写识别或图像的自动分类对象,因此它们能不依赖于文字描述而被检索到,但这些技术尚未用于降低图像噪声。
这些特征会依赖于在邻域内图像对于小波滤波器的响应,例如Meyer滤波器或Gabor滤波器、拉普拉斯金字塔与高斯金字塔,或者现有技术中的任意其他线性滤波器。这些滤波器可以对在邻域内的具有特定方向和/或特定尺度的结构非常敏感。可选地,该滤波器仅应用到该邻域。另外,滤波器可应用于比该邻域更大的区域,例如,在该图像已被过滤后,应用在该邻域内一个或多个像素的灰度值上。这些选项和另外的方式也可应用到这里所描述的任意其他类型的参与对图像像素进行过滤或变换的特征。
附加地或可选地,该特征值取决于图像对高斯滤波器或其他平滑滤波器在不同尺寸参数σ1与σ2之间的响应的差异。在这样的两个滤波器之间的差异倾向于选择在σ1与σ2之间的中间尺度内结构,但不取决于这些结构的方向,如果该滤波器是各向同性的。以此方式定义的特征是特别有用的,如果该图像具有在许多不同方向上类似结构的话。对于其他图像,这些图像是已知的具有倾向于导向特定方向或者窄范围的方向的结构,采用其他特征测量值可以是有利的,这些特征测量值取决于这些结构的导向。
在本发明的一些实施例中,特征可取决于在邻域内图像对非线性变换(例如形态多尺度变换或形态算子)的响应。例如,该特征值取决于正被检查的像素的灰度值,或者在邻域的特定像素,在采用特定尺度参数将非线性多尺度变换应用到图像之后。可选地,该特征值取决于像素的灰度值,采用两个或更多不同尺度参数,例如,对于两个不同尺度参数的像素的灰度值的差异。形态多尺度变换的例子包括形态小波和形态金字塔,例如描述在以下文献:E.H.Adelson,C.H.Anderson,J.R.Bergen,P.J.Burt,and J.M.Ogden,″Pyramid Methods in Image Processing,″RCAEngineer29,no.6,Nov.-Dec.1984,pp.33-41,或者″Nonlinear Multiresolution SignalDecomposition Schemes—Part I:Morphological Pyramids,″John Goutsias,and Henk J.A.M.Heijmans.IEEE Transactions On Image Processing,vol.9,No.11,November 2000。
附加地或可选地,该特征可取决于在应用形态算子后的像素的灰度值。将形态算子应用到图像以增强或抽取特定结构。形态算子的一个例子是顶帽变换,在输入图像与它的由结构元素打开的形态之间的差。这样一个算子将揭示在暗背景上的光亮细节,带有控制所检测特征的尺寸的结构原件的尺寸。可定义一个类似的算子,它从白色背景中抽取暗结构。
关于形状匹配与图像变形的文献包括宽范围的技术,可被用于使图像的形状特征化,以及任意这些方法可被用于定义邻域的特征,在应用如前所述的形态变换或形态算子之前或之后。例子包括:Earth Mover’s距离,介绍于Rubner等人的文献;用于图像变形的Kantorovich-Wasserstein度量,例如Haker等人的文献的描述;由Osada等人定义的形状签名;以及由Belongie等人定义的形状匹配的度量;所有参考文件都在上面引用过。
在本发明的一些实施例中,对于正被检查的像素i与检索像素j的相应特征值是以由几何变换所改变的邻域来计算的,从一个邻域到另一个邻域。例如,两个邻域是处于不同的相对方向、尺度或两者。另外,其中一个邻域可以使相对于另一个邻域的镜像。例如,如果用于找到像素i的特征值的算法使用了在该邻域的像素的灰度值,该像素与像素i在+x方向有一特定距离,则像素j的特征值是采用替代的像素来计算的,该像素与像素j在+y方向(90度旋转)有相同距离,或者在–x方向(反射)有相同距离,或者在+x方向(尺度改变)有两倍距离,或者在+y方向(旋转加尺度改变)有两倍距离,等等。可选地,旋转的和/或缩放的邻域的灰度值是在计算该特征值之前插入的,尤其是如果旋转角度不是90度的整数倍(在像素的笛卡尔网格的情形下)或者缩放因子不是整数。可选地,这些像素是安排在一个三角形或者六角形网格内,或者更复杂的模格内。
采用以这种方式定义的特征将是尤为有用的:如果无需旋转和/或缩放和/或反射而定义相同的特征,以及用于各种不同的旋转角度和/或缩放因子,相同的权重给出作为结果的特征。这可导致一个距离测度,它取决于在图像中的结构的方向和/或尺度,至少对于一些方向和/或尺度的范围。这将是有优势的:如果该图像包括带有不同方向或尺度的类似结构,或者这些图像是互相的镜像图像。
非线性分辨率降低
在日常实践中,通过放射科医生检查的CT图像通常具有在轴方向上约2.5mm-5mm的轴分辨率,切片是沿着该轴方向而获得的。在本说明书的内容中,切片也被称为2D图像,3D高分辨率医学图像是由这些2D图像组成的。换言之,3D高分辨率图像可由多个2D图像来表示。现代CT机的天然轴分辨率是约0.6mm-1.5mm。因此,对于一个三维(3D)高分辨率医学图像,由该CT机获得的数据是远高于由临床医生观察到的3D低分辨率图像。该高分辨率数据通常不在日常实践中被利用,因为它包含高的噪声水平,而且它需要较长的阅读时间。在外科日常观察的每个低分辨率(厚的)切片包括从原始图像中许多薄切片整合得到的数据。
去噪方法与相关联的图像处理方法可选地采用高分辨率数据,该数据通常是从CT机获得的,该数据包含比由外科医生常规观察到的低分辨率数据更多的信息。在去噪程序完成后,产生较低分辨率图像,也就是,对应于比天然CT图像更厚的切片,可由放射科医生来观察。
在临床轴分辨率与原始轴分辨率之间的差异是表现在以下两个方面:第一,更丰富的数据组被用作输入到图像处理算法。由于在薄切片数据中具有高信息内容,该算法有潜在可能产生更多精确的结果。第二,对已处理的图像的下采样是可选地以非线性方式进行的,以产生低分辨率图像,同时保持在原始图像中存在的精细的细节。
A.整数和分数的线性下采样
假设我们希望通过一个整数参数n减少在CT研究中的切片的数量,通过采用简单平均方式采集一组n个切片,从中构造一个切片。也就是:
J ( x , y , z ) = Σ k = 0 n - 1 I ( x , y , z · n + k ) ,
其中,I是CT数据源,而J是下采样的CT数据。
这个方法是明确定义只用于整数数量的切片减少为一个(通过一个整数参数来下采样),并图解地表示在图7中,该图是通过一个整数参数下采样的切片的图。
为简单化,权重显示在图7中,实践中,在标准化之前,这些权重总是标准化的,以致它们的总和等于1。在本实施例中,原始图像的每个切片是乘以1/4。
在实践中,由于矩形窗的频率区域属性,即使对于整数下采样,也很少采用矩形窗。
让我们考虑这样的情形,我们需要通过切片的分数来降低切片的数量。我们能采用分数权重将这些切片进行分配,正如图7B所示,该图是采用分数权重的下采样的切片的图。
所采用的掩膜实际上是滤波器,在本例中,该滤波器是长度为31/3像素(切片)的1D矩形窗口函数,中心在像素11/6、41/2、7 5/6等。也就是,11/6+31/3·z,其中z∈Z+。在图7B中的权重是假设采用线性插补来进行采样。
如上所述,矩形窗口函数具有许多缺陷。其他窗口函数,包括汉恩、海明和高斯窗口函数普遍用于图像处理实践中。人们可通过以下方式构造一个特殊的滤波器Fz:通过使连续窗口函数离散在特殊取样的网格上,然后收集相应的加权平均——
J(x,y,z)=∑kI(x,y,k)·Fz(k)。
然而,在图像处理实践中用于图像下采样的普通方法是:首先将线性滤波器应用于数据,然后对该数据在较低分辨率的网格点上进行点采样。第一个阶段通常被称为“预过滤”数据。该预过滤器是一个低带通滤波器,它被设计为粗采样网格的表示能力之上的抑制频率内容,以便避免在再采样的图像中的锯齿假象。
让我们假设:我们有n个切片,我们想要下采样为m<n个切片。根据采样率r=n/m,合适的预滤波器是这样构造以致在r之上的频率是抑制的。
设F为预滤波器,低分辨率图像是通过沿着3D图像的z维应用该预滤波器来计算的,而对该3D图像的采样(采样一些插补),在z·r区间,其中z∈Z。
J ( x , y , z ) = [ I * F ] ( x , y , z &CenterDot; r + r 2 )
这两种方法都提供了相同的结果——低分辨率图像的每个切片是高分辨率切片的线性组合,且是采用特殊过滤器Fz来构造的。在后面的例子中,它是离散过滤器F的两个移位实例的组合。
B.非线性下采样
首先,高分辨率粗CT图像的每个切片Ik是被分解为拉普拉斯金字塔的带通成分:
I k &RightArrow; L k 0 , L k 1 , L k 2 , . . .
其中,金字塔
Figure BDA00002670283000303
的顶层是Ik的低带通过滤版本,而
Figure BDA00002670283000304
是相同图像的带通过滤版本。拉普拉斯金字塔也被指示作为已转换的数据组。
下采样的图像是采用下面所述的拉普拉斯金字塔分解技术来合成的。对于低分辨率切片Jz,线性滤波器Fz(k)是根据相应的线性下采样方法和它的线性滤波器来构造的。
然后,Jz的拉普拉斯金字塔已构造。该金字塔的顶层是根据非均一(也称为非矩形)权重向量Fz(k)采用相应的权重组通过平均所有的顶层图像
Figure BDA00002670283000305
来合成的。也就是,
Figure BDA00002670283000306
通过采用非均一加权向量对转换的数据组的每一部分的每个像素加权所产生的数据是被称为加权的转换的数据组。
所述金字塔的其他层是通过在每个像素位置从其中一个金字塔取值来形成的。所选择的值是具有在所有金字塔之中的最高权重的绝对值的一个值,其中,每个金字塔是根据Fz(k)来加权的。可选地,如果所选择的值是比特定极限T更高,通过所有金字塔取加权平均,以便避免在该图像的粗边界附近出现夸张的振荡。可选地,可采用软极限代替硬极限。
km l ( x , y ) = arg max k ( abs ( L k l ( x , y ) &CenterDot; F z ( k ) ) )
J z l ( x , y ) = L km l ( x , y ) l ( x , y ) , if L km l ( x , y ) l ( x , y ) < T &Sigma; k L k l ( x , y ) &CenterDot; F z ( k ) otherwise
最后,低分辨率CT切片Jz是通过重构金字塔结构来合成的。
J z 0 , J z 1 , J z 2 , . . . &RightArrow; J z
线性和非线性下采样的例子是显示在图8A至图8C中。图8A是一系列图像,它们是从一个3D高分辨率医学图像中的输入切片。图8B是以线性下采样方法产生的3D低分辨率图像。图8C是以非线性下采样方法产生的3D低分辨率图像。
如果Fz(k)是一个矩形窗口函数,可获得n个切片的多个(整个数量)的智能融合为一个图像,如上所述。
C.去噪方案
为了获得通常可由放射科医生来观察的低分辨率图像,在对由CT机获得的高分辨率图像应用去噪算法之后采取分辨率降低方法。在现有技术中,这样的分辨率降低方法通常是这样进行:首先通过低带通滤波器过滤该图像;然后对图像进行子采样为想要的分辨率。
这里所采用的术语“噪声降低”和“降低噪声”是与“去噪”同义的。这些术语不是指噪声必然完全消失,仅是指噪声在量值上的降低。
上面所述的任意数量可选地是大于或小于10%、20%、50%,或者是2、3、5或10的因数,或者更大、更小或者中间数量。以豪森菲尔德单位表示的任意密度可以是大于或小于10、20、30、50、70、100、150或200HU,或者更大、更小或者中间数量。
这里所采用的术语“平均”可以是指:平均值、中位数、代表值或者加权平均。
这里所采用的术语“约”是指±10%。
术语“包括”、“包含”、“具有”以及它们的连接词是指“包括但不限于”。这个词组包含术语“包括有”和“主要含有”。
词组“主要含有”意指其成分或方法可包括附加成分和/或步骤,但这些附加成分和/或步骤不会从本质上改变所要求的成分或方法的基础特征和创新特征。
这里所用的“一个”和“所述的”包括复数,除非另有明确说明。例如“一个化合物”或“至少一个化合物”可包括复数个化合物,包括它们的混合物。
这里所用的词“示例的”是指“用作一个例子或例证”。任何描述为“示范的”实施例并不是必须作为比其它实施例更好或更有优势的例子,和/或排除从其它实施例中引入的特征。
这里所用的词“可选地”是指“以一些实施例来提供,不以其它实施例来提供”。本发明的任何特定实施例可包括多个“可选的”特征,除非这些特征发生冲突。
通过本申请,本发明的不同实施例可在一个范围内展示。需要明确的是,在该范围内的描述仅是为了方便和简短,而不应视为对本发明的范围固定的限制。相应地,描述的范围应当被视为具有特别揭示了所有可能的子范围以及在该范围内的个别数值。例如,例如从1至6的范围的描述应当被视为也揭示了例如1至3、1至4、1至5、2至4、2至6、3至6等子范围,以及在该范围内的单个数字,例如,1、2、3、4、5和6。不论范围有多宽,都可这样来应用。
无论一个数字范围在这里如何指示,它是指包括在该指定范围内任意引用的数字(部分或整体)。这里所用的词组“归类”或在第一指示数字与第二指示数值“之间归类”以及将第一指示数字“归类”到,或从第一指示数字“归类”到第二指示数字是可交换使用的,它们的意思是指包括第一和第二指示数字以及在两个数字之间的部分或全部数字。
需要明确的是,为表述清晰,本发明的特定特征是描述在分开的实施例中的,但也可在单个实施例中结合提供。相反,为表述清晰,本发明的不同特征是描述在单个实施例中的,但也可分开提供或者在任意合适的子组合中提供,或者在本发明的任意其它描述的实施例中以合适的方式提供。描述在不同实施例中的特定特征不应视为那些实施例的必要特征,除非该实施例缺乏那些特征元素就无法操作。
虽然本发明已经结合特定实施例进行了描述,对于本领域技术人员显然还有许多变动、修饰和改变。相应地,所有这些变动、修饰和改变都将落入本发明的精神与所附的权利要求的范围之内。
在本说明书中提及的所有出版物、专利和专利申请都在这里整体引入作为参考,即使每个个别的出版物、专利或专利申请时特殊的和单独指示为引入作为参考。另外,在本申请中的任意参考文献的引用或鉴定都不应认为该参考文献是本发明的现有技术。以至于所用的章节标题,它们也不应视为必要的限制。

Claims (11)

1.一种用于从三维(3D)高分辨率医学图像产生3D低分辨率图像的方法,其特征在于:所述3D高分辨率医学图像是多个2D图像组成,这些2D图像在轴向方向上具有轴向分辨率,所述多个2D图像是沿着所述轴向方向而获得的,所述方法包括以下步骤:
a.采用可逆带通分解技术分解多个2D图像中的每个图像,以产生多个转换的数据组,每个转换的数据组对应于所述多个2D图像中的一个图像;
b.采用非均一权重向量对转换的数据组内的每一部分中的每个像素进行加权,以产生多个加权的转换的数据组;
c.以非线性方式将对于每个部分的所述多个加权的转换的数据组从每个部分合并成为单独的新的转换的数据组;以及
d.采用可逆带通分解技术的逆向技术从所述单独的新的转换的数据组的每个数据组中产生3D低分辨率图像,该图像在所述轴向方向上具有第一分辨率,所述第一分辨率是小于或等于所述轴向分辨率。
2.根据权利要求1所述的方法,其特征在于:所述3D高分辨率医学图像是CT图像。
3.根据权利要求1所述的方法,其特征在于:所述3D高分辨率医学图像是MRI图像。
4.根据权利要求1所述的方法,其特征在于:所述可逆带通分解技术是拉普拉斯金字塔分解技术。
5.根据权利要求1所述的方法,其特征在于:所述可逆带通分解技术是小波变换分解技术。
6.根据权利要求1所述的方法,其特征在于:所述多个加权的转换的数据组的合并包括:基于加权的像素的期望值,在所述部分的转换的数据组内从对应的位置选择像素。
7.根据权利要求6所述的方法,其特征在于:所述选择是选择具有最高的加权的绝对值的像素。
8.根据权利要求6所述的方法,其特征在于:如果所选择的像素具有比特定极限更大的加权值,在所述部分的转换的数据组内留在对应位置的所有像素的加权平均值是被用于代替所选择的像素值。
9.根据权利要求8所述的方法,其特征在于:所述特定的极限是硬极限。
10.根据权利要求8所述的方法,其特征在于:所述特定的极限是软极限。
11.根据权利要求1所述的方法,其特征在于:所述3D高分辨率医学图像在第二方向上的至少一个第二分辨率是等于所述3D低分辨率图像在对应方向上的对应分辨率,且所述3D高分辨率医学图像的对比度是等于所述3D低分辨率医学图像的对应的对比度。
CN2011800322137A 2010-06-30 2011-06-30 医学图像的非线性分辨率降低方法 Pending CN103069432A (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US35984510P 2010-06-30 2010-06-30
US61/359,845 2010-06-30
PCT/IB2011/052879 WO2012001648A2 (en) 2010-06-30 2011-06-30 Non-linear resolution reduction for medical imagery

Publications (1)

Publication Number Publication Date
CN103069432A true CN103069432A (zh) 2013-04-24

Family

ID=45402491

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2011800322137A Pending CN103069432A (zh) 2010-06-30 2011-06-30 医学图像的非线性分辨率降低方法

Country Status (5)

Country Link
US (1) US20130202177A1 (zh)
EP (1) EP2588374A2 (zh)
JP (1) JP2014505491A (zh)
CN (1) CN103069432A (zh)
WO (1) WO2012001648A2 (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105518477A (zh) * 2013-09-05 2016-04-20 皇家飞利浦有限公司 使用空间自适应正则化以用于图像重建的mri
CN106464798A (zh) * 2014-04-11 2017-02-22 富士胶片株式会社 图像处理装置、摄像装置、图像处理方法以及程序
CN109044356A (zh) * 2018-09-11 2018-12-21 即智数字科技(苏州)有限公司 一种共享式mri三维医疗影像平台
CN109690620A (zh) * 2016-09-12 2019-04-26 松下知识产权经营株式会社 三维模型生成装置以及三维模型生成方法
CN109784191A (zh) * 2018-12-20 2019-05-21 华南理工大学 一种基于商图像的多任务人脸光照编辑方法

Families Citing this family (43)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20120017521A (ko) * 2010-08-19 2012-02-29 삼성전자주식회사 영상 고정밀화를 위한 고속 처리 필터링 장치 및 방법
KR101669819B1 (ko) * 2010-11-11 2016-10-27 삼성전자주식회사 깊이 영상의 고정밀 복원을 위한 필터링 장치 및 방법
US8855394B2 (en) * 2011-07-01 2014-10-07 Carestream Health, Inc. Methods and apparatus for texture based filter fusion for CBCT system and cone-beam image reconstruction
WO2013027138A1 (en) * 2011-08-19 2013-02-28 Koninklijke Philips Electronics N.V. Frequency dependent combination of x-ray images of different modalities
CA2752370C (en) 2011-09-16 2022-07-12 Mcgill University Segmentation of structures for state determination
US9740912B2 (en) * 2011-12-20 2017-08-22 Definiens Ag Evaluation of co-registered images of differently stained tissue slices
EP2631871B1 (en) * 2012-02-27 2015-07-01 ST-Ericsson SA Virtual image generation
JP2013176468A (ja) * 2012-02-28 2013-09-09 Canon Inc 情報処理装置、情報処理方法
CN102789633B (zh) * 2012-07-02 2015-09-02 河海大学常州校区 基于k-svd和局部线性嵌套的图像降噪系统和方法
US9639915B2 (en) * 2012-08-08 2017-05-02 Samsung Electronics Co., Ltd. Image processing method and apparatus
WO2014097065A1 (en) * 2012-12-21 2014-06-26 Koninklijke Philips N.V. Image processing apparatus and method for filtering an image
US9256965B2 (en) * 2013-01-30 2016-02-09 Impac Medical Systems, Inc. Method and apparatus for generating a derived image using images of different types
JP6283875B2 (ja) * 2013-09-05 2018-02-28 キヤノンメディカルシステムズ株式会社 医用画像処理装置、x線診断装置およびx線コンピュータ断層撮影装置
CN104517263B (zh) * 2013-09-30 2019-06-14 Ge医疗系统环球技术有限公司 减少计算机断层扫描图像重构中伪像的方法和装置
US9189834B2 (en) * 2013-11-14 2015-11-17 Adobe Systems Incorporated Adaptive denoising with internal and external patches
US9286540B2 (en) 2013-11-20 2016-03-15 Adobe Systems Incorporated Fast dense patch search and quantization
WO2015083050A1 (en) 2013-12-04 2015-06-11 Koninklijke Philips N.V. Image data processing
US9183648B2 (en) * 2014-03-04 2015-11-10 Ivan Bajic Method and system for high-resolution transforms of frequency-space and image/audio/video-space data
US9767540B2 (en) 2014-05-16 2017-09-19 Adobe Systems Incorporated Patch partitions and image processing
CN105335947B (zh) * 2014-05-26 2019-03-01 富士通株式会社 图像去噪方法和图像去噪装置
FR3026211B1 (fr) 2014-09-19 2017-12-08 Univ Aix Marseille Procede d'identification de l'anisotropie de la texture d'une image numerique
US10542961B2 (en) 2015-06-15 2020-01-28 The Research Foundation For The State University Of New York System and method for infrasonic cardiac monitoring
JP6402282B1 (ja) 2015-09-16 2018-10-10 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. オブジェクト用のx線画像装置
US9799126B2 (en) * 2015-10-02 2017-10-24 Toshiba Medical Systems Corporation Apparatus and method for robust non-local means filtering of tomographic images
US10559073B2 (en) * 2016-03-23 2020-02-11 Intel Corporation Motion adaptive stream processing for temporal noise reduction
US11222404B2 (en) * 2016-03-25 2022-01-11 Koninklijke Philips N.V. Image reconstruction
EP3244368A1 (en) * 2016-05-13 2017-11-15 Stichting Katholieke Universiteit Noise reduction in image data
JP6774813B2 (ja) * 2016-08-12 2020-10-28 日本電子株式会社 画像処理装置、画像処理方法、および分析装置
WO2018231757A1 (en) * 2017-06-12 2018-12-20 The Research Foundation For The State University Of New York System, method and computer-accessible medium for ultralow dose computed tomography image reconstruction
US11317966B2 (en) 2017-07-19 2022-05-03 Biosense Webster (Israel) Ltd. Impedance-based position tracking performance using scattered interpolant
CA3018334A1 (en) * 2017-09-21 2019-03-21 Royal Bank Of Canada Device and method for assessing quality of visualizations of multidimensional data
CN109035137B (zh) * 2018-07-27 2022-11-25 重庆邮电大学 一种基于最优传输理论的多模态医学图像融合方法
GB201817238D0 (en) * 2018-10-23 2018-12-05 Univ Sheffield Medical image processing
CN111325854B (zh) * 2018-12-17 2023-10-24 三菱重工业株式会社 形状模型修正装置及形状模型修正方法以及存储介质
DE102018222595A1 (de) * 2018-12-20 2020-06-25 Siemens Healthcare Gmbh Verfahren zur Bildbearbeitung eines Bilddatensatzes eines Patienten, medizinische Bildgebungseinrichtung, Computerprogramm und elektronisch lesbarer Datenträger
JP7302988B2 (ja) * 2019-03-07 2023-07-04 富士フイルムヘルスケア株式会社 医用撮像装置、医用画像処理装置、及び、医用画像処理プログラム
US11288775B2 (en) * 2019-11-27 2022-03-29 GE Precision Healthcare LLC Methods and systems for parametric noise modulation in x-ray imaging
CN112991165B (zh) * 2019-12-13 2023-07-14 深圳市中兴微电子技术有限公司 一种图像的处理方法及装置
CN111127300B (zh) * 2020-03-26 2020-07-07 眸芯科技(上海)有限公司 基于图像金字塔分解的去噪方法、装置及系统
CN112241759B (zh) * 2020-07-22 2023-06-09 西安交通大学 一种用于多尺度几何分析的金字塔分解方法及系统
CN113077384B (zh) * 2021-03-10 2022-04-29 中山大学 一种数据空间分辨率提高方法、装置、介质及终端设备
CN113129235A (zh) * 2021-04-22 2021-07-16 深圳市深图医学影像设备有限公司 一种医学图像噪声抑制算法
CN115908157A (zh) * 2021-09-30 2023-04-04 想象技术有限公司 渲染3d场景的图像

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07299053A (ja) * 1994-04-29 1995-11-14 Arch Dev Corp コンピュータ診断支援方法
US7206101B2 (en) * 2001-11-21 2007-04-17 Ge Medical Systems Global Technology Company, Llc Computationally efficient noise reduction filter
US8045770B2 (en) * 2003-03-24 2011-10-25 Cornell Research Foundation, Inc. System and method for three-dimensional image rendering and analysis
US7542622B1 (en) * 2003-06-02 2009-06-02 The Trustees Of Columbia University In The City Of New York Spatio-temporal treatment of noisy images using brushlets
US7356174B2 (en) * 2004-05-07 2008-04-08 General Electric Company Contraband detection system and method using variance data
US7593561B2 (en) * 2005-01-04 2009-09-22 Carestream Health, Inc. Computer-aided detection of microcalcification clusters
WO2006114003A1 (en) * 2005-04-27 2006-11-02 The Governors Of The University Of Alberta A method and system for automatic detection and segmentation of tumors and associated edema (swelling) in magnetic resonance (mri) images
US8078255B2 (en) * 2006-03-29 2011-12-13 University Of Georgia Research Foundation, Inc. Virtual surgical systems and methods
US8155409B2 (en) * 2008-04-17 2012-04-10 Ruprecht-Karls-Universitat Wave field microscope with sub-wavelength resolution and methods for processing microscopic images to detect objects with sub-wavelength dimensions
US8953856B2 (en) * 2008-11-25 2015-02-10 Algotec Systems Ltd. Method and system for registering a medical image

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105518477A (zh) * 2013-09-05 2016-04-20 皇家飞利浦有限公司 使用空间自适应正则化以用于图像重建的mri
CN105518477B (zh) * 2013-09-05 2019-08-20 皇家飞利浦有限公司 使用空间自适应正则化以用于图像重建的mri
CN106464798A (zh) * 2014-04-11 2017-02-22 富士胶片株式会社 图像处理装置、摄像装置、图像处理方法以及程序
CN106464798B (zh) * 2014-04-11 2018-01-05 富士胶片株式会社 图像处理装置、摄像装置、图像处理方法以及记录介质
CN109690620A (zh) * 2016-09-12 2019-04-26 松下知识产权经营株式会社 三维模型生成装置以及三维模型生成方法
CN109690620B (zh) * 2016-09-12 2023-05-30 松下知识产权经营株式会社 三维模型生成装置以及三维模型生成方法
CN109044356A (zh) * 2018-09-11 2018-12-21 即智数字科技(苏州)有限公司 一种共享式mri三维医疗影像平台
CN109784191A (zh) * 2018-12-20 2019-05-21 华南理工大学 一种基于商图像的多任务人脸光照编辑方法

Also Published As

Publication number Publication date
WO2012001648A2 (en) 2012-01-05
WO2012001648A3 (en) 2013-01-03
JP2014505491A (ja) 2014-03-06
US20130202177A1 (en) 2013-08-08
EP2588374A2 (en) 2013-05-08

Similar Documents

Publication Publication Date Title
CN102203826B (zh) 医学图像的降噪
CN103069432A (zh) 医学图像的非线性分辨率降低方法
EP2238742B1 (en) Noise reduction of images
JP5038642B2 (ja) ボリュメトリック画像強調システム及び方法
US9262827B2 (en) Lung, lobe, and fissure imaging systems and methods
Yokota et al. Dynamic PET image reconstruction using nonnegative matrix factorization incorporated with deep image prior
US5825909A (en) Automated method and system for image segmentation in digital radiographic images
US6463167B1 (en) Adaptive filtering
US5757953A (en) Automated method and system for region decomposition in digital radiographic images
Rousseau Brain hallucination
Bal et al. An efficient method for PET image denoising by combining multi-scale transform and non-local means
WO2011114243A1 (en) Functional image data enhancement and/or enhancer
Gou et al. Gradient regularized convolutional neural networks for low-dose CT image enhancement
Behrenbruch et al. Image filtering techniques for medical image post-processing: an overview
Silvoster M et al. Efficient segmentation of lumbar intervertebral disc from MR images
Wang et al. Remote sensing image magnification study based on the adaptive mixture diffusion model
Zheng et al. SR-CycleGAN: super-resolution of clinical CT to micro-CT level with multi-modality super-resolution loss
CN115641487B (zh) 一种基于中子和x射线的多级判定融合方法和系统
Koutsouri et al. Image contrast enhancement through regional application of partitioned iterated function systems
Parveen Inactivated Reduced CT Denoising Network Using Image pre-processing Data
Brankov Mesh modeling, reconstruction and spatio-temporal processing of medical images
Li Low-Dose CT Image Denoising Using Deep Learning Methods
Alenezi Novel Methods for Improved Fusion of Medical Images
Madhav et al. A Modern Survey: On Various Approaches of a Natural Scene from Multi Exposure Images
Luo Wavelet-Based Image Registration and Segmentation

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C02 Deemed withdrawal of patent application after publication (patent law 2001)
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20130424