轮廓提取试验

我就是想要类似：https://www.pyimagesearch.com/2015/11/02/watershed-opencv/  这个的效果，将重叠的轮廓分割开成单独的轮廓。

复杂粘连轮廓的处理，要将粘连的轮廓分开，公司的大神用的是轮廓的缺陷点检测，但他的这个算法我看了下，有的凹点会检测不到，所以有的粘连的地方仍然分不开：

像这里有6块石头粘连在一起，有的凹陷点就检测不出来。。。

我看了下谷歌上有人用的：

这个我想试试。

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage import morphology,feature

#创建两个带有重叠圆的图像
x, y = np.indices((80, 80))
x1, y1, x2, y2 = 28, 28, 44, 52
r1, r2 = 16, 20
mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
image = np.logical_or(mask_circle1, mask_circle2)

#现在我们用分水岭算法分离两个圆
distance = ndi.distance_transform_edt(image) #距离变换
local_maxi =feature.peak_local_max(distance, indices=False, footprint=np.ones((3, 3)),
                            labels=image)   #寻找峰值

markers = ndi.label(local_maxi)[0] #初始标记点
labels =morphology.watershed(-distance, markers, mask=image) #基于距离变换的分水岭算法

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 8))
axes = axes.ravel()
ax0, ax1, ax2, ax3 = axes

ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
ax0.set_title("Original")
#ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
#ax1.set_title("Distance")
ax1.imshow(distance, cmap=plt.cm.gray, interpolation='nearest')
ax1.set_title("Distance")
ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
ax2.set_title("Markers")
ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
ax3.set_title("Segmented")

for ax in axes:
    ax.axis('off')

fig.tight_layout()
plt.show()

但我看了下：distance_transform_edt()这个函数的意思：

 Notes
    -----
    The euclidean distance transform gives values of the euclidean
    distance::

                    n
      y_i = sqrt(sum (x[i]-b[i])**2)
                    i

    where b[i] is the background point (value 0) with the smallest
    Euclidean distance to input points x[i], and n is the
    number of dimensions.

    Examples
    --------
    >>> from scipy import ndimage
    >>> a = np.array(([0,1,1,1,1],
    ...               [0,0,1,1,1],
    ...               [0,1,1,1,1],
    ...               [0,1,1,1,0],
    ...               [0,1,1,0,0]))
    >>> ndimage.distance_transform_edt(a)
    array([[ 0.    ,  1.    ,  1.4142,  2.2361,  3.    ],
           [ 0.    ,  0.    ,  1.    ,  2.    ,  2.    ],
           [ 0.    ,  1.    ,  1.4142,  1.4142,  1.    ],
           [ 0.    ,  1.    ,  1.4142,  1.    ,  0.    ],
           [ 0.    ,  1.    ,  1.    ,  0.    ,  0.    ]])

直接可以用OpenCV中的distanceTransform()代替，效果一致。

这个粘连的，下图就是距离变换与python的第二幅图一致。
但寻找局部峰值的那个函数：feature.peak_local_max()就有点难办了。。。我试过用minMaxLoc找到全局最大max，然后按0.9*max找到几个极大值！但这是不对的，理论上就不对，高中三角函数学过极大值可能比最大值的0.9倍甚至0.5倍都小，它们之间没有什么关系的。。。所以不行。

我用了另一种方法：

void localMaxima(Mat src,const int maxZeroSize,const double minPeakValue,vector<Point> &maxLocpts)
{
    double minVal=0,maxVal=0,diff=0;
    Point minLocPt,maxLocPt;

    minMaxLoc(src,&minVal,&maxVal,&minLocPt,&maxLocPt);//计算最大值

    while(maxVal>minPeakValue)
	{
    	maxLocpts.push_back(maxLocPt);
		Mat rect_part = src(Rect(maxLocPt.x-maxZeroSize,maxLocPt.y-maxZeroSize,2*maxZeroSize,2*maxZeroSize));
		rect_part.setTo(Scalar(0));
		//imwrite("temp_peak.jpg",src);
		minMaxLoc(src,&minVal,&maxVal,&minLocPt,&maxLocPt);//计算极大值
	}
}

然后：

int main()
{
	char src_folder[100];
	char dst_folder[100];
	const size_t img_rows=1456;
	const size_t img_cols=1936;
	int threshold_value=150;
	int element_size=3;//morph element size
	int blurwidth=3; //blur window size
	const int area_min=90; //area<area_min,delete it.
	int area_limitable=5800; //area<area_limitable,need judge
	const int area_max=15000;
	int rect_radius_contour=50; //around the contour :neighbors
	Mat element=getStructuringElement(MORPH_RECT,Size(2*element_size+1,2*element_size+1),
					Point(element_size,element_size));

	Mat morph_img(img_rows,img_cols,CV_8UC1);

	morph_img=imread("/home/jumper/Ore_projects/ore_granule/tempBinary.jpg",0);
	threshold(morph_img,morph_img,1,128,CV_THRESH_BINARY);
	Mat distance_img(img_rows,img_cols,CV_32FC1);
	distanceTransform(morph_img,distance_img,CV_DIST_L2,3);
	//imwrite("distance.jpg",distance_img);


	vector<Point> peakpts;
	localMaxima(distance_img,70,80,peakpts);
//	Mat peakimgs(img_rows,img_cols,CV_8UC1,Scalar(0));
//	cout<<"peak num: "<<peakpts.size()<<endl;
//	for(int index=0;index<=peakpts.size();index++)
//	{
//		circle(peakimgs,peakpts[index],5,CV_RGB(255,255,255),2);
//	}
//	imwrite("peakimg2.jpg",peakimgs);


	//marker...
	Mat marker(img_rows,img_cols,CV_32SC1,Scalar(0));
	for(int index=0;index!=peakpts.size();index++)
	{
		int row=peakpts[index].y;
		int col=peakpts[index].x;
		Mat rect_part = marker(Rect(col-70,row-70,2*70,2*70));
		rect_part.setTo(short(index+1));
		//marker.ptr<short>(row)[col]=index+1;
	}
//	imwrite("marker.jpg",marker);

	Mat src_img(img_rows,img_cols,CV_8UC3);
	src_img=imread("/home/jumper/Ore_projects/ore_granule/22db_1.bmp");
	Mat water_src(img_rows,img_cols,CV_8UC3);
	src_img.copyTo(water_src,morph_img);
//	imwrite("water_src.jpg",water_src);
	watershed(water_src,marker);

	Mat contrast_img(img_rows,img_cols,CV_8UC1,Scalar(0));
	for(int r=0;r!=marker.rows;r++)
	{
		for(int c=0;c!=marker.cols;c++)
		{
			if(marker.ptr<short>(r)[c]==-1)
			{
				circle(contrast_img,Point(c,r),5,CV_RGB(255,255,255),2);
			}
		}
	}
	imwrite("contrast_img.jpg",contrast_img);
	//marker.convertTo(marker,CV_8UC1);
	//imwrite("watershed.jpg",marker);

	return 0;
}

我的确得到了那两个极大值点：

我怕就2个点作为初始标记点太寒酸，我还特意将这两个点附近一块地方标记为1和2：显示画出来就是这样！

然后就可以运行分水岭分割了：这是原图部分:

but I got nothing!!!

没有正确分割 没有得到我想要的！！！？？？？

后来我又用了传说中的grabCut()算法，通过给它赋予可能的前景以及确定的前景和确定的背景，让它找到真正的前景，，，然而。。。不过可能是我操作的问题咯。

最后我用了另一种方法，既然我找到了盆地的聚水点（不知道是不是这个名词），那我找这个聚水点距离整个轮廓最近的两条边作为它的rect：

这还是不完善的版本，因为图像的长宽限制没加进去。。。在几处地方都还没做限制，因为我还只是试试效果而已。。。

#include<opencv2/opencv.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<iostream>
#include<fstream>
#include <vector>
using namespace cv;
using namespace std;

void localMaxima(Mat src,const int maxZeroSize,const double minPeakValue,vector<Point> &maxLocpts);

int main()
{
	char src_folder[100];
	char dst_folder[100];
	const size_t img_rows=1456;
	const size_t img_cols=1936;
	int blurwidth=3;

	Mat background_img(img_rows,img_cols,CV_8UC3);
	background_img=imread("/home/jumper/Ore_projects/ore_granule/imgs/1204-backgrond/back.bmp");



	Mat src_img(img_rows,img_cols,CV_8UC3);
	for(int ind=1;ind!=6;ind++)
	{
		cout<<"proceed "<<ind<<endl;
		sprintf(src_folder,"/home/jumper/Ore_projects/ore_granule/imgs/1204/22db_%d.bmp",ind);
		src_img=imread(src_folder);
		src_img=background_img-src_img;
		cvtColor(src_img,src_img,CV_BGR2HSV);
		vector<Mat> hsvmats(3);
		split(src_img,hsvmats);
		Mat vimg=hsvmats[2];
		threshold(vimg,vimg,110,255,CV_THRESH_BINARY);
		medianBlur(vimg,vimg,blurwidth);
//		sprintf(dst_folder,"/home/jumper/Ore_projects/ore_granule/1204result/%d.jpg",ind);
//		imwrite(dst_folder,vimg);

		Mat contours_img(img_rows,img_cols,CV_8UC1,Scalar(0));
		vector<vector<Point> > contours;
		findContours( vimg, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE );
		cv::Rect rectangle;
		Mat contours_img2(img_rows,img_cols,CV_8UC1,Scalar(0));
		for (int i = 0; i< contours.size(); i++)
		{
			contours_img2.setTo(Scalar(0));
			vector<Point> pts;
			vector<Point> result;
			double area=cv::contourArea(contours[i]);
			if(area<490)
			{
				continue;
			}
			cv::drawContours(contours_img,contours,i,Scalar(255),CV_FILLED);
			cv::drawContours(contours_img2,contours,i,Scalar(255),CV_FILLED);


			rectangle=cv::boundingRect(contours[i]);
			if(area<44000)
			{
				cv::rectangle(contours_img,rectangle,Scalar(255),2);
				continue;
			}


			Mat distance_img(rectangle.height,rectangle.width,CV_32FC1);
			distanceTransform(contours_img2,distance_img,CV_DIST_L2,3);

			vector<Point> peakpts;
			localMaxima(distance_img,70,80,peakpts);
			if(peakpts.size()==0)
			{
				continue;
			}
			Mat peakimgs(img_rows,img_cols,CV_8UC1,Scalar(0));
			cout<<"peak num: "<<peakpts.size()<<endl;
			for(int index=0;index<=peakpts.size();index++)
			{
				circle(peakimgs,peakpts[index],5,Scalar(255),2);
			}

			int col_left=rectangle.x;
			int col_right=rectangle.x+rectangle.width;
			int row_up=rectangle.y;
			int row_bottom=rectangle.y+rectangle.height;
			for(int index=0;index!=peakpts.size();index++)
			{
				int row=peakpts[index].y;
				int col=peakpts[index].x;

				int col_left_distance=col-col_left;
				int col_right_distance=col_right-col;
				int row_up_distance=row-row_up;
				int row_bottom_distance=row_bottom-row;
				int rect_height=min(row_up_distance,row_bottom_distance);
				int rect_width=min(col_left_distance,col_right_distance);
				cv::Point start_pt(col-rect_width,row-rect_height);
				cv::Point end_pt(col+rect_width,row+rect_height);
				cv::rectangle(contours_img,start_pt,end_pt,Scalar(255),2);
				circle(contours_img,peakpts[index],5,Scalar(0),2);
			}
		}
		sprintf(dst_folder,"/home/jumper/Ore_projects/ore_granule/1204result/%d_rect.jpg",ind);
		imwrite(dst_folder,contours_img);

	}
	return 0;
}

void localMaxima(Mat src,const int maxZeroSize,const double minPeakValue,vector<Point> &maxLocpts)
{
    double minVal=0,maxVal=0,diff=0;
    Point minLocPt,maxLocPt;

    minMaxLoc(src,&minVal,&maxVal,&minLocPt,&maxLocPt);//计算最大值

    while(maxVal>minPeakValue)
	{
    	maxLocpts.push_back(maxLocPt);
		Mat rect_part = src(Rect(maxLocPt.x-maxZeroSize,maxLocPt.y-maxZeroSize,2*maxZeroSize,2*maxZeroSize));
		rect_part.setTo(Scalar(0));
		//imwrite("temp_peak.jpg",src);
		minMaxLoc(src,&minVal,&maxVal,&minLocPt,&maxLocPt);//计算极大值
	}
}

黑色的圈是peak点：

这是结果。

但可能阈值等还要调整，效果还不够好：

这样子不够准确而且有的会有丢失。。。

我再想想。。。

路漫漫。。。

今天调整了下值：

not enough...

不稳妥。。。

修改后：

但如果有大量粘连的情况还是不够好，当大量粘连时我的质心点会有点多，也就是这个距离变换还不够。。。还要不断改善。

今天是2019.2.18，看到了一种改善的分水岭效果：

同时发现分水岭来切割粘连图像比较适合于目标差不多大小的情况，如果大小相差悬殊，则不适合，分割不开！

看这幅图虽然右边不是最后的结果，但也是结果的原型，可以看到有的分割开的都是因为目标相差不大，而分割不开的基本都是

要么粘连特别紧密，要么粘连没有那么紧密但是大小相差较大

调节抑制局部极小值后：还是有分不开的很多。

bw0 = imread('G:\src.png');
%bw=imdilate(bw0,ones(3));
bw=im2bw(bw0);
imshow(bw,[]);
D = -bwdist(~bw);
figure;
imshow(D,[]);
mask = imextendedmin(D,2);
figure;
imshow(mask,[]);
imshowpair(bw,mask,'blend');
D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
bw3 = bw;
bw3(Ld2 == 0) = 0;
figure;
imshow(bw3);

等等，我发现一个地方我错了，对于函数imextendedmin(D,4).我以为是对距离变换后的图像D，对每一个连通域，取其数值整数大于前4个的置1.结果发现不是！！

imextendedmin(D,4).的结果是这样，看我用笔画出来的这个连通域，这个连通域的数值是：

看数值，这个连通域经过这个函数后应该大于4、5、6、7的值位置都置1，也就是这个大的连通域应该被分为2个小的置1的块。但经过这个函数后的结果是：

结果只有1个小块，所以也就是我用笔画出来的那样，所以后面用分水岭也是分不开这个粘连体的。

所以我要改造一下这个函数，让它按我以为的那样实现。得到类似各个石头的核心区域mask。

1：然后我看了一下函数D2=imimposemin(D,mask);的源码，这个函数的意思是：强制最小。D是距离变换结果，令fm=mask；且fm核心区域即为1的地方全部置为-Inf即负无穷大；为0的地方置为正无穷大Inf；

2：距离变换图像D的最大值max-最小值min的绝对值为range，如果range不等于0，则h=range*0.001；否则为0.1；

3：修改距离变换的结果fp1=D+h 即给距离变换的值都加上常数h。注意Matlab的距离变换结果是负数D的值都是负的。如果转opencv要注意；

4：g=min（fp1，fm）；所以此时图像g中要么是负无穷大代表的每个目标的核心区域，要么就是别的数值代表的背景以及目标的边缘包含粘连边缘；

5：对fm取反得到fm2，对g取反得到g2；g2=1-g；

6：对fm2，g2，以8连通域做形态学重建J=imreconstruct(fm2，g2,8)；这个函数我追源码追不下去了，只看到J根g2其实差不多，但也有差别，这个差别就是形态学重建导致的，而形态学重建究竟是怎么实现的 我不知道！？

7：Ld2=1-J； 对Ld2进行分水岭 。Ld3=watershed(Ld2)；即Ld3==0的位置就是各个粘连目标的分界线；其它都是正整数表示不同的区域。

8：对于最原始的粘连图像bw(Ld3==0)=0；这样就分开了粘连目标。

但opencv的实现我还在写。

我发现另一个更好的算法，而且可并行。中间是原版效果，最后一个是现在的效果：

 

 参考代码放在：overlap-segmentation.zip-其它文档类资源-CSDN下载