写在前面

再过2个小时，就是2020年了，回顾2019，有得也必有失，只感叹时间过得真的很快，像我这样的贫困人口也要被消灭了，哈哈。2019最后这2个小时就用一篇博客作为结尾吧。

昨天CSDN上一位同学说他们要做一个豆子识别的任务，但是豆子之间的粘连比较严重，如下图，不利于豆子的识别，问我有什么方法可以解决，当然深度学习完全可以做到完美分割与识别，但是他的需求只是预处理，传统的图像处理就可以，所以这里只讨论传统方法。

看二值图，确实有些地方粘连挺严重，如红色框展示的地方。如果用形态学腐蚀，开运算，可以解决占粘连，但是会严重影响豆子本身的尺寸信息，所以这个问题还是很有趣的，值得探索一下。

 灰度方向统计滤波

观察二值图，我最直接的感受就是粘连处相对豆子本身是比较细小的，假设二值图豆子及粘连处为1，那么就存在一个特点：粘连处总有一个方向为1的个数很少。所以我就分别统计二值图中一个滤波窗口内4个方向为1的个数，然后以 最小值/最大值 作为滤波结果，就叫它：灰度方向统计滤波吧。

这样的好处就是豆子本身最小值/最大值的比值相对很大，甚至为1，而粘连处则一定会很小；不足之处就是像二值图中红色框中粘连过大，甚至和小的豆子的尺寸快一样了，那这样的处理就不那么如意了。

我最后的结果是这样的：

可以看到绝大部分的粘连已经被消除，然后红色框中粘连过大的地方依然没有消除。所以还必须施加某种约束以优化效果。

 梯度方向统计滤波 

再次观察二值图，发现粘连很大的地方其内部的梯度是为0的，这就恰好就是我需要的结果，但是需要考虑的是：豆子内部的梯度也是为0的。所以如果直接计算梯度，结果只是提取了边缘，相当于Canny边缘检测的第二步，这不是我想要的结果。

   那如何保证只抑制粘连处，而不抑制豆子呢（尤其是小的豆子）？ 

     观察二值图，可以发现，粘连处的边缘远远少于豆子的边缘。所以又会有一个特征：在一定滤波尺寸内，粘连处满足一定梯度方向要求的像素一定也远远少于豆子（因为只有边缘存在梯度）。  

所以我只要统计滤波窗口内满足一定梯度方向要求的像素，然后将它们梯度幅值的均方或平均作为滤波结果即可，就叫它：梯度方向统计滤波吧。

上面的图是一篇论文的，我就借鉴一下吧，将滤波窗划分为4个象限，然后。。。，唉，实在不想打字画图了，直接简单地写在纸上吧：

请忽略字，太丑，赶紧写完回去睡觉了。

结果还是可以的：既抑制了粘连部位，又保住了豆子，如红色框处。

结果 

那就把两者的结果的优点结合一下吧，两个滤波结果相乘：完美解决！

代码 

LIG.h:

#ifndef LIG_H
#define LIG_H
/*************************************************************************************
/编译会出现重定义的报错，运行时请将属性管理器连接器中的强制文件输出改为【Enabled (/FORCE)】
*************************************************************************************/
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
 
class LIG{
 
public:
 
	cv::Mat Calc_Imap(cv::Mat src, int wsize);   //src必须为8位无符号整型图像
	cv::Mat Calc_Gmap(cv::Mat src, int wsize);   
	void Calc_LIG(cv::Mat src, cv::Mat &Imap, cv::Mat &Gmap, cv::Mat &IGmap,int wsize, double k);
 
private:
	inline cv::Mat get_filter0(int s);
	inline cv::Mat get_filter45(int s);
	inline cv::Mat get_filter90(int s);
	inline cv::Mat get_filter135(int s);
	void Nozerosfilter(cv::Mat src, cv::Mat &dst, cv::Mat &filter);
 
	inline double Mean_surround(cv::Mat &src);
	inline cv::Mat get_NW_data(cv::Mat &src, int ROI_center);
	inline cv::Mat get_NE_data(cv::Mat &src, int ROI_center);
	inline cv::Mat get_SW_data(cv::Mat &src, int ROI_center);
	inline cv::Mat get_SE_data(cv::Mat &src, int ROI_center);
	void Calc_gradient(cv::Mat img, cv::Mat_<double> &Gx, cv::Mat_<double> &Gy);
 
};
 
 
inline cv::Mat LIG::get_filter45(int s){
	cv::Mat dst = cv::Mat::zeros(s, s, CV_8U);
	for (int i = 0; i < s; ++i){
		dst.at<uchar>(i, s - i - 1) = 1;
	}
	return dst;
}
 
inline cv::Mat LIG::get_filter135(int s){
	cv::Mat dst = cv::Mat::zeros(s, s, CV_8U);
	for (int i = 0; i < s; ++i){
		dst.at<uchar>(i, i) = 1;
	}
	return dst;
}
 
inline cv::Mat LIG::get_filter0(int s){
	cv::Mat dst = cv::Mat::zeros(s, s, CV_8U);
	for (int i = 0; i < s; ++i){
		dst.at<uchar>((s - 1) / 2, i) = 1;
	}
	return dst;
}
 
inline cv::Mat LIG::get_filter90(int s){
	cv::Mat dst = cv::Mat::zeros(s, s, CV_8U);
	for (int i = 0; i < s; ++i){
		dst.at<uchar>(i, (s - 1) / 2) = 1;
	}
	return dst;
}
 
void LIG::Nozerosfilter(cv::Mat src, cv::Mat &dst, cv::Mat &filter){
	filter2D(src, dst, -1, filter, cv::Point(-1, -1), 0);
}
 
 
inline double LIG::Mean_surround(cv::Mat &src){
	double sum = 0.0, mean = 0.0;
	size_t mid = (src.total()) / 2;
	cv::MatIterator_<double> it = src.begin<double>();
	cv::MatIterator_<double> it_end = src.end<double>();
	cv::MatIterator_<double> it_mid = it + mid;
	for (; it != it_end; ++it){
		if (it != it_mid)
			sum = sum + *it;
	}
	mean = sum / (src.total() - 1);
	return mean;
}
 
inline cv::Mat LIG::get_NW_data(cv::Mat &src, int ROI_center){
	cv::Rect area(0, 0, ROI_center + 1, ROI_center + 1);
	cv::Mat	NW_data = src(area);
	return NW_data;
}
 
inline cv::Mat LIG::get_NE_data(cv::Mat &src, int ROI_center){
	cv::Rect area(ROI_center, 0, ROI_center + 1, ROI_center + 1);
	cv::Mat	NE_data = src(area);
	return NE_data;
}
 
inline cv::Mat LIG::get_SW_data(cv::Mat &src, int ROI_center){
	cv::Rect area(0, ROI_center, ROI_center + 1, ROI_center + 1);
	cv::Mat	SW_data = src(area);
	return SW_data;
}
 
inline cv::Mat LIG::get_SE_data(cv::Mat &src, int ROI_center){
	cv::Rect area(ROI_center, ROI_center, ROI_center + 1, ROI_center + 1);
	cv::Mat	SE_data = src(area);
	return SE_data;
}
 
void LIG::Calc_gradient(cv::Mat img, cv::Mat_<double> &Gx, cv::Mat_<double> &Gy){
	cv::Mat_<double> dImg;
	img.convertTo(dImg, CV_64F);
	int rows = img.rows;
	int cols = img.cols;
 
	cv::Mat_<double> xTopVec = dImg.col(1) - dImg.col(0);
	cv::Mat_<double> xBotVec = dImg.col(cols - 1) - dImg.col(cols - 2);
	cv::Mat_<double> xForwMat = dImg(cv::Range(0, rows), cv::Range(0, cols - 2));
	cv::Mat_<double> xBackMat = dImg(cv::Range(0, rows), cv::Range(2, cols));
	cv::Mat_<double> centGx = (xBackMat - xForwMat) / 2;
	cv::Mat_<double> tmpGx = cv::Mat::zeros(rows, cols, CV_64F);
	for (int i = 1; i < cols - 1; i++) {
		centGx.col(i - 1).copyTo(tmpGx.col(i));
	}
	xTopVec.copyTo(tmpGx.col(0));
	xBotVec.copyTo(tmpGx.col(cols - 1));
 
	cv::Mat_<double> yTopArr = dImg.row(1) - dImg.row(0);
	cv::Mat_<double> yBotArr = dImg.row(rows - 1) - dImg.row(rows - 2);
	cv::Mat_<double> yForwMat = dImg(cv::Range(0, rows - 2), cv::Range(0, cols));
	cv::Mat_<double> yBackMat = dImg(cv::Range(2, rows), cv::Range(0, cols));
	cv::Mat_<double> centGy = (yBackMat - yForwMat) / 2;
	cv::Mat_<double> tmpGy = cv::Mat::zeros(rows, cols, CV_64F);
	for (int i = 1; i < rows - 1; i++) {
		centGy.row(i - 1).copyTo(tmpGy.row(i));
	}
	yTopArr.copyTo(tmpGy.row(0));
	yBotArr.copyTo(tmpGy.row(rows - 1));
 
	Gx = tmpGx;
	Gy = tmpGy;
}
 
#endif

LIG.cpp: 

#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
#include "LIG.h"
 
cv::Mat LIG::Calc_Imap(cv::Mat src, int s){
 
	if (s % 2 == 0){
		fprintf(stderr, "Please enter odd size!");
		exit(-1);
	}
	if (src.channels() > 1) cvtColor(src, src, CV_RGB2GRAY);
	if (src.type() != CV_8U)
		src.convertTo(src, CV_8U, 1, 0);
 
	cv::Mat filter0, filter45, filter135, filter90,  //四个方向滤波器
		dst0 = cv::Mat::zeros(src.rows, src.cols, CV_8U),
		dst45 = cv::Mat::zeros(src.rows, src.cols, CV_8U),
		dst135 =cv::Mat::zeros(src.rows, src.cols, CV_8U),
		dst90 = cv::Mat::zeros(src.rows, src.cols, CV_8U);
 
	//获取四个方向滤波核
	filter0 = get_filter0(s);
	filter45 = get_filter45(s);
	filter90 = get_filter90(s);
	filter135 = get_filter135(s);
 
	//四个方向分别滤波
	Nozerosfilter(src, dst0, filter0);
	Nozerosfilter(src, dst45, filter45);
	Nozerosfilter(src, dst90, filter90);
	Nozerosfilter(src, dst135, filter135);
 
	std::vector<cv::Mat> filtedst;
	filtedst.push_back(dst0);
	filtedst.push_back(dst45);
	filtedst.push_back(dst90);
	filtedst.push_back(dst135);
	cv::Mat min = cv::Mat::zeros(src.rows, src.cols, CV_8U);//四个方向中为1的最小个数
	cv::Mat max = cv::Mat::zeros(src.rows, src.cols, CV_8U);//四个方向中为1的最大个数
	double num_Min, num_Max;
	double DK[4] = {};
	for (int i = 0; i < src.rows; ++i){
		for (int j = 0; j < src.cols; ++j){
			for (int k = 0; k < 4; ++k){
				uint pix = filtedst[k].at<uchar>(i, j);
				DK[k] = pix;
			}
			num_Max = *(std::max_element)(DK, DK + 4);
			num_Min = *(std::min_element)(DK, DK + 4);
			min.at<uchar>(i, j) = num_Min;
			max.at<uchar>(i, j) = num_Max;
		}
	}
 
	//最小值/最大值
	min.convertTo(min, CV_64F, 1, 0);
	max.convertTo(max, CV_64F, 1, 0);
	cv::Mat ration = cv::Mat::zeros(src.rows, src.cols, CV_64F);
	ration = min / max;
	
	return ration;
}
 
cv::Mat LIG::Calc_Gmap(cv::Mat src, int wsize){
 
	if (wsize % 2 == 0){
		fprintf(stderr, "Please enter odd size!");
		exit(-1);
	}
	if (src.channels() > 1) cvtColor(src, src, CV_RGB2GRAY);
	if (src.type() == CV_8U)
		src.convertTo(src, CV_64F, 1.0, 0);
 
	int Border = (wsize - 1) / 2;
	cv::Mat Newsrc;
	cv::copyMakeBorder(src, Newsrc, Border, Border, Border, Border, cv::BORDER_REFLECT_101);//以边缘为轴，对称
	cv::Mat dst = cv::Mat::zeros(src.size(), src.type());
 
	//Gradient characteristics
	cv::Mat_<double > img_gx, img_gy, img_gvalue;
	Calc_gradient(src, img_gx, img_gy);
	cv::magnitude(img_gx, img_gy, img_gvalue);
	cv::copyMakeBorder(img_gx, img_gx, Border, Border, Border, Border, cv::BORDER_REFLECT_101);
	cv::copyMakeBorder(img_gy, img_gy, Border, Border, Border, Border, cv::BORDER_REFLECT_101);
	cv::copyMakeBorder(img_gvalue, img_gvalue, Border, Border, Border, Border, cv::BORDER_REFLECT_101);
	int ROI_center = Border;
	int nss = 0;
 
	for (int r = Border; r < src.rows + Border; ++r){
		for (int c = Border; c < src.cols + Border; ++c){
			cv::Rect area(c - Border, r - Border, wsize, wsize);
			cv::Mat Slide_Window_gx = img_gx(area);
			cv::Mat Slide_Window_gy = img_gy(area);
			cv::Mat Slide_Window_gvalue = img_gvalue(area);
 
			//First quadrant - NE
			cv::Mat GX1 = get_NE_data(Slide_Window_gx, ROI_center);
			cv::Mat GY1 = get_NE_data(Slide_Window_gy, ROI_center);
			cv::Mat Gvalue1 = get_NE_data(Slide_Window_gvalue, ROI_center);
			int num = 0;
			double gm1 = 0.0;
			double G1 = 0.0;
			for (int i = 0; i < GX1.rows; ++i){
				for (int j = 0; j < GX1.cols; ++j){
					double x = GX1.at<double>(i, j);
					double y = GY1.at<double>(i, j);
					double value = Gvalue1.at<double>(i, j);
					if (x <= 0 && y >= 0){
						gm1 = gm1 + value*value;
						++num;
					}
				}
			}
			if (num>0)
				G1 = gm1 / num;
 
			//The second quadrant - NW
			cv::Mat GX2 = get_NW_data(Slide_Window_gx, ROI_center);
			cv::Mat GY2 = get_NW_data(Slide_Window_gy, ROI_center);
			cv::Mat Gvalue2 = get_NW_data(Slide_Window_gvalue, ROI_center);
			num = 0;
			double gm2 = 0.0;
			double G2 = 0.0;
			for (int i = 0; i < GX2.rows; ++i){
				for (int j = 0; j < GX2.cols; ++j){
					double x = GX2.at<double>(i, j);
					double y = GY2.at<double>(i, j);
					double value = Gvalue2.at<double>(i, j);
					if (x >= 0 && y >= 0){
						gm2 = gm2 + value*value;
						++num;
					}
				}
			}
			if (num > 0)
				G2 = gm2 / num;
 
			//The third quadrant -SW
			cv::Mat GX3 = get_SW_data(Slide_Window_gx, ROI_center);
			cv::Mat GY3 = get_SW_data(Slide_Window_gy, ROI_center);
			cv::Mat Gvalue3 = get_SW_data(Slide_Window_gvalue, ROI_center);
			num = 0;
			double gm3 = 0.0;
			double G3 = 0.0;
			for (int i = 0; i < GX3.rows; ++i){
				for (int j = 0; j < GX3.cols; ++j){
					double x = GX3.at<double>(i, j);
					double y = GY3.at<double>(i, j);
					double value = Gvalue3.at<double>(i, j);
					if (x >= 0 && y <= 0){
						gm3 = gm3 + value*value;
						++num;
					}
				}
			}
			if (num>0)
				G3 = gm3 / num;
 
			//The fourth quadrant -SE
			cv::Mat GX4 = get_SE_data(Slide_Window_gx, ROI_center);
			cv::Mat GY4 = get_SE_data(Slide_Window_gy, ROI_center);
			cv::Mat Gvalue4 = get_SE_data(Slide_Window_gvalue, ROI_center);
			num = 0;
			double gm4 = 0.0;
			double G4 = 0.0;
			for (int i = 0; i < GX4.rows; ++i){
				for (int j = 0; j < GX4.cols; ++j){
					double x = GX4.at<double>(i, j);
					double y = GY4.at<double>(i, j);
					double value = Gvalue4.at<double>(i, j);
					if (x <= 0 && y <= 0){
						gm4 = gm4 + value*value;
						++num;
					}
				}
			}
			if (num>0)
				G4 = gm4 / num;
 
			//Ration
			double G[4] = { G1, G2, G3, G4 };
			double max_G = *(std::max_element)(G, G + 4);
			double min_G = *(std::min_element)(G, G + 4);
			//double gradient_ration = min_G / (max_G + 0.000000001);
			double gradient_ration =  (G1 + G2 + G3 + G4);
 
			dst.at<double>(r - Border, c - Border) = gradient_ration;
 
		}
	}
	src.convertTo(src, CV_8U, 1, 0);
	cv::normalize(dst, dst, 0, 1.0, cv::NORM_MINMAX); //If the image is to be displayed, it should be normalized
	return dst;
}
 
 
void LIG::Calc_LIG(cv::Mat src, cv::Mat &Imap, cv::Mat &Gmap, cv::Mat &IGmap, int wsize, double K){
	Imap = Calc_Imap(src, wsize);
	Gmap = Calc_Gmap(src, wsize);
	IGmap = Imap.mul(Gmap);
}