CN108846802A - 一种去除染色体g显带中期灰度图像噪声方法 - Google Patents

Abstract

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

Description

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

技术领域：

本发明涉及一种去除染色体G显带中期灰度图像噪声的方法，属于图像处理领域。

背景技术：

对染色体G显带中期灰度图像进行核型分析时，一幅干净且带纹清晰的图像有利于医务人员进行诊断分析。然而实际情况由于在细胞中期进行制片过程中，染色体图像难免会混入噪声。对大量染色体图像观察发现，染色体图像中常伴有团块噪声、细胞碎片噪声、组织液留下形成的不规则絮状或颗粒状噪声、甚至还有不同细胞染色体溅入对方染色体簇的情况。目前在相关生殖遗传专科医院中，对染色体图像进行噪声去除仍然依赖医生借助相关图像处理软件进行手动操作。

当前图像预处理方法涉及到的经典算法已经是相当成熟，然而在解决实际问题中任何一种算法都无法独立胜任。对染色体G显带中期灰度图像进行噪声去除时，阈值分割是常用的算法，考虑到算法的自适应性及运算复杂度问题，本发明方法首先采用大津法(Otsu法)[N.Otsu,A threshold selection method from gray-level histograms,IEEETransactions on Systems,Man and Cybernetics 9(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_mean²

步骤九、根据伽马校正系数γ校正灰度图像m_imgGray，得到图像m_imgGam；

步骤十、定义任意灰度图像的灰度直方图H：

假设任意一幅灰度图像Img，其数据类型为8bits整数型，即灰度值变化为从0到255，那么图像Img的灰度直方图是一个一维向量，其定义如下：

其中N_p表示图像Img中灰度值为p的像素总数；

对伽马校正系数γ校正的图像m_imgGam进行7×7中值滤波，并统计m_imgGam的灰度直方图H_G；

步骤十一、对校正图像m_imgGam的直方图H_G前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，开运算的结构元素是半径为3的圆盘；

步骤十四、对步骤十三的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，腐蚀运算的结构元素是半径为5的圆盘；

步骤十六、对步骤十五的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...maxlabel2

sum(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′的连通区域坐标集合；

步骤十九、按步骤十八，计算连通区域排序后的面积变化率拐点n_T，拐点n_T判断准则为：

S′(n_T)＞S′(n_T-1)+S′(n_T-2)+S′(n_T-3)

具体在计算符合拐点判断准则条件的n_T时，确定计算区间为n_T∈[round(maxlabel2*0.6),maxlabel2]，round()为四舍五入运算符；当计算得到符合条件的n_T后，依照步骤十八的排序结果将标签号介于区间[n_T,maxlabel2]的连通区域判定为大面积噪声区域，当没有符合条件的n_T时，说明当前处理的染色体图像不存在大面积噪声区域；

步骤二十、根据步骤十八更新后的L2(n′)和步骤十九的拐点判断准则定义大面积噪声区域坐标集合E2(n_D)：

E2(n_D)＝L2(n_D)

其中n_D∈[n_T,maxlabel2]，由此可见，E2(n_D)是步骤十八L2(n′)的子集；

步骤二十一、将步骤十四中的大面积噪声区域去除；

步骤二十二、运用聚类算法对步骤十四的连通区域集合进一步清洗，将不符合聚类条件的连通区域视为噪声区域进行去除，符合条件的连通区域则为最终有效的染色体连通区域，得到最终只保留了染色体区域而去除了噪声区域的图像。

进一步的改进，所述步骤九中，根据伽马校正系数γ校正灰度图像m_imgGray的计算方法如下：

a)创建一个伽马校正查询数组LUT，用于表示输入灰度值p与校正后灰度值LUT(p)之间的映射关系：

其中round()为四舍五入运算符，LUT的维度是256×1，即为256行1列的数组，p代表区间[0,255]之间的任意整数；

b)对灰度图像m_imgGray的前景区域校正结果m_imgGam为：

进一步的改进，所述步骤二十一中，将步骤十四中的大面积噪声区域去除的方法为通过判断噪声区域集合E2(n_D)与L1(m)是否存在交集，若存在标签号为m的连通区域符合下式，则将m标记的连通区域去除；

m＝1,2...maxlabel1,n_D＝n_T,n_T+1...maxlabel2}

其中表示空集。

进一步的改进，所述步骤二十二包括如下步骤：

(a)将步骤十四中标签号为m的连通区域坐标集合L1(m)映射到灰度图像m_imgGray上，统计标签号为m的连通区域的灰度直方图H_m，将H_m中元素平均分成m^*个区间，将相邻区间元素进行合并，最后作归一化处理形成m^*维度的向量：

其中为向量H_m的第m^*个元素；

(b)计算步骤十四中标签号为m的连通区域坐标集合L1(m)的形心坐标D_m，D_m为2维向量：

其中x_m和y_m分别表示标签号为m的连通区域的形心x坐标和y坐标，i，j正如步骤十四所指分别为图像的行、列坐标；

(c)遍历L1(m)中所有连通区域，对任意两个连通区域进行聚类判断，将符合如下聚类条件的连通区域视为同一类别：

||H_m1-H_m2||₂＞T5,m1,m2∈[1,maxlabel1]

||D_m1-D_m2||₂＜T6，m1,m2∈[1,maxlabel1]

其中T5为直方图相似度阈值，T6为欧氏距离阈值；

(d)去除不符合聚类条件的连通区域，得到最终只保留了染色体区域而去除了噪声区域的图像，得到最终只保留了染色体区域而去除了噪声区域的图像。

进一步的改进，所述阈值T1的取值范围为3.5～4.5，阈值T2的取值为10.0，阈值T5的取值为0.78，T6的取值范围一般为150～300。

进一步的改进，所述T1取值为4.0，T6取值为150。

聚类算法解释如下：

假设给定K组向量v₁、v₂、v₃…v_K，其中v₁＝[v₁₁,v₁₂,v₁₃...v_1M]，v_K＝[v_K1,v_K2,v_K3...v_KM]，其余向量表达式同理。那么本专利采用的聚类算法实现步骤如下：

步骤一、初始化预定义变量：类编号计数class＝1，当前第class类的向量集为空

步骤二、访问当前第class类的向量集是否为空，若为空，则从给定K组向量中任意拿出一个向量放入第class类的向量集中，若不为空，则执行步骤三；

步骤三、从第class类的向量集中取出第1个向量作为基准比较向量v₀；

步骤四、逐个计算给定的K组向量与v₀的欧氏距离，例如向量v_k与v₀之间的欧氏距离值为d＝||v₀-v_k||₂，其中k∈[1,K]。若向量v_k与v₀满足聚类条件则说明这两个向量属于同一个类别，并将v_k放入第class类的向量集中(这里的聚类条件可根据实际情况定义，例如||v₀-v_k||₂＜T^*，T^*为距离阈值)；

步骤五、依次从第class类的向量集中取出第2个、第3个…向量作为基准比较向量v₀，重复步骤四；

步骤六、当第class类的向量集中所有元素都用作过基准向量后，此时属于同一个class类别的向量集就创建完成。接着整理给定的K组向量，将已经赋值为类别号的向量从给定K组向量中剔除，然后class的数值自增1，初始化新的第class类的向量集为空并对剩余的K组(此时K的数值会减小)向量执行步骤二。

如上所述步骤旨在将一幅染色体G显带中期灰度图像通过一系列经典图像处理算法的操作以达到正确标记染色体区域，剔除大面积的团块噪声、颗粒噪声及背景噪声区域的目的。一般情况下，由于染色体本身的姿态多样性，以及制片过程的复杂性，一幅染色体图像中除了包含各类型噪声外，染色体自身也会形成类似两体粘连、多体粘连、重叠以及交叉等多种形态，这就导致难以寻找一种普适性较好的图像处理方法既能达到剔除噪声的目的，又能很好保留染色体区域。上述步骤十三和十五都运用了数学形态学方法对染色体二值化图像进行处理，其中步骤十三的开运算一方面能够剔除颗粒噪声区域，另一方面也能将轻度粘连的染色体进行分割。这个步骤中结构元素的尺寸选择很大程度会影响步骤十四的染色体区域标记，因为较大尺寸的结构元素会损失染色体信息。例如将一条完整的染色体区域分断，又如将21和22号染色体的随体剔除，而较小的结构元素在去除颗粒噪声区域方面有限。步骤十五运用的腐蚀运算主要是为了确定大面积噪声区域，由于染色体的粘连性，为了避免删除多体粘连的染色体区域，这个步骤的结构元素尺寸选择既要考虑能将多体粘连尽量分开，又要避免腐蚀掉过多染色体区域。除此之外，染色体带纹强度的深浅不一很容易在步骤五和十一的大津法阈值分割中将完整染色体分断，因此中值滤波主要是为了尽量模糊染色体带纹细节，以避免这种情况。

总之，对于形态多样，噪声丰富的染色体G显带中期灰度图像，本发明方法能够自适应地完成染色体核型分析前的去除噪声功能。通过MATLAB平台编程对5000幅染色体G显带中期灰度图像进行处理并统计将染色体灰度图像全部噪声剔除的准确率为86.38％，实验结果说明本发明方法基本能够降低目前医生的手动操作压力，提高相应的工作效率。

上述的中值过滤、步长、像素尺寸、圆盘半径等均是试验得出的最佳竖直，也可根据实际情况进行调整变更。

附图说明：

附图1：一幅原始染色体G显带中期灰度图像，其背景比较均匀；

附图2：对附图1进行3×3中值滤波后的图像；

附图3：附图2的背景掩模区域图像；

附图4：对附图2进行大津法阈值分割的图像；

附图5：对附图4进行膨胀运算后的二值化图像，并定义其为染色体前景区域；

附图6：附图2的染色体前景区域图像；

附图7：对附图6进行伽马校正后的图像；

附图8：对附图7进行7×7中值滤波后的图像；

附图9：对附图8的前景区域进行大津法阈值分割后并作开运算的二值化图像；

附图10：对附图9中大津法分割后的图像进行腐蚀运算的二值化图像；

附图11：从附图9中删除面积较大的噪声区域后的图像；

附图12：通过聚类算法进一步去除噪声区域后染色体图像；

附图13：一幅原始染色体G显带中期灰度图像，其背景分布不均匀；

附图14：对附图13进行3×3中值滤波后的图像；

附图15：从附图14提取的背景的图像；

附图16：对附图14进行背景校正后的结果；

附图17：附图16的前景区域图像；

附图18：对附图17进行伽马校正、中值滤波后的结果；

附图19：对附图18进行大津法分割以及开运算操作的结果；

附图20：去除大面积噪声、颗粒噪声等噪声区域后的染色体图像。

具体实施方式：

先结合具体实施例及附图，来进一步阐述本发明。

实施例一：针对背景比较均匀的染色体G显带中期灰度图像

(1)读入一副原始染色体G显带中期灰度图像m_imgSrc，如附图1；

(2)对m_imgSrc进行3×3中值滤波得到m_imgGray，如附图2；

(3)提取m_imgGray的背景掩模区域，如附图3，白色区域就标记着背景；

(4)计算m_imgGray的背景掩模区域的灰度分布标准差std_bk＝3.477，设置阈值T2＝10，由于std_bk<T2，所以背景区域基本均匀无需进行校正；

(5)对灰度图像m_imgGray进行大津法二值化分割得到m_imgBW，结果如附图4；

(6)对步骤(5)m_imgBW进行膨胀运算得到染色体前景区域m_imgFore，如附图5，其对应的前景区域灰度图像如附图6；

(7)对附图6进行伽马校正如附图7，并作7×7中值滤波如附图8；

(8)对附图8染色体前景灰度图像进行大津法分割，并用开运算剔除部分颗粒噪声，如附图9；

(9)对步骤(8)大津法分割后的结果进行腐蚀运算，如附图10，自适应计算面积较大的连通区域；

(10)从附图9剔除面积较大连通区域如附图11；

(11)通过聚类算法进一步剔除噪声区域后的染色体图像如附图(12)，到此就完成了对一幅染色体G显带中期灰度图像m_imgSrc的噪声去除。

实施例二：针对背景不均匀的染色体G显带中期灰度图像

(1)读入一副原始染色体G显带中期灰度图像m_imgSrc，如附图13；

(2)对m_imgSrc作3×3中值滤波得m_imgGray，如附图14；

(3)提取m_imgGray背景掩模区域并计算其灰度分布标准差std_bk＝30.96，设置阈值T2＝10，由于std_bk>T2，所以m_imgGray需要校正；

(4)提取m_imgGray的背景图像m_imgBKgray如附图15；

(5)对m_imgGray进行背景校正，如附图16；

(6)运用大津法分割、膨胀运算提取m_imgGray的前景区域，如附图17，并做伽马校正如附图18；

(7)对附图18的前景区域进行大津法分割，并用开运算剔除部分颗粒噪声，如附图19；

(8)剔除大面积噪声区域后再通过聚类算法进一步剔除噪声的结果如附图20，到此就完成了对一幅染色体G显带中期灰度图像m_imgSrc的噪声去除。

Claims (7)

1.一种去除染色体G显带中期灰度图像噪声方法，其特征在于，包括如下步骤：

步骤一、读入一幅原始染色体G显带中期灰度图像m_imgSrc，其中任意一幅图像像素灰度值的索引方式为：以图像m_imgSrc为例，其坐标(i,j)的像素灰度值为m_imgSrc(i,j)，i表示图像行坐标，j表示图像列坐标；

步骤二、对原始图像m_imgSrc进行R₁×R₂中值滤波得到图像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进行分块操作，每一子块窗口区域的像素尺寸为C₁×C₂，即为一个尺寸为C₁×C₂步长为S₁的滑窗；

二)分别计算所有子块窗口内灰度分布的标准差，将标准差小于阈值T1的子块窗口区域归类为背景区域；

步骤四、计算灰度图像m_imgGray的背景掩模区域m_imgBKmask的标准差std_bk，当std_bk大于阈值T2时，判断图像m_imgGray的背景具有不均匀性，此时对m_imgGray进行背景校正：

a)初始化背景图像m_imgBKgray为m_imgGray；

b)在背景图像m_imgBKgray上放置一个像素尺寸为C₃×C₄,步长为S₂个像素的滑窗；

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进行R₃×R₄中值滤波即完成灰度图像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_mean²

步骤九、根据伽马校正系数γ校正灰度图像m_imgGray，得到图像m_imgGam；

步骤十、定义任意灰度图像的灰度直方图H：

假设任意一幅灰度图像Img，其数据类型为8bits整数型，即灰度值变化为从0到255，那么图像Img的灰度直方图是一个一维向量，其定义如下：

其中N_p表示图像Img中灰度值为p的像素总数；

对伽马校正系数γ校正的图像m_imgGam进行R₅×R₆中值滤波，并统计m_imgGam的灰度直方图H_G；

步骤十一、对校正图像m_imgGam的直方图H_G前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，开运算的结构元素是半径为B₁的圆盘；

步骤十四、对步骤十三的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，腐蚀运算的结构元素是半径为B₂的圆盘；

步骤十六、对步骤十五的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...maxlabel2

sum(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′的连通区域坐标集合；

步骤十九、按步骤十八，计算连通区域排序后的面积变化率拐点n_T，拐点n_T判断准则为：

S′(n_T)＞S′(n_T-1)+S′(n_T-2)+S′(n_T-3)

具体在计算符合拐点判断准则条件的n_T时，确定计算区间为

n_T∈[round(maxlabel2*0.6),maxlabel2]，round()为四舍五入运算符；当计算得到符合条件的n_T后，依照步骤十八的排序结果将标签号介于区间[n_T,maxlabel2]的连通区域判定为大面积噪声区域，当没有符合条件的n_T时，说明当前处理的染色体图像不存在大面积噪声区域；

步骤二十、根据步骤十八更新后的L2(n′)和步骤十九的拐点判断准则定义大面积噪声区域坐标集合E2(n_D)：

E2(n_D)＝L2(n_D)

其中n_D∈[n_T,maxlabel2]，由此可见，E2(n_D)是步骤十八L2(n′)的子集；

步骤二十一、将步骤十四中的大面积噪声区域去除；

步骤二十二、运用聚类算法对步骤十四的连通区域集合进一步清洗，将不符合聚类条件的连通区域视为噪声区域进行去除，符合条件的连通区域则为最终有效的染色体连通区域，得到最终只保留了染色体区域而去除了噪声区域的图像。

2.如权利要求1所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，所述步骤九中，根据伽马校正系数γ校正灰度图像m_imgGray的计算方法如下：

a)创建一个伽马校正查询数组LUT，用于表示输入灰度值p与校正后灰度值LUT(p)之间的映射关系：

其中round()为四舍五入运算符，LUT的维度是256×1，即为256行1列的数组，p代表区间[0,255]之间的任意整数；

b)对灰度图像m_imgGray的前景区域校正结果m_imgGam为：

3.如权利要求1所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，所述步骤二十一中，将步骤十四中的大面积噪声区域去除的方法为通过判断噪声区域集合E2(n_D)与L1(m)是否存在交集，若存在标签号为m的连通区域符合下式，则将m标记的连通区域去除；

其中表示空集。

4.如权利要求1所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，所述步骤二十二包括如下步骤：

(a)将步骤十四中标签号为m的连通区域坐标集合L1(m)映射到灰度图像m_imgGray上，统计标签号为m的连通区域的灰度直方图H_m，将H_m中元素平均分成m^*个区间，将相邻区间元素进行合并，最后作归一化处理形成m^*维度的向量：

其中为向量H_m的第m^*个元素；

(b)计算步骤十四中标签号为m的连通区域坐标集合L1(m)的形心坐标D_m，D_m为2维向量：

其中x_m和y_m分别表示标签号为m的连通区域的形心x坐标和y坐标，i，j正如步骤十四所指分别为图像的行、列坐标；

(c)遍历L1(m)中所有连通区域，对任意两个连通区域进行聚类判断，将符合如下聚类条件的连通区域视为同一类别：

||H_m1-H_m2||₂＞T5,m1,m2∈[1,maxlabel1]

||D_m1-D_m2||₂＜T6,m1,m2∈[1,maxlabel1]

其中T5为直方图相似度阈值，T6为欧氏距离阈值；

(d)去除不符合聚类条件的连通区域，得到最终只保留了染色体区域而去除了噪声区域的图像，得到最终只保留了染色体区域而去除了噪声区域的图像。

5.如权利要求1所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，所述阈值T1的取值范围为3.5～4.5，阈值T2的取值为10.0，阈值T5的取值为0.78，T6的取值范围一般为150～300。

6.如权利要求5所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，所述T1取值为4.0，T6取值为150。

7.如权利要求5所述的去除染色体G显带中期灰度图像噪声方法，其特征在于，R₁和R₂等于3，C₁和C₂均等于50，C₃和C₄均等于32，S₁等于50，S₂等于8，R₃和R₄均等于5，R₅和R₆均等于7，B₁等于3，B₂等于5。