CN108852385A - 一种x射线造影方法及基于x射线造影的动态阅片方法 - Google Patents
一种x射线造影方法及基于x射线造影的动态阅片方法 Download PDFInfo
- Publication number
- CN108852385A CN108852385A CN201810203889.1A CN201810203889A CN108852385A CN 108852385 A CN108852385 A CN 108852385A CN 201810203889 A CN201810203889 A CN 201810203889A CN 108852385 A CN108852385 A CN 108852385A
- Authority
- CN
- China
- Prior art keywords
- frequency
- signal
- ray imaging
- contrast agent
- angiographic image
- 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
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 65
- 238000000034 method Methods 0.000 title claims abstract description 41
- 238000001914 filtration Methods 0.000 claims abstract description 16
- 230000009466 transformation Effects 0.000 claims abstract description 6
- 238000004458 analytical method Methods 0.000 claims abstract description 4
- 239000002872 contrast media Substances 0.000 claims description 66
- 210000004204 blood vessel Anatomy 0.000 claims description 53
- 238000002601 radiography Methods 0.000 claims description 21
- 230000010412 perfusion Effects 0.000 claims description 19
- 230000004907 flux Effects 0.000 claims description 9
- 238000010521 absorption reaction Methods 0.000 claims description 5
- 150000002496 iodine Chemical class 0.000 claims description 5
- XQZXYNRDCRIARQ-LURJTMIESA-N iopamidol Chemical compound C[C@H](O)C(=O)NC1=C(I)C(C(=O)NC(CO)CO)=C(I)C(C(=O)NC(CO)CO)=C1I XQZXYNRDCRIARQ-LURJTMIESA-N 0.000 claims description 5
- 150000003839 salts Chemical class 0.000 claims description 5
- AMDBBAQNWSUWGN-UHFFFAOYSA-N Ioversol Chemical compound OCCN(C(=O)CO)C1=C(I)C(C(=O)NCC(O)CO)=C(I)C(C(=O)NCC(O)CO)=C1I AMDBBAQNWSUWGN-UHFFFAOYSA-N 0.000 claims description 3
- 229960004359 iodixanol Drugs 0.000 claims description 3
- NBQNWMBBSKPBAY-UHFFFAOYSA-N iodixanol Chemical compound IC=1C(C(=O)NCC(O)CO)=C(I)C(C(=O)NCC(O)CO)=C(I)C=1N(C(=O)C)CC(O)CN(C(C)=O)C1=C(I)C(C(=O)NCC(O)CO)=C(I)C(C(=O)NCC(O)CO)=C1I NBQNWMBBSKPBAY-UHFFFAOYSA-N 0.000 claims description 3
- 229960002603 iopromide Drugs 0.000 claims description 3
- DGAIEPBNLOQYER-UHFFFAOYSA-N iopromide Chemical compound COCC(=O)NC1=C(I)C(C(=O)NCC(O)CO)=C(I)C(C(=O)N(C)CC(O)CO)=C1I DGAIEPBNLOQYER-UHFFFAOYSA-N 0.000 claims description 3
- 229960004537 ioversol Drugs 0.000 claims description 3
- 229960001025 iohexol Drugs 0.000 claims description 2
- NTHXOOBQLCIOLC-UHFFFAOYSA-N iohexol Chemical compound OCC(O)CN(C(=O)C)C1=C(I)C(C(=O)NCC(O)CO)=C(I)C(C(=O)NCC(O)CO)=C1I NTHXOOBQLCIOLC-UHFFFAOYSA-N 0.000 claims description 2
- 150000001412 amines Chemical class 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 11
- 210000001519 tissue Anatomy 0.000 description 20
- 238000005516 engineering process Methods 0.000 description 15
- 238000002583 angiography Methods 0.000 description 11
- 239000000523 sample Substances 0.000 description 9
- 230000005855 radiation Effects 0.000 description 6
- 230000029058 respiratory gaseous exchange Effects 0.000 description 5
- 230000002123 temporal effect Effects 0.000 description 5
- 210000000269 carotid artery external Anatomy 0.000 description 4
- 230000000737 periodic effect Effects 0.000 description 4
- 239000012472 biological sample Substances 0.000 description 3
- 238000006073 displacement reaction Methods 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 230000002708 enhancing effect Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000003672 processing method Methods 0.000 description 3
- 230000011218 segmentation Effects 0.000 description 3
- 230000005469 synchrotron radiation Effects 0.000 description 3
- 230000002792 vascular Effects 0.000 description 3
- 239000008280 blood Substances 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 210000000988 bone and bone Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 239000003795 chemical substances by application Substances 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- MIKKOBKEXMRYFQ-WZTVWXICSA-N meglumine amidotrizoate Chemical compound C[NH2+]C[C@H](O)[C@@H](O)[C@H](O)[C@H](O)CO.CC(=O)NC1=C(I)C(NC(C)=O)=C(I)C(C([O-])=O)=C1I MIKKOBKEXMRYFQ-WZTVWXICSA-N 0.000 description 2
- 230000008520 organization Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000000241 respiratory effect Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- ZCYVEMRRCGMTRW-UHFFFAOYSA-N 7553-56-2 Chemical compound [I] ZCYVEMRRCGMTRW-UHFFFAOYSA-N 0.000 description 1
- 238000012935 Averaging Methods 0.000 description 1
- 241001465754 Metazoa Species 0.000 description 1
- 241000699666 Mus <mouse, genus> Species 0.000 description 1
- 241000699670 Mus sp. Species 0.000 description 1
- 238000005267 amalgamation Methods 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 229960000074 biopharmaceutical Drugs 0.000 description 1
- 210000004556 brain Anatomy 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000036461 convulsion Effects 0.000 description 1
- 238000002586 coronary angiography Methods 0.000 description 1
- 210000004351 coronary vessel Anatomy 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000011049 filling Methods 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 229910052740 iodine Inorganic materials 0.000 description 1
- 239000011630 iodine Substances 0.000 description 1
- 230000002262 irrigation Effects 0.000 description 1
- 238000003973 irrigation Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
- 238000002560 therapeutic procedure Methods 0.000 description 1
- 210000005166 vasculature Anatomy 0.000 description 1
- 238000010626 work up procedure Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/48—Diagnostic techniques
- A61B6/481—Diagnostic techniques involving the use of contrast agents
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
- A61B6/504—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of blood vessels, e.g. by angiography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
- A61B6/5264—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to motion
-
- G—PHYSICS
- G02—OPTICS
- G02B—OPTICAL ELEMENTS, SYSTEMS OR APPARATUS
- G02B27/00—Optical systems or apparatus not provided for by any of the groups G02B1/00 - G02B26/00, G02B30/00
- G02B27/02—Viewing or reading apparatus
- G02B27/022—Viewing apparatus
- G02B27/023—Viewing apparatus for viewing X-ray images using image converters, e.g. radioscopes
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Physics & Mathematics (AREA)
- Optics & Photonics (AREA)
- Radiology & Medical Imaging (AREA)
- Molecular Biology (AREA)
- High Energy & Nuclear Physics (AREA)
- Veterinary Medicine (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Biophysics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Vascular Medicine (AREA)
- Dentistry (AREA)
- Oral & Maxillofacial Surgery (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明提供一种X射线造影方法，包括：原始的造影图像序列；在原始的造影图像序列的其中一帧上选取两个特征点，提取这两个特征点的时序信号并将其傅里叶变换为频域信号；对两个特征点的频域信号进行特征分析，选取这两个特征点在相同频率下的最大的至少两个模值差；遍历所有像素点，以0Hz和其中一个特征频率为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度；输出以调制深度为造影成像参量的图像。此外，还提供了基于该X射线造影方法的动态阅片方法。本发明的X射线造影方法及基于X射线造影的动态阅片方法实现了低X射线剂量下的高信噪比成像，省去减影和帧间配准，并实现图像融合的效果。
Description
技术领域
本发明涉及一种基于多帧医学影像的X射线血管造影方法，尤其涉及一种基于时空滤波的血管造影数据处理方法及动态阅片方法。
背景技术
X射线的血管造影方法是从最终图像中移除背景结构的X射线成像方法，主要分为两种。一是K边减影血管造影方法，此方法是在略低于和略高于造影剂的K边能量下，分别得到样品的X射线影像图。然后将这两幅图像做差，利用造影剂在K边能量上对X射线的强吸收效应可得到只含有造影剂信息的血管图像；二是单能时间减影技术，在略高于造影剂K边的能量上，采集包含背景结构的掩膜图像。然后，注射造影剂，随后采集含有造影剂的血管和灌注组织的图像。从含有所注射的血管和灌注组织的图像中减去掩膜图像。这样便得到了只含有造影剂信息的血管图像；目前，一般常用单能时间减影技术针对血管病情的检查诊断、血管位置的测量、介入治疗的手术规划以及生物医学等方面研究，通过放射医学影像设备及其后续数据处理方法实现多帧血管造影图像。但是，由于人(或其他活体生物)不可避免的呼吸、心跳、以及其他部分组织的运动、血管对造影剂的应激反应和外部设备的振动等，直接减影往往会带来难以接受的人工伪影。
专利文件CN 106999129 A、CN 101822545 B、CN 102103745 A、CN 102663709 A、CN 1897033 A、CN 101799918 B中分别涉及数字减影血管造影技术的设备、方法的软件实现、以及后续的血管增强、血管分割、运动伪迹的消除、不同区域血管造影的图像融合等。
其中，CN106999129 A中方法通过额外采集不含造影剂的完整呼吸周期图像序列和心跳周期图像序列，在随后的减影操作中分别通过呼吸相位和心跳相位来匹配蒙片与盈片分别减去呼吸伪影和心跳伪影以实现减影图像质量的提升；CN 101822545 B中方法属于一种图像配准技术，使得数字减影操作蒙片与盈片的图像配准更加精确，从而去除运动伪影。CN106999129 A、专利CN 101822545 B中方法是解决运动伪迹的消除问题，但操作复杂繁琐，对原始数据质量要求较高。
CN102103745 A中用计算机软件的方法实现了数字减影操作，可实现减影图像的动态阅片效果，但本质上对数字减影成像技术未做提升。
CN 102663709 A中方法是一种血管增强算法，利用经验模态分解X射线冠脉造影图像为一系列本征模态函数，利用噪声的分布规律去除背景噪声，进而选择一些特定层本征模态函数进行加权构造冠状动脉图像，再利用基于Hessian矩阵特征值的血管测度函数增强图像中的血管结构，以改善造影图像的视觉效果；CN 1897033 A针对数字减影图像不能完全消除人体组织引起的噪声信号，以及背景和目标混杂在一起的情况，提出一种基于血管直径的先验知识以及重叠分块技术的局部阈值血管分割方法，以提取出目标血管的方法。CN 102663709 A涉及血管增强，CN 1897033 A涉及血管分割，都是基于数字减影成像操作后进行的，以期进一步提高血管图像质量。
针对造影剂在血管中分布分散不连续的情况，CN 101799918B涉及图像融合，是基于数字减影成像操作后进行的。CN 101799918 B中的方法针对在数字减影成像，单幅减影图像只能代表局部血管影像状况，提出一种通过动态模糊技术和Ridgelet变换理论的数字减影图像的自适应融合方法。使用该专利文件中提出的数据处理方法，可以一定程度的提高成像质量。
综上，上述这些专利文件中公开的X射线造影方法基于传统的数字减影成像技术，往往难以解决以下问题：
1)传统的基于数字减影成像技术的X射线造影方法为造影剂减影成像，对造影效果不理想的图像一般采用待测结构增强的算法来提高其成像质量，算法过于复杂且提升效果有限，并且这种算法强烈依赖于高辐射剂量，由于放射剂量低或造影剂浓度过低可能导致的原始造影图像序列中造影效果不理想，进而影响调制深度，降低图像的分辨率。
2)传统的基于数字减影成像技术的X射线造影方法需要减影操作以得到血管图像，其关键的一步为找到合适的蒙片与盈片。但由于造影过程中样品不可避免的运动导致蒙片与盈片间存在位移，则会极大的影响成像的信噪比。一般的采用帧间配准的方法解决此问题，有时还需要引入入侵性的外部特征点以完成帧间匹配。这样会极大的增加算法的复杂度，且对图像质量提升有限。
3)传统的基于数字减影成像技术的X射线造影方法，选取的盈片中造影剂无法完全灌注到所有血管中而是部分灌注。所以，单幅造影图像只能代表局部血管影像状况。医生观察多个造影图像序列时费时费力且容易漏诊。因此将每一幅数字减影血管造影图像融合到完整的图像中，使之包含血管整体结构，对医疗辅助诊断具有重要意义。但图像融合技术不但操作复杂，甚至会引入新的噪声。
4)传统的基于数字减影成像技术的X射线造影方法，对任意一张盈片来说无法全部找到合适的蒙片，且盈片中血管有较大的位移。所以无法实施良好的数字减影血管造影过程的动态阅片播放。
发明内容
本发明提出一种基于时空滤波的X射线造影方法及基于X射线造影的动态阅片方法，以实现低X射线剂量下的高信噪比成像，省去减影和帧间配准，并实现图像融合的效果。
为了实现上述目的，本发明提供了一种X射线造影方法，包括：步骤S1：以采样频率Fs获取N帧的原始的造影图像序列，作为当前的造影图像序列；步骤S2：在原始的造影图像序列的其中一帧上选取血管特征点和背景特征点，提取这两个特征点的时序信号并将其时序信号傅里叶变换为频域信号；步骤S3：对血管特征点和背景特征点的频域信号进行特征分析，依次选取这两个特征点在相同频率下的最大的至少两个模值差，第i大的模值差所对应的频率作为特征频率xi*Fs/N Hz；步骤S4：在当前的造影图像序列中遍历造影区域的所有像素点，以0Hz和其中一个特征频率xi*Fs/N Hz为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度MDi；步骤S5：输出以调制深度MDi为造影成像参量的造影图像；其中xi为正整数。
所述步骤S1包括：对造影对象灌注造影剂，在X光光源的照射下，以造影剂未灌注之前为起始帧至以造影剂灌注完毕为结束帧，采集原始的造影图像序列。
所述X光的能量高于所述造影剂的吸收边，光子通量为109—1010photonssecond‐ 1mm‐2。
所述造影剂为非离子型碘盐类造影剂，该非离子型碘盐类造影剂选自碘必乐、碘海醇、碘普罗胺、碘佛醇、碘克沙醇、泛影葡胺中的一种。
所述造影剂的用量为100—500微升，浓度为150—350 280mg·ml‐1。
所述步骤S3中所选取的模值差为2‐4个，且所述特征频率的个数等于模值差的个数。
所述每个像素点的调制深度MDi通过下式计算：
其中，S0(x,y,t)、Sxi(x,y,t)分别是对每个像素点的时序信号做傅里叶带通滤波处理后获得的直流信号和对应xi*Fs/N Hz频率的交流信号；
Mean(abs(Sxi(x,y,t)))为对应xi*Fs/N Hz频率的交流信号的绝对值的平均值；
Mean(abs(S0(x,y,t)))为直流信号的绝对值的平均值。
所述步骤S4还包括：在遍历造影区域的所有像素点之前，先对其每一帧图像做去噪处理，该去噪处理包括中值去噪、高斯去噪或维纳去噪。
进一步地，本发明还提供了一种基于X射线造影的动态阅片方法，包括：
步骤S1’：采用权利要求1所述X射线造影方法，得到特征频率xi*Fs/N Hz，并得到以调制深度MDi为造影成像参量的造影图像，作为原始的造影图像序列的前N帧所对应的造影图像，取j＝1；步骤S2’：选取原始的造影图像序列的前N‐j帧作为当前的造影图像序列，并以步骤S1’得到的特征频率计算所述当前的造影图像序列所对应的特征频率xi*Fs/(N‐j)Hz；步骤S3’：在当前的造影图像序列中遍历造影区域的所有像素点，以0Hz和步骤S2’得到的特征频率xi*Fs/(N‐j)Hz为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度MDi；得到以调制深度MDi为造影成像参量的造影图像，作为原始的造影图像序列的前N‐j帧所对应的造影图像；步骤S4’：当N‐j为n时，依次播放所述前n、n+1…、N‐1、N‐2、N帧所对应的造影图像；否则将j+1更新为j并回到步骤S2’。
当i为1时，所述以调制深度MD1为造影成像参量的造影图像为血管造影图像。
与现有的基于数字减影技术的X射线造影方法相比，本发明的基于时空滤波的X射线造影方法具有以下特点和有益效果。
1)本发明采用的成像参量为调制深度即造影剂灌注的特征频率信号与0频直流信号的比值，算法较为简单，并且即使由于放射剂量低或造影剂浓度过低导致的原始造影图像序列中造影效果不理想，但其调制深度几乎不受影响，从而可以实现低辐射剂量成像，和传统的数字减影造影方法相比有更高的分辨率。根据该特性，本发明中的血管造影方法不仅适合于同步辐射X射线光源和医院放射医学影像设备，对实验室X射线光管光源也同样适用。
2)本发明中的血管造影方法中并不需要减影操作，而是在频域中将造影剂灌注的特征信号与生物运动噪声的特征信号相互区别，使得所得到的血管造影图像中只含有血管信息，而背景组织的运动信息可单独成像。因此，也就不需要帧间配准，大大减小了算法的复杂度，并且使得本发明中的血管造影方法成像信噪比更高，自适应地实现运动伪迹的消除，更适合应用在具有明显运动伪影的X射线血管造影成像中。
3)本发明中的血管造影方法成像参量为调制深度MDi，其充分利用了整个造影剂灌注序列的信息，因此无需再用图像融合技术便能得到高质量的血管整体结构。针对造影剂在血管中分布分散不连续的情况，可自动实现图像融合的效果。
4)本发明可提供造影剂灌注过程的血管造影影像序列，尤其适合应用在具有明显运动伪影的X射线血管造影成像中。本发明的血管造影影像序列的动态阅片播放，其动态画面连续、清晰、血管轨迹稳定不漂移。因为血管平衡位置处的信号更倾向于阶跃信号占主导，非平衡位置处信号更倾向于类周期信号占主导。这在频域中容易区分处理，经过频域处理并不会干扰造影剂灌注信息的提取。从而，即使样品中血管也有周期性的位移，图像中的血管也稳定在血管的平衡位置，而不需要寻找相应蒙片。
附图说明
图1为根据本发明的一个实施例的原始的造影图像序列的其中一帧图像。
图2A-图2C为如图1所示的血管特征点v的对应0Hz、1*Fs/N Hz、3*Fs/NHz的时序信号图。
图2D-图2F为如图1所示的b点对应0Hz、1*Fs/N Hz、3*Fs/N Hz的时序信号图。
图3A‐3B为如图1所示的原始的造影图像序列的以MD1、MD2为成像参量的造影图像。
图4为传统的时间减影血管造影图像。
具体实施方式
以下结合具体实施例，对本发明做进一步说明。应理解，以下实施例仅用于说明本发明而非用于限制本发明的范围。
根据本发明的一个实施例，本发明的X射线造影方法基于时空滤波，具体包括：
步骤S1：原始造影数据的采集。与传统的X射线数字减影血管造影方法类似，对造影对象灌注造影剂，在X光光源的照射下，以造影剂未灌注之前为起始帧至以造影剂灌注完毕为结束帧，采集数十至数百帧的原始的造影图像序列。
在本实施例中，造影对象为实验小鼠，在动物准备过程中，在颈外动脉(ECA)与颈总动脉(CCA)的分叉点上***一个血管造影管，造影区域为小鼠脑部；造影剂为180微升的浓度为280mgIml‐1的碘必乐(Iopamiro)造影剂。其中，造影剂用量根据成像的生物样品、成像区域、X射线光子通量等因素还可在100—500微升范围内具体调节，造影剂浓度可根据生物样品耐受度、成像区域、X射线光子通量等因素还可在150—350mgI·ml‐1范围内具体调节。碘必乐(Iopamiro)造影剂为一种非离子型碘盐类造影剂，也可以由其他常用的非离子型碘盐类造影剂，如碘海醇、碘普罗胺、碘佛醇、碘克沙醇、泛影葡胺等中的一种来替代。血管造影管注入颈外动脉(ECA)，尽管造影剂本身起着提升对比度的作用，但由于生物样品骨骼或其他较厚的软组织的影响，使造影剂和生物样品组织交织叠加在一起在无减影操作下图像对比度较低，或者由于造影剂浓度较低、X射线光子通量较低、探测器探测效率低下等因素的影响，使得采集到的造影图像序列对比度相对较低；灌注速率由保定兰格恒流泵有限公司(Longerpump)生产的LSP01‐1A型微型注射泵控制，注射率为133.3μls‐1。X光光源可以是同步辐射X射线、实验室X光机、医院放射医学影像设备，本实施例中X光光源为在上海同步辐射光源BL13W实验站的实验室光源；其中X光的能量由所使用的造影剂类型决定，一般略高于造影剂的吸收边即可，在本实施例中X光的能量为33.3keV；光子通量根据X光装置的运行状态、成像可接受的辐射剂量大小等因素决定，且在可接受的辐射剂量范围内提高光子通量可相应提高成像质量。光子通量在109—1010photonssecond‐1mm‐2范围内均可，在本实施例中光子通量为2.38×1010photonssecond‐1mm‐2。探测器为PCOX射线CCD(像素大小为6.5×6.5μm，视场大小为13×13mm，PCO‐TECH Inc，德国)；探测器曝光时间为10毫秒、采集帧率即采样频率Fs为100fps。
此外，原始数据的获得也可以采用从传统的时间减影血管造影序列中截取对应造影剂灌注的那一部分序列。
对于步骤S1得到的造影图像序列，由于造影剂对X射线的吸收作用，血管部分与背景组织部分有不一样的时序变化特点。背景组织在理想状态下是一系列的静态信号，但实际操作中由于不可避免的呼吸、心跳、以及其他部分组织的运动、血管对造影剂的应激反应和外部设备的振动等原因，其在时序上是一组类周期震荡信号或是其他类型的杂乱信号。相应的血管区域，由于造影剂的灌注，理想状态下应是随着造影剂的填充情况以及浓度变化情况而变化的一组信号。实际操作中，其是一组类周期信号叠加在一个阶跃信号之上。由此，对应血管区域的信号其与背景组织信号相比的区别仅仅在于，血管信号多出一个阶跃信号，此阶跃信号对应着有造影剂灌注的血管信息。
步骤S2：在原始的造影图像序列的其中一帧上选取两个特征点，提取这两个特征点的时序信号并将其时序信号傅里叶变换为频域信号。
如图1所示为原始的造影图像序列的其中一帧。在原始的造影图像序列的这一帧上选取如下两个特征点，其一选取明显位于血管上的血管区域像素点作为血管特征点v。其二选取具有运动代表性的背景组织像素点作为背景特征点b。提取出该两特征点的时序信号，可以发现背景特征点b处的信号是周期性信号和静态直流信号的叠加。血管特征v处的信号则类似于静态直流信号和部分周期信号叠加在一个阶跃信号之上。随后对提取出的两时序信号采用快速傅里叶变换(FFT)的方法进行傅里叶频域变换。
由此得到的频域信号，设采样频率为Fs；采样个数与帧数相同，设为N；则频域中0Hz频率对应着直流(DC)信号，1*Fs/N Hz、2*Fs/N Hz、3*Fs/N Hz、…、i*Fs/N Hz、…(N‐1)*Fs/N Hz的频率对应着高频交流(AC)信号。血管特征点v在频率i*Fs/N Hz下对应的模值为Avi，即血管特征点v对应的模值包括Av1、Av2、Av3、…Avi…AvN‐1；背景特征点b在频率i*Fs/NHz下对应的模值为Abi，即背景特征点b对应的模值包括Ab1、Ab2、Ab3、…Abi…AbN‐1。
步骤S3：对血管特征点v、背景特征点b的频域信号进行特征分析，两特征点b、v在相同频率下的模值差随着频率的变化而变化，从大到小依次选取2个模值差，其所对应的频率作为特征频率。
在本实施例中，仅仅选取了模值差最大的2个频率作为特征频率，这两个特征频率x1*Fs/N Hz＝1*Fs/N Hz、x2*Fs/N Hz＝3*Fs/N Hz，其中，最大的模值差所对应的特征频率由x1*Fs/N Hz表示，大小为1*Fs/N Hz，对应着血管；次最大的模值差所对应的特征频率由x2*Fs/N Hz表示，大小为3*Fs/N Hz，对应着运动的背景组织。
其中，模值差与X射线强度、探测器曝光时间、造影剂剂量、样品类型厚度等因素有关。我们在实践中，在样品可接受的辐射剂量下，尽量采集图像质量较好的原始的造影图像序列，一般情况下，模值具有明显的可区分的差值，大约为5‐20倍的模值。因此，通常也可以依次取模值差大于模值的5倍的2‐4个模值差所分别对应的频率作为特征频率(第i大的模值差所对应的频率作为特征频率xi*Fs/N Hz)，而不仅仅局限于2个，以轻易区分血管信号和背景组织信号。最大的模值差对应的特征频率x1*Fs/N Hz必为造影剂灌注信号的最佳频率，且该频率下血管特征点v对应的模值必定大于背景特征点b对应的模值；而其他特征频率往往对应着运动噪声或非最佳的血管信号或其两者的叠加，血管特征点v对应的模值和背景特征点b对应的模值差相应较小且大小关系不定。
大多数情况下，造影剂为浓度恒定的连续的一次灌注，因此在原始的造影图像序列中，血管特征点v的信号对应着阶跃信号，特征频率为1*Fs/N Hz。若造影剂灌注时浓度变化大或者时断时续，或者在某些局部的特殊的血管中有特殊的灌注特点，则造影剂灌注信号非阶跃信号，血管特征频率也可能是其他频率。
步骤S4：在取i＝1或2的情况下，在原始的造影图像序列中遍历造影区域的所有像素点，以0Hz和其中一个特征频率为基准对每个像素点的时序信号做傅里叶带通滤波处理，计算每个像素点的造影成像参量，即调制深度(MDi)。
其中，若原始的造影图像序列的图像中有许多杂散噪声，可以在遍历造影区域的所有像素点之前，即在计算每一点的造影成像参量之前，先对其每一帧图像做去噪处理，例如中值去噪、高斯去噪、维纳去噪。
以任一像素点s(x,y)为例，计算该点的造影成像参量。首先，取出该像素点的时序信号，以0Hz、x1*Fs/N Hz、x2*Fs/N Hz为基准对该时序信号做傅里叶带通滤波处理。带通滤波之后的信号为S0(x,y,t)、Sx1(x,y,t)、Sx2(x,y,t)。其中，S0(x,y,t)为直流信号，Sx1(x,y,t)、Sx2(x,y,t)为相应频率的交流信号。例如，在本实施例中，对血管特征点v和背景特征点b的时序信号，以0Hz、x1*Fs/N Hz、x2*Fs/N Hz为基准对该时序信号做傅里叶带通滤波处理。带通滤波之后的v点信号为Sv0(x,y,t)、Svx1(x,y,t)、Svx2(x,y,t)，其中，如图2a所示Sv0(x,y,t)为直流信号；如图2b、2c所示，相应频率的交流信号为Svx1(x,y,t)、Svx2(x,y,t)。带通滤波之后的b点信号为Sb0(x,y,t)、Sbx1(x,y,t)、Sbx2(x,y,t)，其中，如图2d所示，Sb0(x,y,t)为直流信号。如图2e、2f所示，相应频率的交流信号为Sbx1(x,y,t)、Sbx2(x,y,t)。交流信号对应着造影剂的灌注信号或背景组织的运动噪声信号。
随后，计算该像素点的输出造影成像参量为调制深度(MDi)，公式如下：
Mean(abs(Sxi(x,y,t)))为对应xi*Fs/N Hz频率的交流信号的绝对值的平均值。Mean(abs(S0(x,y,t)))为直流信号的绝对值的平均值。分母Mean(abs(S0(x,y,t)))的作用为对成像的不同区域做归一化处理。这是因为由于入射X射线的不均一性，以及造影图像不同区域的密度、厚度、元素分布的不同。不同区域的像素点有较大差异的像素值分布。导致不同区域的血管，有不同的造影效果。不利于整体的血管造影成像。因此，以直流信号作归一化，可去除该因素的影响。
此外，当步骤S3得到的特征频率的总个数不为2时，该步骤S4中的i可以取1到特征频率的总个数之间的任意一个数，而不局限于取i＝1或2。
步骤S5：输出以MD1、MD2为成像参量的造影图像。
如图3a、3b所示为在本实施例中以MD1、MD2为成像参量的造影图像。图3a中主要包含血管信息，对应着最大的模值差对应的特征频率即x1*Fs/NHz。这部分信号对应着血管的灌注信号，血管区域像素中该组分信号强度是其背景组织的10倍。图3b主要为运动的背景组织，对应着次最大的模值差对应的特征频率即x2*Fs/N Hz，该部分信号主要为呼吸运动。本例成功的区分了血管信号和运动噪声。
一般地，造影图像中血管区域的像素，由于造影剂的灌注，其时序信号理想情况下为阶跃信号，经过频域滤波后主要特征信号分布在最大的模值差所对应的特征频率x1*Fs/N Hz上，而相应的，造影图像中的背景区域其时序信号的主要特征分布在其他频率上，此频率和血管灌注频率有一定的区别。本发明中的方法可以很轻易的区分出造影剂信号和背景组织的运动噪声，从而达到抑制噪声的效果，提高成像的分辨率。因此，在这些造影图像中，一般性的MD1为血管造影图像，MD2为对应频率的高频运动组织图像，例如受呼吸影响的运动组织图像、受肌肉抽搐影响的运动组织图像等。这说明本发明中的方法可以很轻易的区分出传统数字减影技术无法区分开的造影剂信号和背景组织的运动噪声。作为对照，本例给出了传统的时间减影血管造影图像如图4所示，其中造影剂信号和背景组织的运动噪声远远不如本发明中的清晰，并且也无法区分开来。
此外，本发明还提供了一种基于X射线造影的动态阅片方法，包括：
步骤S1’：取i＝1，采用上文所述的X射线造影方法，得到特征频率x1*Fs/NHz，并得到以MD1为造影成像参量的血管造影图像，作为原始的造影图像序列的前N帧所对应的血管造影图像，并取j＝1；
步骤S2’：选取原始的造影图像序列的前N‐j帧作为当前的造影图像序列，并以步骤S1’得到的特征频率计算该当前的造影图像序列所对应的特征频率xi*Fs/(N‐j)Hz；特征频率由x1*Fs/N Hz相应地替换为x1*Fs/(N‐j)Hz，是因为对应造影剂信号的阶跃信号，在不同长度的时序信号中其特征频率也相应改变。
步骤S3’：重复上文所述的X射线造影方法的步骤S4和步骤S5，具体包括：在当前的造影图像序列中遍历造影区域的所有像素点，以0Hz和当前步骤S2’得到的特征频率x1*Fs/(N‐j)Hz为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度MD1；得到以调制深度MD1为造影成像参量的血管造影图像，作为原始的造影图像序列的前N‐j帧所对应的血管造影图像；
步骤S4’：当N‐j为n时，依次播放前n、n+1…、N‐1、N‐2、N帧所对应的造影图像；否则将j+1更新为j并回到步骤S2’。其中n为造影剂恰好灌注到造影区域的影像帧数。
这样就可实现造影过程的动态阅片播放，动态画面连续、清晰、血管轨迹稳定不抖动。
以上所述的，仅为本发明的较佳实施例，并非用以限定本发明的范围，本发明的上述实施例还可以做出各种变化。即凡是依据本发明申请的权利要求书及说明书内容所作的简单、等效变化与修饰，皆落入本发明专利的权利要求保护范围。本发明未详尽描述的均为常规技术内容。
Claims (10)
1.一种X射线造影方法，包括：
步骤S1：以采样频率Fs获取N帧的原始的造影图像序列，作为当前的造影图像序列；
步骤S2：在原始的造影图像序列的其中一帧上选取血管特征点和背景特征点，提取这两个特征点的时序信号并将其时序信号傅里叶变换为频域信号；
步骤S3：对血管特征点和背景特征点的频域信号进行特征分析，依次选取这两个特征点在相同频率下的最大的至少两个模值差，第i大的模值差所对应的频率作为特征频率xi*Fs/N Hz；
步骤S4：在当前的造影图像序列中遍历造影区域的所有像素点，以0Hz和其中一个特征频率xi*Fs/N Hz为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度MDi；
步骤S5：输出以调制深度MDi为造影成像参量的造影图像；
其中xi为正整数。
2.根据权利要求1所述的X射线造影方法，其特征在于，所述步骤S1包括：对造影对象灌注造影剂，在X光光源的照射下，以造影剂未灌注之前为起始帧至以造影剂灌注完毕为结束帧，采集原始的造影图像序列。
3.根据权利要求2所述的X射线造影方法，其特征在于，所述X光的能量高于所述造影剂的吸收边，光子通量为109—1010photons second‐1mm‐2。
4.根据权利要求2所述的X射线造影方法，其特征在于，所述造影剂为非离子型碘盐类造影剂，该非离子型碘盐类造影剂选自碘必乐、碘海醇、碘普罗胺、碘佛醇、碘克沙醇、泛影葡胺中的一种。
5.根据权利要求2所述的X射线造影方法，其特征在于，所述造影剂的用量为100—500微升，浓度为150—350mgI·ml‐1。
6.根据权利要求1所述的X射线造影方法，其特征在于，所述步骤S3中所选取的模值差为2‐4个，且所述特征频率的个数等于模值差的个数。
7.根据权利要求1所述的X射线造影方法，其特征在于，所述每个像素点的调制深度MDi通过下式计算：
其中，S0(x,y,t)、Sxi(x,y,t)分别是对每个像素点的时序信号做傅里叶带通滤波处理后获得的直流信号和对应xi*Fs/N Hz频率的交流信号；
Mean(abs(Sxi(x,y,t)))为对应xi*Fs/N Hz频率的交流信号的绝对值的平均值；
Mean(abs(S0(x,y,t)))为直流信号的绝对值的平均值。
8.根据权利要求1所述的X射线造影方法，其特征在于，所述步骤S4还包括：在遍历造影区域的所有像素点之前，先对其每一帧图像做去噪处理，该去噪处理包括中值去噪、高斯去噪或维纳去噪。
9.一种基于X射线造影的动态阅片方法，包括：
步骤S1’：采用权利要求1所述X射线造影方法，得到特征频率xi*Fs/N Hz，并得到以调制深度MDi为造影成像参量的造影图像，作为原始的造影图像序列的前N帧所对应的造影图像，取j＝1；
步骤S2’：选取原始的造影图像序列的前N‐j帧作为当前的造影图像序列，并以步骤S1’得到的特征频率计算所述当前的造影图像序列所对应的特征频率xi*Fs/(N‐j)Hz；
步骤S3’：在当前的造影图像序列中遍历造影区域的所有像素点，以0Hz和步骤S2’得到的特征频率xi*Fs/(N‐j)Hz为基准对每个像素点的时序信号做傅里叶带通滤波处理，并通过得到的信号计算每个像素点的调制深度MDi；得到以调制深度MDi为造影成像参量的造影图像，作为原始的造影图像序列的前N‐j帧所对应的造影图像；
步骤S4’：当N‐j为n时，依次播放所述前n、n+1…、N‐1、N‐2、N帧所对应的造影图像；否则将j+1更新为j并回到步骤S2’。
10.根据权利要求9所述的基于X射线造影的动态阅片方法，其特征在于，当i为1时，所述以调制深度MD1为造影成像参量的造影图像为血管造影图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810203889.1A CN108852385B (zh) | 2018-03-13 | 2018-03-13 | 一种x射线造影方法及基于x射线造影的动态阅片方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810203889.1A CN108852385B (zh) | 2018-03-13 | 2018-03-13 | 一种x射线造影方法及基于x射线造影的动态阅片方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108852385A true CN108852385A (zh) | 2018-11-23 |
| CN108852385B CN108852385B (zh) | 2022-03-04 |
Family
ID=64326063
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810203889.1A Active CN108852385B (zh) | 2018-03-13 | 2018-03-13 | 一种x射线造影方法及基于x射线造影的动态阅片方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108852385B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1141156A (zh) * | 1995-05-22 | 1997-01-29 | 株式会社岛津制作所 | 数字血管造影设备 |
| US6005917A (en) * | 1995-04-05 | 1999-12-21 | Andersson; Mats | Velocity adaptive filtered angiography |
| US20030013953A1 (en) * | 1998-04-10 | 2003-01-16 | Mistretta Charles A. | Time resolved computed tomography angiography |
| CN101052992A (zh) * | 2004-07-27 | 2007-10-10 | 杜尔牙科器械两合公司 | 用于提高射线照片上不同结构的感觉能力的方法和设备 |
| CN101236647A (zh) * | 2007-12-07 | 2008-08-06 | 华中科技大学 | 一种融合上下文信息的数字血管造影图像增强方法 |
| CN101999885A (zh) * | 2010-12-21 | 2011-04-06 | 中国人民解放军国防科学技术大学 | 一种自动分离动静脉血管的内源光学成像方法 |
| CN104517301A (zh) * | 2014-12-30 | 2015-04-15 | 华中科技大学 | 多参数模型指导的迭代提取血管造影图像运动参数的方法 |
-
2018
- 2018-03-13 CN CN201810203889.1A patent/CN108852385B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6005917A (en) * | 1995-04-05 | 1999-12-21 | Andersson; Mats | Velocity adaptive filtered angiography |
| CN1141156A (zh) * | 1995-05-22 | 1997-01-29 | 株式会社岛津制作所 | 数字血管造影设备 |
| US20030013953A1 (en) * | 1998-04-10 | 2003-01-16 | Mistretta Charles A. | Time resolved computed tomography angiography |
| CN101052992A (zh) * | 2004-07-27 | 2007-10-10 | 杜尔牙科器械两合公司 | 用于提高射线照片上不同结构的感觉能力的方法和设备 |
| CN101236647A (zh) * | 2007-12-07 | 2008-08-06 | 华中科技大学 | 一种融合上下文信息的数字血管造影图像增强方法 |
| CN101999885A (zh) * | 2010-12-21 | 2011-04-06 | 中国人民解放军国防科学技术大学 | 一种自动分离动静脉血管的内源光学成像方法 |
| CN104517301A (zh) * | 2014-12-30 | 2015-04-15 | 华中科技大学 | 多参数模型指导的迭代提取血管造影图像运动参数的方法 |
Non-Patent Citations (1)
| Title |
|---|
| TIMUR AKSOY AND GOZDE UNAL: "Template-based CTA to x-ray angio rigid registration of coronary arteries in frequency domain with automatic x-ray segmentation", 《MEDICAL PHYSICS》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108852385B (zh) | 2022-03-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Halliburton et al. | SCCT guidelines on radiation dose and dose-optimization strategies in cardiovascular CT | |
| US7840255B2 (en) | X-ray CT apparatus and myocardial perfusion image generating system | |
| JP5921132B2 (ja) | 医用画像処理システム | |
| US8498463B2 (en) | Mask construction for cardiac subtraction | |
| Robson et al. | Correction of respiratory and cardiac motion in cardiac PET/MR using MR-based motion modeling | |
| US20100067765A1 (en) | List Mode-Based Respiratory and Cardiac Gating in Positron Emission Tomography | |
| US20130261445A1 (en) | Representation of blood vessels and tissue in the heart | |
| CN106562797A (zh) | 单次曝光数字减影血管造影成像***与方法 | |
| CN114423348A (zh) | 降低灌注成像噪声和辐射或造影剂剂量的基于k空间方法 | |
| Ramirez‐Giraldo et al. | A strategy to decrease partial scan reconstruction artifacts in myocardial perfusion CT: Phantom and in vivo evaluation | |
| Koivumäki et al. | An integrated bioimpedance—ECG gating technique for respiratory and cardiac motion compensation in cardiac PET | |
| Yan et al. | Image quality of automatic coronary CT angiography reconstruction for patients with HR≥ 75 bpm using an AI-assisted 16-cm z-coverage CT scanner | |
| Liu et al. | Renal perfusion and hemodynamics: accurate in vivo determination at CT with a 10-fold decrease in radiation dose and HYPR noise reduction | |
| CN108852385A (zh) | 一种x射线造影方法及基于x射线造影的动态阅片方法 | |
| Paul et al. | Effect of contrast material on radiation dose in an adult cardiac dual-energy CT using retrospective ECG-gating | |
| Lehmann et al. | Angle‐independent measure of motion for image‐based gating in 3D coronary angiography | |
| Liu et al. | Effect of deep learning image reconstruction with high-definition standard scan mode on image quality of coronary stents and arteries | |
| Liang et al. | Dynamic chest image analysis: model-based perfusion analysis in dynamic pulmonary imaging | |
| US20240138795A1 (en) | Systems and methods for pulmonary perfusion analysis using dynamic radiography | |
| Luo et al. | Digital subtraction angiography image segmentation based on multiscale Hessian matrix applied to medical diagnosis and clinical nursing of coronary stenting patients | |
| Xiong et al. | Evaluation the effect of wide-body detector CT under free breathing combined with cardiac motion correction technology on CCTA image quality | |
| RU2644928C1 (ru) | Способ определения кадров, соответствующих границам фаз кровообращения, при проведении ангиографического исследования (варианты) | |
| Sharma | X-Ray Particle Image Velocimetry in 3D-Printed Phantoms Using High-Speed Angiography | |
| Abou-Issa et al. | Effect of low tube kV on radiation dose and image quality in retrospective ECG-gated coronary CT angiography | |
| Cai et al. | Dynamic cone beam CT angiography of carotid and cerebral arteries using canine model |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |