CN110428380A - 紫色土土壤图像阴影增强方法 - Google Patents

紫色土土壤图像阴影增强方法 Download PDF

Info

Publication number
CN110428380A
CN110428380A CN201910705747.XA CN201910705747A CN110428380A CN 110428380 A CN110428380 A CN 110428380A CN 201910705747 A CN201910705747 A CN 201910705747A CN 110428380 A CN110428380 A CN 110428380A
Authority
CN
China
Prior art keywords
value
soil
pixel
region
shadow
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
CN201910705747.XA
Other languages
English (en)
Other versions
CN110428380B (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.)
Chongqing Normal University
Original Assignee
Chongqing Normal University
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 Chongqing Normal University filed Critical Chongqing Normal University
Priority to CN201910705747.XA priority Critical patent/CN110428380B/zh
Publication of CN110428380A publication Critical patent/CN110428380A/zh
Application granted granted Critical
Publication of CN110428380B publication Critical patent/CN110428380B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/40—Image enhancement or restoration using histogram techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/90—Dynamic range modification of images or parts thereof
    • G06T5/94—Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/90—Determination of colour characteristics

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Analysis (AREA)

Abstract

本发明提供一种紫色土土壤图像阴影增强方法，包括步骤：S1：将紫色土土壤的原始图像划分为非阴影区域、半阴影区域和全阴影区域；S2至S9：将阴影区域亮度进行初始增强，S10：对阴影区域亮度增强后的紫色土土壤的初始图像进行色温校正；S11：对色温校正后的紫色土土壤图像进行局部对比度校正；本发明通过建立紫色土土壤图像的光照模型，对分割出紫色土土壤图像的半阴影区域和全阴影区域，进行亮度增强，以方便后续对紫色土土种的准确识别。

Description

紫色土土壤图像阴影增强方法
技术领域
本发明涉及紫色土图像识别的技术领域，具体涉及一种紫色土土壤图像阴影增强方法。
背景技术
机器视觉识别土壤在农业生产中有着重要的实用价值。在农业生产中，识别土壤是非常重要的。由于土壤分类系统复杂，只有极少数科研院所的专家才能准确辨识土壤，生产一线的农业技术人员要完全准确的把当地土壤辨识清楚也是非常困难的。基于土壤分类系统的土壤类型识别是一个农业生产实践中的共性难题。人工智能技术的发展，使机器视觉识别土壤成为可能。机器视觉辨识土壤是对野外自然条件下拍摄的具有复杂背景的土壤图像进行识别，提高对紫色土土壤图像的分割提取精度。
但是提取出紫色土图像受到直接光照和间接光照的影响，往往存在阴影，阴影的存在，对后续对紫色土土属和土种的识别造成影响，使识别出的紫色土土属和土种过程时间开销大，识别出的结果并不准确，甚至相差甚远。
发明内容
有鉴于此，本发明的目的是提供一种紫色土土壤图像阴影增强方法，通过建立紫色土土壤图像的光照模型，对分割出紫色土土壤图像的半阴影区域和全阴影区域，进行亮度增强，以方便后续对紫色土土属和土种的准确识别。
本发明提供一种紫色土土壤图像阴影增强方法，包括步骤：
S1：将紫色土土壤的原始图像划分为非阴影区域、半阴影区域和全阴影区域；
S2：将紫色土土壤的原始图像转换为HSI颜色模型的H域值矩阵H0；其中，矩阵H0对应于原始图像相同位置的像素的H域值；
将紫色土土壤的原始图像转换为HSI颜色模型的I域值矩阵I0；其中，矩阵I0对应于原始图像相同位置的像素的I域值；
S3：按1o整化统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图Hist_l、Hist_f和Hist_p；其中，Hist_l为非阴影区域中像素点H分量的直方图，Hist_f为全影区域中像素点H分量的直方图，Hist_p为半影区域中像素点H分量的直方图；其中，按1°整化统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图是指将紫色土土壤的原始图像的H分量等角度划分为360等分，每个等分对应一个灰度级；
S4：初始化的标记向量Hist；其中，标记向量Hist初始化为360维的0向量，标记向量Hist中的每一维分别对应为按1°整化得到各等分对应的灰度级的标签；
S5：依次遍历向量Hist_l、Hist_f和Hist_p中的各个灰度级，若灰度级对应的向量Hist_l、Hist_f和Hist_p的值都大于零，则更新该灰度级对应的标记向量中的标签值为1；
S6：重复步骤S5，直到所有的向量Hist_l、Hist_f和Hist_p中的各个灰度级都被遍历，得到更新后的标记向量HistⅠ；
S7：计算非阴影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μl；计算全影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μf；计算半影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μp；
S8：根据μl、μf和μp，计算直接光照和间接光照的比值r、半影区整体的直接光照的遮蔽系数α和全局直接光照强度Id；
S9：根据r、α和Id，计算半影区域和全影区域各像素点增强后的亮度值，并将矩阵I0中各个像素点对应位置的元素更新为增强后的亮度，得到阴影增强后的紫色土土壤图像；
S10：对阴影区域亮度增强后的紫色土土壤的初始图像进行色温校正，得到色温校正后的紫色土土壤图像；
S11：对色温校正后的紫色土土壤图像进行局部对比度校正，得到局部对比度校正后的紫色土土壤图像，即得到最终的阴影区域亮度增强后的紫色土土壤图像。
进一步，所述步骤S9中，半影区域和全影区域各像素点增强后的亮度值的计算公式为：
其中，Ip non_shadow(x,y)为增强后的半影区像素点(x,y)的亮度值，Ip(x,y)为增强前半影区像素点(x,y)的亮度值，If non_shadow(x,y)为增强后的全影区像素点(x,y)的亮度值，If(x,y)为增强前全影区像素点(x,y)的亮度值，αi为像素点(x,y)的直接光照遮蔽系数，r为直接光照与间接光照的比值。
进一步，所述像素点(x,y)的直接光照遮蔽系数αi的计算公式为：
其中，μM×M表示像素点(x,y)的M×M邻域内像素点亮度算术平均值；
所述Id的计算公式为：
所述α的计算公式为：
进一步，所述r的计算公式为：
进一步，所述步骤S10包括步骤：
S101：将最终的阴影区域亮度增强后的紫色土土壤的初始图像转换为Lab颜色模型的a分量矩阵A；其中，矩阵A对应于原始图像相同位置的像素的a分量值；
将最终的阴影区域亮度增强后的紫色土土壤图像转换为Lab颜色模型的b分量矩阵B；其中，矩阵B对应于原始图像相同位置的像素的b分量值；
S102：分别统计得到阴影区域a分量和b分量的统计直方图Hist_as和Hist_bs，其中阴影区域由全阴影区域和半阴影区域组成；分别统计得到非阴影区域a分量和b分量的统计直方图Hist_al和Hist_bl；
S103：分别对阴影区域a分量和b分量的统计直方图进行均衡化，得到阴影区域直方图向量Equal_as和Equal_bs；
分别对非阴影区域a分量和b分量的统计直方图进行均衡化，得到非阴影区直方图向量Equal_al和Equal_bl；
S104：按照单映射规则，分别寻找向量Equal_as中各个灰度级在Equal_al中最接近的灰度级，并将Equal_as中各个灰度级对应的a分量的值ai分别映射到Equal_al中与Equal_as中各个灰度级最接近的灰度级对应矩阵A中的元素，得到新的a分量矩阵Anew；
按照单映射规则，分别寻找向量Equal_bs中各个灰度级在Equal_bl中最接近的灰度级，并将Equal_bs中各个灰度级对应的b分量的值bi分别映射到Equal_bl中与Equal_bs中各个灰度级最接近的灰度级对应矩阵B中的元素，得到新的b分量矩阵Bnew，即得到色温校正后的紫色土土壤图像。
进一步，步骤S104中单映射规则为：
所述Equal_as、Equal_al、Equal_bs、Equal_bl的计算公式如下：
其中，asj为阴影区域中第j个灰度级的a分量，m为阴影区域中灰度级的总数量，nsj为阴影区域中a分量值为asj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(asj)为阴影区域中a分量值为asj的概率；
其中，ali为非阴影区域中第i个灰度级的a分量，n为非阴影区域中灰度级的总数量，kli为非阴影区域中a分量值为ali的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(ali)为阴影区域中a分量值为ali的概率；
其中，bsj为阴影区域中第j个灰度级的b分量，m为阴影区域中灰度级的总数量，nsj'为阴影区域中b分量值为bsj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(bsj)为阴影区域中b分量值为bsj的概率；
其中，bli为非阴影区域中第i个灰度级的b分量，n为非阴影区域中灰度级的总数量，kli'为非阴影区域中b分量值为bli的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(bli)为阴影区域中b分量值为bli的概率。
进一步，所述步骤S11包括步骤：
S111：将最终的阴影区域亮度增强后的紫色土土壤图像转换为HSI颜色模型的I域值矩阵I2；其中，矩阵I2对应于原始图像相同位置的像素的I域值；
S112：计算色温增强后的紫色土土壤图像中每个像素点的局部对比度；
S113：根据局部对比度对矩阵I2中的所有元素进行初步对比度校正，得到进行初步对比度校正后的矩阵I3；
S114：初始化与矩阵I3大小相同的异常标记矩阵Markerror；
S115：选取任一未被遍历的矩阵I3中的元素作为当前遍历元素，判断当前遍历元素值是否大于非阴影区域的像素点的算术平均值，若是，则判断当前遍历元素为亮度初始对比度校正后的异常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置1；若否，则判断当前遍历元素为亮度初始对比度校正后的正常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置0；
S116：重复步骤S115，直到矩阵I3中所有的元素都被遍历，得到更新后的异常标记矩阵MarkerrorI；
S117：随机选取异常标记矩阵MarkerrorI中未被补充校正遍历且值为1的元素，作为当前补充校正遍历元素，计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值，并将矩阵I3中该元素值更新为该补充校正值；
S118：判断异常标记矩阵MarkerrorI中是否还有未被遍历且值为1的元素，若有，则返回步骤S117，若无，则结束对异常标记矩阵MarkerrorI的补充校正遍历，得到补充校正后的I域值矩阵I4，即得到局部对比度校正后的紫色土土壤图像，也就是得到最终的阴影区域亮度增强后的紫色土土壤图像。
进一步，所述步骤S112中局部对比度的计算公式为：
其中，Cg×g(x,y)表示色温增强后的紫色土土壤图像中像素点(x,y)的局部对比度，Iold(x,y)表示紫色土土壤的原始图像中像素点(x,y)的亮度值；μold g×g(x,y)表示紫色土土壤的原始图像中像素点(x,y)的g×g邻域的g×g个像素点亮度的算术平均值；
所述步骤S113中根据局部对比度对矩阵I2中的所有元素进行初步对比度校正的公式为：
Inew(x,y)＝Cg×g(x,y)×μnew g×g(x,y) (12)
其中，Inew(x,y)表示初步对比度校正后的矩阵I3中元素(x,y)的值；μnew g×g(x,y)表示矩阵I2中元素(x,y)的g×g邻域的g×g个元素值的算术平均值。
进一步，所述步骤S117中计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值的公式为：
其中，INEW(x,y)表示矩阵I3中元素(x,y)的补充校正值，表示紫色土土壤的原始图像中非阴影区域的所有像素点亮度的算术平均值；表示紫色土土壤的原始图像中阴影区域的所有像素点亮度的算术平均值；表示阴影区域亮度增强后的紫色土土壤的初始图像中阴影区域的所有像素点亮度的算术平均值；Ierror(x,y)表示矩阵I3中当前补充校正遍历元素对应的元素值。
本发明的有益效果：本发明通过建立紫色土土壤图像的光照模型，对分割出紫色土土壤图像的半阴影区域和全阴影区域，进行亮度增强，以方便后续对紫色土土属和土种的准确识别。
附图说明
下面结合附图和实施例对本发明作进一步描述：
图1为本发明的方法流程图
图2为本发明的实验效果对比图。
具体实施方式
如图1所示，本发明提供的一种紫色土土壤图像阴影增强方法，包括步骤：
S1：将紫色土土壤的原始图像划分为非阴影区域、半阴影区域和全阴影区域；彩色图像中的像素点的亮度由两种光照决定，在图像的非阴影区域像素，其亮度主要由直接光照和间接光照决定；在图像阴影区域的像素点，其亮度主要由部分直接光照和间接光照决定。通过对阴影区域中像素点所受直接光照的量进行判断，将阴影区域像素分为全影区域和半影区域。当像素点所受直接光照为0时，该像素点即可视为全影区域；当像素点受到部分直接光照时，该像素点即可视为半影区域。在半影区域中的像素点由于受到部分直接光照，因此，半影区像素点的亮度是介于非阴影区域和全影区域的。将全影区域像素和半影区域像素视为两类，通过像素点亮度求取分割阈值，实现阴影区域像素的分类，方便后续进一步对紫色土土壤的土种和土属的识别。
具体地，根据图像的光照原理，一幅图像可以视为由反射分量与照度分量构成，反射分量代表了图像中像素点的内在属性，而照度分量可以视为图像所受的外部光照。在自然光照条件下拍摄的土壤图像中，每个像素点的照度分量可以看作是受直接光照和间接光照的影响，其中直接光照可以视为太阳光的照射，间接光照可以视为物体表面光线的漫反射。因此，构成了图像的基本光照模型如下：
I(x,y)＝L(x,y)×R(x,y) (14)
其中，I(x,y)表示原始图像，R(x,y)表示反射图像，L(x,y)表示照度图像中。将照度分量分解为直接照度分量和间接照度分量，得到下式：
I(x,y)＝[Ld(x,y)+Le(x,y)]×R(x,y) (15)
其中，Ld(x,y)和Le(x,y)分别为直接照度图像和间接照度图像。
在自然条件下拍摄的图像中，处于图像非阴影区域的像素，其所受光照为直接光照和间接光照，对于处于图像全影区域的像素，其所受光照仅为间接光照。即为：
Il(x,y)＝[Ld(x,y)+Le(x,y)]×R(x,y) (16)
Is(x,y)＝Le(x,y)×R(x,y) (17)
对于一幅自然条件下拍摄的土壤图像，可以视为光照充足图像与阴影图像的组合，得到自然条件下图像光照模型：
I(x,y)＝ki×Il(x,y)+(1-ki)×Is(x,y) (18)
其中Il(x,y)表示光照充足图像，Is(x,y)表示阴影图像，ki表示遮蔽系数。将式(16)和式(17)带入式(18)，得到：
上式即为图像的基本光照模型，其中当ki＝1时，表示当前像素点位于图像的非阴影区域；当ki＝0时，表示当前像素点位于图像的全影区域；当0＜ki＜1时，表示当前像素点位于图像的半影区域。在本实施例中，将紫色土土壤的原始图像划分为非阴影区域、半阴影区域和全阴影区域可以采用现有的划分算法，也可以采用如下划分方法。
紫色土土壤的原始图像的非阴影区域、半阴影区域和全阴影区域的划分的方法，包括步骤：
S1-1：将紫色土土壤的原始图像分割为阴影区域和非阴影区域；本实施例中，所述紫色土土壤的原始图像是指只包含有紫色土土壤的彩色图像，可以通过现有的紫色土图像分割提取方法对采集到的包含有紫色土的图像进行紫色土图像的分割提取得到，在此不赘述。
S1-2：消除离散的阴影区域，得到保留的阴影区域；
S1-3：将保留的阴影区域分割为全影区域和半影区域。通过上述方法，可以精确分割出紫色土土壤图像的非阴影区域、半阴影区域和全阴影区域，方便后续对紫色土土属和土种的准确识别。
进一步，所述步骤S1-1包括步骤：
S1-11：将原始图像转换为HSI颜色模型的I域值矩阵；初始化与I域值矩阵大小相同的标记矩阵Mark；其中，I域值矩阵的元素值对应原始图像相同位置的像素的I域值，标记矩阵Mark初始化为零矩阵；
S1-12：计算I域值矩阵的整体均值μI，将μI作为分割非阴影区域和阴影区域的初始分类阈值T0；将T0作为当前分类阈值；
S1-13：将I域值矩阵中元素值小于当前分类阈值的元素划分为当前阴影区域，将I域值矩阵中元素值大于或等于当前分类阈值的元素划分为当前非阴影区域；
S1-14：分别计算当前阴影区域的算术平均值和当前非阴影区域的元素值的算术平均值；
S1-15：在当前阴影区域的元素值的算术平均值与当前非阴影区域的算术平均值的之间范围内，获得新分类阈值Tnew；即μ1＜Tnew＜μ2，其中，μ1和μ2分别为当前阴影区域的所有元素值的算术平均值和当前非阴影区域的所有元素值的算术平均值；
S1-16：判断新分类阈值Tnew是否等于当前分类阈值；若是，停止迭代，将最后一次迭代中的当前阴影区域和当前非阴影区域对应的原始图像中的像素区域分别作为步骤S1-1最终分割得到的阴影区域和非阴影区域，进入步骤S1-17；若否，则将当前分类阈值更新为新分类阈值，返回步骤S1-13；
S1-17：更新标记矩阵Mark中的元素值，得到更新后的标记矩阵MarkⅠ，具体包括：将标记矩阵Mark中与最终分割得到的阴影区域位置相同的元素值置3，将标记矩阵Mark中与最终分割得到的非阴影区域位置相同的元素值置0，得到更新后的标记矩阵MarkⅠ。紫色土土壤图像的阴影区域相比正常光照区域具有平均亮度低的特点，因此，阴影区域的像素满足像素的亮度分量低于图像的平均亮度。根据阴影区域与正常光照区域的亮度特点，可以采用阈值分割的方法提取土壤图像中的阴影区域，而图像的整体亮度均值反应了图像的光照强度，因此，将图像的平均亮度视为阴影区域的分割阈值。当图像阴影区域的像素点较少时，利用图像的平均亮度作为紫色土土壤图像的分割阈值就能简单有效的分割提取图像的阴影区域，但对于图像中存在大量阴影区域时，图像的整体亮度均值较低，所得分割阈值不能准确有效的分割阴影区域。因此，上述方法中，采用类间方差迭代更新算法求取图像阴影区域的分割阈值，能够很好地适应各种阴影区域占比情况的紫色土土壤图像的阴影区域分割提取要求，准确地将紫色土土壤图像中的非阴影区域和阴影区域分割检测出来。
进一步，所述步骤S1-15中新分类阈值Tnew的获得公式为：
其中，表示使得f(X)取得最大值所对应变量X的值；和分别为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的两类元素的类间方差和类内方差；新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的两类元素具体为：新分类阈值Tnew作为分类阈值，重新将I域值矩阵中元素值小于Tnew和元素值大于或等于Tnew的元素划分为两类元素。通过上述计算公式，能够得到适应能力强的新分类阈值，以适应包含不同占比的阴影区域的紫色土图像的分割阴影区域和非阴影区域的需求。
进一步，所述和的计算公式分别为：
其中，μ为I域值矩阵中所有元素值的的算术平均值；μ1'表示由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新阴影区的所有元素值的算术平均值，μ2'表示由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新非阴影区的所有元素值的算术平均值；新分类阈值Tnew作为分类阈值将I域值矩阵中元素值小于Tnew的元素组成的区域划分为新阴影区，将元素值大于或等于Tnew的元素组成的元素划分为新非阴影区；σ1表示由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新阴影区的所有元素值的标准差，σ2表示由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新非阴影区的所有元素值的标准差；w1和w2分别为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新阴影区和新非阴影区所包含的元素数量占I域值矩阵总元素数量的比例；N1为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新阴影区的所有元素数量，N2为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新非阴影区的所有元素数量。
具体地，σ1和σ2的计算公式为：
其中，Ii为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新阴影区中第i个元素值(I域值)，Ij为由新分类阈值Tnew作为分类阈值将I域值矩阵中元素划分得到的新非阴影区中第j个元素值(I域值)。通过上述方法，能够得到适应能力强的新分类阈值，以适应包含不同占比的阴影区域的紫色土图像的分割阴影区域和非阴影区域的需求。
进一步，所述步骤S1-2包括步骤：
S1-21：初始化与原始图像像素大小相同的类标签矩阵lab，用于存储标记矩阵MarkⅠ各个元素的类标签，类标签矩阵lab中各个位置元素值分别表示在标记矩阵MarkⅠ中与类标签矩阵lab位置对应的元素的类标签值；类标签矩阵lab的所有元素初始化为0；所述位置对应是指位置相同；
初始化元素数量阈值T；
S1-22：对类标签矩阵lab进行初始更新，得到初始更新后的类标签矩阵labI；
S1-23：对类标签矩阵lab I进行更新，得到第一次更新后的类标签矩阵labⅡ；
S1-24：对类标签矩阵labⅡ进行更新，得到第二次更新后的类标签矩阵labⅢ；
S1-25：对类标签矩阵labⅢ进行更新，得到第三次更新后的类标签矩阵labⅣ；
S1-26：对类标签矩阵labⅣ进行更新，得到第四次更新后的类标签矩阵labⅤ；步骤S1-23至S1-26是执行对合并类标签的操作，在减少类标签值数量同时，简化并方便了后续对离散阴影区域的判定。
S1-27：统计类标签矩阵labⅤ中各个非零类标签值对应的元素总数量；
S1-28：随机选取未被判断遍历的非零类标签值，判断该非零类标签值对应的元素总数量是否小于T，若是，则判定该非零类标签值对应的原始图像区域为离散阴影区域，将标记矩阵MarkⅠ中离散阴影区域对应的元素值更新为0，若否，则不做处理；
S1-29：重复步骤S1-28直到所有的非零类标签值均被判断遍历，得到更新后的标记矩阵MarkⅡ，其中，原始图像中与标记矩阵MarkⅡ中元素值为3的元素位置相同的元素组成区域为保留的阴影区域。通过步骤S1-1对紫色土土壤图像分割后，分割得到土壤的阴影区域和非阴影区域。在本实施例中，阴影区域的元素不仅可以赋值为3，还可以赋值为其他值，只要能与非阴影区域进行区别即可。由于自然条件下拍摄的土壤图像中存在大量的自然断口，导致分割后阴影区域中存在较多的离散阴影区域。离散阴影区域具有面积小，分布不连续的特点，对离散阴影区域增强将导致增强后土壤图像存在大量的伪边界。同时，由于离散阴影区域面积小，对土壤类型识别的干扰较小，因此，可以将离散阴影区域中离散区域剔除，以方便后续对紫色土土壤的土种和土属的识别，同时，也能提高后续对紫色土土壤的土种和土属的识别精度。
进一步，所述步骤S1-22包括步骤：
S1-221：按照预设的标签初始更新遍历顺序Ⅰ，选择一个未经初始更新遍历的标记矩阵MarkⅠ中的元素，作为当前初始更新遍历元素；其中，预设的标签初始更新遍历顺序Ⅰ为：以标记矩阵MarkⅠ中左上角的元素作为起始遍历元素，从上到下逐行遍历标记矩阵MarkⅠ中的元素；在遍历每一行元素时，按照从左到右的顺序遍历每一行的元素；
S1-222：判断当前初始更新遍历元素值是否为3，若是，则进入步骤S1-223；若否，则进入步骤S1-224；
S1-223：判断标记矩阵MarkⅠ中当前初始更新遍历元素的八连通邻域点中是否存在类标签矩阵lab中与当前初始更新遍历元素的八连通邻域点位置对应的元素值(类标签值)大于0的点，若存在，则将lab矩阵中当前初始更新遍历元素的八连通邻域点中类标签值最小的非零元素值作为当前初始更新遍历元素的类标签值，将类标签矩阵lab当前初始更新遍历元素对应位置的元素值更新为当前初始更新遍历元素的类标签值，进入步骤S1-224；若不存在，则将类标签矩阵lab中值最大的元素值增加1，作为当前初始更新遍历元素的类标签值，将类标签矩阵lab当前初始更新遍历元素对应位置的元素值更新为当前初始更新遍历元素的类标签值，进入步骤S1-224；
S1-224：判断标记矩阵MarkⅠ中是否还有未经初始更新遍历的元素，若有，则返回步骤S1-221，若无，则结束对标记矩阵MarkⅠ的初始更新遍历，得到初始更新后的类标签矩阵lab I。例如：一个300×300像素的原始图像，标记矩阵MarkⅠ也是300×300的大小，在类标签矩阵lab中，当初始更新遍历到第2行第4列的类标签值为3时，第2行第4列的元素的八连通邻域点(位于第1行第3列的元素、位于第1行第4列的元素、位于第1行第5列的元素、位于第2行第3列的元素、位于第2行第5列的元素、位于第3行第3列的元素、位于第3行第4列的元素、位于第3行第5列的元素)中，位于第1行第5列的元素值的类标签值为3，位于第3行第3列的类标签值为5，将位于第1行第5列的类标签值3作为第2行第4列的元素的类标签值存储在类标签矩阵lab中第2行第4列的元素中，即类标签矩阵lab中第2行第4列的元素值更新为3。再例如：一个300×300像素的原始图像，标记矩阵MarkⅠ也是300×300的大小，，在类标签矩阵lab中，当初始更新遍历到第2行第4列的类标签值为3时，第2行第4列的元素的八连通邻域点(位于第1行第3列的元素、位于第1行第4列的元素、位于第1行第5列的元素、位于第2行第3列的元素、位于第2行第5列的元素、位于第3行第3列的元素、位于第3行第4列的元素、位于第3行第5列的元素)的类标签值均为0，此时，类标签矩阵lab中最大的元素值(类标签值)为15，则将类标签矩阵lab中第2行第4列的元素置为16。通过上述方法，用类标签将阴影区域的各个像素点进行分类。
进一步，所述步骤S1-23具体包括：
S1-231：按照预设的更新遍历顺序Ⅱ，选择类标签矩阵lab I中一个非零元素作为当前遍历元素；其中，更新遍历顺序Ⅱ为：以类标签矩阵lab I中左上角的元素作为起始遍历元素，从上至下逐行遍历；在遍历每一行元素时，按照从左到右的顺序逐个遍历元素；
S1-232：判断类标签矩阵lab I中当前遍历元素的八连通邻域点中是否存在满足条件1的点：若存在，则将类标签矩阵lab I中当前遍历元素值更新为当前遍历元素的八连通邻域点中值非零的最小值，然后进入步骤S1-233；若不存在，则进入步骤S1-233；其中，条件1为：当前遍历元素的八连通邻域点的元素值小于当前遍历元素值，且值非零；
S1-233：判断类标签矩阵lab I中是否还有未被更新遍历的元素，若存在，则返回步骤S1-231，若不存在，则对类标签矩阵lab I进行第一次更新结束，得到第一次更新后类标签矩阵labⅡ；
所述步骤S1-24具体包括：
S1-241：按照预设的更新遍历顺序Ⅱ，选择类标签矩阵labⅡ中一个非零元素作为当前遍历元素；其中，更新遍历顺序Ⅱ为：以类标签矩阵labⅡ中右上角的元素作为起始遍历元素，从上至下逐行遍历；在遍历一行时，按照从右到左的顺序逐个遍历元素；
S1-242：判断类标签矩阵labⅡ中当前遍历元素的八连通邻域点中是否存在满足条件2的点：若存在，则将类标签矩阵labⅡ中当前遍历元素值更新为其八连通邻域点中值非零的最小值，然后进入步骤S1-243；若不存在，则进入步骤S1-243；其中，条件2为：当前遍历元素的八连通邻域点的元素值小于当前遍历元素值，且值非零；
S1-243：判断类标签矩阵labⅡ中是否还有未被更新遍历的元素，若存在，则返回步骤S1-241，若不存在，则对类标签矩阵labⅡ进行第二次更新结束，得到第二次更新后类标签矩阵labⅢ；
所述步骤S1-25具体包括：
S1-251：按照预设的更新遍历顺序Ⅲ，选择类标签矩阵labⅢ中一个非零元素作为当前遍历元素；其中，更新遍历顺序Ⅲ为：以类标签矩阵labⅢ中左下角的元素作为起始遍历元素，从下至上逐行遍历；在遍历一行时，按照从左到右的顺序逐个遍历元素；
S1-252：判断类标签矩阵labⅢ中当前遍历元素的八连通邻域点中是否存在满足条件3的点：若存在，则将类标签矩阵labⅢ中当前遍历元素值更新为其八连通邻域点中值非零的最小值，然后进入步骤S1-253；若不存在，则进入步骤S1-253；其中，条件3为：当前遍历元素的八连通邻域点的元素值小于当前遍历元素值，且值非零；
S1-253：判断类标签矩阵labⅢ中是否还有未被更新遍历的元素，若存在，则返回步骤S1-251，若不存在，则对类标签矩阵labⅢ进行第三次更新结束，得到第三次更新后类标签矩阵labⅣ；
所述步骤S1-26具体包括：
S1-261：按照预设的更新遍历顺序Ⅳ，选择类标签矩阵labⅣ中一个非零元素作为当前遍历元素；其中，更新遍历顺序Ⅳ为：以类标签矩阵labⅣ中右下角的元素作为起始遍历元素，从下至上逐行遍历；在遍历一行时，按照从右到左的顺序逐个遍历元素；
S1-262：判断类标签矩阵labⅣ中当前遍历元素的八连通邻域点中是否存在满足条件4的点：若存在，则将类标签矩阵labⅣ中当前遍历元素值更新为其八连通邻域点中值非零的最小值，然后进入步骤S1-263；若不存在，则进入步骤S1-263；其中，条件4为：当前遍历元素的八连通邻域点的元素值小于当前遍历元素值，且值非零；
S1-263：判断类标签矩阵labⅣ中是否还有未被更新遍历的元素，若存在，则返回步骤S1-261，若不存在，则对类标签矩阵labⅣ进行第四次更新结束，得到第四次更新后类标签矩阵labⅤ。通过上述方法，用类标签将阴影区域的各个像素点进行分类，然后再合并相互连通(相互连通是指相互为对方的八连通邻域点)的像素，形成连通区域，再判断各连通区域(各连通区域用类标签进行区分)的大小是否小于元素数量阈值T，若小于则可判定该连通区域是离散阴影区域，这样可以精确地将离散阴影区域识别并去除。例如：将初始化元素数量阈值T为20，也就是连通区域的元素总数量少于20，则将该连通区域作为离散阴影区域并去除。
进一步，所述步骤S1-3包括步骤：
S1-31：计算步骤S1-11得到的I域值矩阵中保留的阴影区域对应位置的元素值的算术平均值μs，将μs作为分割全影区域和半影区域的初始分类阈值T0’；将T0’作为当前分类阈值；
S1-32：将I域值矩阵中的保留的阴影区域对应位置区域中，元素值小于当前分类阈值的元素划分为当前全影区域，将I域值矩阵中元素值大于或等于当前分类阈值的元素划分为当前半阴影区域；
S1-33：分别计算当前半阴影区域的算术平均值和当前全阴影区域的元素值的算术平均值；
S1-34：在当前全阴影区域的元素值的算术平均值与当前半阴影区域的算术平均值的之间范围内，获得新分类阈值Tnew’，即μ全＜Tnew'＜μ半；其中，μ全和μ半分别为当前全阴影区域的所有元素值的算术平均值和当前半阴影区域的所有元素值的算术平均值；
S1-35：判断新分类阈值Tnew’是否等于当前分类阈值；若是，停止迭代，将最后一次迭代中的当前半阴影区域和当前全阴影区域对应的原始图像中的像素区域分别作为步骤S3最终分割得到的半阴影区域和全阴影区域；若否，则将当前分类阈值更新为新分类阈值，返回步骤S32。彩色图像中的像素点的亮度由两种光照决定，在图像的非阴影区域像素，其亮度主要由直接光照和间接光照决定；在图像阴影区域的像素点，其亮度主要由部分直接光照和间接光照决定。通过对阴影区域中像素点所受直接光照的量进行判断，将阴影区域像素分为全影区域和半影区域。当像素点所受直接光照为0时，该像素点即可视为全影区域；当像素点受到部分直接光照时，该像素点即可视为半影区域。在半影区域中的像素点由于受到部分直接光照，因此，半影区像素点的亮度是介于非阴影区域和全影区域的。将全影区域像素和半影区域像素视为两类，通过像素点亮度求取分割阈值，实现阴影区域像素的分类，方便后续进一步对紫色土土壤的土种和土属的识别。
具体地，根据图像的光照原理，一幅图像可以视为由反射分量与照度分量构成，反射分量代表了图像中像素点的内在属性，而照度分量可以视为图像所受的外部光照。在自然光照条件下拍摄的土壤图像中，每个像素点的照度分量可以看作是受直接光照和间接光照的影响，其中直接光照可以视为太阳光的照射，间接光照可以视为物体表面光线的漫反射。因此，构成了图像的基本光照模型如下：
I(x,y)＝L(x,y)×R(x,y) (14)
其中，I(x,y)表示原始图像，R(x,y)表示反射图像，L(x,y)表示照度图像。将照度分量分解为直接照度分量和间接照度分量，得到下式：
I(x,y)＝[Ld(x,y)+Le(x,y)]×R(x,y) (15)
其中，Ld(x,y)和Le(x,y)分别为直接照度图像和间接照度图像。
进一步，所述步骤S1-34中新分类阈值Tnew’的获得公式为：
其中，表示使得f(X)取得最大值所对应变量X的值；μ全为当前全阴影区域的元素值的算术平均值，μ半为当前半阴影区域的算术平均值，和分别为由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域元素划分得到的两类元素的类间方差和类内方差；新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域元素划分得到的两类元素具体为：新分类阈值Tnew’作为分类阈值，重新将I域值矩阵中保留的阴影区域中元素值小于Tnew’和元素值大于或等于Tnew’的元素划分为两类元素。
进一步，所述和的计算公式分别为：
其中，μ'为I域值矩阵中保留的阴影区域中所有元素值的算术平均值；μ全'表示由新分类阈值Tnew’作为分类阈值将I域值矩阵中元素划分得到的全阴影区的所有元素值的算术平均值，μ半'表示由新分类阈值Tnew’作为分类阈值将I域值矩阵中元素划分得到的半阴影区的所有元素值的算术平均值；新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素值小于Tnew’的元素组成的区域划分为新全阴影区，将元素值大于或等于Tnew’元素组成的元素划分为新半阴影区；σ全表示由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素划分得到的新全阴影区的所有元素值的标准差，σ半表示由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素划分得到的新半阴影区的所有元素值的标准差；w全和w半分别为由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素划分得到的新全阴影区和新半阴影区所包含的元素数量占I域值矩阵中保留的阴影区域总元素数量的比例； N全为由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素划分得到的新全阴影区的所有元素数量，N半为由新分类阈值Tnew’作为分类阈值将I域值矩阵中保留的阴影区域中元素划分得到的新半阴影区的所有元素数量。由于w全+w半＝1，则Tnew'始终满足μ全＜Tnew'＜μ半。通过上述计算方法，可以得到适应能力强的新分类阈值，可以适应各种全阴影区域和半阴影区域占比的阴影区域，以精确划分全阴影区域和半阴影区域，方便后续对紫色土土壤的土种和土属的准确识别。
回到本文方法，在执行了步骤S1后，还要执行如下步骤：
S2：将紫色土土壤的原始图像转换为HSI颜色模型的H域值矩阵H0；其中，矩阵H0对应于原始图像相同位置的像素的H域值；
将紫色土土壤的原始图像转换为HSI颜色模型的I域值矩阵I0；其中，矩阵I0对应于原始图像相同位置的像素的I域值；
S3：按1°整化统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图Hist_l、Hist_f和Hist_p；其中，Hist_l为非阴影区域中像素点H分量的直方图，Hist_f为全影区域中像素点H分量的直方图，Hist_p为半影区域中像素点H分量的直方图；其中，按1°整化，统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图是指将紫色土土壤的原始图像的H分量等角度划分为360等分，每个等分对应一个灰度级；本实施例中，H分量的取值范围一共对应360°，按1°整化将紫色土土壤的原始图像的H分量等角度划分为360个等分，每个等分对应1°。本实施例中，[0°，1°)对应一个等分，[1°，2°)对应一个等分，[2°，3°)对应一个等分，依次类推按1°整化得到360个等分。本实施例中，统计得到图像区域的像素点H分量的直方图为现有技术在此不赘述。
S4：初始化的标记向量Hist；其中，标记向量Hist初始化为360维的0向量，标记向量Hist中的每一维分别对应为按1°整化得到各等分对应的灰度级的标签；本实施例中，按1°整化，将紫色土土壤的原始图像的H分量等角度划分为360个等分(也就是360个色度(灰度级))，故设置360维的标签向量。
S5：依次遍历向量Hist_l、Hist_f和Hist_p中的各个灰度级，若灰度级对应的向量Hist_l、Hist_f和Hist_p的值都大于零，则更新该灰度级对应的标记向量中的标签值为1；
S6：重复步骤S5，直到所有的向量Hist_l、Hist_f和Hist_p中的各个灰度级都被遍历，得到更新后的标记向量HistⅠ；
S7：计算非阴影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μl；计算全影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μf；计算半影区域中标记向量HistⅠ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μp；例如：灰度级为2和69两个灰度级对应的标记向量HistⅠ中标签值为1；非阴影区域中灰度级为2的像素点有60个，灰度级为69的像素点有40个，则μl为这60+40＝100个像素点亮度的算术平均值，其中，像素点亮度可以直接从矩阵I0中直接读取这些像素点对应位置的元素值得到；同理，μf和μp也采用相同的方法获得。
S8：根据μl、μf和μp，计算直接光照和间接光照的比值r、半影区整体的直接光照的遮蔽系数α和全局直接光照强度Id；
S9：根据r、α和Id，计算半影区域和全影区域各像素点增强后的亮度值，并将矩阵I0中各个像素点对应位置的元素更新为增强后的亮度。其中，所述增强是指亮度增强，得到阴影增强后的紫色土土壤图像。通过上述方法，对分割出紫色土土壤图像的半阴影区域和全阴影区域，进行亮度增强，以方便后续对紫色土土属和土种的准确识别；
S10：对阴影区域亮度增强后的紫色土土壤的初始图像进行色温校正，得到色温校正后的紫色土土壤图像；
S11：对色温校正后的紫色土土壤图像进行局部对比度校正，得到局部对比度校正后的紫色土土壤图像，即得到最终的阴影区域亮度增强后的紫色土土壤图像。
进一步，所述步骤S9中，半影区域和全影区域各像素点增强后的亮度值的计算公式为：
其中，Ip non_shadow(x,y)为增强后的半影区像素点(x,y)的亮度值，Ip(x,y)为增强前半影区像素点(x,y)的亮度值，If non_shadow(x,y)为增强后的全影区像素点(x,y)的亮度值，If(x,y)为增强前全影区像素点(x,y)的亮度值，αi为像素点(x,y)的直接光照遮蔽系数，r为直接光照与间接光照的比值。
进一步，所述像素点(x,y)的直接光照遮蔽系数αi的计算公式为：
其中，μM×M表示像素点(x,y)的M×M邻域内像素点亮度算术平均值；例如：求像素点(6,6)的3×3邻域内像素点亮度算术平均值μ3×3，即求像素点(5,5)、(5,6)、(5,7)、(6,5)、(6,6)、(6,7)、(7,5)、(7,6)和(7,7)九个像素点亮度的算术平均值，所述像素点亮度可直接矩阵I0中直接读取这些像素点对应位置的元素值得到。
所述Id的计算公式为：
所述α的计算公式为：
进一步，所述r的计算公式为：
在紫色土土壤图像中，相邻土壤像素所受光照强度是相似的，因此，将像素的相邻点平均光照强度作为当前点的像素光照强度，来增强阴影区域的亮度。通过上述的计算公式，能够考虑到紫色土土壤的原始图像的实际光照情况，来增强紫色土土壤图像中阴影区域的亮度，以避免失真，提高了增强阴影区域亮度的精确度。
上述公式的推导过程如下：
令紫色土土壤图像中非阴影区域像素点集合为U，半影区像素点集合为P，全影区像素点集合为F。分别统计得到图像非阴影区域、全影区和半影区中像素点H分量的直方图Hist_l、Hist_f和Hist_p，分别计算直方图具有相同灰度级i时亮度的算术平均值：算术平均值的计算公式为现有技术，在此不赘述。
根据彩色图像的光照模型[(14)至(19)式]，半影区像素点可视为接受部分直接光照，全影区像素点可视为无直接光照，因此，根据像素点所处位置可建立像素点的光照模型如下：
将上式中非阴影区域亮度均值μl分别与半影区亮度均值μp、全影区亮度均值μf相除，得到直接光照Ld、间接光照Le与亮度均值的关系：
令直接光照与间接光照的比值可得：
其中，α表示半影区整体的直接光照的遮蔽系数。
根据半影区像素点亮度均值μp与非阴影区域像素点亮度均值μl，可以估算出整幅图像的直接光照强度Id：
由光照模型可知，半影区像素点亮度增强后与增强前的比例为：
全影区像素点亮度增强后和增强前的比例为：
其中，Ip non_shadow(x,y)为增强后的半影区像素点亮度值，Ip(x,y)为增强前半影区像素亮度值，If non_shadow(x,y)为增强后的全影区像素点亮度值，If(x,y)为增强前全影区像素点亮度值，αi表示像素点(x,y)的直接光照遮蔽系数。
进一步化简得到阴影增强公式：
令得到：
同理，得到全影区增强公式：
对于(27)中像素点直接光照遮蔽系数αi，计算公式如下：
进一步，所述步骤S10包括步骤：
S101：将最终的阴影区域亮度增强后的紫色土土壤的初始图像转换为Lab颜色模型的a分量矩阵A；其中，矩阵A对应于原始图像相同位置的像素的a分量值；
将最终的阴影区域亮度增强后的紫色土土壤图像转换为Lab颜色模型的b分量矩阵B；其中，矩阵B对应于原始图像相同位置的像素的b分量值；
S102：分别统计得到阴影区域a分量和b分量的统计直方图Hist_as和Hist_bs，其中阴影区域由全阴影区域和半阴影区域组成；分别统计得到非阴影区域a分量和b分量的统计直方图Hist_al和Hist_bl；
S103：分别对阴影区域a分量和b分量的统计直方图进行均衡化，得到阴影区域直方图向量Equal_as和Equal_bs；
分别对非阴影区域a分量和b分量的统计直方图进行均衡化，得到非阴影区直方图向量Equal_al和Equal_bl；
S104：按照单映射规则，分别寻找向量Equal_as中各个灰度级在Equal_al中最接近的灰度级，并将Equal_as中各个灰度级对应的a分量的值ai分别映射到Equal_al中与Equal_as中各个灰度级最接近的灰度级对应矩阵A中的元素，得到新的a分量矩阵Anew；本实施例中，步骤S104中的映射是指替换；将Equal_as中各个灰度级对应的a分量的值ai映射到Equal_al中与Equal_as中最接近的灰度级是指将Equal_as中各个灰度级对应的a分量的值ai替换给Equal_al中与Equal_as中最接近的灰度级；例如：Equal_as中灰度级i在Equal_al中最接近灰度级j，则将灰度级i对应的a分量的值ai映射到灰度级j对应矩阵A中的元素(使灰度级j对应矩阵A中的元素值为ai)；
按照单映射规则，分别寻找向量Equal_bs中各个灰度级在Equal_bl中最接近的灰度级，并将Equal_bs中各个灰度级对应的b分量的值bi分别映射到Equal_bl中与Equal_bs中各个灰度级最接近的灰度级对应矩阵B中的元素，得到新的b分量矩阵Bnew，即得到色温校正后的紫色土土壤图像。将Equal_bs中各个灰度级对应的b分量的值bi分别映射到Equal_bl中与Equal_bs中各个灰度级最接近的灰度级对应矩阵B中的元素是指将Equal_bl中各个灰度级对应的b分量的值bi替换给Equal_bl中与Equal_bs中最接近的灰度级；例如：Equal_bs中灰度级i在Equal_bl中最接近灰度级j，则将灰度级i对应的b分量的值bi映射到灰度级j对应矩阵A中的元素(使灰度级j对应矩阵B中的元素值为bi)；
经亮度增强后的图像阴影区域，在色温上与非阴影区域存在较大的差别，导致增强后图像失真。因此，在Lab颜色空间中，以非阴影区域a分量与b分量的统计直方图为基础，对增强后阴影区域a分量与b分量进行调整，使得增强后土壤图像色温保持一致。
进一步，步骤S104中单映射规则为：
所述Equal_as、Equal_al、Equal_bs、Equal_bl的计算公式如下：
其中，asj为阴影区域中第j个灰度级的a分量，m为阴影区域中灰度级的总数量，nsj为阴影区域中a分量值为asj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(asj)为阴影区域中a分量值为asj的概率；
其中，ali为非阴影区域中第i个灰度级的a分量，n为非阴影区域中灰度级的总数量，kli为非阴影区域中a分量值为ali的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(ali)为阴影区域中a分量值为ali的概率；
其中，bsj为阴影区域中第j个灰度级的b分量，m为阴影区域中灰度级的总数量，nsj'为阴影区域中b分量值为bsj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(bsj)为阴影区域中b分量值为bsj的概率；
其中，bli为非阴影区域中第i个灰度级的b分量，n为非阴影区域中灰度级的总数量，kli'为非阴影区域中b分量值为bli的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(bli)为阴影区域中b分量值为bli的概率。通过上述单映射规则能够获得量Equal_as中各个灰度级在Equal_al中最接近的灰度级和向量Equal_bs中各个灰度级在Equal_bl中最接近的灰度级，方便对阴影区域亮度增强后的紫色土土壤的初始图像进行色温校正，以避免增强阴影区域的亮度后紫色土土壤图像失真。
进一步，所述步骤S11包括步骤：
S111：将最终的阴影区域亮度增强后的紫色土土壤图像转换为HSI颜色模型的I域值矩阵I2；其中，矩阵I2对应于原始图像相同位置的像素的I域值；
S112：计算色温增强后的紫色土土壤图像中每个像素点的局部对比度；
S113：根据局部对比度对矩阵I2中的所有元素进行初步对比度校正，得到进行初步对比度校正后的矩阵I3；
S114：初始化与矩阵I3大小相同的异常标记矩阵Markerror；
S115：选取任一未被遍历的矩阵I3中的元素作为当前遍历元素，判断当前遍历元素值是否大于非阴影区域的像素点的算术平均值，若是，则判断当前遍历元素为亮度初始对比度校正后的异常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置1；若否，则判断当前遍历元素为亮度初始对比度校正后的正常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置0；初步对比度校正后的像素点可能存在过增强的情况，即增强后阴影区域像素点亮度值大于非阴影区域像素点的平均亮度，可认为该点为亮度异常点。
S116：重复步骤S115，直到矩阵I3中所有的元素都被遍历，得到更新后的异常标记矩阵MarkerrorI；
S117：随机选取异常标记矩阵MarkerrorI中未被补充校正遍历且值为1的元素，作为当前补充校正遍历元素，计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值，并将矩阵I3中该元素值更新为该补充校正值；
S118：判断异常标记矩阵MarkerrorI中是否还有未被遍历且值为1的元素，若有，则返回步骤S117，若无，则结束对异常标记矩阵MarkerrorI的补充校正遍历，得到补充校正后的I域值矩阵I4，即得到局部对比度校正后的紫色土土壤图像，也就是得到最终的阴影区域亮度增强后的紫色土土壤图像。对土壤图像阴影区域像素点亮度进行增强并校正色温后，在图像的局部区域内，像素点亮度会增强不同的倍数，而图像的阴影区域的原始梯度会随之产生变化，导致图像阴影区域的对比度改变。因此，在图像阴影区域增强后，需要按照阴影区域的原始对比度进行土壤图像阴影区域进行校正。
进一步，所述步骤S112中局部对比度的计算公式为：
其中，Cg×g(x,y)表示色温增强后的紫色土土壤图像中像素点(x,y)的局部对比度，Iold(x,y)表示紫色土土壤的原始图像中像素点(x,y)的亮度值；μold g×g(x,y)表示紫色土土壤的原始图像中像素点(x,y)的g×g邻域的g×g个像素点亮度的算术平均值；在土壤图像阴影区域中，像素点的邻域亮度均值随图像增强而产生了变化，为保证增强后像素点的亮度依然满足原始梯度关系，增强后的局部对比度应当保持不变，因此通过(9)式计算局部对比度。
所述步骤S113中根据局部对比度对矩阵I2中的所有元素进行初步对比度校正的公式为：
Inew(x,y)＝Cg×g(x,y)×μnew g×g(x,y) (12)
其中，Inew(x,y)表示初步对比度校正后的矩阵I3中元素(x,y)的值；μnew g×g(x,y)表示矩阵I2中元素(x,y)的g×g邻域的g×g个元素值的算术平均值。
进一步，所述步骤S117中计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值的公式为：
其中，INEW(x,y)表示矩阵I3中元素(x,y)的补充校正值，表示紫色土土壤的原始图像中非阴影区域的所有像素点亮度的算术平均值；表示紫色土土壤的原始图像中阴影区域的所有像素点亮度的算术平均值；表示阴影区域亮度增强后的紫色土土壤的初始图像中阴影区域的所有像素点亮度的算术平均值；Ierror(x,y)表示矩阵I3中当前补充校正遍历元素对应的元素值。由于增强后图像中非阴影区域的亮度均值未发生改变，因此，计算异常点的原始亮度相对于非阴影区域的亮度均值以及阴影区域原亮度均值的比值，通过比值对增强后异常亮度点进行校正，校正公式为(11)式。
在野外自然环境下用土锹锹出耕层0～20cm左右的紫色土，拍摄其无锹痕的(心土)自然断口图像(心土能最大限度地保持紫色土自然颜色和原状土壤结构)100张，通过手动分割提取无背景杂质的土壤图像区域，作为算法实验图像样本。
随机选取4张紫色土土壤区域分割结果图像作为一组实验对象，共计25组紫色土土壤图像。选取文献1中的算法为对比算法1；选取文献2中的算法为对比算法2。其中，文献1：Anoopa S,Dhanya V,Kizhakkethottam J J.Shadow Detection and Removal UsingTri-Class Based Thresholding and Shadow Matting Technique[J].ProcediaTechnology,2016,24:1358-1365；文献2：Silva G F,Carneiro G B,Doth R,et al.Nearreal-time shadow detection and removal in aerial motion imagery application[J].ISPRS Journal of Photogrammetry and Remote Sensing,2017,11:005。仿真实验环境为：CPU为Intel(R)Core(TM)i5 8300H，2.30GHz；内存为8GB，显卡为GTX1050Ti，操作系统为Windows 10专业版，VC++2015与OpenCV3.4。对所有25组实验图像进行实验，所得实验结果是类似的。图2为本发明的实验效果对比图，其中，图2第一行的图是一组实验对象的四个原始图像，图2第二行的图是第一行的图经对比算法1增强阴影区域亮度后的图，图2第三行的图是第一行的图经对比算法2增强阴影区域亮度后的图，图2第四行的图是第一行的图经本文算法中步骤S1至S9增强阴影区域亮度后的图；图2第五行的图是第四行的图经本文算法中步骤S10色温校正后的图；图2第六行的图是第五行的图经本文算法中步骤S11局部对比度校正后的图。
实验的数据结果如表1和表2所示。
表1半影区和全影区亮度均值与标准差
表2阴影增强算法增强效果及时间花销
实验结果显示对比算法1中基于三类阈值和阴影抠图技术的阴影检测与去除算法对土壤图像增强效果不佳，平均增强倍数为1.75，增强后阴影区域亮度值依然远低于非阴影区域亮度值；对比算法2中航空运动图像中的近实时阴影检测与去除算法增强土壤图像阴影区域时，平均增强倍数为3.31，增强后阴影区域亮度均值大于非阴影区域，导致图像整体失真严重。本文阴影区域增强算法的平均增强倍数为2.62，增强后阴影区域亮度均值与非阴影区域接近。
本文阴影增强算法平均时间花销为6.49s，相比于对比算法时间花销较大，是因其统计相同色度像素点在阴影区域与非阴影区域的亮度值，通过亮度比值求取直接亮度强度以及遮蔽系数，并对增强后图像进行色温校正以及局部对比度校正，从而增加了算法的时间花销。但通过本文算法增强后土壤图像，阴影区域亮度与非阴影区域近似，同时失真较小。
综合表1和表2实验结果数据，本文阴影增强算法增强后图像阴影区域亮度与非阴影区域一致，增强结果符合人类视觉特性。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (9)

1.一种紫色土土壤图像阴影增强方法，其特征在于：包括步骤：
S1：将紫色土土壤的原始图像划分为非阴影区域、半阴影区域和全阴影区域；
S2：将紫色土土壤的原始图像转换为HSI颜色模型的H域值矩阵H0；其中，矩阵H0对应于原始图像相同位置的像素的H域值；
将紫色土土壤的原始图像转换为HSI颜色模型的I域值矩阵I0；其中，矩阵I0对应于原始图像相同位置的像素的I域值；
S3：按1°整化统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图Hist_l、Hist_f和Hist_p；其中，Hist_l为非阴影区域中像素点H分量的直方图，Hist_f为全影区域中像素点H分量的直方图，Hist_p为半影区域中像素点H分量的直方图；其中，按1°整化统计得到非阴影区域、全影区域和半影区域中像素点H分量的直方图是指将紫色土土壤的原始图像的H分量等角度划分为360等分，每个等分对应一个灰度级；
S4：初始化的标记向量Hist；其中，标记向量Hist初始化为360维的0向量，标记向量Hist中的每一维分别对应为按1°整化得到各等分对应的灰度级的标签；
S5：依次遍历向量Hist_l、Hist_f和Hist_p中的各个灰度级，若灰度级对应的向量Hist_l、Hist_f和Hist_p的值都大于零，则更新该灰度级对应的标记向量中的标签值为1；
S6：重复步骤S5，直到所有的向量Hist_l、Hist_f和Hist_p中的各个灰度级都被遍历，得到更新后的标记向量Hist Ⅰ；
S7：计算非阴影区域中标记向量Hist Ⅰ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μl；计算全影区域中标记向量Hist Ⅰ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μf；计算半影区域中标记向量Hist Ⅰ中标签值为1的所有灰度级对应的所有像素点亮度的算术平均值μp；
S8：根据μl、μf和μp，计算直接光照和间接光照的比值r、半影区整体的直接光照的遮蔽系数α和全局直接光照强度Id；
S9：根据r、α和Id，计算半影区域和全影区域各像素点增强后的亮度值，并将矩阵I0中各个像素点对应位置的元素更新为增强后的亮度，得到阴影增强后的紫色土土壤图像；
S10：对阴影区域亮度增强后的紫色土土壤的初始图像进行色温校正，得到色温校正后的紫色土土壤图像；
S11：对色温校正后的紫色土土壤图像进行局部对比度校正，得到局部对比度校正后的紫色土土壤图像，即得到最终的阴影区域亮度增强后的紫色土土壤图像。
2.根据权利要求1所述紫色土土壤图像阴影增强方法，其特征在于：所述步骤S9中，半影区域和全影区域各像素点增强后的亮度值的计算公式为：
其中，Ip non_shadow(x,y)为增强后的半影区像素点(x,y)的亮度值，Ip(x,y)为增强前半影区像素点(x,y)的亮度值，If non_shadow(x,y)为增强后的全影区像素点(x,y)的亮度值，If(x,y)为增强前全影区像素点(x,y)的亮度值，αi为像素点(x,y)的直接光照遮蔽系数，r为直接光照与间接光照的比值。
3.根据权利要求2所述紫色土土壤图像阴影增强方法，其特征在于：所述像素点(x,y)的直接光照遮蔽系数αi的计算公式为：
其中，μM×M表示像素点(x,y)的M×M邻域内像素点亮度算术平均值；
所述Id的计算公式为：
所述α的计算公式为：
4.根据权利要求3所述紫色土土壤图像阴影增强方法，其特征在于：所述r的计算公式为：
5.根据权利要求1所述紫色土土壤图像阴影增强方法，其特征在于：所述步骤S10包括步骤：
S101：将最终的阴影区域亮度增强后的紫色土土壤的初始图像转换为Lab颜色模型的a分量矩阵A；其中，矩阵A对应于原始图像相同位置的像素的a分量值；
将最终的阴影区域亮度增强后的紫色土土壤图像转换为Lab颜色模型的b分量矩阵B；其中，矩阵B对应于原始图像相同位置的像素的b分量值；
S102：分别统计得到阴影区域a分量和b分量的统计直方图Hist_as和Hist_bs，其中阴影区域由全阴影区域和半阴影区域组成；分别统计得到非阴影区域a分量和b分量的统计直方图Hist_al和Hist_bl；
S103：分别对阴影区域a分量和b分量的统计直方图进行均衡化，得到阴影区域直方图向量Equal_as和Equal_bs；
分别对非阴影区域a分量和b分量的统计直方图进行均衡化，得到非阴影区直方图向量Equal_al和Equal_bl；
S104：按照单映射规则，分别寻找向量Equal_as中各个灰度级在Equal_al中最接近的灰度级，并将Equal_as中各个灰度级对应的a分量的值ai分别映射到Equal_al中与Equal_as中各个灰度级最接近的灰度级对应矩阵A中的元素，得到新的a分量矩阵Anew；
按照单映射规则，分别寻找向量Equal_bs中各个灰度级在Equal_bl中最接近的灰度级，并将Equal_bs中各个灰度级对应的b分量的值bi分别映射到Equal_bl中与Equal_bs中各个灰度级最接近的灰度级对应矩阵B中的元素，得到新的b分量矩阵Bnew，即得到色温校正后的紫色土土壤图像。
6.根据权利要求5所述紫色土土壤图像阴影增强方法，其特征在于：步骤S104中单映射规则为：
所述Equal_as、Equal_al、Equal_bs、Equal_bl的计算公式如下：
其中，asj为阴影区域中第j个灰度级的a分量，m为阴影区域中灰度级的总数量，nsj为阴影区域中a分量值为asj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(asj)为阴影区域中a分量值为asj的概率；
其中，ali为非阴影区域中第i个灰度级的a分量，n为非阴影区域中灰度级的总数量，kli为非阴影区域中a分量值为ali的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(ali)为阴影区域中a分量值为ali的概率；
其中，bsj为阴影区域中第j个灰度级的b分量，m为阴影区域中灰度级的总数量，nsj'为阴影区域中b分量值为bsj的像素点的总数量；Ns为阴影区域中像素点的总数量；Ps(bsj)为阴影区域中b分量值为bsj的概率；
其中，bli为非阴影区域中第i个灰度级的b分量，n为非阴影区域中灰度级的总数量，kli'为非阴影区域中b分量值为bli的像素点的总数量；Kl为非阴影区域中像素点的总数量；Pl(bli)为阴影区域中b分量值为bli的概率。
7.根据权利要求1所述紫色土土壤图像阴影增强方法，其特征在于：所述步骤S11包括步骤：
S111：将最终的阴影区域亮度增强后的紫色土土壤图像转换为HSI颜色模型的I域值矩阵I2；其中，矩阵I2对应于原始图像相同位置的像素的I域值；
S112：计算色温增强后的紫色土土壤图像中每个像素点的局部对比度；
S113：根据局部对比度对矩阵I2中的所有元素进行初步对比度校正，得到进行初步对比度校正后的矩阵I3；
S114：初始化与矩阵I3大小相同的异常标记矩阵Markerror；
S115：选取任一未被遍历的矩阵I3中的元素作为当前遍历元素，判断当前遍历元素值是否大于非阴影区域的像素点的算术平均值，若是，则判断当前遍历元素为亮度初始对比度校正后的异常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置1；若否，则判断当前遍历元素为亮度初始对比度校正后的正常点，将异常标记矩阵Markerror中当前遍历元素对应位置的元素置0；
S116：重复步骤S115，直到矩阵I3中所有的元素都被遍历，得到更新后的异常标记矩阵MarkerrorI；
S117：随机选取异常标记矩阵MarkerrorI中未被补充校正遍历且值为1的元素，作为当前补充校正遍历元素，计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值，并将矩阵I3中该元素值更新为该补充校正值；
S118：判断异常标记矩阵MarkerrorI中是否还有未被遍历且值为1的元素，若有，则返回步骤S117，若无，则结束对异常标记矩阵MarkerrorI的补充校正遍历，得到补充校正后的I域值矩阵I4，即得到局部对比度校正后的紫色土土壤图像，也就是得到最终的阴影区域亮度增强后的紫色土土壤图像。
8.根据权利要求7所述紫色土土壤图像阴影增强方法，其特征在于：所述步骤S112中局部对比度的计算公式为：
其中，Cg×g(x,y)表示色温增强后的紫色土土壤图像中像素点(x,y)的局部对比度，Iold(x,y)表示紫色土土壤的原始图像中像素点(x,y)的亮度值；μold g×g(x,y)表示紫色土土壤的原始图像中像素点(x,y)的g×g邻域的g×g个像素点亮度的算术平均值；
所述步骤S113中根据局部对比度对矩阵I2中的所有元素进行初步对比度校正的公式为：
Inew(x,y)＝Cg×g(x,y)×μnew g×g(x,y) (12)
其中，Inew(x,y)表示初步对比度校正后的矩阵I3中元素(x,y)的值；μnew g×g(x,y)表示矩阵I2中元素(x,y)的g×g邻域的g×g个元素值的算术平均值。
9.根据权利要求8所述紫色土土壤图像阴影增强方法，其特征在于：所述步骤S117中计算矩阵I3中当前补充校正遍历元素对应位置的元素的补充校正值的公式为：
其中，INEW(x,y)表示矩阵I3中元素(x,y)的补充校正值，表示紫色土土壤的原始图像中非阴影区域的所有像素点亮度的算术平均值；表示紫色土土壤的原始图像中阴影区域的所有像素点亮度的算术平均值；表示阴影区域亮度增强后的紫色土土壤的初始图像中阴影区域的所有像素点亮度的算术平均值；Ierror(x,y)表示矩阵I3中当前补充校正遍历元素对应的元素值。
CN201910705747.XA 2019-08-01 2019-08-01 紫色土土壤图像阴影增强方法 Active CN110428380B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910705747.XA CN110428380B (zh) 2019-08-01 2019-08-01 紫色土土壤图像阴影增强方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910705747.XA CN110428380B (zh) 2019-08-01 2019-08-01 紫色土土壤图像阴影增强方法

Publications (2)

Publication Number Publication Date
CN110428380A true CN110428380A (zh) 2019-11-08
CN110428380B CN110428380B (zh) 2022-02-18

Family

ID=68413647

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910705747.XA Active CN110428380B (zh) 2019-08-01 2019-08-01 紫色土土壤图像阴影增强方法

Country Status (1)

Country Link
CN (1) CN110428380B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111667432A (zh) * 2020-06-09 2020-09-15 中国电子科技集团公司第五十四研究所 一种基于物理模型的遥感影像阴影去除方法
CN111754501A (zh) * 2020-06-30 2020-10-09 重庆师范大学 基于fcm算法的自适应土壤图像阴影检测方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106023113A (zh) * 2016-05-27 2016-10-12 哈尔滨工业大学 基于非局域稀疏的卫星高分图像阴影区域恢复方法
CN109064411A (zh) * 2018-06-13 2018-12-21 长安大学 一种基于光照补偿的路面图像去阴影方法
CN109785329A (zh) * 2018-10-29 2019-05-21 重庆师范大学 基于改进slic算法的紫色土图像分割提取方法
WO2019104230A1 (en) * 2017-11-21 2019-05-31 The Charles Stark Draper Laboratory, Inc. System and method for imaging and localization of contrast-enhanced features in the presence of accumulating contrast agent in a body and amplification of features in ultrasound-generated images
US20190220959A1 (en) * 2015-12-16 2019-07-18 Dropbox, Inc. Enhancing a digital image

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20190220959A1 (en) * 2015-12-16 2019-07-18 Dropbox, Inc. Enhancing a digital image
CN106023113A (zh) * 2016-05-27 2016-10-12 哈尔滨工业大学 基于非局域稀疏的卫星高分图像阴影区域恢复方法
WO2019104230A1 (en) * 2017-11-21 2019-05-31 The Charles Stark Draper Laboratory, Inc. System and method for imaging and localization of contrast-enhanced features in the presence of accumulating contrast agent in a body and amplification of features in ultrasound-generated images
CN109064411A (zh) * 2018-06-13 2018-12-21 长安大学 一种基于光照补偿的路面图像去阴影方法
CN109785329A (zh) * 2018-10-29 2019-05-21 重庆师范大学 基于改进slic算法的紫色土图像分割提取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HUSSAIN K: "" A histogram specification technique for dark image enhancement using a local transformation method"", 《IPSJ TRANSACTIONS ON COMPUTER VISION & APPLICATIONS》 *
赵展: ""高分辨率遥感影像阴影提取与恢复研究"", 《测绘地理信息》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111667432A (zh) * 2020-06-09 2020-09-15 中国电子科技集团公司第五十四研究所 一种基于物理模型的遥感影像阴影去除方法
CN111667432B (zh) * 2020-06-09 2022-08-02 中国电子科技集团公司第五十四研究所 一种基于物理模型的遥感影像阴影去除方法
CN111754501A (zh) * 2020-06-30 2020-10-09 重庆师范大学 基于fcm算法的自适应土壤图像阴影检测方法

Also Published As

Publication number Publication date
CN110428380B (zh) 2022-02-18

Similar Documents

Publication Publication Date Title
CN108573276B (zh) 一种基于高分辨率遥感影像的变化检测方法
CN103971126B (zh) 一种交通标志识别方法和装置
CN103035013B (zh) 一种基于多特征融合的精确运动阴影检测方法
CN103186904B (zh) 图片轮廓提取方法及装置
CN108319973A (zh) 一种树上柑橘果实检测方法
CN103020970B (zh) 玉米果穗图像籽粒分割方法
CN107273896A (zh) 一种基于图像识别的车牌检测识别方法
CN107578035A (zh) 基于超像素‑多色彩空间的人体轮廓提取方法
CN108027248A (zh) 具有基于特征的定位和导航的工业车辆
CN110276386A (zh) 一种基于机器视觉的苹果分级方法及系统
CN104598908A (zh) 一种农作物叶部病害识别方法
CN109657632A (zh) 一种车道线检测识别方法
CN104102928B (zh) 一种基于纹理基元的遥感图像分类方法
CN109636784A (zh) 基于最大邻域和超像素分割的图像显著性目标检测方法
CN105095892A (zh) 基于图像处理的学生文档管理系统
CN108537239A (zh) 一种图像显著性目标检测的方法
CN108765336A (zh) 基于暗亮原色先验与自适应参数优化的图像去雾方法
CN107423735B (zh) 一种利用水平梯度和饱和度的车牌定位方法
CN106971397B (zh) 基于改进jseg算法的城市高分辨率遥感影像分割方法
CN109448015A (zh) 基于显著图融合的图像协同分割方法
CN110472598A (zh) 基于区域特征的svm机采棉花含杂图像分割方法及系统
CN104217440B (zh) 一种从遥感图像中提取建成区的方法
CN102208099A (zh) 一种抗光照变化的作物彩色图像分割方法
CN110687121B (zh) 一种陶瓦智能在线检测与自动分级方法和系统
CN109754440A (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