canny 算法的实现、SSIM评价图片的相似性

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
struct position
{
	int x,y;//行号与列号
};
//滞后阈值化的递归
void search(Mat pixelValue,float lowThreshold,float highThreshold,vector<struct position>&weakEdge,int n)
{
	if (n== 0)
	{
		for (int i = 0; i < weakEdge.size(); i++)
		{
			int x = weakEdge[i].x;
			int y = weakEdge[i].y;
			if (pixelValue.at<float>(x, y) != 255)
			{
				pixelValue.at<float>(x, y) = 0;
			}
		}
	}
	else
	{
		for (int k = 0; k < weakEdge.size(); k++)
		{
			int i = weakEdge[k].x;
			int j = weakEdge[k].y;
			if (pixelValue.at<float>(i - 1, j - 1) >= highThreshold || pixelValue.at<float>(i - 1, j) >= highThreshold || pixelValue.at<float>(i - 1, j + 1) >= highThreshold || pixelValue.at<float>(i, j - 1) >= highThreshold || pixelValue.at<float>(i, j) >= highThreshold || pixelValue.at<float>(i, j + 1) >= highThreshold || pixelValue.at<float>(i + 1, j - 1) >= highThreshold || pixelValue.at<float>(i + 1, j) >= highThreshold || pixelValue.at<float>(i + 1, j + 1) >= highThreshold)
			{
				pixelValue.at<float>(i, j) = 255;
			}
		}
		search(pixelValue, lowThreshold, highThreshold, weakEdge, n - 1);
	}
}
void doubleThreshold(Mat &pixelValue,float lowThreshold,float highThreshold)
{
	vector<struct position> weakEdge;
	for (int i = 1; i < pixelValue.rows - 1; i++)
	{
		for (int j = 1; j < pixelValue.cols - 1; j++)
		{
			if (pixelValue.at<float>(i, j) <= lowThreshold)
			{
			   pixelValue.at<float>(i, j) = 0;
			}
			else if(pixelValue.at<float>(i, j) >= highThreshold)
			{
			    pixelValue.at<float>(i, j) = 255;
			}
			else
			{
				if (pixelValue.at<float>(i - 1, j - 1) >= highThreshold || pixelValue.at<float>(i - 1, j) >= highThreshold || pixelValue.at<float>(i - 1, j + 1) >= highThreshold || pixelValue.at<float>(i, j - 1) >= highThreshold|| pixelValue.at<float>(i, j) >= highThreshold || pixelValue.at<float>(i, j + 1) >= highThreshold || pixelValue.at<float>(i + 1, j - 1) >= highThreshold || pixelValue.at<float>(i + 1, j) >= highThreshold || pixelValue.at<float>(i + 1, j + 1) >= highThreshold)
				{
					pixelValue.at<float>(i, j) = 255;
				}
				else
				{
					struct position p;
					p.x = i;
					p.y = j;
					weakEdge.push_back(p);
				}
			}		
		}
	}
	int n = 50;//迭代次数
	search(pixelValue, lowThreshold, highThreshold, weakEdge,n);
}
Mat myCanny(Mat srcImage)
{
	//高斯平滑滤波
	Mat grayImage;
	GaussianBlur(srcImage, grayImage, Size(5, 5), 0,0);
	//计算像素点的梯度强度与方向,采用Sobel算子
	Mat edge = Mat::zeros(grayImage.rows, grayImage.cols,CV_32FC1);
	Mat pixelValue = edge.clone();
	Mat edge_x = edge.clone();
	Mat edge_y = edge.clone();
	Mat edge_theta = edge.clone();
	float sobel_x, sobel_y, sobel,theta;
	for (int i = 1; i < grayImage.rows - 1; i++)
	{
		for (int j = 1; j < grayImage.cols - 1; j++)
		{
			//计算x方向
			sobel_x = grayImage.at<uchar>(i + 1, j - 1) + 2 * grayImage.at<uchar>(i + 1, j) + grayImage.at<uchar>(i + 1, j + 1) - grayImage.at<uchar>(i - 1, j - 1) - 2 * grayImage.at<uchar>(i - 1, j) - grayImage.at<uchar>(i - 1, j + 1);
			//计算y方向
			sobel_y = grayImage.at<uchar>(i - 1, j + 1) + 2 * grayImage.at<uchar>(i, j + 1) + grayImage.at<uchar>(i + 1, j + 1) - grayImage.at<uchar>(i - 1, j - 1) - 2 * grayImage.at<uchar>(i, j - 1) - grayImage.at<uchar>(i + 1, j - 1);
			sobel = sqrt(sobel_x*sobel_x + sobel_y * sobel_y);
			//判断角度的范围
			theta = atan2f(sobel_y,sobel_x) * 180 / CV_PI;
			if (theta >= 67.5&&theta < 112.5 || theta >= -112.5&&theta <= -67.5)
			{
				theta = 90;
			}
			else if (theta >= -22.5&&theta < 22.5 || theta >=157.5&&theta <= 180||theta>=-180&&theta<-157.5)
			{
				theta = 0;
			}
			else if (theta>=22.5&&theta<67.5||theta>=-157.5&&theta<-112.5)
			{
				theta = 45;
			}
			else
			{
				theta = 135;
			}
			edge_x.at<float>(i, j) = abs(sobel_x);
			edge_y.at<float>(i, j) = abs(sobel_y);
			edge_theta.at<float>(i, j) = theta;
			edge.at<float>(i, j) =sobel;
		}
	}
	//边缘非极大值抑制
	for (int i = 1; i < grayImage.rows - 1; i++)
	{
		for (int j = 1; j < grayImage.cols - 1; j++)
		{
			if (edge_theta.at<float>(i,j)== 0&&edge.at<float>(i,j)==max(edge.at<float>(i,j),max(edge.at<float>(i-1,j), edge.at<float>(i+1, j))))
			{
				pixelValue.at<float>(i, j) = edge.at<float>(i, j);
			}
			else if (edge_theta.at<float>(i, j) == 90 && edge.at<float>(i, j) == max(edge.at<float>(i, j), max(edge.at<float>(i, j-1), edge.at<float>(i, j+1))))
			{
				pixelValue.at<float>(i, j) = edge.at<float>(i, j);
			}
			else if (edge_theta.at<float>(i, j) == 45 && edge.at<float>(i, j) == max(edge.at<float>(i, j), max(edge.at<float>(i-1, j - 1), edge.at<float>(i+1, j + 1))))
			{
				pixelValue.at<float>(i, j) = edge.at<float>(i, j);
			}
			else if (edge_theta.at<float>(i, j) == 135 && edge.at<float>(i, j) == max(edge.at<float>(i, j), max(edge.at<float>(i - 1, j + 1), edge.at<float>(i + 1, j - 1))))
			{
				pixelValue.at<float>(i, j) = edge.at<float>(i, j);
			}
		}
	}
	//滞后阈值化
	float lowThreshold, highThreshold;
	lowThreshold = 9;
	highThreshold = 64;
	doubleThreshold(pixelValue, lowThreshold, highThreshold);
	pixelValue.convertTo(pixelValue, CV_8UC1);
	return pixelValue;
}
double calculateSSIM(Mat opencv_dstImage,Mat my_dstImage)
{
	int r = opencv_dstImage.rows;
	int c = opencv_dstImage.cols;
	double SSIM ,meanValue1, meanValue2, variance1, variance2, covariance,c1,c2;
	meanValue1 = 0;
	meanValue2 = 0;
	variance1 = 0;
	variance2 = 0;
	covariance = 0;
	c1 = 6.5025;
	c2 = 58.5225;
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < c; j++)
		{
			meanValue1 += double(opencv_dstImage.at<uchar>(i, j));
			meanValue2 += double(my_dstImage.at<uchar>(i, j));
		}
	}
	meanValue1 = meanValue1 / (r*c);
	meanValue2 = meanValue2 / (r*c);
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < c; j++)
		{
			variance1 += pow(double(opencv_dstImage.at<uchar>(i, j))-meanValue1,2);
			variance2 += pow(double(my_dstImage.at<uchar>(i, j))-meanValue2,2);
			covariance += (double(opencv_dstImage.at<uchar>(i, j)) - meanValue1)*(double(my_dstImage.at<uchar>(i, j)) - meanValue2);
		}
	}
	variance1 = variance1 / (r*c - 1);
	variance2 = variance2 / (r*c - 1);
	covariance = covariance / (r*c - 1);
	SSIM = (2 * meanValue1*meanValue2 + c1)*(2 * covariance + c2) / ((pow(meanValue1,2)+pow(meanValue2,2)+c1)*(variance1+variance2+c2));
	return SSIM;
}
int main()
{
	Mat grayImage = imread("D:\\figure\\house.jpg",0);
	Mat my_dstImage=myCanny(grayImage);
	Mat opencv_dstImage;
	GaussianBlur(grayImage,grayImage,Size(3,3),1,1);
	Canny(grayImage, opencv_dstImage, 10, 100);
	double SSIM = calculateSSIM(opencv_dstImage, my_dstImage);
	cout << " SSIM " << SSIM << endl;
	imshow("opencv Canny", opencv_dstImage);
	imshow("my Canny", my_dstImage);
	waitKey(0);
	return 0;
}