一种去除染色体G显带中期灰度图像噪声方法技术

本发明专利技术提出了一种去除染色体G显带中期灰度图像噪声方法。该方法提出一种局部阈值分割方案，并将经典图像处理算法进行巧妙整合与优化，形成一套具有自适应性且能自动化地处理染色体图像的方法。对于5000幅染色体G显带中期灰度图像数据集，对比本发明专利技术方法处理结果和医务人员手动处理结果得知本发明专利技术方法自动处理准确率达到86.38％。该方法能够较大程度减轻医务人员阅片压力，提高染色体核型分析工作效率。

【技术实现步骤摘要】
一种去除染色体G显带中期灰度图像噪声方法
：本专利技术涉及一种去除染色体G显带中期灰度图像噪声的方法，属于图像处理领域。
技术介绍
：对染色体G显带中期灰度图像进行核型分析时，一幅干净且带纹清晰的图像有利于医务人员进行诊断分析。然而实际情况由于在细胞中期进行制片过程中，染色体图像难免会混入噪声。对大量染色体图像观察发现，染色体图像中常伴有团块噪声、细胞碎片噪声、组织液留下形成的不规则絮状或颗粒状噪声、甚至还有不同细胞染色体溅入对方染色体簇的情况。目前在相关生殖遗传专科医院中，对染色体图像进行噪声去除仍然依赖医生借助相关图像处理软件进行手动操作。当前图像预处理方法涉及到的经典算法已经是相当成熟，然而在解决实际问题中任何一种算法都无法独立胜任。对染色体G显带中期灰度图像进行噪声去除时，阈值分割是常用的算法，考虑到算法的自适应性及运算复杂度问题，本专利技术方法首先采用大津法(Otsu法)[N.Otsu,Athresholdselectionmethodfromgray-levelhistograms,IEEETransactionsonSystems,ManandCybernetics9(1)(1979)62-66]进行图像阈值计算。实验中发现制片后的部分染色体图像因光照不均匀使得图像背景明暗不一，这必然会影响自适应阈值的计算。因此阈值分割前当判断染色体图像的背景不均匀性达到一定阈值条件时，本专利技术方法就对染色体图像进行背景校正，尽量使染色体图像的背景明暗分布均匀。此外染色体区域与背景的对比度较弱时也会影响阈值分割结果，为了避免这种情况，本专利技术提出自适应计算染色体区域对比度，进而可以自适应修正gamma校正的系数达到增强染色体区域对比度的目的。其次本专利技术方法采用了数学形态学中的腐蚀、膨胀和开闭运算对阈值分割后的二值化图像进一步分析其几何形状特征，从而达到剔除大量团块噪声和颗粒噪声等区域。最后利用连通区域标记算法提取出每个独立连通的染色体以及剩余残留噪声区域，并对其执行聚类算法达到进一步剔除噪声的目的。针对染色体灰度图像，本专利技术方法具有两大创新点，其一所提方法并未将大津法直接运用到染色体图像的全局区域，而是结合数学形态学的膨胀运算提取出染色体所在的局部区域。其二所提方法综合和优化了多种经典图像处理算法基本达到了自适应去除噪声的目的。名词解释聚类算法：将具有相同或相似性特征的向量标记为同一个类别的算法。
技术实现思路
：由于染色体中期灰度图像存在各类噪声区域，本专利技术目的是提出一种去除染色体G显带中期灰度图像噪声的方法。该方法能够自适应地去除染色体中期灰度图像噪声、提取染色体有效区域，为接下来的核型分析奠定基础。本专利技术通过以下方案实现：一种去除染色体G显带中期灰度图像噪声方法，包括如下步骤：步骤一、读入一幅原始染色体G显带中期灰度图像m_imgSrc，数据类型为8bits整数型，图像分辨率为950×760，即行像素数目R＝760，列像素数目C＝950，其中任意一幅图像像素灰度值的索引方式为：以图像m_imgSrc为例，其坐标(i,j)的像素灰度值为m_imgSrc(i,j)，i表示图像行坐标，j表示图像列坐标；步骤二、对原始图像m_imgSrc进行3×3中值滤波得到图像m_imgGray，如下式；m_imgGray(i,j)＝med{m_imgSrc(i+s,j+t)|s,t∈{-1,0,1}}其中med{}为中值滤波运算符，表示对med{}运算符中的元素进行排序后取中间数值，s和t表示像素坐标(i,j)的邻域范围，表示图像m_imgGray在坐标(i,j)的灰度值；步骤三、提取步骤二图像m_imgGray的背景掩模区域m_imgBKmask：一)将灰度图像m_imgGray进行分块操作，每一子块窗口区域的像素尺寸为50×50，即为一个尺寸为50×50步长为50的滑窗；二)分别计算所有子块窗口内灰度分布的标准差，将标准差小于阈值T1的子块窗口区域归类为背景区域；步骤四、计算灰度图像m_imgGray的背景掩模区域m_imgBKmask的标准差std_bk，当std_bk大于阈值T2时，判断图像m_imgGray的背景具有不均匀性，此时对m_imgGray进行背景校正：a)初始化背景图像m_imgBKgray为m_imgGray；b)在背景图像m_imgBKgray上放置一个尺寸为32×32步长为8个像素的滑窗；c)对m_imgBKgray从左到右、从上至下移动滑窗，并计算每次滑窗内图像灰度分布的标准差σ、平均值μ、最大值maxval；d)将滑窗内的所有像素灰度值重新赋为min{maxval,μ+3*σ}，其中min{}为计算最小值运算符，表示计算min{}运算符中的元素最小值；解释：min{maxval,μ+3*σ}的计算结果为maxval和μ+3*σ两个元素的最小那个。e)重复步骤c)和d)，完成背景图像m_imgBKgray的构建；f)计算背景校正时的差分图像m_imgDiff，即m_imgGray减去背景图像m_imgBKgray，计算结果为m_imgDiff(i，j)＝uint8(255-|m_imgGray(i，j)-m_imgBKgray(i，j)|)其中uint8()为8bits整数型强制转化运算符，i，j的取值范围为i∈[1,R],j∈[1,C]；g)对m_imgDiff进行5×5中值滤波即完成灰度图像m_imgGray的背景校正；步骤五、对灰度图像m_imgGray运用大津法进行自适应阈值分割，分割结果为一幅二值化图像m_imgBW，即m_imgBW的灰度值只有0或1，分割阈值为T3，计算公式如下：其中i，j的取值范围为i∈[1,R],j∈[1,C]；步骤六、对步骤五的m_imgBW进行形态学膨胀运算得到一幅新的二值化图像m_imgFore，定义灰度图像m_imgGray的前景区域的坐标集合为：F＝{(i,j)|m_imgFore(i,j)＝1,i＝1,2...R,j＝1,2...C}其中膨胀运算的结构元素是半径为9的圆盘，F表示灰度图像m_imgGray的前景区域的坐标集合；步骤七、计算灰度图像m_imgGray的前景区域的灰度平均值gray_mean：步骤八、根据步骤七计算的灰度图像m_imgGray的前景区域灰度平均值定义伽马校正系数γ：γ＝10·gray_mean2步骤九、根据伽马校正系数γ校正灰度图像m_imgGray，得到图像m_imgGam；步骤十、定义任意灰度图像的灰度直方图H：假设任意一幅灰度图像Img，其数据类型为8bits整数型，即灰度值变化为从0到255，那么图像Img的灰度直方图是一个一维向量，其定义如下：其中Np表示图像Img中灰度值为p的像素总数；对伽马校正系数γ校正的图像m_imgGam进行7×7中值滤波，并统计m_imgGam的灰度直方图HG；步骤十一、对校正图像m_imgGam的直方图HG前255个灰阶分布统计运用大津法计算一个自适应阈值T4，并用该阈值对m_imgGam进行二值化分割，同时更新步骤五的m_imgBW，结果如下：步骤十二、更新步骤六灰度图像m_imgGray的前景区域坐标集合：F＝{(i,j)|m_imgBW(i,j)＝1,i＝1,2...R,本文档来自技高网...

【技术保护点】
1.一种去除染色体G显带中期灰度图像噪声方法，其特征在于，包括如下步骤：步骤一、读入一幅原始染色体G显带中期灰度图像m_imgSrc，其中任意一幅图像像素灰度值的索引方式为：以图像m_imgSrc为例，其坐标(i,j)的像素灰度值为m_imgSrc(i,j)，i表示图像行坐标，j表示图像列坐标；步骤二、对原始图像m_imgSrc进行R1×R2中值滤波得到图像m_imgGray，如下式；m_imgGray(i,j)＝med{m_imgSrc(i+s,j+t)|s,t∈{‑1,0,1}}其中med{}为中值滤波运算符，表示对med{}运算符中的元素进行排序后取中间数值，s和t表示像素坐标(i,j)的邻域范围，表示图像m_imgGray在坐标(i,j)的灰度值；步骤三、提取步骤二图像m_imgGray的背景掩模区域m_imgBKmask：一)将灰度图像m_imgGray进行分块操作，每一子块窗口区域的像素尺寸为C1×C2，即为一个尺寸为C1×C2步长为S1的滑窗；二)分别计算所有子块窗口内灰度分布的标准差，将标准差小于阈值T1的子块窗口区域归类为背景区域；步骤四、计算灰度图像m_imgGray的背景掩模区域m_imgBKmask的标准差std_bk，当std_bk大于阈值T2时，判断图像m_imgGray的背景具有不均匀性，此时对m_imgGray进行背景校正：a)初始化背景图像m_imgBKgray为m_imgGray；b)在背景图像m_imgBKgray上放置一个像素尺寸为C3×C4,步长为S2个像素的滑窗；c)对m_imgBKgray从左到右、从上至下移动滑窗，并计算每次滑窗内图像灰度分布的标准差σ、平均值μ、最大值maxval；d)将滑窗内的所有像素灰度值重新赋为min{maxval,μ+3*σ}，其中min{}为计算最小值运算符，表示计算min{}运算符中的元素最小值；e)重复步骤c)和d)，完成背景图像m_imgBKgray的构建；f)计算背景校正时的差分图像m_imgDiff，即m_imgGray减去背景图像m_imgBKgray，计算结果为m_imgDiff(i,j)＝uint8(255‑|m_imgGray(i,j)‑m_imgBKgray(i,j))其中uint8()为8bits整数型强制转化运算符，i，j的取值范围为i∈[1,R],j∈[1,C]；g)对m_imgDiff进行R3×R4中值滤波即完成灰度图像m_imgGray的背景校正；步骤五、对灰度图像m_imgGray运用大津法进行自适应阈值分割，分割结果为一幅二值化图像m_imgBW，即m_imgBW的灰度值只有0或1，分割阈值为T3，计算公式如下：...

【技术特征摘要】
1.一种去除染色体G显带中期灰度图像噪声方法，其特征在于，包括如下步骤：步骤一、读入一幅原始染色体G显带中期灰度图像m_imgSrc，其中任意一幅图像像素灰度值的索引方式为：以图像m_imgSrc为例，其坐标(i,j)的像素灰度值为m_imgSrc(i,j)，i表示图像行坐标，j表示图像列坐标；步骤二、对原始图像m_imgSrc进行R1×R2中值滤波得到图像m_imgGray，如下式；m_imgGray(i,j)＝med{m_imgSrc(i+s,j+t)|s,t∈{-1,0,1}}其中med{}为中值滤波运算符，表示对med{}运算符中的元素进行排序后取中间数值，s和t表示像素坐标(i,j)的邻域范围，表示图像m_imgGray在坐标(i,j)的灰度值；步骤三、提取步骤二图像m_imgGray的背景掩模区域m_imgBKmask：一)将灰度图像m_imgGray进行分块操作，每一子块窗口区域的像素尺寸为C1×C2，即为一个尺寸为C1×C2步长为S1的滑窗；二)分别计算所有子块窗口内灰度分布的标准差，将标准差小于阈值T1的子块窗口区域归类为背景区域；步骤四、计算灰度图像m_imgGray的背景掩模区域m_imgBKmask的标准差std_bk，当std_bk大于阈值T2时，判断图像m_imgGray的背景具有不均匀性，此时对m_imgGray进行背景校正：a)初始化背景图像m_imgBKgray为m_imgGray；b)在背景图像m_imgBKgray上放置一个像素尺寸为C3×C4,步长为S2个像素的滑窗；c)对m_imgBKgray从左到右、从上至下移动滑窗，并计算每次滑窗内图像灰度分布的标准差σ、平均值μ、最大值maxval；d)将滑窗内的所有像素灰度值重新赋为min{maxval,μ+3*σ}，其中min{}为计算最小值运算符，表示计算min{}运算符中的元素最小值；e)重复步骤c)和d)，完成背景图像m_imgBKgray的构建；f)计算背景校正时的差分图像m_imgDiff，即m_imgGray减去背景图像m_imgBKgray，计算结果为m_imgDiff(i,j)＝uint8(255-|m_imgGray(i,j)-m_imgBKgray(i,j))其中uint8()为8bits整数型强制转化运算符，i，j的取值范围为i∈[1,R],j∈[1,C]；g)对m_imgDiff进行R3×R4中值滤波即完成灰度图像m_imgGray的背景校正；步骤五、对灰度图像m_imgGray运用大津法进行自适应阈值分割，分割结果为一幅二值化图像m_imgBW，即m_imgBW的灰度值只有0或1，分割阈值为T3，计算公式如下：其中i，j的取值范围为i∈[1,R],j∈[1,C]；步骤六、对步骤五的m_imgBW进行形态学膨胀运算得到一幅新的二值化图像m_imgFore，定义灰度图像m_imgGray的前景区域的坐标集合为：F＝{(i,j)|m_imgFore(i,j)＝1,i＝1,2...R,j＝1,2...C}其中膨胀运算的结构元素是半径为9的圆盘，F表示灰度图像m_imgGray的前景区域的坐标集合；步骤七、计算灰度图像m_imgGray的前景区域的灰度平均值gray_mean：步骤八、根据步骤七计算的灰度图像m_imgGray的前景区域灰度平均值定义伽马校正系数γ：γ＝10·gray_mean2步骤九、根据伽马校正系数γ校正灰度图像m_imgGray，得到图像m_imgGam；步骤十、定义任意灰度图像的灰度直方图H：假设任意一幅灰度图像Img，其数据类型为8bits整数型，即灰度值变化为从0到255，那么图像Img的灰度直方图是一个一维向量，其定义如下：其中Np表示图像Img中灰度值为p的像素总数；对伽马校正系数γ校正的图像m_imgGam进行R5×R6中值滤波，并统计m_imgGam的灰度直方图HG；步骤十一、对校正图像m_imgGam的直方图HG前255个灰阶分布统计运用大津法计算一个自适应阈值T4，并用该阈值对m_imgGam进行二值化分割，同时更新步骤五的m_imgBW，结果如下：步骤十二、更新步骤六灰度图像m_imgGray的前景区域坐标集合：F＝{(i,j)|m_imgBW(i,j)＝1,i＝1,2...R,j＝1,2...C}m_imgBW表示步骤十一的二值化图像；步骤十三、对步骤十一的二值化图像m_imgBW进行形态学开运算，结果为m_imgOpen，开运算的结构元素是半径为B1的圆盘；步骤十四、对步骤十三的m_imgOpen进行连通区域标记，即将每个独立的、具有8邻接关系的连通区域赋值为一个标签号，标记后的图像为m_imgOpenLabel，标签号为m的连通区域坐标集合L1(m)定义如下：L1(m)＝{(i,j)|m_imgOpenLabel(i,j)＝m}其中i＝1,2...R,j＝1,2...C,m＝1,2...maxlabel1，maxlabel1为m_imgOpenLabel中最大标签号，m表示区间[1,maxlabel1]中的所有整数；步骤十五、对步骤十一的二值化图像m_imgBW进行形态学腐蚀运算，结果为m_imgErode，腐蚀运算的结构元素是半径为B2的圆盘；步骤十六、对步骤十五的m_imgErode进行连通区域...

【专利技术属性】
技术研发人员：丰生日，蔡自兴，卢光琇，蔡昱峰，林戈，穆阳，谭跃球，李仪，
申请(专利权)人：湖南省自兴人工智能研究院，湖南霖德医疗产业有限公司，
类型：发明
国别省市：湖南,43