CN104463815B - Dsa图像生成方法及系统 - Google Patents
Dsa图像生成方法及系统 Download PDFInfo
- Publication number
- CN104463815B CN104463815B CN201410681600.9A CN201410681600A CN104463815B CN 104463815 B CN104463815 B CN 104463815B CN 201410681600 A CN201410681600 A CN 201410681600A CN 104463815 B CN104463815 B CN 104463815B
- Authority
- CN
- China
- Prior art keywords
- image
- mrow
- dsa
- msub
- mapping
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000013507 mapping Methods 0.000 claims abstract description 75
- 238000012545 processing Methods 0.000 claims abstract description 47
- 230000009466 transformation Effects 0.000 claims abstract description 21
- 238000004364 calculation method Methods 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 11
- 230000001186 cumulative effect Effects 0.000 claims description 5
- 238000005315 distribution function Methods 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 5
- 230000002708 enhancing effect Effects 0.000 abstract description 4
- 238000002601 radiography Methods 0.000 abstract 2
- 230000009467 reduction Effects 0.000 description 6
- 238000002583 angiography Methods 0.000 description 5
- 239000002872 contrast media Substances 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 238000003384 imaging method Methods 0.000 description 4
- 210000000988 bone and bone Anatomy 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 210000004204 blood vessel Anatomy 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000008570 general process Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- 230000002792 vascular Effects 0.000 description 1
Landscapes
- Image Processing (AREA)
Abstract
本发明提供一种DSA图像生成方法及系统，其中的方法包括根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像；将锐化图像的灰度值映射到指定灰度范围，形成与造影片位深相同的映射图像；将映射图像进行局部直方图均衡处理，形成细节突出的增强图像，并将细节突出的增强图像以一定的迭加权重迭加到映射图像上，生成DSA图像；其中，迭加权重可调。利用上述发明能够在灰度映射前增强差值图像，使图像的对比度更强，有效信息更突出。
Description
技术领域
本发明涉及二维图像处理技术领域，更为具体地，涉及一种数字减影血管造影过程中应用的DSA图像生成方法及系统。
背景技术
数字减影血管造影(Digital Subtraction Angiography，DSA)技术，是通过软导管将造影剂注射到感兴趣的血管部位进行X射线成像，获得血管造影片；然后通过计算机把血管造影片上的骨与软组织的影像消除，仅在造影片上突出血管的一种摄影技术。这种技术应用计算机程序进行两次成像完成，在注入造影剂之前，首先进行第一次成像，并用计算机将成像后获得的血管造影图像转换成数字信号储存起来。注入造影剂后，再次成像并转换成数字信号。将注入造影剂前后两次成像的数字信号相减，消除相同的信号，得到一幅只有造影剂的血管图像。
在实际应用中，数字减影血管造影的一般处理流程为：
(1)将蒙片和造影片分别取对数进行相减；
(2)将对数差值图像从浮点重新映射到指定的灰度范围；
(3)采用传统方法对获取的映射结果增强及去噪。
但是，传统的DSA减影处理方法只将蒙片和造影片分别取对数做差，再将差值图像按某一系数或某种规则映射到原图的灰度范围内，这样不仅增强效果不理想，还会丢失图像的部分信息，而且处理后的图像细节也不清晰，进而会影响医生诊断。
发明内容
鉴于上述问题，本发明的目的是提供一种DSA图像生成方法及系统，以解决现有数字减影血管造影的增强效果不理想，处理后的图像细节不清晰，存在丢失部分有用信息的问题。
根据本发明的一个方面，提供了一种方法，包括根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像；将锐化图像的灰度值映射到指定灰度范围，形成与造影片位深相同的映射图像；将映射图像进行局部直方图均衡处理，形成细节突出的增强图像，并将细节突出的增强图像以一定的迭加权重迭加到映射图像上，生成DSA图像；其中，迭加权重可调。
其中，对差值图像进行锐化处理，形成锐化图像的计算公式为：
Isharp(n,m)＝Isub(n,m)+μz(n,m)
其中，Isharp(n,m)为锐化图像，Isub(n,m)为差值图像，μ为增强因子，z(n,m)为对差值图像Isub(n,m)进行滤波后的结果，(n,m)表示二维图像的行像素值与列像素值。
其中，z(n,m)的计算公式为：
其中，Sd＝{Isub(i′,j′),|Isub(i′,j′)-Isub(i,j)|≤d}，即在以像素点(i,j)为中心的窗口(n′,m′)中，与像素点(i,j)的距离不大于d的像素点的集合，M为像素点集合Sd的像素个数，d为窗口参数，窗口(n′,m′)∈(n,m)。
其中，将映射图像进行局部直方图均衡处理，形成细节突出的增强图像的过程包括根据映射图像中的每个像素点，获取以每个像素点为中心的R区域内的直方图，并根据直方图获取细节突出的增强图像；其中，
对于映射图像中的每个像素点(i,j)，根据公式：
计算以象素点(i,j)为中心的矩形区域内的直方图；
其中，矩形区域为R：W*W，其中，W＝2w+1，w为步长，p(rk)为灰度级为rk的发生频率估计，rk为第k个灰度级，nk为图像中灰度级为rk的像素个数，其中k＝0,1,…,2L-1，L为蒙片的位深；
根据灰度级为rk的发生频率估计p(rk)计算累积分布函数：
然后，对像素点(i,j)做相应变换，根据公式：
Ieq(i,j)＝(2L-1)Pw(Iunsign(i,j))得到均衡后的图形数据；
其中，Ieq为细节突出的增强图像，Iunsign为映射图像。
其中，将细节突出的增强图像迭加到映射图像上，生成DSA图像的计算公式为：
IDSA＝α*Iunsign+(1-α)*Ieq
其中，α取0.8～0.9，为可以调整的迭加权重，IDSA为DSA图像，Iunsign为映射图像，Ieq为细节突出的增强图像。
根据本发明的另一方面，提供了一种DSA图像生成系统，包括锐化图像形成单元，用于根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像；映射图像形成单元，用于将锐化图像的灰度值映射到指定灰度范围，形成与造影片位深相同的映射图像；DSA图像形成单元，用于将映射图像进行局部直方图均衡处理，形成细节突出的增强图像，并将细节突出的增强图像以一定的迭加权重迭加到映射图像上，生成DSA图像；其中，迭加权重可调。
利用上述根据本发明的DSA图像生成方法及系统，用户可以结合临床经验，自行调节迭加权重参数，设定经验值，得到细节更加清晰的减影图像，使用户获得最佳体验。同时可以避免进口价格昂贵的设备，仅在常规X光机上嵌入DSA系统即可增加减影功能，提高我国中小医院DSA系统的普及率。
为了实现上述以及相关目的，本发明的一个或多个方面包括后面将详细说明并在权利要求中特别指出的特征。下面的说明以及附图详细说明了本发明的某些示例性方面。然而，这些方面指示的仅仅是可使用本发明的原理的各种方式中的一些方式。此外，本发明旨在包括所有这些方面以及它们的等同物。
附图说明
通过参考以下结合附图的说明及权利要求书的内容，并且随着对本发明的更全面理解，本发明的其它目的及结果将更加明白及易于理解。在附图中：
图1为根据本发明实施例的DSA图像生成方法的流程图；
图2为根据本发明实施例的DSA图像生成方法的造影片图像示意图；
图3为根据图2中的造影片生成的映射图像示意图；
图4-1为根据图3中的映射图像调整步长生成的DSA图像一；
图4-2为根据图3中的映射图像调整步长生成的DSA图像二；
图4-3为根据图3中的映射图像调整步长生成的DSA图像三；
图4-4为根据图3中的映射图像调整步长生成的DSA图像四；
图5-1为根据图3中的映射图像调整迭加权重生成的DSA图像一；
图5-2为根据图3中的映射图像调整迭加权重生成的DSA图像二；
图5-3为根据图3中的映射图像调整迭加权重生成的DSA图像三；
图6为根据本发明实施例的DSA图像生成方法的流程图；
图7为根据本发明实施例的DSA图像生成系统的方框示意图。
在所有附图中相同的标号指示相似或相应的特征或功能。
具体实施方式
在下面的描述中，出于说明的目的，为了提供对一个或多个实施例的全面理解，阐述了许多具体细节。然而，很明显，也可以在没有这些具体细节的情况下实现这些实施例。
由于目前数字减影血管造影的处理效果不理想，存在丢失部分有用信息的问题，本发明在传统处理方法的基础上，首先对图像进行锐化增强处理，增强差值图像的边界信息及对比度，提高边界提取的准确性，然后通过调节局部直方图均衡中的窗口大小以及迭加权重来选择想要达到的图像增强效果，增强整体图像数据处理的灵活性，易于测试调节，提高中小医院DSA系统的普及率。
对于下述根据本发明提供的DSA图像生成方法及系统，二维图像均使用I表示(例如，IMask表示蒙片)，其中图像I表示对应二维图像的灰度值，对图像进行的各种运算即对各图像的灰度值进行的运算，以下不做具体区分。
以下将结合附图对本发明的具体实施例进行详细描述。
图1示出了根据本发明实施例的DSA图像生成方法的流程。
如图1所示，本发明提供的DSA图像生成方法包括：
S110：根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像。
具体地，首先确定用于对比的基准图像作为蒙片，并对该蒙片的灰度值作对数变换，生成相应的蒙片对数数据，然后对造影片的灰度值作对数变换，生成相应的造影片对数数据，最后将获取的造影片对数数据与蒙片对数数据相减形成差值图像数据，并通过差值图像数据获取差值图像，将该差值图像用锐化方式增强，形成锐化图像，增强图像对比度，使差值图像的边缘更清晰。
在本发明的一个具体实施方式中，确定用于对比的基准图像作为蒙片IMask，并对该蒙片的灰度值进行对数变换，生成蒙片对数数据IMask′。考虑到DSA图像的灰度最小值为0，对灰度值作对数变换的公式修改为：IMask′＝log(IMask+1)；依次遍历其他不同帧的造影片ILive，并分别对其灰度值作对数变换，生成造影片对数数据ILive′＝log(ILive+1)，将其分别与蒙片对数数据IMask′相减，得到多个抑制骨骼区域的差值图像Isub＝ILive′-IMask′；然后将各差值图像Isub用锐化方式增强，锐化图像的计算方法如下：
Isharp(n,m)＝Isub(n,m)+μz(n,m)
其中，Isharp(n,m)为锐化图像，Isub(n,m)为差值图像，μ为增强因子，z(n,m)为对差值图像Isub(n,m)进行滤波后的结果，(n,m)表示二维图像的行像素值与列像素值。
其中，z(n,m)的计算公式为：
其中，Sd＝{Isub(i′,j′),|Isub(i′,j′)-Isub(i,j)|≤d}，即在窗口(n′,m′)中，以像素点(i,j)为中心与像素点(i,j)的距离不大于d的像素点的集合，M为像素点集合Sd的像素个数，d为窗口参数，通过调整d的数值，可以调整计算量及计算速度等，窗口(n′,m′)∈(n,m)。
该步骤得到锐化图像Isharp，使变换前后的差值图像的边缘更清晰，增强差值图像的对比度。
S120：将锐化图像的灰度值映射到指定灰度范围，形成与造影片位深相同的映射图像。
具体地，考虑到将图像对数变换相减后会大大缩小差值图像灰度的动态范围，因此需要进行相应的灰度映射来扩大图像的动态范围。在本发明中，利用线性映射方法将锐化图像Isharp的灰度值映射到指定灰度范围[0，2L-1]，其中L为蒙片的位深，这样能够得到与原造影片位深相等的映射图像Iunsign。
S130：将映射图像进行局部直方图均衡处理，形成细节突出的增强图像，并将细节突出的增强图像以一定的迭加权重迭加到映射图像上，形成DSA图像；其中，迭加权重可调。
具体地，对映射图像Iunsign做局部直方图均衡处理，提取出图像中更为突出的纹理信息，主要包括根据映射图像中的每个像素点，获取以每个像素点为中心的R区域内的直方图，并根据直方图获取细节突出的增强图像；其中，
(1)设矩形区域为R：W*W，其中，W＝2w+1，w为步长，对于映射图像Iunsign中的每个像素点(i,j)，根据公式：
计算以像素点(i,j)为中心的矩形区域内的直方图；
其中，p(rk)为灰度级为rk的发生频率估计，rk为第k个灰度级，nk为图像中灰度级为rk的像素个数，其中k＝0,1,…,2L-1，L为蒙片的位深。
(2)根据灰度级为rk的发生频率估计p(rk)计算累积分布函数：
然后，对像素点(i,j)做相应变换，根据公式：
Ieq(i,j)＝(2L-1)Pw(Iunsign(i,j))得到均衡后的图形数据；
其中，Ieq为细节突出的增强图像，Iunsign为映射图像。
该步骤中，可调的参数只有用来控制计算窗口大小的W，小窗口能够增强图像的细节，而大窗口能够使图像的整体轮廓更加清晰。在实际应用中，根据具体需求相应设置窗口的大小。
通过局部直方图均衡处理获得细节突出的增强图像时，会引入较多噪声，因此，在获取到细节突出的增强图像Ieq之后，将细节突出的增强图像Ieq以一定的迭加权重迭加到映射后的映射图像Iunsign上，即可得到DSA图像IDSA，其计算公式为：
IDSA＝α*Iunsign+(1-α)*Ieq
其中，IDSA为DSA图像，α取0.8～0.9，为可以调整的迭加权重，通过调整迭加权重α来调节图像细节增强所占的比重，Iunsign为映射图像。
需要说明的是，在本发明提供的DSA图像生成方法中，α的取值在0.8～0.9范围内具有较好的图像视觉效果。其中，在α值小于0.8时，细节突出的增强图像所占比重会较大，获取的DSA图像会较粗糙，不易于阅片；在α值大于0.9时，DSA图像的细节不能有效的突出，对比度不强，失去图像增强的意义。参数α的设置较灵活，可以设定默认值，也可以设定为用户自行调节的选项，用户可以结合临床经验，自行调节该参数，设定经验值。本发明融入用户的临床经验，提供更为灵活的调节选项，使用户获得最佳体验。
具体地，作为示例，图2示出了根据本发明实施例的DSA图像生成方法的造影片ILive。
在该造影片的基础上，通过其与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像，将锐化图像的灰度值映射到指定灰度范围，形成与所述造影片位深相同的映射图像；其中，图3示出了根据图2中的造影片生成的映射图像。
在该映射图像Iunsign的基础上，对其做局部直方图均衡处理，获取细节突出的增强图像，并将细节突出的增强图像以一定的迭加权重迭加到映射图像上，形成DSA图像。具体通过以下实例进行说明：
实施例一
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝5，α＝0.9，矩形区域R为：(5*2+1)*(5*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.9)的迭加权重迭加至映射图像上，获取DSA图像；其中，图4-1示出了根据图3中的映射图像调整步长生成的DSA图像一，该实施例中的DSA图像如图4-1所示。
实施例二
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝10，α＝0.9，矩形区域R为：(10*2+1)*(10*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.9)的迭加权重迭加至映射图像上，获取DSA图像；其中，图4-2示出了根据图3中的映射图像调整步长生成的DSA图像二，该实施例中的DSA图像如图4-2所示。
实施例三
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝15，α＝0.9，矩形区域R为：(15*2+1)*(15*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.9)的迭加权重迭加至映射图像上，获取DSA图像；其中，图4-3示出了根据图3中的映射图像调整步长生成的DSA图像三，该实施例中的DSA图像如图4-3所示。
实施例四
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝30，α＝0.9，矩形区域R为：(30*2+1)*(30*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.9)的迭加权重迭加至映射图像上，获取DSA图像；其中，图4-4示出了根据图3中的映射图像调整步长生成的DSA图像四，该实施例中的DSA图像如图4-4所示
实施例五
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝15，α＝0.95，矩形区域R为：(15*2+1)*(15*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.95)的迭加权重迭加至映射图像上，获取DSA图像；其中，图5-1示出了根据图3中的映射图像调整迭加权重生成的DSA图像一，该实施例中的DSA图像如图5-1所示。
实施例六
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝15，α＝0.8，矩形区域R为：(15*2+1)*(15*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.8)的迭加权重迭加至映射图像上，获取DSA图像；其中，图5-2示出了根据图3中的映射图像调整迭加权重生成的DSA图像二，该实施例中的DSA图像如图5-2所示。
实施例七
在对图3中的映射图像进行局部直方图均衡处理、迭加细节突出的增强图像的过程中，设置步长w＝15，α＝0.75，矩形区域R为：(15*2+1)*(15*2+1)，通过对映射图形的局部直方图均衡处理及将细节突出的增强图像以(1-0.75)的迭加权重迭加至映射图像上，获取DSA图像；其中，图5-3示出了根据图3中的映射图像调整迭加权重生成的DSA图像三，该实施例中的DSA图像如图5-3所示。
由上述各实施例可以看出，在步长w＝5时，DSA图像引入了明显的噪声；步长w＝30时，计算量会增大，图像处理效率低，而步长w＝15时噪声稍低，计算量也较小；此外，在α值小于0.8时，细节突出的增强图像所占比重会较大，获取的DSA图像会较粗糙，不易于阅片；在α值大于0.9时，DSA图像的细节不能有效的突出，对比度不强。可知，在本发明提供的DSA图像生成方法中，α的取值在0.8～0.9范围内具有较好的图像视觉效果。
需要说明的是，在获取到DSA图像后，可以将DSA图像进行降噪处理，再次调高图像质量。
具体地，将得到的DSA图像IDSA进行降噪处理，具体计算如下：
该降噪处理，计算量小，在提高图像质量的同时能够提高图像的处理速度。
为了更详细的描述本发明的DSA图像生成方法，图6示出了根据本发明另一实施例的DSA图像生成方法流程。如图6所示，本实施例的DSA图像生成方法包括：
分别输入Mask图像(步骤S601)和Live图像(步骤S603)；
根据输入的Mask图像和Live图像分别计算Mask图像灰度值的自然对数(步骤S602)和Live图像灰度值的自然对数(步骤S604)；
S605：将Live图像的自然对数与Mask图像的自然对数做差，形成差值图像数据；
S606：将差值图像进行锐化增强，形成锐化图像；
S607：将锐化图像灰度值映射到指定范围，形成映射图像；
S608：对映射图像做局部直方图均衡，形成细节突出的增强图像；
S609：将细节突出的增强图像与Mask图像按不同的迭加权重做和，形成DSA图像；
S610：对DSA图像进行去噪处理，得到最终的减影图像。
S611：输出减影图像。
与上述DSA图像生成方法对应，本发明还提供一种DSA图像生成系统，具体地，图7示出了根据本发明实施例的DSA图像生成系统的逻辑结构。
如图7所示，本发明提供的DSA图像生成系统700，包括锐化图像形成单元710、映射图像形成单元720和DSA图像形成单元730。
其中，锐化图像形成单元710用于根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对差值图像进行锐化处理，形成锐化图像；映射图像形成单元720，用于将锐化图像的灰度值映射到指定灰度范围，形成与造影片位深相同的映射图像；DSA图像形成单元730，用于将映射图像进行局部直方图均衡处理，形成细节突出的增强图像，并将细节突出的增强图像迭加到映射图像上，生成DSA图像。
具体地，在锐化图像形成单元710中，首先确定用于对比的基准图像作为蒙片IMask，并对该蒙片的灰度值进行对数变换，生成蒙片对数数据IMask′。考虑到DSA图像的灰度最小值为0，作对数变换的公式修改为：IMask′＝log(IMask+1)；依次遍历其他不同帧的造影片ILive，并分别对其灰度值作对数变换，生成造影片对数数据ILive′＝log(ILive+1)，分别将其与蒙片对数数据IMask′相减，得到各抑制骨骼区域的差值图像Isub＝ILive′-IMask′(以下将以一幅造影片的情况进行描述，对多幅造影片执行相同操作即可)；
然后，将差值图像Isub用锐化方式增强，计算方法如下：
Isharp(n,m)＝Isub(n,m)+μz(n,m)
其中，Isharp(n,m)为锐化图像，Isub(n,m)为差值图像，μ为增强因子，z(n,m)为对差值图像Isub(n,m)进行滤波后的结果。其中，z(n,m)的计算公式为：
其中，Sd＝{Isub(i′,j′),|Isub(i′,j′)-Isub(i,j)|≤d}，即在窗口(n′,m′)中，以像素点(i,j)为中心与像素点(i,j)的距离不大于d的像素点的集合，M为像素点集合Sd的像素个数，(n,m)表示二维图像的行像素值与列像素值。
锐化图像形成单元710得到的锐化图像Isharp(n,m)，使变换前后的差值图像的边缘更清晰，增强差值图像的对比度。
考虑到将图像对数变换相减后会大大缩小差值图像灰度的动态范围，因此需要进行相应的灰度映射来扩大图像的动态范围。在映射图像形成单元720中，利用线性映射方法将锐化图像形成单元710得到的锐化图像Isharp(n,m)的灰度值映射到指定灰度范围[0，2L-1]，其中L为蒙片的位深，这样能够得到与原造影片位深相同的映射图像Iunsign。
在DSA图像形成单元730中，对映射后的映射图像Iunsign做局部直方图均衡处理，提取出图像中更为突出的纹理信息，主要包括根据映射图像中的每个像素点，获取以每个像素点为中心的R区域内的直方图，并根据直方图获取细节突出的增强图像，其中，
(1)设矩形区域为R：W*W，其中，W＝2w+1，w为步长，对于映射图像Iunsign中的每个像素点(i,j)，计算以像素点(i,j)为中心的矩形区域内的直方图的公式为：
其中，p(rk)为灰度级为rk的发生频率估计，rk为第k个灰度级，nk为图像中灰度级为rk的像素个数，其中k＝0,1,…,2L-1，L为蒙片的位深。
(2)根据所述灰度级为rk的发生频率估计p(rk)计算累积分布函数：
然后，对像素点(i,j)做相应变换，根据公式：
Ieq(i,j)＝(2L-1)Pw(Iunsign(i,j))得到均衡后的图形数据；
其中，Ieq为细节突出的增强图像，Iunsign为映射图像，L为蒙片的位深。
该单元中，可调的参数只有用来控制计算窗口大小的W，小窗口能够增强图像的细节，而大窗口能够使图像的整体轮廓更加清晰。在实际应用中，根据具体需求响应设置窗口的大小。
在获取到细节突出的增强图像Ieq之后，将细节突出的增强图像Ieq以一定的迭加权重迭加到映射后的映射图像Iunsign上，即可得到增强细节的DSA图像IDSA，其计算公式为：
IDSA＝α*Iunsign+(1-α)*Ieq
其中，IDSA为DSA图像，α为可以调整的迭加权重，通过调整迭加权重α来调节图像细节增强所占的比重。
将DSA图像进行降噪处理，形成最终减影图像。
将得到的DSA图像IDSA进行降噪处理，具体计算如下：
该降噪处理，计算量小，能够提高图像的处理速度。
需要说明的是，本发明提供的DSA图像生成方法，α的取值在0.8～0.9范围内具有较好的图像视觉效果。其中，迭加权重α的设置较灵活，可以设定默认值，也可以根据临床经验自行进行调节，图像处理灵活，使用户获得最佳体验。
本发明提供的DSA图像生成方法及系统，首先对差值图像进行锐化增强，增强差值图像的边界信息及其对比度，提高边界提取的准确性，然后对图像进行局部直方图均衡处理，提取图像细节，捕捉图像的有效信息，以得到更为清晰的减影图像，同时，通过调节局部直方图均衡中的窗口大小来选择想要达到的增强效果，小窗口能增强图像的细节、大窗口能增强图像的轮廓，并通过调整迭加权重α来调节细节增强所占的比重，增强整体图像数据处理的灵活性，易于测试调节。因此，根据本发明提供的DSA图像生成方法及系统能够避免进口设备价格昂贵的问题，仅在常规X光机上嵌入DSA系统即可实现减影功能，提高我国中小医院DSA系统的普及率。
如上参照附图以示例的方式描述根据本发明的DSA图像生成方法及系统。但是，本领域技术人员应当理解，对于上述本发明所提出的DSA图像生成方法及系统，还可以在不脱离本发明内容的基础上做出各种改进。因此，本发明的保护范围应当由所附的权利要求书的内容确定。
Claims (10)
1.一种DSA图像生成方法，包括：
根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对所述差值图像进行锐化处理，形成锐化图像；
将所述锐化图像的灰度值映射到指定灰度范围[0，2L-1]，形成与所述造影片位深相同的映射图像，L为所述蒙片的位深；
将所述映射图像进行局部直方图均衡处理，根据所述映射图像中的每个像素点，获取以每个像素点为中心的R区域内的直方图，并根据所述直方图获取细节突出的增强图像，并将所述细节突出的增强图像以一定的迭加权重迭加到所述映射图像上，生成DSA图像；其中，所述迭加权重可调，所述R区域为：W*W，其中，W＝2w+1，w为步长。
2.如权利要求1所述的DSA图像生成方法，其中，对所述差值图像进行锐化处理，形成锐化图像的计算公式为：
Isharp(n,m)＝Isub(n,m)+μz(n,m)
其中，Isharp(n,m)为锐化图像，Isub(n,m)为差值图像，μ为增强因子，z(n,m)为对差值图像Isub(n,m)进行滤波后的结果，(n,m)表示二维图像的行像素值与列像素值。
3.如权利要求2所述的DSA图像生成方法，其中，所述z(n,m)的计算公式为：
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>I</mi>
<mrow>
<mi>s</mi>
<mi>u</mi>
<mi>b</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mi>M</mi>
</mfrac>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>d</mi>
</msub>
</mrow>
</munder>
<msub>
<mi>I</mi>
<mrow>
<mi>s</mi>
<mi>u</mi>
<mi>b</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</mrow>
其中，Sd＝{Isub(i′,j′),|Isub(i′,j′)-Isub(i,j)|≤d}，即在以像素点(i,j)为中心的窗口(n′,m′)中，与像素点(i,j)的距离不大于d的像素点的集合，M为像素点集合Sd的像素个数，d为窗口参数，窗口(n′,m′)∈(n,m)。
4.如权利要求1所述的DSA图像生成方法，其中，
根据所述映射图像中的每个像素点，获取以所述每个像素点为中心的R区域内的直方图，并根据所述直方图获取细节突出的增强图像的过程包括：
对于所述映射图像中的每个像素点(i,j)，根据公式：
计算以象素点(i,j)为中心的矩形区域内的直方图；
其中，矩形区域为R：W*W，其中，W＝2w+1，w为步长，p(rk)为灰度级为rk的发生频率估计，rk为第k个灰度级，nk为图像中灰度级为rk的像素个数，其中k＝0,1,…,2L-1；
根据所述灰度级为rk的发生频率估计p(rk)计算累积分布函数：
<mrow>
<msub>
<mi>p</mi>
<mi>w</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>r</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>t</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mi>k</mi>
</munderover>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>r</mi>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
然后，对像素点(i,j)做相应变换，根据公式：
Ieq(i,j)＝(2L-1)Pw(Iunsign(i,j))得到均衡后的图形数据；
其中，Ieq为细节突出的增强图像，Iunsign为映射图像。
5.如权利要求1所述的DSA图像生成方法，其中，将所述细节突出的增强图像迭加到所述映射图像上，生成DSA图像的计算公式为：
IDSA＝α*Iunsign+(1-α)*Ieq
其中，α取0.8～0.9，为可以调整的迭加权重，IDSA为DSA图像，Iunsign为映射图像，Ieq为细节突出的增强图像。
6.一种DSA图像生成系统，包括：
锐化图像形成单元，用于根据造影片与蒙片的灰度值对数变换后相减的差值图像数据，获取差值图像，并对所述差值图像进行锐化处理，形成锐化图像；
映射图像形成单元，用于将所述锐化图像的灰度值映射到指定灰度范围[0，2L-1]，形成与所述造影片位深相同的映射图像，L为所述蒙片的位深；
DSA图像形成单元，用于将所述映射图像进行局部直方图均衡处理，根据所述映射图像中的每个像素点，获取以每个像素点为中心的R区域内的直方图，并根据所述直方图获取细节突出的增强图像，并将所述细节突出的增强图像以一定的迭加权重迭加到所述映射图像上，生成DSA图像；其中，所述迭加权重可调，所述R区域为：W*W，其中，W＝2w+1，w为步长。
7.如权利要求6所述的DSA图像生成系统，其中，
在所述锐化图像形成单元中，对所述差值图像进行锐化处理，形成锐化图像的计算公式为：
Isharp(n,m)＝Isub(n,m)+μz(n,m)
其中，Isharp(n,m)为锐化图像，Isub(n,m)为差值图像，μ为增强因子，z(n,m)为对差值图像Isub(n,m)进行滤波后的结果，(n,m)表示二维图像的行像素值与列像素值。
8.如权利要求7所述的DSA图像生成系统，其中，
所述z(n,m)的计算公式为：
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>I</mi>
<mrow>
<mi>s</mi>
<mi>u</mi>
<mi>b</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mi>M</mi>
</mfrac>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>d</mi>
</msub>
</mrow>
</munder>
<msub>
<mi>I</mi>
<mrow>
<mi>s</mi>
<mi>u</mi>
<mi>b</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</mrow>
其中，Sd＝{Isub(i′,j′),|Isub(i′,j′)-Isub(i,j)|≤d}，即在以像素点(i,j)为中心的窗口(n′,m′)中，与像素点(i,j)的距离不大于d的像素点的集合，M为像素点集合Sd的像素个数，d为窗口参数，窗口(n′,m′)∈(n,m)。
9.如权利要求6所述的DSA图像生成系统，其中，
根据所述映射图像中的每个像素点，获取以所述每个像素点为中心的R区域内的直方图，并根据所述直方图获取细节突出的增强图像的过程包括：
对于所述映射图像中的每个像素点(i,j)，根据公式：
计算以像素点(i,j)为中心的矩形区域内的直方图；
其中，矩形区域为R：W*W，其中，W＝2w+1，w为步长，p(rk)为灰度级为rk的发生频率估计，rk为第k个灰度级，nk为图像中灰度级为rk的像素个数，其中k＝0,1,…,2L-1；
根据所述灰度级为rk的发生频率估计p(rk)计算累积分布函数：
<mrow>
<msub>
<mi>p</mi>
<mi>w</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>r</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>t</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mi>k</mi>
</munderover>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>r</mi>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
然后，对像素点(i,j)做相应变换，根据公式：
Ieq(i,j)＝(2L-1)Pw(Iunsign(i,j))得到均衡后的图形数据；
其中，Ieq为细节突出的增强图像，Iunsign为映射图像。
10.如权利要求6所述的DSA图像生成系统，其中，
在所述DSA图像形成单元中，将所述细节突出的增强图像迭加到所述映射图像上，生成DSA图像的计算公式为：
IDSA＝α*Iunsign+(1-α)*Ieq
其中，α取0.8～0.9，为可以调整的迭加权重，IDSA为DSA图像，Iunsign为映射图像，Ieq为细节突出的增强图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410681600.9A CN104463815B (zh) | 2014-11-24 | 2014-11-24 | Dsa图像生成方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410681600.9A CN104463815B (zh) | 2014-11-24 | 2014-11-24 | Dsa图像生成方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104463815A CN104463815A (zh) | 2015-03-25 |
| CN104463815B true CN104463815B (zh) | 2017-08-29 |
Family
ID=52909803
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410681600.9A Active CN104463815B (zh) | 2014-11-24 | 2014-11-24 | Dsa图像生成方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104463815B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106238347A (zh) * | 2016-07-27 | 2016-12-21 | 合肥高晶光电科技有限公司 | 一种ccd色选机的图像处理方法 |
| CN107633512B (zh) * | 2017-09-18 | 2020-04-17 | 浙江莱达信息技术有限公司 | 一种医学x光影像滤波处理的二级修正方法 |
| CN111613120B (zh) * | 2020-04-08 | 2023-03-14 | 宁波创导三维医疗科技有限公司 | 一种介入手术操作造影成像效果模拟系统 |
| CN111862074B (zh) * | 2020-07-30 | 2023-11-24 | 国网湖南省电力有限公司 | 一种电缆阻水缓冲层缺陷识别方法及装置 |
| CN111915538B (zh) * | 2020-08-19 | 2024-03-19 | 南京普爱医疗设备股份有限公司 | 一种用于数字血管减影的图像增强方法及系统 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080051648A1 (en) * | 2006-08-25 | 2008-02-28 | Suri Jasjit S | Medical image enhancement system |
| CN101853497A (zh) * | 2010-02-25 | 2010-10-06 | 杭州海康威视软件有限公司 | 一种图像增强方法和装置 |
| CN104156931B (zh) * | 2014-09-04 | 2017-03-29 | 成都金盘电子科大多媒体技术有限公司 | 一种数字减影血管造影方法 |
-
2014
- 2014-11-24 CN CN201410681600.9A patent/CN104463815B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN104463815A (zh) | 2015-03-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9349198B2 (en) | Robust artifact reduction in image reconstruction | |
| US10643319B2 (en) | Apparatus and method for context-oriented blending of reconstructed images | |
| CN104463815B (zh) | Dsa图像生成方法及系统 | |
| Wang et al. | Metal artifact reduction in CT using fusion based prior image | |
| US9002134B2 (en) | Multi-scale image normalization and enhancement | |
| US10867375B2 (en) | Forecasting images for image processing | |
| EP2579207A2 (en) | Method of noise reduction in digital x-ray frames series | |
| CN106530236B (zh) | 一种医学图像处理方法及系统 | |
| CN105741241B (zh) | 基于合成增强图像的肿瘤区域图像增强方法及系统 | |
| JPS59229670A (ja) | 空間分散フイルタを用いる多重測定雑音低減装置 | |
| CN105321155A (zh) | 一种cbct图像环形伪影消除方法 | |
| JP6071444B2 (ja) | 画像処理装置及びその作動方法、プログラム | |
| CN107680057A (zh) | 超声图像增强的方法及装置 | |
| CN111402150B (zh) | 一种ct图像金属伪影去除方法及装置 | |
| JP2019524356A (ja) | 異なる反復から抽出された特徴画像を使用する特徴ベースの画像処理 | |
| JP2019088735A (ja) | 医用画像処理装置および医用画像処理方法 | |
| WO2011114243A1 (en) | Functional image data enhancement and/or enhancer | |
| CN116018611A (zh) | 使用深度卷积网络的噪声抑制 | |
| Johari et al. | Metal artifact suppression in dental cone beam computed tomography images using image processing techniques | |
| CN113205461B (zh) | 一种低剂量ct影像去噪模型训练方法、去噪方法及装置 | |
| US11771390B2 (en) | Method and device for determining the contour of anatomical structures in a digital X-ray-based fluoroscopic image | |
| EP3658031B1 (en) | Motion compensated cardiac valve reconstruction | |
| US9153016B2 (en) | Method and system for enhancing contrast of spatially-localized phenomena in mammography image | |
| CN110506294B (zh) | 数字x射线图像中具有低信息内容的区域的探测 | |
| JP2005012771A (ja) | 放射線撮像画像のコントラスト及び輝度を設定するための方法及び装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |