#include "pch.h"
#include <iostream> 

#include <opencv2/opencv.hpp>  
using namespace std;
using namespace cv;

// 拉普拉斯锐化函数
void LaplacianSharpDeal(const Mat &src, Mat &dst) {
 if (!src.data)return;
 for (int i = 0; i < src.rows; ++i)
  for (int j = 0; j < src.cols; ++j) {
   float a;
   if (i > 1 && i < src.rows - 1 && j > 1 && j < src.cols - 1) {
    a = 5 * (float)src.at<uchar>(i, j) - (float)src.at<uchar>(i - 1, j) - (float)src.at<uchar>(i, j - 1) -
     (float)src.at<uchar>(i, j + 1) - (float)src.at<uchar>(i + 1, j);
   }
   else {//边缘赋值
    a = src.at<uchar>(i, j);
   }
   if (a > 255 || a < 0) {
    dst.at<uchar>(i, j) = src.at<uchar>(i, j);
   }
   else {
    dst.at<uchar>(i, j) = a;
   }
  }
}

// 基本全局阈值方法函数
int BasicGlobalThreshold(Mat src, float oldValue)
{ 
 int cols = src.cols;
 int rows = src.rows;
 float G1 = 0;
 float G2 = 0;
 float g1 = 0;
 float g2 = 0;
 float u1 = 0;
 float u2 = 0;
 float T0 = 0;
 // 计算灰度直方图分布，统计像素数和频率
 for (int i = 0; i < rows; i++)
 {
  for (int j = 0; j < cols; j++)
  {
   if (src.at<uchar>(i, j) > oldValue)
   {
    G1 += src.at<uchar>(i, j);
    g1 += 1;
   }
   else
   {
    G2 += src.at<uchar>(i, j);
    g2 += 1;
   }
  }
 }
 u1 = G1 / g1;
 u2 = G2 / g2;
 T0 = (u1 + u2) / 2;
 std::cout << T0 << std::endl;
 if (abs(oldValue - T0) < 0.1) {
  return T0;
 }
 else
 {
  BasicGlobalThreshold(src, T0);
 }
}

// Otsu方法函数
int Otsu(Mat src)
{
 int cols = src.cols;
 int rows = src.rows;
 int nPixelNum = cols * rows;
 // 初始化
 int pixelNum[256];
 double probability[256];
 for (int i = 0; i < 256; i++)
 {
  pixelNum[i] = 0;
  probability[i] = 0.0;
 }
 // 统计像素数和频率
 for (int i = 0; i < rows; i++)
 {
  for (int j = 0; j < cols; j++)
  {
   pixelNum[src.at<uchar>(i, j)]++;
  }
 }
 for (int i = 0; i < 256; i++)
 {
  probability[i] = (double)0.1*pixelNum[i] / nPixelNum;
 }
 // 计算
 int Threshold = 0;          // 最佳阈值
 double MaxDelta = 0.0;      // 最大类间方差
 double Mean_0 = 0.0;        // 左边平均值
 double Mean_1 = 0.0;        // 右边平均值
 double Delta = 0.0;         // 类间方差
 double Mean_0_temp = 0.0;   // 左边平均值中间值
 double Mean_1_temp = 0.0;   // 右边平均值中间值
 double Probability_0 = 0.0;       // 左边频率值
 double Probability_1 = 0.0;       // 右边频率值
 for (int j = 0; j < 256; j++)
 {
  for (int i = 0; i < 256; i++)
  {
   if (i < j)// 前半部分
   {
    Probability_0 += probability[i];
    Mean_0_temp += i * probability[i];
   }
   else      // 后半部分
   {
    Probability_1 += probability[i];
    Mean_1_temp += i * probability[i];
   }
  }
  // 计算平均值
  // Mean_0_teamp计算的是前半部分的灰度值的总和除以总像素数，
  // 所以要除以前半部分的频率才是前半部分的平均值，后半部分同样
  Mean_0 = Mean_0_temp / Probability_0;
  Mean_1 = Mean_1_temp / Probability_1;
  Delta = (double)(Probability_0 * Probability_1 * pow((Mean_0 - Mean_1), 2));
  if (Delta > MaxDelta)
  {
   MaxDelta = Delta;
   Threshold = j;
  }
  // 相关参数归零
  Probability_0 = 0.0;
  Probability_1 = 0.0;
  Mean_0_temp = 0.0;
  Mean_1_temp = 0.0;
  Mean_0 = 0.0;
  Mean_1 = 0.0;
  Delta = 0.0;
 }
 return Threshold;
}

void main() {
 Mat image = imread("A1.bmp", 0);
 Mat image1,image2;
 Mat image3(image.size(), image.type());
 Mat image4(image.size(), image.type());

 std::cout << "基本全局阈值方法" << std::endl;

 int OstuThreshold1 = BasicGlobalThreshold(image, 0.01);
 int OstuThreshold2 = Otsu(image);

 std::cout << "Otsu方法" << std::endl;
 std::cout << OstuThreshold2 << std::endl;
 threshold(image, image1, OstuThreshold1, 255, CV_THRESH_OTSU);
 threshold(image, image2, OstuThreshold2, 255, CV_THRESH_OTSU);

 LaplacianSharpDeal(image2, image3);
 LaplacianSharpDeal(image, image4);
 
 imshow("基本全局阈值方法", image1);
 imshow("Otsu方法", image2);
 imshow("先阈值分割后边缘检测", image3);
 imshow("直接对图像进行边缘检测", image4);
 waitKey();
}