#include <iostream>
#include<opencv2/opencv.hpp>
#include<opencv2/photo.hpp>
using namespace cv;
using namespace std;
void inpaint_test()
{
    Mat img=imread("D:/Qt/MyImage/jingsai_1.jpg",0);
    imshow("original image",img);
    Mat mask=Mat::zeros(img.rows,img.cols,img.type());
    mask(Rect(377,228,(438-377),(269-228)))=1;
    Mat img_c;
    cv::inpaint(img,mask,img_c,3,INPAINT_TELEA);

    imshow("clear image",img_c);
    waitKey();
}
Mat getHist(const Mat &img)//img必須是單通道圖像，img的概率直方圖
{
    if(img.channels()!=1)
    {
        cout<<"channels of input image !=1"<<endl;
        exit(1) ;
    }
    int L=256;//灰度級(jí)數(shù)
     Mat hist(1,L,CV_64F,Scalar::all(0));
     int v=0;
    for(int i=0;i<img.rows;i++)
        for(int j=0;j<img.cols;j++)
        {
            v=img.ptr<uchar>(i)[j];
            hist.ptr<double>(0)[v]+=1;
        }
    return hist;
}
//計(jì)算累積概率分布函數(shù)，Mat &pi是概率分布，本函數(shù)可以指定某一灰度范圍內(nèi)(k1~k2)的累積概率
double getLocAccHist(const Mat&pi,int k1,int k2)
{
    double P=0;
    for(int k=k1;k<=k2;k++)
    {
        P +=pi.ptr<double>(0)[k];
    }
    return P;
}

Mat getMeanPerGrayLevel(const Mat& pi)//計(jì)算每個(gè)灰度級(jí)上的像素均值
{
    int L=256;
    Mat ipi(pi.size(),pi.type(),Scalar::all(0));
//    m1k.ptr<double>(0)[0]=pi.ptr<double>(0)[0]*0;
    for(int k=1;k<L;k++)
    {
        ipi.ptr<double>(0)[k]=ipi.ptr<double>(0)[(k-1)]+pi.ptr<double>(0)[k]*k;
    }
    return ipi;
}
/*計(jì)算某一灰度級(jí)范圍內(nèi)的所有像素灰度均值,下面是參數(shù)說(shuō)明
k1:灰度級(jí)小端；k2：灰度級(jí)大端；
hist_img:原圖像的直方圖
*/
double getMeanvalueOfC(int k1,int k2,const Mat&hist_img)
{
    int pixTotalNum=0;
    double pixTotalValue=0;
    double mk1k2=0;
    for(int k=k1;k<=k2;k++)
    {
        pixTotalNum +=hist_img.ptr<double>(0)[k];//計(jì)算k1與k2之間所有像素的個(gè)數(shù)
        //計(jì)算k1與k2之間所有像素的灰度值之和
        pixTotalValue +=k*hist_img.ptr<double>(0)[k];
    }
    if(pixTotalNum!=0)
        mk1k2=pixTotalValue/pixTotalNum;
    return mk1k2;
}

void getOstu2Threshvalue(const Mat&img,int &k11,int &k22)
{
    Mat hist_img=getHist(img);
    Mat pi=hist_img/(img.rows*img.cols);//1.灰度級(jí)概率分布
    //3.全局均值、全局方差
    Scalar means,stddev;
    meanStdDev(img,means,stddev);
    double mG=means[0];
    //5.計(jì)算m1,第一類(lèi)C1中的像素灰度均值
    int L=256;

    double SigmaB2=0;
    for(int k1=1;k1<L-2;k1++)
    {
        for(int k2=k1+1;k2<L-1;k2++)
        {
            //6.計(jì)算mk,第一類(lèi)C2中的像素灰度均值
            double m1=getMeanvalueOfC(0,k1,hist_img);
            double m2=getMeanvalueOfC(k1+1,k2,hist_img);
            double m3=getMeanvalueOfC(k2+1,L-1,hist_img);
            //
            double P1=getLocAccHist(pi,0,k1);
            double P2=getLocAccHist(pi,k1+1,k2);
            double P3=getLocAccHist(pi,k2+1,L-1);
            //7.計(jì)算類(lèi)間方差SigmaB2
            double temp=saturate_cast<double>(P1*(m1-mG)*(m1-mG)
                                              +P2*(m2-mG)*(m2-mG)
                                              +P3*(m3-mG)*(m3-mG));
            if(SigmaB2<temp)
            {
                SigmaB2=temp;
                k11=k1;
                k22=k2;
            }
        }
    }
}
Mat  myThreshHold(const Mat&img,int k1,int k2,int v1=0,int v2=90,int v3=255)
{
    Mat dst(img.size(),img.type(),Scalar::all(0));
    for(int i=0;i<img.rows;i++)
        for(int j=0;j<img.cols;j++)
        {
            int v=img.ptr<uchar>(i)[j];
            if(v<=k1)
                dst.ptr<uchar>(i)[j]=v1;
            else if(v<=k2)
                dst.ptr<uchar>(i)[j]=v2;
            else
                dst.ptr<uchar>(i)[j]=v3;
        }
    return dst;
}
int main()
{
    Mat img=imread("D:/Qt/MyImage/ch10/Fig1043(a)(yeast_USC).tif",0);
//    GaussianBlur(img,img,Size(7,7),0);
//    Mat kernel=(Mat_<int>(3,3)<<-1,-1,-1,
//                                -1,9,-1,
//                                -1,-1,-1);
//    filter2D(img,img,CV_8U,kernel);
    imshow("original image",img);
    int k11,k22;
    getOstu2Threshvalue(img,k11,k22);
    Mat dst=myThreshHold(img,k11,k22);

    imshow("binary image",dst);
    waitKey();

    return 0;
}