CN115100164A - 一种基于边缘检测的超声ct渡越时间自动提取方法 - Google Patents

一种基于边缘检测的超声ct渡越时间自动提取方法 Download PDF

Info

Publication number
CN115100164A
CN115100164A CN202210784649.1A CN202210784649A CN115100164A CN 115100164 A CN115100164 A CN 115100164A CN 202210784649 A CN202210784649 A CN 202210784649A CN 115100164 A CN115100164 A CN 115100164A
Authority
CN
China
Prior art keywords
channel
edge detection
transit time
transition time
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.)
Granted
Application number
CN202210784649.1A
Other languages
English (en)
Other versions
CN115100164B (zh
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.)
Wuhan Weishi Medical Science Image Co ltd
Huazhong University of Science and Technology
Original Assignee
Wuhan Weishi Medical Science Image Co ltd
Huazhong University of Science and Technology
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 Wuhan Weishi Medical Science Image Co ltd, Huazhong University of Science and Technology filed Critical Wuhan Weishi Medical Science Image Co ltd
Priority to CN202210784649.1A priority Critical patent/CN115100164B/zh
Publication of CN115100164A publication Critical patent/CN115100164A/zh
Application granted granted Critical
Publication of CN115100164B publication Critical patent/CN115100164B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • 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
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding
    • 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/10132Ultrasound image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明属于超声断层成像领域,具体公开了一种基于边缘检测的超声CT渡越时间自动提取方法,包括:(1)带通滤波;(2)边缘检测算子与像素卷积,得到通道信号集图像;(3)计算各通道二值化阈值;(4)二值化处理,得到二值化的通道信号集图像;(5)获取渡越时间,基于二值化的通道信号集图像,取各通道第一个非零值作为该通道的渡越时间。本发明通过对处理方法的整体流程设计、关键阈值的自适应取值规则等进行改进,基于图像边缘检测理论,利用各通道采用自适应阈值,并可优选配合设定检测范围,以及对渡越时间差异矩阵进行滤波处理,能够避免人为设定阈值参数,弱化噪声的影响,解决低信噪比信号渡越时间提取不准确的技术问题。

Description

一种基于边缘检测的超声CT渡越时间自动提取方法
技术领域
本发明属于超声断层成像领域,更具体地,涉及一种基于边缘检测的超声CT渡越时间自动提取方法,该方法能够对超声断层声速成像模式中透射信号的渡越时间自动提取。
背景技术
超声断层成像(超声CT)是指通过超声探头对物体发射超声波并接收反射数据或者透射数据,利用这些数据重建出超声断层图像,以便观测物体内部的三维信息。超声检测具有价格低廉、对人体无害等优点,随着探头加工工艺和计算机高性能运算的快速发展,超声断层成像技术近些年来又再次成为了医学成像研究热点。
超声CT成像有两种成像模式,反射式成像和透射式成像。由于采集多方位的反射信息,超声CT的反射图像具有更高的图像分辨率,以便辅助医生看到更微小的病变组织。而通过透射数据可以重建出声速、衰减系数等功能参数,属于功能像的领域。有研究指出,在病变早期,病变组织先有功能参数的变化,后有结构变化。因此,透射式成像对病变的早期成像和诊断具有重要意义,能更早地诊断出病变。
超声CT声速成像根据从透射波提取的渡越时间求解声速参数,因而渡越时间是成像过程中重要的输入参数。渡越时间是信号到达的时间(即信号的起跳时间),当前渡越时间提取主要难点在于信号信噪比低,难以确定信号的到达时间。现有的渡越时间提取方法有幅度法、赤池信息标准法、互相关法等。赤池信息标准法易受开关切换噪声扰动。幅度法实施过程中的幅度需要人为设定,因而存在不确定性。互相关法由于实际人体的非线性效应而存在误差。
边缘检测算法作为一种图像处理算法,现有的图像边缘检测算法采用单一的全局阈值进行二值化,未考虑各通道信号强度存在的差异,易检测到沿超声探头内壁传播的信号起跳时间。
发明内容
针对现有技术的以上缺陷或改进需求,本发明的目的在于提供一种基于边缘检测的超声CT渡越时间自动提取方法,其中通过对处理方法的整体流程设计、关键阈值的自适应取值规则等进行改进,基于图像边缘检测理论,利用各通道采用自适应阈值,并可优选配合设定检测范围,以及对渡越时间差异矩阵进行滤波处理,能够避免人为设定阈值参数,弱化噪声的影响,由此解决低信噪比信号渡越时间提取不准确的技术问题。
为实现上述目的,按照本发明的一个方面,提供了一种基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,包括以下步骤:
(1)带通滤波:根据探头的截止频率,确定带宽范围,并对采集到的超声断层透射信号进行带通滤波处理;
(2)边缘检测算子与像素卷积:根据步骤(1)得到的带通滤波处理后的信号,绘制单个阵元发射、其他接收阵元接收的初始通道信号集图像;该初始通道信号集图像是以接收阵元序号和采样时间分别作为坐标轴,图中各像素点的灰度值则对应相应接收到的信号幅度;
接着,采用边缘检测算子与初始通道信号集图像进行卷积,得到边缘检测后的通道信号集图像;
(3)计算各通道二值化阈值:基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以预先选定的阈值T0为界,将该通道对应的像素点划分成两部分,并就这两部分中的每一部分分别计算它们的平均灰度值;记这两部分得到的平均灰度值分别为A_mean和B_mean,计算T1=(A_mean+B_mean)/2,若abs(T1-T0)<epsilon,则将T1作为该通道的二值化分割阈值;否则,将T1赋值给T0,并重复该步骤;其中,epsilon为取值预先设定的正数;
对各个通道进行处理,即可得到各通道的二值化分割阈值;
(4)二值化处理:基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以步骤(3)得到的该通道的二值化分割阈值为界,将该通道对应的像素点划分成两部分,其中,对于像素点灰度值大于或大于等于阈值的这一部分,将这些像素点灰度值设置为预先选定的非零值;对于像素点灰度值小于等于或小于阈值的另一部分,将这些像素点灰度值设置为0;
对各个通道进行处理,即可得到二值化的通道信号集图像;
(5)获取渡越时间:基于步骤(4)得到的二值化的通道信号集图像,对各通道像素灰度进行检索,结合各采样点对应的采样时间先后,取各通道第一个非零值作为该通道的渡越时间;如此即可得到各通道的渡越时间,完成超声透射信号的渡越时间自动提取;所述渡越时间能够用于声速图像重建。
作为本发明的进一步优选,该方法在步骤(3)开始前,还包括:
计算各通道二值化阈值检测范围:基于步骤(2)得到的边缘检测后的通道信号集图像,根据当前待处理数据采集过程中的温度,获取该温度下的耦合剂声速;接着,利用探头内相应发射阵元与接收阵元的位置距离,结合所述耦合剂声速,计算各通道的理论渡越时间;然后,基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以该通道对应的理论渡越时间前后各m个采样点所涵盖的范围,作为该通道二值化阈值的检测范围,并将该范围外的各像素点灰度值设置为0;其中,m为预先选定的正整数;
如此,即可得到各通道的二值化阈值检测范围;
相应的,所述步骤(3)中,所述预先选定的阈值T0即为相应通道检测范围内对应的各像素点的灰度平均值。
作为本发明的进一步优选,m为50到200。
作为本发明的进一步优选,所述步骤(5)中,在得到各通道的渡越时间之后,还包括步骤:
(S1)获取渡越时间矩阵:
基于所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为待测信号的渡越时间矩阵;
(S2)获取渡越时间差异矩阵:
采集纯耦合剂条件下超声断层透射信号,记为参考信号;然后基于该参考信号,重复步骤(1)至步骤(5),并基于参考信号条件下所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为参考信号的渡越时间矩阵;
将所述步骤(S1)得到的待测信号的渡越时间矩阵与所述参考信号的渡越时间矩阵相减后,即可得到渡越时间差异矩阵;
(S3)渡越时间差异矩阵滤波:
建立a×b大小的模板,使该模板的中心遍历渡越时间差异矩阵中的各个元素,在遍历过程中,以模板大小范围内各元素的中值或均值作为相应模板中心元素对应的更新值,并在遍历结束后完成对渡越时间差异矩阵的更新值赋值,即可得到滤波后的渡越时间差异矩阵;其中,a、b均为预先设定的正整数;
所述滤波后的渡越时间差异矩阵中的各个元素,即对应相应发射阵元、相应接收阵元的待测信号与参考信号的渡越时间差值。
作为本发明的进一步优选,所述步骤(S3)中,所述重复遍历若干次具体是重复遍历1~5次。
作为本发明的进一步优选,所述探头具体为环形探头;
所述步骤(1)中,采集超声断层透射信号,具体是以环形探头上排除以发射阵元为中心的、不超过1/4圆周范围的其他阵元作为接收阵元,所对应得到的超声断层透射信号。
按照本发明的另一方面,本发明提供了一种基于边缘检测的超声CT渡越时间自动提取系统,其特征在于,包括:
带通滤波功能模块,用于:根据探头的截止频率确定的带宽范围,对采集到的超声断层透射信号进行带通滤波处理;
边缘检测算子与像素卷积功能模块,用于:根据通过带通滤波功能模块得到的带通滤波处理后的信号,绘制单个阵元发射、其他接收阵元接收的初始通道信号集图像,并采用边缘检测算子与初始通道信号集图像进行卷积,得到边缘检测后的通道信号集图像;其中,所述初始通道信号集图像是以接收阵元序号和采样时间分别作为坐标轴,图中各像素点的灰度值则对应相应接收到的信号幅度;
各通道二值化阈值计算功能模块,用于:基于通过边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以预先选定的阈值T0为界,将该通道对应的像素点划分成两部分,并就这两部分中的每一部分分别计算它们的平均灰度值;记这两部分得到的平均灰度值分别为A_mean和B_mean,计算T1=(A_mean+B_mean)/2,若abs(T1-T0)<epsilon,则将T1作为该通道的二值化分割阈值;否则,将T1赋值给T0,并重复该步骤;其中,epsilon为取值预先设定的正数;对各个通道进行处理,即可得到各通道的二值化分割阈值;
二值化处理功能模块,用于:基于边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以通过各通道二值化阈值计算功能模块得到的该通道的二值化分割阈值为界,将该通道对应的像素点划分成两部分,其中,对于像素点灰度值大于或大于等于阈值的这一部分,将这些像素点灰度值设置为预先选定的非零值;对于像素点灰度值小于等于或小于阈值的另一部分,将这些像素点灰度值设置为0;对各个通道进行处理,即可得到二值化的通道信号集图像;
获取渡越时间功能模块,用于:基于通过二值化处理功能模块得到的二值化的通道信号集图像,对各通道像素灰度进行检索,结合各采样点对应的采样时间先后,取各通道第一个非零值作为该通道的渡越时间;如此即可得到各通道的渡越时间,完成超声透射信号的渡越时间自动提取;所述渡越时间能够用于声速图像重建。
作为本发明的进一步优选,在所述边缘检测算子与像素卷积功能模块与所述各通道二值化阈值计算功能模块之间,还包括:
各通道二值化阈值检测范围计算功能模块,用于:基于通过边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,根据当前待处理数据采集过程中的温度,获取该温度下的耦合剂声速;接着,利用探头内相应发射阵元与接收阵元的位置距离,结合所述耦合剂声速,计算各通道的理论渡越时间;然后,基于边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以该通道对应的理论渡越时间前后各m个采样点所涵盖的范围,作为该通道二值化阈值的检测范围,并将该范围外的各像素点灰度值设置为0;其中,m为预先选定的正整数;如此,即可得到各通道的二值化阈值检测范围;
相应的,对于所述各通道二值化阈值计算功能模块,所述预先选定的阈值T0即为相应通道检测范围内对应的各像素点的灰度平均值。
作为本发明的进一步优选,所述获取渡越时间功能模块之后,还包括以下子功能模块:
获取渡越时间矩阵子功能模块,用于:基于所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为待测信号的渡越时间矩阵;
获取渡越时间差异矩阵子功能模块,用于:采集纯耦合剂条件下超声断层透射信号,记为参考信号;然后基于该参考信号,利用所述带通滤波功能模块、所述边缘检测算子与像素卷积功能模块、所述各通道二值化阈值计算功能模块、所述二值化处理功能模块和所述获取渡越时间功能模块,并基于参考信号条件下所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为参考信号的渡越时间矩阵;将获取渡越时间矩阵子功能模块得到的待测信号的渡越时间矩阵与所述参考信号的渡越时间矩阵相减后,即可得到渡越时间差异矩阵;
渡越时间差异矩阵滤波子功能模块,用于:建立a×b大小的模板,使该模板的中心遍历渡越时间差异矩阵中的各个元素,在遍历过程中,以模板大小范围内各元素的中值或均值作为相应模板中心元素对应的更新值,并在遍历结束后完成对渡越时间差异矩阵的更新值赋值,即可得到滤波后的渡越时间差异矩阵;其中,a、b均为预先设定的正整数;
所述滤波后的渡越时间差异矩阵中的各个元素,即对应相应发射阵元、相应接收阵元的待测信号与参考信号的渡越时间差值。
通过本发明所构思的以上技术方案,与现有技术相比,首先通过边缘检测算子对滤波后的通道信号集图像进行卷积处理,通过迭代法计算各通道内的阈值实现边缘检测,比传统的单一阈值边缘检测方法能更好地提取到低信噪比通道的渡越时间;并且,基于本发明还可进一步配合计算参考信号的理论渡越时间作为各通道阈值检测范围,以及通过对渡越时间矩阵滤波来消除异常值,能够进一步弱化噪声的影响,得到的渡越时间更加准确。
与现有边缘检测算法相似,本发明通过边缘检测模板与像素的卷积能弱化噪声的影响,进一步的,与现有边缘检测算法不同的是,本发明不需要人为设定阈值。本发明基于图像的边缘检测理论,在通道信号集图像的边缘检测过程中引入了通道自适应阈值,考虑了各通道幅度差异,在各通道计算阈值相比于全局阈值能更准确地提取渡越时间,实现了超声透射信号的渡越时间自动提取。
而进一步的,利用设定检测范围,以及对渡越时间差异矩阵进行滤波处理等优选处理,能够进一步提升处理效率及准确性。
综上,本发明基于边缘检测的超声断层成像中渡越时间自动提取方法,避免了人为设定阈值参数,减弱了噪声的影响,实现了低信噪比信号的渡越时间自动提取。
附图说明
图1是实际环形探头超声CT系统的一例来自乳腺癌患者的通道信号集图像。
图2是图1中数据经sobel算子边缘检测后。
图3是在图2之后采用迭代法实现的单一阈值sobel算子边缘检测图像。
图4是图3中迭代法单一阈值sobel算子边缘检测后检测到的渡越时间。
图5是在图2之后采用迭代法实现的自适应通道内阈值sobel算子边缘检测图像。
图6是图5中迭代法自适应通道内阈值sobel算子边缘检测后检测到的渡越时间。
图7是图1中数据限定检测范围后采用迭代法自适应通道内阈值sobel算子边缘检测图像。
图8是图7中限定检测范围和采用迭代法自适应通道内阈值sobel算子边缘检测后检测到的渡越时间。
图9是本发明实施例方法提取到的乳腺癌患者数据的渡越时间差异矩阵(矩阵中的各个元素,即,相应发射阵元、相应接收阵元的待测信号与参考信号的渡越时间差值,已进行了灰度图显示)。
图10是对图1超声CT数据采用本发明实施例方法提取渡越时间并重建的声速图。
图11是本发明实施例方法所对应的框图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。此外,下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明基于边缘检测的超声CT渡越时间自动提取方法,在具体实现时,可以包括以下步骤:
1)带通滤波
首先对系统采集到的超声断层信号进行带通滤波处理。带宽范围为频率f1~f2(f1和f2为上下截止频率,为探头出厂参数),设计通带范围为f1~f2的带通滤波器对信号进行滤波。带通滤波器可选取FIR滤波器、IIR滤波器等。
2)边缘检测算子与像素卷积
基于步骤1)得到的带通滤波处理后的信号,绘制单个阵元发射而其他阵元接收的通道信号集图像,作为待处理的图像。与现有技术相似,具体绘制方法可以是:x轴为接收通道序号,y轴为采样时间,不同信号幅度体现为亮度的深浅。采用边缘检测算子与通道信号集图像像素进行卷积,得到边缘检测后的灰度图像。此灰度图像并非二值图像。边缘检测算子可选取Sobel算子、Roberts算子、Prewitt算子等。
3)计算各通道二值化阈值检测范围
根据当前待处理数据采集过程中的温度,通过查表获取耦合剂的声速。耦合剂一般为水或其他与人体参数相近的材料。根据探头出厂时提供的各阵元位置和耦合剂声速,计算当前阵元发射而其他阵元接收信号时各通道的理论渡越时间。与现有技术相同,理论渡越时间计算方法具体为:探头出厂的阵元距离/水的声速,其中,水的声速可采用声速仪进行测量。该理论渡越时间前后各m个采样点范围内作为各通道二值化阈值的检测范围。m为自然数,其取值可预先设定,例如可选取为50到200。m个采样点的时间=(m-1)/采样频率。
4)计算各通道二值化阈值
(a)对于某一通道(通道对应接收阵元序号),选择一个初始阈值T0,例如可取该通道内所有像素的均值;
(b)根据阈值T0将该通道分为A和B两部分。分别求出A和B的平均灰度值A_mean和B_mean;
(c)计算T1=(A_mean+B_mean)/2,比较abs(T1-T0)<epsilon(epsilon为一个较小的正数,具体取值可预先设定;abs表示绝对值函数,下同)是否成立。若成立,则停止迭代,此时的T1就是该通道内分割阈值;否则,将T1赋值给T0(T0=T1,此时T0为处在更新过程的阈值),继续从第(b)步开始。
5)二值化
根据第4)步计算的各通道阈值Tx,对各通道检测范围内的像素进行二值化,即小于Tx的像素灰度置为0,大于等于Tx的像素灰度置为1。对各通道检测范围外的像素值置0。
6)获取渡越时间
通过对各通道像素灰度进行检索,取各通道第一个非零值(与步骤5)的二值化规则相对应)作为该通道的渡越时间。
7)获取渡越时间矩阵
对所有发射阵元进行第1)到第6)步骤的处理,得到所有发射阵元的渡越时间矩阵。
8)获取渡越时间差异矩阵
对参考信号(无人体,纯耦合剂状态下采集的信号)和人体信号(将人体器官置于耦合剂当中,进行信号采集)进行第1)到第7)步骤的处理,并将人体信号的渡越时间矩阵减去参考信号的渡越时间矩阵,得到渡越时间差异矩阵。
9)渡越时间差异矩阵滤波
对渡越时间差异矩阵进行滤波处理。处理方式为:
(a)计算某一矩阵元素3×3范围内的渡越时间差异的中值(或均值)。模板大小可不限定于3×3。
(b)若该元素的渡越时间差异数值与该中值(或均值)的差值大于该中值(或均值)的两倍以上,则将该元素的渡越时间差异数值置为该中值(或均值)。
(c)对各矩阵元素逐一检索,进行步骤(a)(b)的处理。
(d)对渡越时间差异矩阵进行n次步骤(a)(b)(c)的处理,n可为1到5。
10)图像重建
可进一步基于步骤9)得到的渡越时间差异矩阵,对其进行声速图像重建,从而得到声速图,完成超声断层成像输出。
以下为具体实施例:
实施例1
本实施例对于环形探头采集到的人体乳腺信号进行了实例演示。
该系统环形探头由2048阵元组成,阵元中心频率为1.5-4MHz。如图1所示,为1号阵元发射,其他阵元接收的通道信号集图像绘制。这里去除了发射阵元左右各256阵元。基于本发明,超声断层成像中的透射波渡越时间自动提取方法具体处理过程如下:
首先对系统采集到的超声断层信号进行带通滤波处理。对信号进行频谱分析,得到主要的带宽范围为1.5-4MHz,如图2所示,为经过1.5-4MHz通带范围的带通滤波后的通道信号集图像,采用的滤波器为FIR滤波器。
绘制单个阵元发射而其他阵元接收的通道信号集图像,作为待处理的图像。采用边缘检测Sobel算子与像素进行卷积,得到边缘检测后的灰度图像。图3为迭代法单一阈值的边缘检测效果,使用的边缘检测算子为Sobel算子(所基于的迭代法单一阈值的边缘检测,指的是迭代计算阈值是在整幅图像进行的,即阈值T0将整幅图像分为A和B两部分,而非针对单一通道分别进行迭代计算阈值;针对单一通道分别计算阈值,即阈值T0将单个通道分为A和B两部分)。图4为对应图3检测到的渡越时间。
(a)对于某一通道(通道对应接收阵元序号),选择一个初始阈值T0,如取该通道内所有像素的均值作为T0;(b)根据阈值T0将该通道分为A和B两部分。分别求出A和B的平均灰度值A_mean和B_mean;(c)计算T1=(A_mean+B_mean)/2,比较abs(T1-T0)<epsilon(epsilon为一个较小值,本实施例中epsilon取10-5)是否成立。若成立,则停止迭代,此时的T1就是该通道内分割阈值;否则,将T1赋值给T0(T0=T1,此时T0为处在更新过程的阈值),继续从第(b)步开始。图5为迭代法通道内自适应阈值的边缘检测效果,图6为对应图5检测到的渡越时间。通过对比图6与图4,可以发现部分图4中检测到异常值(3750)的位置,检测到正常的数值。
根据当前待处理数据采集过程中的温度,通过查表获取耦合剂(此实例中为水)的声速。根据探头出厂时提供的各阵元位置和水的声速,计算当前阵元发射而其他阵元接收信号时各通道的理论渡越时间。该理论渡越时间前后各150个采样点范围内作为各通道二值化阈值的检测范围。图7为对阈值检测范围进行限定后使用迭代法通道内自适应阈值的边缘检测效果。通过对各通道像素灰度进行检索,取各通道第一个非零值作为该通道的渡越时间。图8为对应图7检测到的渡越时间。图8即为实施步骤1)到6)的效果。可以发现,图8相对于图6,进一步减少了异常值(3750)的数量。图8相比于图4,可见,通过进一步依靠限定检测范围和通道内自适应阈值,逐步减少了渡越时间异常值的数量。
对所有发射阵元进行上述步骤的处理,得到所有发射阵元的渡越时间矩阵。对参考信号(无人体,纯耦合剂状态下采集的信号)和人体信号(将人体器官置于耦合剂当中,进行信号采集)进行上述步骤的处理,并将人体信号的渡越时间矩阵减去参考信号的渡越时间矩阵,得到渡越时间差异矩阵。对渡越时间差异矩阵进行滤波处理。处理方式为:
(a)计算某一矩阵元素3×3范围内的渡越时间差异的中值(或均值)。
(b)若该元素的渡越时间差异数值与该中值(或均值)的差值大于该中值(或均值)的两倍以上,则将该元素的渡越时间差异数值置为该中值(或均值)。
(c)对各矩阵元素逐一检索,进行步骤(a)(b)的处理。
(d)对渡越时间差异矩阵进行2次步骤(a)(b)(c)的处理。
得到渡越时间差异矩阵如图9所示。对图9进行声速图像重建,得到声速图如图10所示。
可见,基于本发明方法,首先通过边缘检测算子对滤波后的通道信号集图像进行卷积处理,接着计算参考信号的理论渡越时间作为各通道阈值检测范围,然后再通过迭代法计算各通道内的阈值实现渡越时间检测,最后通过对渡越时间差异矩阵滤波来消除异常值,效果最优。当然,考虑到图6所对应的效果仍然优于图4,也可以只采用迭代法自适应通道内阈值这一改进,仍能取得优于现有技术的处理效果。
本领域的技术人员容易理解,以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,包括以下步骤:
(1)带通滤波:根据探头的截止频率,确定带宽范围,并对采集到的超声断层透射信号进行带通滤波处理;
(2)边缘检测算子与像素卷积:根据步骤(1)得到的带通滤波处理后的信号,绘制单个阵元发射、其他接收阵元接收的初始通道信号集图像;该初始通道信号集图像是以接收阵元序号和采样时间分别作为坐标轴,图中各像素点的灰度值则对应相应接收到的信号幅度;
接着,采用边缘检测算子与初始通道信号集图像进行卷积,得到边缘检测后的通道信号集图像;
(3)计算各通道二值化阈值:基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以预先选定的阈值T0为界,将该通道对应的像素点划分成两部分,并就这两部分中的每一部分分别计算它们的平均灰度值;记这两部分得到的平均灰度值分别为A_mean和B_mean,计算T1=(A_mean+B_mean)/2,若abs(T1-T0)<epsilon,则将T1作为该通道的二值化分割阈值;否则,将T1赋值给T0,并重复该步骤;其中,epsilon为取值预先设定的正数;
对各个通道进行处理,即可得到各通道的二值化分割阈值;
(4)二值化处理:基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以步骤(3)得到的该通道的二值化分割阈值为界,将该通道对应的像素点划分成两部分,其中,对于像素点灰度值大于或大于等于阈值的这一部分,将这些像素点灰度值设置为预先选定的非零值;对于像素点灰度值小于等于或小于阈值的另一部分,将这些像素点灰度值设置为0;
对各个通道进行处理,即可得到二值化的通道信号集图像;
(5)获取渡越时间:基于步骤(4)得到的二值化的通道信号集图像,对各通道像素灰度进行检索,结合各采样点对应的采样时间先后,取各通道第一个非零值作为该通道的渡越时间;如此即可得到各通道的渡越时间,完成超声透射信号的渡越时间自动提取;所述渡越时间能够用于声速图像重建。
2.如权利要求1所述基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,该方法在步骤(3)开始前,还包括:
计算各通道二值化阈值检测范围:基于步骤(2)得到的边缘检测后的通道信号集图像,根据当前待处理数据采集过程中的温度,获取该温度下的耦合剂声速;接着,利用探头内相应发射阵元与接收阵元的位置距离,结合所述耦合剂声速,计算各通道的理论渡越时间;然后,基于步骤(2)得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以该通道对应的理论渡越时间前后各m个采样点所涵盖的范围,作为该通道二值化阈值的检测范围,并将该范围外的各像素点灰度值设置为0;其中,m为预先选定的正整数;
如此,即可得到各通道的二值化阈值检测范围;
相应的,所述步骤(3)中,所述预先选定的阈值T0即为相应通道检测范围内对应的各像素点的灰度平均值。
3.如权利要求2所述基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,m为50到200。
4.如权利要求1所述基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,所述步骤(5)中,在得到各通道的渡越时间之后,还包括步骤:
(S1)获取渡越时间矩阵:
基于所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为待测信号的渡越时间矩阵;
(S2)获取渡越时间差异矩阵:
采集纯耦合剂条件下超声断层透射信号,记为参考信号;然后基于该参考信号,重复步骤(1)至步骤(5),并基于参考信号条件下所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为参考信号的渡越时间矩阵;
将所述步骤(S1)得到的待测信号的渡越时间矩阵与所述参考信号的渡越时间矩阵相减后,即可得到渡越时间差异矩阵;
(S3)渡越时间差异矩阵滤波:
建立a×b大小的模板,使该模板的中心遍历渡越时间差异矩阵中的各个元素,在遍历过程中,以模板大小范围内各元素的中值或均值作为相应模板中心元素对应的更新值,并在遍历结束后完成对渡越时间差异矩阵的更新值赋值,即可得到滤波后的渡越时间差异矩阵;其中,a、b均为预先设定的正整数;
所述滤波后的渡越时间差异矩阵中的各个元素,即对应相应发射阵元、相应接收阵元的待测信号与参考信号的渡越时间差值。
5.如权利要求1所述基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,所述步骤(S3)中,所述重复遍历若干次具体是重复遍历1~5次。
6.如权利要求1所述基于边缘检测的超声CT渡越时间自动提取方法,其特征在于,所述探头具体为环形探头;
所述步骤(1)中,采集超声断层透射信号,具体是以环形探头上排除以发射阵元为中心的、不超过1/4圆周范围的其他阵元作为接收阵元,所对应得到的超声断层透射信号。
7.一种基于边缘检测的超声CT渡越时间自动提取系统,其特征在于,包括:
带通滤波功能模块,用于:根据探头的截止频率确定的带宽范围,对采集到的超声断层透射信号进行带通滤波处理;
边缘检测算子与像素卷积功能模块,用于:根据通过带通滤波功能模块得到的带通滤波处理后的信号,绘制单个阵元发射、其他接收阵元接收的初始通道信号集图像,并采用边缘检测算子与初始通道信号集图像进行卷积,得到边缘检测后的通道信号集图像;其中,所述初始通道信号集图像是以接收阵元序号和采样时间分别作为坐标轴,图中各像素点的灰度值则对应相应接收到的信号幅度;
各通道二值化阈值计算功能模块,用于:基于通过边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以预先选定的阈值T0为界,将该通道对应的像素点划分成两部分,并就这两部分中的每一部分分别计算它们的平均灰度值;记这两部分得到的平均灰度值分别为A_mean和B_mean,计算T1=(A_mean+B_mean)/2,若abs(T1-T0)<epsilon,则将T1作为该通道的二值化分割阈值;否则,将T1赋值给T0,并重复该步骤;其中,epsilon为取值预先设定的正数;对各个通道进行处理,即可得到各通道的二值化分割阈值;
二值化处理功能模块,用于:基于边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以通过各通道二值化阈值计算功能模块得到的该通道的二值化分割阈值为界,将该通道对应的像素点划分成两部分,其中,对于像素点灰度值大于或大于等于阈值的这一部分,将这些像素点灰度值设置为预先选定的非零值;对于像素点灰度值小于等于或小于阈值的另一部分,将这些像素点灰度值设置为0;对各个通道进行处理,即可得到二值化的通道信号集图像;
获取渡越时间功能模块,用于:基于通过二值化处理功能模块得到的二值化的通道信号集图像,对各通道像素灰度进行检索,结合各采样点对应的采样时间先后,取各通道第一个非零值作为该通道的渡越时间;如此即可得到各通道的渡越时间,完成超声透射信号的渡越时间自动提取;所述渡越时间能够用于声速图像重建。
8.如权利要求7所述基于边缘检测的超声CT渡越时间自动提取系统,其特征在于,在所述边缘检测算子与像素卷积功能模块与所述各通道二值化阈值计算功能模块之间,还包括:
各通道二值化阈值检测范围计算功能模块,用于:基于通过边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,根据当前待处理数据采集过程中的温度,获取该温度下的耦合剂声速;接着,利用探头内相应发射阵元与接收阵元的位置距离,结合所述耦合剂声速,计算各通道的理论渡越时间;然后,基于边缘检测算子与像素卷积功能模块得到的边缘检测后的通道信号集图像,针对某一接收阵元通道,以该通道对应的理论渡越时间前后各m个采样点所涵盖的范围,作为该通道二值化阈值的检测范围,并将该范围外的各像素点灰度值设置为0;其中,m为预先选定的正整数;如此,即可得到各通道的二值化阈值检测范围;
相应的,对于所述各通道二值化阈值计算功能模块,所述预先选定的阈值T0即为相应通道检测范围内对应的各像素点的灰度平均值。
9.如权利要求7所述基于边缘检测的超声CT渡越时间自动提取系统,其特征在于,所述获取渡越时间功能模块之后,还包括以下子功能模块:
获取渡越时间矩阵子功能模块,用于:基于所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为待测信号的渡越时间矩阵;
获取渡越时间差异矩阵子功能模块,用于:采集纯耦合剂条件下超声断层透射信号,记为参考信号;然后基于该参考信号,利用所述带通滤波功能模块、所述边缘检测算子与像素卷积功能模块、所述各通道二值化阈值计算功能模块、所述二值化处理功能模块和所述获取渡越时间功能模块,并基于参考信号条件下所有发射阵元的各通道渡越时间,得到所有发射阵元的渡越时间矩阵,作为参考信号的渡越时间矩阵;将获取渡越时间矩阵子功能模块得到的待测信号的渡越时间矩阵与所述参考信号的渡越时间矩阵相减后,即可得到渡越时间差异矩阵;
渡越时间差异矩阵滤波子功能模块,用于:建立a×b大小的模板,使该模板的中心遍历渡越时间差异矩阵中的各个元素,在遍历过程中,以模板大小范围内各元素的中值或均值作为相应模板中心元素对应的更新值,并在遍历结束后完成对渡越时间差异矩阵的更新值赋值,即可得到滤波后的渡越时间差异矩阵;其中,a、b均为预先设定的正整数;
所述滤波后的渡越时间差异矩阵中的各个元素,即对应相应发射阵元、相应接收阵元的待测信号与参考信号的渡越时间差值。
CN202210784649.1A 2022-06-29 2022-06-29 一种基于边缘检测的超声ct渡越时间自动提取方法 Active CN115100164B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210784649.1A CN115100164B (zh) 2022-06-29 2022-06-29 一种基于边缘检测的超声ct渡越时间自动提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210784649.1A CN115100164B (zh) 2022-06-29 2022-06-29 一种基于边缘检测的超声ct渡越时间自动提取方法

Publications (2)

Publication Number Publication Date
CN115100164A true CN115100164A (zh) 2022-09-23
CN115100164B CN115100164B (zh) 2024-09-13

Family

ID=83297273

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210784649.1A Active CN115100164B (zh) 2022-06-29 2022-06-29 一种基于边缘检测的超声ct渡越时间自动提取方法

Country Status (1)

Country Link
CN (1) CN115100164B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115670504A (zh) * 2022-10-24 2023-02-03 浙江衡玖医疗器械有限责任公司 一种三维超声断层成像系统原始信号质量判断方法
CN118379315A (zh) * 2024-04-23 2024-07-23 盐城工学院 一种基于FPGA的8方向Sobel边缘检测系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170115423A1 (en) * 2015-10-26 2017-04-27 Schlumberger Technology Corporation Downhole Caliper Using Multiple Acoustic Transducers
CN110051387A (zh) * 2019-04-11 2019-07-26 华中科技大学 一种基于射线理论的超声ct图像重建方法及系统
CN110575202A (zh) * 2019-08-27 2019-12-17 华中科技大学 一种基于费马原理的超声ct图像重建方法及系统
CN110772281A (zh) * 2019-10-23 2020-02-11 哈尔滨工业大学(深圳) 基于改进射线追踪法的超声ct成像系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170115423A1 (en) * 2015-10-26 2017-04-27 Schlumberger Technology Corporation Downhole Caliper Using Multiple Acoustic Transducers
CN110051387A (zh) * 2019-04-11 2019-07-26 华中科技大学 一种基于射线理论的超声ct图像重建方法及系统
CN110575202A (zh) * 2019-08-27 2019-12-17 华中科技大学 一种基于费马原理的超声ct图像重建方法及系统
CN110772281A (zh) * 2019-10-23 2020-02-11 哈尔滨工业大学(深圳) 基于改进射线追踪法的超声ct成像系统

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115670504A (zh) * 2022-10-24 2023-02-03 浙江衡玖医疗器械有限责任公司 一种三维超声断层成像系统原始信号质量判断方法
CN115670504B (zh) * 2022-10-24 2024-01-09 浙江衡玖医疗器械有限责任公司 一种三维超声断层成像系统原始信号质量判断方法
CN118379315A (zh) * 2024-04-23 2024-07-23 盐城工学院 一种基于FPGA的8方向Sobel边缘检测系统

Also Published As

Publication number Publication date
CN115100164B (zh) 2024-09-13

Similar Documents

Publication Publication Date Title
CN115100164B (zh) 一种基于边缘检测的超声ct渡越时间自动提取方法
CN110325119B (zh) 卵巢卵泡计数和大小确定
US8094893B2 (en) Segmentation tool for identifying flow regions in an image system
US20060184021A1 (en) Method of improving the quality of a three-dimensional ultrasound doppler image
EP2016905B1 (en) Ultrasound diagnostic apparatus
US20210236095A1 (en) Ultrasound ct image reconstruction method and system based on ray theory
US20200015786A1 (en) Method for generating an enhanced image of a volume of tissue
CN111035411B (zh) 一种基于螺旋扫描的超声断层三维成像方法及系统
CN117056642B (zh) 一种微波数据处理方法及系统
WO1995020357A1 (en) Method for reducing speckle in vessel imaging
US20110002518A1 (en) Method and system for processing ultrasound data
CN105631879A (zh) 一种基于线型阵列的超声层析成像系统及方法
CN110706236B (zh) 血管图像的三维重建方法及装置
Gemmeke et al. Hardware setup for the next generation of 3D ultrasound computer tomography
CN114972567B (zh) 一种基于波动方程的医学超声ct多参数图像重建方法
CN111354006A (zh) 超声图像中目标组织的描迹方法及装置
US20220370037A1 (en) Vascular tissue characterization devices, systems, and methods
JP2000126178A (ja) 立体表面形状定量化方法、及びこれを応用した悪性腫瘍自動識別方法
JP7349202B2 (ja) 超音波信号処理方法、装置、機器及び記憶媒体
JP2000126181A (ja) 腫瘍抽出処理方法
WO2021018103A1 (zh) 超声信号处理方法、装置、设备及存储介质
US6358206B1 (en) Ultrasound process for the determination of the location of a parietal surface in a tissue and of the absolute radius of an artery, and ultrasound apparatus for carrying out such process
CN115644811A (zh) 一种基于光学相干断层成像系统的快速的投影图像重建方法
CN111829956B (zh) 基于超声结构分层引导的光声内窥定量层析成像方法及其系统
CN111260606B (zh) 诊断装置和诊断方法

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