轮廓提取试验

来源：互联网发布：铁血工造知乎编辑：程序博客网时间：2024/05/24 06:27

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

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

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

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

这个我想试试。

import numpy as npimport matplotlib.pyplot as pltfrom scipy import ndimage as ndifrom skimage import morphology,feature#创建两个带有重叠圆的图像x, y = np.indices((80, 80))x1, y1, x2, y2 = 28, 28, 44, 52r1, r2 = 16, 20mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2image = 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 = axesax0.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 sizeint blurwidth=3; //blur window sizeconst int area_min=90; //area<area_min,delete it.int area_limitable=5800; //area<area_limitable,need judgeconst int area_max=15000;int rect_radius_contour=50; //around the contour :neighborsMat 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点：

这是结果。