【OpenCV】第十七章: 图像分割与提取_bgdmodel = np.zeros

第十七章: 图像分割与提取

我们在图像处理中，经常会需要从图像中将前景对象作为目标图像分割或者提取出来，比如监控视频中的车辆、行人等提取出来。
而实现图像分割可以用：形态学变换、阈值算法、图像金字塔、图像轮廓、边缘检测等方法实现。但是本章介绍使用分水岭算法及GrabCut算法对图像进行分割和提取

一、分水岭算法

极好的参考资料：    图像分割的经典算法：分水岭算法 - 知乎

• 算法原理
  分水岭算法的启发思路是：把一幅灰度图像看成地理上的地形表面，每个像素的灰度值代表高度。灰度值大的区域看成山丘，灰度值小的区域看成凹地。假如开始下雨，凹地首先被雨水填上，如果雨水一直下直到下到地平面(假设地平面的灰度值是100，小于100的都是凹地，大于100的都是山丘)，此时灰度值小于100的都变成黑色了，大于100的像素组成的图案就是一幅灰度图的分水岭线，其实也就是用阈值找到图像的轮廓。找到轮廓后，假设雨继续下，此时我们要在轮廓和轮廓之间筑坝防止水互相注入，然后雨继续下，每个轮廓又不断注水，被水淹的地方就变成黑色，然后每个轮廓区域就又形成自己的轮廓，其实就是找到每个轮廓的轮廓，就实现了图像的分割。

  或者有的教材说：灰度图像中灰度值较大的像素连成的线可以看做山脊，也就是分水岭。其中的水就是用于二值化的gray threshold level，二值化阈值可以理解为水平面，比水平面低的区域会被淹没，刚开始用水填充每个孤立的山谷(局部最小值)。当水平面上升到一定高度时，水就会溢出当前山谷，可以通过在分水岭上修大坝，从而避免两个山谷的水汇集，这样图像就被分成2个像素集，一个是被水淹没的山谷像素集，一个是分水岭线像素集。最终这些大坝形成的线就对整个图像进行了分区，实现对图像的分割。
  

   但是，用上面的分水岭算法进行图像分割时，由于噪声点或其它因素的干扰，可能会得到密密麻麻的小区域，即图像被分得太细（over-segmented，过度分割），这因为图像中有非常多的局部极小值点，每个点都会自成一个小区域。
  解决方法有：
  1、对图像进行高斯平滑操作，抹除很多小的最小值，这些小分区就会合并。
  2、不从最小值开始增长，可以将相对较高的灰度值像素作为起始点（需要用户手动标记），从标记处开始进行淹没，则很多小区域都会被合并为一个区域，这被称为基于图像标记(mark)的分水岭算法。 下面三个图分别是原图，分水岭过分割的图以及基于标记的分水岭算法得到的图：
  
  其中标记的每个点就相当于分水岭中的注水点，从这些点开始注水使得水平面上升，但是如上图所示，图像中需要分割的区域太多了，手动标记太麻烦，我们可是使用距离转换的方法进行标记，OpenCV中就是使用的这种方法。

• 距离变换函数：cv2.distanceTransform()
  该函数是计算二值图像内任意点到最近背景点的距离，也就是计算图像内非零值点像素点到其最近的零值点像素点的距离。就是计算二值图像中所有像素点距离其最近的零值像素点的距离。如果某个像素点本身就是零值像素点，那么这个像素点计算出来的距离就是0。

  所以，这个计算结果反映的是各个像素点与图像背景(图像背景就是像素值为0的区域)之间的距离。如果一个像素点是前景对象的质心或者中心，这个像素点就距离零值点最远，那么计算出来的结果就最大；如果一个像素点是前景对象的边缘点，那这个像素就离0值点较近，其计算结果就较小。

  所以，对上述的距离结果进行阈值化处理后，就可以得到图像前景对象的中心、骨架等信息，还能细化前景对象的轮廓，就可以准确的获得前景图像。

  API：dst = cv2.distanceTransform(img, distanceType, maskSize)
  　　img：是8为单通道二值图像
  　　distanceType：距离计算方法。
  　　　　　　　cv2.DIST_USER：用户自定义距离
  　　　　　　　cv2.DIST_L1：街道距离，distance=|x1-x2|+|y1-y2|
  　　　　　　　cv2.DIST_L2：欧式距离
  　　　　　　　cv2.DIST_C：distance=max(|x1-x2|,|y1-y2|)
  　　　　　　　cv2.DIST_L12，cv2.DIST_FAIR，cv2.DIST_WELSCH，cv2.DIST_HUBER等7中距离计算方式。
  　　dst：函数返回值是一个图像类型CV_32F的浮点数，尺寸和img相同。
  
  说明：这个函数是用来确定前景对象的。我们之前确定前景对象，我们是用边缘检测以及基于边缘检测的各种轮廓来提取前景对象，也或者用形态学变化来提取的。但是如果前景对象是紧挨着的或者是重叠遮挡的着的，我们还想把逐个逐个的前景对象的轮廓都分别提取出来，前面的方法就效果大打折扣了，我们就需要本函数的方法来提取。 本函数用的方法效果非常好！因为如果有些前景对象本身就是粘连在一起的，用轮廓提取就根本做不到逐个对象一个个提取出来，而虽然开闭运算可以提取出来，但是我们画轮廓的时候也希望单个轮廓和单个轮廓之间紧挨而行，此时开闭运算就做不到了，就得要用本函数来实现！

• 标注函数：cv2.connectedComponents()
  通过距离变换函数，计算出图像中各个区域是前景还是背景还是未知区域后，我们就开始对这些区域进行标注了。
  API：retval, labels = cv2.connectedComponents(img)
  　　　img：8位单通道待标注图像
  　　　retval：返回的是标注数量
  　　　labels：是标注结果图像
  函数cv2.connectedComponents()在标注图像时，会将背景标注为0，将其他对象从1开始整数往后标，表示不同的前景区域。
  我们可以将标注区域理解为进行分水岭算法分割的"种子"区域。
• 分水岭算法函数：markers = cv2.watershed(image, labels)
  image：输入图像，必须是8位三通道的图像。
  labels：是32位单通道的标注结果，尺寸和image相同，也就是上面的标注函数的返回值。
  但是要说明的是，分水岭函数中的参数labels中，0代表未知区域。所以我们还要把上面标注函数的返回值labels进行调整：
  把cv2.connectedComponents()的返回值labels = labels+1 ，还得 laels[未知区域]=0 进行调整后才能传入cv2.watershed()里面。
  markers：函数返回值中，每个像素要么被置为初期的"种子值"，要么被设置为-1，表示边界。

  #例17.1 分水岭算法图像分割实例  
  import cv2
  import numpy as np
  import matplotlib.pyplot as plt
   
  img = cv2.imread(r'C:\Users\25584\Desktop\water_coins.jpg') #img.shape返回：(312, 252, 3)
  img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  #img_gray.shape返回：(312, 252)
   
  #------------Otsu阈值处理,变成二值图像--------------
  t, otsu = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)   #t返回162.0， otsu.shape返回(312, 252)
   
  #------------形态学的开运算，就是先腐蚀erode后膨胀dilate,目的一是去噪，二是先把前景对象重叠的部分分开，方便后面计数或者画每个对象的轮廓-------------
  img_opening = cv2.morphologyEx(otsu, cv2.MORPH_OPEN, kernel=np.ones((3,3), np.uint8), iterations=2)  #A 这还是一个二值图像
   
  #-------计算距离,确定前景对象--------------------
  dist = cv2.distanceTransform(img_opening, cv2.DIST_L2, 5)   #float32的浮点型数组，dist.shape返回(312, 252),dist是一个灰度图像
  th, sure_fg = cv2.threshold(dist, 0.5*dist.max(), 255, cv2.THRESH_BINARY) #把dist阈值化处理，变成一个0和255的二值图像，此时就是我们要的确定前景
  sure_fg = np.uint8(sure_fg)
   
  #-----计算确定背景、计算未知区域------------------
  sure_bg = cv2.dilate(img_opening, kernel=np.ones((3,3), np.uint8), iterations=3) #对前景膨胀来确定背景
  unknown = cv2.subtract(sure_bg, sure_fg)  #确定背景图-确定前景图，生成未知区域图
   
  #------标注确定前景图,调整标注规则---------------------------
  ret, labels = cv2.connectedComponents(sure_fg)  #有24个硬币，ret返回是25, labels是一个形状是(312, 252)的int32的数组
  labels = labels+1  #把背景标为1，前景对象依次为2，3，，，26
  labels[unknown==255]=0  #0代表未知区域
   
  #------------使用分水岭算法对图像进行分割---------------
  img1 = img.copy()
  markers = cv2.watershed(img1,labels)
  img1[markers==-1]=[0,255,0]
   
  #可视化：
  plt.figure(figsize=(12,6))
  plt.subplot(251), plt.imshow(img[:,:,::-1])   #原图 
  plt.subplot(252), plt.imshow(img_gray, cmap='gray') #灰度图
  plt.subplot(253), plt.imshow(otsu, cmap='gray')  #otsu阈值处理后的二值图
  plt.subplot(254), plt.imshow(img_opening, cmap='gray')  #开运算去噪后的图像
  plt.subplot(255), plt.imshow(dist, cmap='gray')   #距离图像
  plt.subplot(256), plt.imshow(sure_fg, cmap='gray')  #确定前景
  plt.subplot(257), plt.imshow(sure_bg, cmap='gray')  #确定背景
  plt.subplot(258), plt.imshow(unknown, cmap='gray')  #确定未知区域图
  plt.subplot(259), plt.imshow(labels, cmap='gray')  #标注图
  plt.subplot(2,5,10), plt.imshow(img1[:,:,::-1])  #分割结果
  plt.show()
  

  

  A：我们这里主要是去噪目的的，因为这张图的前景对象粘连在一起，用开运算把各个对象分开以后画出来的轮廓就也是分开的，我们希望画出来的轮廓该埃还得挨，因为硬币的轮廓本来就是挨着的，所以我们要用其他方法寻找硬币的轮廓！用cv2.distanceTransform()来计算硬币像素点到其最近的零值像素点的距离，根据这个距离结果我们去判断硬币的真正轮廓，也就是寻找图像的确定前景和确定背景。

  二、交互式前景提取：GrabCut算法

  交互式前景提取方法是：先用矩形框将要提取的前景对象框出来，也就是先指出前景区域的大致位置范围，然后算法不断迭代分割直到达到最好的效果。但是有时这种做法效果不理想，就需要用户干预提取，怎么干预？制作提取掩膜，也就是模板，就是制作一幅和原始图像大小一样的任意一幅图像，这副图像里面的白色标注代表要提取的前景区域，黑色标注表示的是背景区域。然后将标注图像作为掩膜，让算法继续迭代提取前景。

• GrabCut算法原理：
  1、将前景所在的大致位置使用矩形框标注出来。由于矩形框框起来的有前景还部分背景，所以矩形框区域实际上是未确定区域，而矩形框以外是确定背景区域。
  2、根据矩形框外的确定背景的数据来分析框内区域中的前景和背景。
  3、用高斯混合模型(Gaussians Mixture Model, GMM)对前景和背景建模。GMM模型根据确定背景的像素点数据进行建模，创建确定背景的像素分布。然后判断对框内像素与已知分类像素(前景和背景)之间的关系进行分类。
  4、根据像素分布情况生成一幅图，图中的节点就是各个像素点。除了像素点外，还有两个节点：前景节点和背景节点。所有的前景像素都和前景节点相连，所有的背景像素都和背景节点相连。每个像素连接到前景节点或背景节点的边的权重由像素是前景或背景的概率来决定。
  5、图中的每个像素除了与前景点或背景节点相连外，彼此之间还存在着连接。两个像素连接的边的权重值由它们的相似性决定，两个像素的颜色越接近，边的权重值越大。
  6、完成节点连接后，需要解决的问题变成了一幅连通的图。该图上根据各自边的权重关系进行切割，将不同的点划分为前景节点和背景节点。
  7、不断重复上述过程，直至分类收敛为止。

• API：mask, bgdModel, fgdModel = cv2.grabCut(img, mask, rect, bgdModel, fgdModel, iterCount[, mode])
  img：输入图像，要求是8位3通道的图像。
  mask：掩膜图像，要求是8位单通道图像，该参数用于确定前景区域、背景区域和不确定区域，可以设置为4种形式：
  　　　cv2.GC_BGD：表示确定背景，也可以用0表示
  　　　cv2.GC_FGD：表示确定前景，也可以用1表示
  　　　cv2.GC_PR_BGD：表示可能的背景，也可以用2表示
  　　　cv2.GC_PR_FGD：表示可能的前景，也可以用3表示
  rect：指包含前景对象的区域，格式是(x,y,w,h)，该区域外都被默认为确定背景
  bgdModel：算法内部使用的数组，只需要创建大小为(1,65)的numpy.float64数组。
  fgdModel：算法内部使用的数组，只需要创建大小为(1,65)的numpy.float64数组。
  iterCount：表示迭代次数
  mode：表示迭代模式，有4种模式：
  　　　cv2.GC_INIT_WITH_RECT：使用矩形模板
  　　　cv2.GC_INIT_WITH_MASK：使用自定义模板，所有roi区域外的像素都会自动被处理为背景
  　　　cv2.GC_EVAL：修复模式
  　　　cv2.GC_EVAL_FREEZE_MODEL：使用固定模式

• 交互式前景提取的具体操作步骤：
  1、加载要提取前景的图片img，可以是3通道的彩图。
  2、制作掩膜图像(有的地方也叫模板)mask。掩膜图像mask的尺寸应该和img大小相同，但是必须是二通道。可以先将mask的值全部置为0，就是先假设都是确定背景，在后续再根据需要修改部分值。
  3、指定前景对象的大致位置rect, rect的格式是(x,y,w,h),xy表示方框的左上角像素点的坐标值，w表示框的宽，h表示框的高。
  4、创建bgdModel, fgdModel对象。
  5、用cv2.grabCut(img, mask, rect, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_RECT)算法迭代，迭代完毕生成一个新的mask和bgdModel、fgdModel。其中bgdModel和fgdModel我们暂时用不到，我们现在需要的是新生成的mask，这个mask的尺寸大小还是和img一样，也是uint8类型，但是里面的像素值只有0、1、2、3，四种数值。其中0表示确定背景、1表示确定前景、2表示可能背景、3表示可能前景。所以，我们要把mask里面的==0和==2的像素点都置为0，把==1和==3的像素点都置为1。这样就是我们的前景对象提取的掩膜图像了。
  6、原图img和掩膜图像mask进行按位与运算，提取出前景。

  #例17.2 仅使用矩形模板，让grabcut算法迭代，提取前景  
  import cv2
  import numpy as np
  import matplotlib.pyplot as plt
   
  img = cv2.imread(r'C:\Users\25584\Desktop\lenacolor.png')   #img.shape返回：(512, 512, 3)
   
  mask = np.zeros(img.shape[:2], np.uint8)   #掩膜值都设为0
  bgdModel = np.zeros((1,65), np.float64)
  fgdModel = np.zeros((1,65), np.float64)
  rect = (50, 50, 420, 500)   #矩形框
   
  img0 = img.copy()   #把矩形框画出来看看
  cv2.rectangle(img0, (50, 50), (420, 500), (0,255,0), 3)  
   
  mask1, bgd1, fgd1 = cv2.grabCut(img, mask, rect, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_RECT)  #mask11.shape是(512, 512)，np.unique(mask11.ravel())返回array([0, 2, 3])
  mask11 = np.where((mask1==0)|(mask1==2), 0,1).astype('uint8')
  img1 = img.copy()
  ogc1 = img1*mask11[:,:, np.newaxis]
   
  #可视化：
  plt.figure(figsize=(16,6))
  plt.subplot(151), plt.imshow(img[:,:,::-1])   #原图 
  plt.subplot(152), plt.imshow(img0[:,:,::-1]) 
  plt.subplot(153), plt.imshow(mask1, cmap='gray')  
  plt.subplot(154), plt.imshow(mask11, cmap='gray')  
  plt.subplot(155), plt.imshow(ogc1[:,:,::-1])  
  plt.show()
  

  
   

  #例17.3 使用标注图像作为模板，迭代提取前景 
  import cv2
  import numpy as np
  import matplotlib.pyplot as plt
   
  img = cv2.imread(r'C:\Users\25584\Desktop\lenacolor.png')   #img.shape返回：(512, 512, 3)
   
  mask = np.zeros(img.shape[:2], np.uint8)   #掩膜值都设为0,后面再迭代
  bgdModel = np.zeros((1,65), np.float64)
  fgdModel = np.zeros((1,65), np.float64)
  rect = (50, 50, 420, 500)   #矩形框
   
  img0 = img.copy()   #把矩形框画出来看看
  cv2.rectangle(img0, (50, 50), (420, 500), (0,255,0), 3)  
   
  mask1, bgd1, fgd1 = cv2.grabCut(img, mask, rect, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_RECT) 
   
  #-----------制作标注图像------------------------------
  mask2 = cv2.imread(r'C:\Users\25584\Desktop\mask.png')   #mask2_2.shape返回：(512, 512, 3)
  mask2_1 = cv2.cvtColor(mask2, cv2.COLOR_BGR2GRAY)
  mask3 = mask1.copy()
  mask3[mask2_1==0]=0
  mask3[mask2_1==255]=1
  #--------------根据标注模板继续迭代---------------------
  mask4, bgd4,fgd4 = cv2.grabCut(img, mask3, None, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_MASK)
  mask5 = np.where((mask4==0)|(mask4==2), 0,1).astype('uint8')  #合并可能区域
  img1 = img.copy()
  ogc = img1*mask5[:,:, np.newaxis]  #原图与模板与运算，提取前景对象
   
  #可视化：
  plt.figure(figsize=(22,8))
  plt.subplot(191), plt.imshow(img[:,:,::-1])   #原图 
  plt.subplot(192), plt.imshow(img0[:,:,::-1]) #矩形框图
  plt.subplot(193), plt.imshow(mask1, cmap='gray')  #mask第一个迭代后的图像mask1
  plt.subplot(194), plt.imshow(mask2[:,:,::-1])  #掩膜图像，原图，彩图
  plt.subplot(195), plt.imshow(mask2_1, cmap='gray')  #掩膜图像，灰度图
  plt.subplot(196), plt.imshow(mask3, cmap='gray')  #mask1和人工标注过的掩膜图像用切片提取标注信息的图片mask3
  plt.subplot(197), plt.imshow(mask4, cmap='gray')  #mask3作为掩膜用算法迭代的结果mask4
  plt.subplot(198), plt.imshow(mask5, cmap='gray') #合并可能区域
  plt.subplot(199), plt.imshow(ogc[:,:,::-1])  #提取图像
  plt.show()
  

  
   

  #例17.4 直接使用模板提取图像前景  
  import cv2
  import numpy as np
  import matplotlib.pyplot as plt
   
  img = cv2.imread(r'C:\Users\25584\Desktop\lenacolor.png') 
  img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
   
  mask = np.zeros(img.shape[:2], np.uint8)   
  bgdModel = np.zeros((1,65), np.float64)
  fgdModel = np.zeros((1,65), np.float64)
  mask[30:512, 50:400]=3  #lena头像的可能区域
  mask[70:300, 150:200]=1  #lena头像的确定区域，如果不设置这个区域，头像的提取不完整
  mask1=mask.copy()
  mask2,bgd,fgd = cv2.grabCut(img, mask1, None, bgdModel, fgdModel, 5, cv2.GC_INIT_WITH_MASK)
  mask3 = np.where((mask2==2)|(mask2==0), 0,1).astype('uint8')
  ogc = img*mask3[:,:,np.newaxis]
   
  #可视化：
  plt.figure(figsize=(12,4))
  plt.subplot(151), plt.imshow(img[:,:,::-1])   #原图 
  plt.subplot(152), plt.imshow(mask,cmap='gray') #可能区域和确定区域
  plt.subplot(153), plt.imshow(mask2,cmap='gray')  #算法迭代后的mask
  plt.subplot(154), plt.imshow(mask3,cmap='gray')  #合并可能区域后的mask
  plt.subplot(155), plt.imshow(ogc[:,:,::-1])  
  plt.show()