[发明专利]一种去除染色体G显带中期灰度图像噪声方法

摘要

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

主权项

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进行连通区域标记，即将每个独立的、具有8邻接关系的连通区域赋值为一个标签号，标记后的图像为m_imgErodeLabel，标签号为n的连通区域坐标集合L2(n)定义如下：L2(n)＝{(i,j)|m_imgErodeLabel(i,j)＝n}其中i＝1,2...R,j＝1,2...C,n＝1,2...maxlabel2，maxlabel2为m_imgErodeLabel中最大标签号，n表示区间[1,maxlabel2]中的所有整数；步骤十七、定义标签号为q的连通区域面积为sum(q)；其中sum(q)表示标签号为q的连通区域坐标集合的元素个数；接着统计步骤十六中L2(n)的每个标签号对应的连通区域面积S(n)：S(n)＝sum(n),n＝1,2...maxlabel2sum(n)＝∑L2(n)其中∑L2(n)表示标签号为n的连通区域坐标集合元素个数，n如步骤十六所指；步骤十八、对步骤十七中每个标签号的连通区域按面积由小到大进行排序，并按照排序顺序修改排序后连通区域的标签号：S′(n′)＝sort(S(n))其中sort()为面积由小到大排序运算符，n＝1,2...maxlabel2，n′为按面积排序后的连通区域标签号；与此同时更新步骤十六的L2(n)为L2(n′)：L2(n′)＝{(i,j)|m_imgErodeLabel(i,j)＝n′}n′即为排序后的标签号，L2(n′)表示标签号为n′的连通区域坐标集合；步骤十九、按步骤十八，计算连通区域排序后的面积变化率拐点nT，拐点nT判断准则为：S′(nT)＞S′(nT‑1)+S′(nT‑2)+S′(nT‑3)具体在计算符合拐点判断准则条件的nT时，确定计算区间为nT∈[round(maxlabel2*0.6),maxlabel2]，round()为四舍五入运算符；当计算得到符合条件的nT后，依照步骤十八的排序结果将标签号介于区间[nT,maxlabel2]的连通区域判定为大面积噪声区域，当没有符合条件的nT时，说明当前处理的染色体图像不存在大面积噪声区域；步骤二十、根据步骤十八更新后的L2(n′)和步骤十九的拐点判断准则定义大面积噪声区域坐标集合E2(nD)：E2(nD)＝L2(nD)其中nD∈[nT,maxlabel2]，由此可见，E2(nD)是步骤十八L2(n′)的子集；步骤二十一、将步骤十四中的大面积噪声区域去除；步骤二十二、运用聚类算法对步骤十四的连通区域集合进一步清洗，将不符合聚类条件的连通区域视为噪声区域进行去除，符合条件的连通区域则为最终有效的染色体连通区域，得到最终只保留了染色体区域而去除了噪声区域的图像。