//otsu算法实现函数
int Otsu(Mat &image)
{
  int width = image.cols;
  int height = image.rows;
  int x = 0, y = 0;
  int pixelCount[256];
  float pixelPro[256];
  int i, j, pixelSum = width * height, threshold = 0;

  uchar* data = (uchar*)image.data;

  //初始化 
  for (i = 0; i < 256; i++)
  {
    pixelCount[i] = 0;
    pixelPro[i] = 0;
  }

  //统计灰度级中每个像素在整幅图像中的个数 
  for (i = y; i < height; i++)
  {
    for (j = x; j<width; j++)
    {
      pixelCount[data[i * image.step + j]]++;
    }
  }


  //计算每个像素在整幅图像中的比例 
  for (i = 0; i < 256; i++)
  {
    pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum);
  }

  //经典ostu算法,得到前景和背景的分割 
  //遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值 
  float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
  for (i = 0; i < 256; i++)
  {
    w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;

    for (j = 0; j < 256; j++)
    {
      if (j <= i) //背景部分 
      {
        //以i为阈值分类，第一类总的概率 
        w0 += pixelPro[j];
        u0tmp += j * pixelPro[j];
      }
      else    //前景部分 
      {
        //以i为阈值分类，第二类总的概率 
        w1 += pixelPro[j];
        u1tmp += j * pixelPro[j];
      }
    }

    u0 = u0tmp / w0;    //第一类的平均灰度 
    u1 = u1tmp / w1;    //第二类的平均灰度 
    u = u0tmp + u1tmp;   //整幅图像的平均灰度 
                //计算类间方差 
    deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);
    //找出最大类间方差以及对应的阈值 
    if (deltaTmp > deltaMax)
    {
      deltaMax = deltaTmp;
      threshold = i;
    }
  }
  //返回最佳阈值; 
  return threshold;
}

int main()
{
  Mat White=imread("white.tif");//读取图像
  int threshold_white = otsu(White);//阈值计算，利用otsu
  cout << "最佳阈值：" << threshold_white << endl;
  Mat thresholded = Mat::zeros(White.size(), White.type());
  threshold(White, thresholded, threshold_white, 255, CV_THRESH_BINARY);//二值化
  vector<vector<Pointcontours;
  vector<Vec4i>hierarchy;
  findContours(thresholded, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);//查找轮廓

  int i = 0;
  int count = 0;
  Point pt[10];//假设有三个连通区域
  Moments moment;//矩
  vector<Point>Center;//创建一个向量保存重心坐标
  for (; i >= 0; i = hierarchy[i][0])//读取每一个轮廓求取重心
  {
    Mat temp(contours.at(i));
    Scalar color(0, 0, 255);
    moment = moments(temp, false);
    if (moment.m00 != 0)//除数不能为0
    {
      pt[i].x = cvRound(moment.m10 / moment.m00);//计算重心横坐标
      pt[i].y = cvRound(moment.m01 / moment.m00);//计算重心纵坐标

    }
      Point p = Point(pt[i].x, pt[i].y);//重心坐标
      circle(White, p, 1, color, 1, 8);//原图画出重心坐标
      count++;//重心点数或者是连通区域数
      Center.push_back(p);//将重心坐标保存到Center向量中
    }
  }
  cout << "重心点个数：" << Center.size() << endl;
  cout << "轮廓数量：" << contours.size() << endl;
  imwrite("Center.tif", White);
}

void getAlikeContours(std::vector<cv::Point> Inputlist, cv::Mat InputImage, std::vector<cv::Point> &Outputlist)
{
 Mat image;
 InputImage.copyTo(image);
 vector<vector<Point> > contours;
 findContours(image, contours, RETR_EXTERNAL, CHAIN_APPROX_NONE);//查找最外层轮廓
 for (int idx = contours.size() - 1; idx >= 0; idx--)
  {
	for (int i = contours[idx].size() - 1; i >= 0; i--)
	{
		if (contours[idx][i].x == 1 || contours[idx][i].y == 1 || contours[idx][i].x == image.cols - 2 || contours[idx][i].y == image.rows - 2)
		{
			swap(contours[idx][i], contours[idx][contours[idx].size() - 1]);
			contours[idx].pop_back();
			
		}
	}
	//可能会存在空的轮廓，把他们删除
	for (int idx = contours.size() - 1; idx >= 0; idx--)
	{
		if (contours[idx].size() == 0) contours.erase(contours.begin() + idx);
	}
 
	while (true)
	{
		if (contours.size() == 0) break;
		if (contours.size() == 1)
		{
			vector<Point> finalList;
			finalList.assign(contours[0].begin(), contours[0].end());
			convexHull(Mat(finalList), Outputlist, true);
			break;
		}
 
		int maxContourIdx = 0;
		int maxContourPtNum = 0;
		for (int index = contours.size() - 1; index >= 0; index--)
		{
			if (contours[index].size() > maxContourPtNum)
			{
				maxContourPtNum = contours[index].size();
				maxContourIdx = index;
			}
		}
		//第二大轮廓
		int secondContourIdx = 0;
		int secondContourPtNum = 0;
		for (int index = contours.size() - 1; index >= 0; index--)
		{
			if (index == maxContourIdx) continue;
			if (contours[index].size() > secondContourPtNum)
			{
				secondContourPtNum = contours[index].size();
				secondContourIdx = index;
			}
		}
		vector<Point> maxlist;
		vector<Point> maxAndseclist;
		vector<Point> maxlistHull;
		vector<Point> maxAndseclistHull;
		maxlist.insert(maxlist.end(), contours[maxContourIdx].begin(), contours[maxContourIdx].end());
		maxAndseclist.insert(maxAndseclist.end(), contours[maxContourIdx].begin(), contours[maxContourIdx].end());
		maxAndseclist.insert(maxAndseclist.end(), contours[secondContourIdx].begin(), contours[secondContourIdx].end());
		convexHull(Mat(maxlist), maxlistHull, true);
		convexHull(Mat(maxAndseclist), maxAndseclistHull, true);
		double maxcontourScore = matchShapes(Inputlist, maxlistHull, CV_CONTOURS_MATCH_I1, 0);
		double maxandseccontourScore = matchShapes(Inputlist, maxAndseclistHull, CV_CONTOURS_MATCH_I1, 0);
		if (maxcontourScore>maxandseccontourScore)
		{
			contours[maxContourIdx].insert(contours[maxContourIdx].end(), contours[secondContourIdx].begin(), contours[secondContourIdx].end());
		}
		contours.erase(contours.begin() + secondContourIdx);
	}
}