#include "stdio.h"  
#include "cv.h"  
#include "highgui.h"  
#include "Math.h"
#define max_corners 4


int Otsu(IplImage* src);

int main(int argc, char*argv[])
{
    IplImage* img = cvLoadImage("C:\\Users\\aoe\\Desktop\\the last\\20160716_3DMeasurement\\1.jpg", 0);
    cvNamedWindow("img", 1);
    cvShowImage("img"
 
, img);

    double sum = 0.0;//求影象灰度平均值
    for (int y = 0; y < img->height; y++)
    {
        for (int x = 0; x < img->width; x++)
        {
            CvScalar cs = { 0.0 };
            cs=cvGet2D(img, y, x);
            sum += cs.val[0];
        }
    }
    double average = 0.0;//as.val[0]是double型別的資料
    average = sum / ((img->height)*(
 
img->width));
    printf("sum=%lf\n", sum);
    printf("average=%lf\n", average);
    IplImage* aver = cvCreateImage(cvGetSize(img), 8, 1);//生成一幅畫素值均為上一幅影象的畫素的均值的影象
    for (int y = 0; y < aver->height; y++)
    {
        for (int x = 0; x < aver->width; x++)
        {
            CvScalar cs = { 0
 
.0 };
            cs.val[0]=average;
            cvSet2D(aver, y, x,cs);
        }
    }
    cvNamedWindow("aver", 1);
    cvShowImage("aver", aver);


    IplImage* next = cvCreateImage(cvGetSize(img), 8, 1);//生成一幅畫素值減去影象畫素的均值的影象
    for (int y = 0; y < next->height; y++)
    {
        for (int x = 0; x < next->width; x++)
        {
            CvScalar cs = {0.0};
            cs = cvGet2D(img, y, x);
            CvScalar ds = { 0.0 };
            if (cs.val[0] > average)
                ds.val[0] = average;
            else
                ds.val[0] = cs.val[0];          
            cvSet2D(next, y, x, ds);
        }
    }
    cvNamedWindow("next", 1);
    cvShowImage("next", next);

    IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    int threshold = Otsu(next);//最大類間方差閾值分割
    printf("threshold = %d\n", threshold);
    cvThreshold(next, dst, threshold, 255, CV_THRESH_BINARY);

    cvNamedWindow("dst", 1);
    cvShowImage("dst", dst);

/*  IplImage* im_median_filter = cvCreateImage(cvGetSize(img), IPL_DEPTH_8U, 1);
    cvSmooth(dst, im_median_filter, CV_MEDIAN, 3);//中值濾波
    cvNamedWindow("im_median_filter", 1);
    cvShowImage("im_median_filter", im_median_filter);*/

    CvRect roi = cvRect(30, 30, 120, 120);//去除複雜背景
    IplImage* img1 = cvCreateImage(cvGetSize(dst), dst->depth, dst->nChannels);
    for (int y = 0; y < img1->height; y++)
    {
        for (int x = 0; x < img1->width; x++)
        {
            CvScalar cs = (255);
            cvSet2D(img1, y, x, cs);
        }
    }
    CvRect roi1 = cvRect(30, 30, 120, 120);
    cvNamedWindow("img1");
    cvShowImage("img1", img1);

    cvSetImageROI(dst, roi);
    cvSetImageROI(img1, roi1);
    cvCopy(dst, img1);
    cvResetImageROI(dst);
    cvResetImageROI(img1);

    cvNamedWindow("result", 1);
    cvShowImage("result", img1);

    IplImage*edge = cvCreateImage(cvGetSize(img), 8, 1);//canny邊緣檢測
    int edgeThresh = 1;
    cvCanny(img1, edge, edgeThresh, edgeThresh * 3, 3);

    cvNamedWindow("canny", 1);
    cvShowImage("canny", edge);
    int count = 0;
    for (int yy = 0; yy < edge->height; yy++)//統計邊緣影象中共有多少個黑色畫素點
    {
        for (int xx = 0; xx < edge->width; xx++)
        {
            //CvScalar ss = (255);
            double ds = cvGet2D(edge, yy, xx).val[0];
            if (ds == 0)
                count++;
        }
    }
    int dianshu_threshold = (176 * 144 - count) / 4;//將白色畫素點數的四分之一作為hough變換的閾值
    IplImage* houghtu = cvCreateImage(cvGetSize(edge), IPL_DEPTH_8U, 1);//hough直線變換
    CvMemStorage*storage = cvCreateMemStorage();
    CvSeq*lines = 0;
    int i, j, k, m, n;
    while (true)//迴圈找出合適的閾值，使檢測到的直線的數量在8-12之間
    {
        lines = cvHoughLines2(edge, storage, CV_HOUGH_STANDARD, 1, CV_PI / 180, dianshu_threshold, 0, 0);
        int line_number = lines->total;
        if (line_number <8)
        {
            dianshu_threshold = dianshu_threshold - 2;
        }
        else if (line_number > 12)
        {
            dianshu_threshold = dianshu_threshold + 1;
        }
        else
        {
            printf("line_number=%d\n", line_number);
            break;
        }
    }

    int A = 1;
    double B = CV_PI / 180;

    while (1)
    {
        for (i = 0; i <lines->total; i++)//將多條非常相像的直線剔除
        {
            for (j = 0; j < lines->total; j++)
            {
                if (j != i)
                {
                    float*line1 = (float*)cvGetSeqElem(lines, i);
                    float*line2 = (float*)cvGetSeqElem(lines, j);
                    float rho1 = line1[0];
                    float threta1 = line1[1];
                    float rho2 = line2[0];
                    float threta2 = line2[1];
                    if (abs(rho1 - rho2) < A && abs(threta1 - threta2) < B)
                        cvSeqRemove(lines, j);
                }
            }
        }
        printf("lines->total=%d\n", lines->total);
        if ((lines->total > 4)&&(A<8))//剔除一圈後如何直線的數量大於4，則改變A和B，繼續刪除相似的直線
        {
            A = A + 1;
            B = B + CV_PI / 180;
        }
        else
        {
            printf("lines->total=%d\n", lines->total);
            break;
        }
    }




    for (k = 0; k < lines->total; k++)//畫出直線
    {
        float*line = (float*)cvGetSeqElem(lines, k);
        float rho = line[0];//r=line[0]
        float threta = line[1];//threta=line[1]
        CvPoint pt1, pt2;
        double a = cos(threta), b = sin(threta);
        double x0 = a*rho;
        double y0 = b*rho;
        pt1.x = cvRound(x0 + 100 * (-b));//定義直線的終點和起點，直線上每一個點應該滿足直線方程r=xcos(threta)+ysin(threta);
        pt1.y = cvRound(y0 + 100 * (a));
        pt2.x = cvRound(x0 - 1200 * (-b));
        pt2.y = cvRound(y0 - 1200 * (a));
        cvLine(houghtu, pt1, pt2, CV_RGB(0, 255, 255), 1, 8);
    }
    const double C = CV_PI / 3;
    const double D = CV_PI / 1.5;
    int num = 0;
    CvPoint arr[12] = { { 0, 0 } };
    for (m = 0; m < lines->total; m++)//畫出直線的交點
    {
        for (n = 0; n < lines->total; n++)
        {
            if (n != m)
            {
                float*Line1 = (float*)cvGetSeqElem(lines, m);
                float*Line2 = (float*)cvGetSeqElem(lines, n);
                float Rho1 = Line1[0];
                float Threta1 = Line1[1];
                float Rho2 = Line2[0];
                float Threta2 = Line2[1];
                if ((abs(Threta1 - Threta2) > C) && (abs(Threta1 - Threta2) <D))
                {
                    double a1 = cos(Threta1), b1 = sin(Threta1);
                    double a2 = cos(Threta2), b2 = sin(Threta2);
                    CvPoint pt;
                    pt.x = (Rho2*b1 - Rho1*b2) / (a2*b1 - a1*b2);//直線交點公式
                    pt.y = (Rho1 - a1*pt.x) / b1;
                    cvCircle(houghtu, pt, 3, CV_RGB(255, 255, 0));
                    arr[num++] = pt;//將點的座標儲存在一個數組中
                }
            }

        }
    }
    printf("num=%d\n", num);
    printf("arr[0].x=%d\n", arr[0].x);
    printf("arr[0].y=%d\n", arr[0].y);
    printf("arr[1].x=%d\n", arr[1].x);
    printf("arr[1].y=%d\n", arr[1].y);
    printf("arr[2].x=%d\n", arr[2].x);
    printf("arr[2].y=%d\n", arr[2].y);
    printf("arr[3].x=%d\n", arr[3].x);
    printf("arr[3].y=%d\n", arr[3].y);
    printf("arr[4].x=%d\n", arr[4].x);
    printf("arr[4].y=%d\n", arr[4].y);
    printf("arr[5].x=%d\n", arr[5].x);
    printf("arr[5].y=%d\n", arr[5].y);
    printf("arr[6].x=%d\n", arr[6].x);
    printf("arr[6].y=%d\n", arr[6].y);
    printf("arr[7].x=%d\n", arr[7].x);
    printf("arr[7].y=%d\n", arr[7].y);
    printf("arr[8].x=%d\n", arr[8].x);
    printf("arr[8].y=%d\n", arr[8].y);
    printf("arr[9].x=%d\n", arr[9].x);
    printf("arr[9].y=%d\n", arr[9].y);
    printf("arr[10].x=%d\n", arr[10].x);
    printf("arr[10].y=%d\n", arr[10].y);
    printf("arr[11].x=%d\n", arr[11].x);
    printf("arr[11].y=%d\n", arr[11].y);

    CvPoint arr1[12] = { { 0, 0 } };//將重複的角點剔除
    int num1 = 0;
    for (int r = 0; r < 12; r++)
    {
        int s = 0;
        for (; s < num1; s++)
        {
            if (abs(arr[r].x - arr1[s].x) <= 2 && abs(arr[r].y - arr1[s].y) <= 2)
                break;
        }
        if (s == num1)
        {
            arr1[num1] = arr[r];
            num1++;
        }

    }


    printf("num1=%d\n", num1);
    printf("arr1[0].x=%d\n", arr1[0].x);
    printf("arr1[0].y=%d\n", arr1[0].y);
    printf("arr1[1].x=%d\n", arr1[1].x);
    printf("arr1[1].y=%d\n", arr1[1].y);
    printf("arr1[2].x=%d\n", arr1[2].x);
    printf("arr1[2].y=%d\n", arr1[2].y);
    printf("arr1[3].x=%d\n", arr1[3].x);
    printf("arr1[3].y=%d\n", arr1[3].y);
    printf("arr1[4].x=%d\n", arr1[4].x);
    printf("arr1[4].y=%d\n", arr1[4].y);
    printf("arr1[5].x=%d\n", arr1[5].x);
    printf("arr1[5].y=%d\n", arr1[5].y);
    printf("arr1[6].x=%d\n", arr1[6].x);
    printf("arr1[6].y=%d\n", arr1[6].y);
    printf("arr1[7].x=%d\n", arr1[7].x);
    printf("arr1[7].y=%d\n", arr1[7].y);
    printf("arr1[8].x=%d\n", arr1[8].x);
    printf("arr1[8].y=%d\n", arr1[8].y);
    printf("arr1[9].x=%d\n", arr1[9].x);
    printf("arr1[9].y=%d\n", arr1[9].y);
    printf("arr1[10].x=%d\n", arr1[10].x);
    printf("arr1[10].y=%d\n", arr1[10].y);
    printf("arr1[11].x=%d\n", arr1[11].x);
    printf("arr1[11].y=%d\n", arr1[11].y);
    for (int w = 0; w < 5; w++)
    {
        CvPoint ps;
        ps = arr1[w];
        cvCircle(img, ps, 3, CV_RGB(255, 0, 0));
    }
    cvNamedWindow("img", 1);
    cvShowImage("img", img);
    cvNamedWindow("houghtu", 1);
    cvShowImage("houghtu", houghtu);
    cvWaitKey(-1);

    cvReleaseImage(&img);
    //cvReleaseImage(&edge);

    return 0;
}


int Otsu(IplImage* src)
{
    int height = src->height;
    int width = src->width;

    //histogram    
    float histogram[256] = { 0 };
    for (int i = 0; i < height; i++)
    {
        unsigned char* p = (unsigned char*)src->imageData + src->widthStep * i;
        for (int j = 0; j < width; j++)
        {
            histogram[*p++]++;
        }
    }
    //normalize histogram    
    int size = height * width;
    for (int i = 0; i < 256; i++)
    {
        histogram[i] = histogram[i] / size;
    }

    //average pixel value    
    float avgValue = 0;
    for (int i = 0; i < 256; i++)
    {
        avgValue += i * histogram[i];  //整幅影象的平均灰度  
    }

    int threshold;
    float maxVariance = 0;
    float w = 0, u = 0;
    for (int i = 0; i < 256; i++)
    {
        w += histogram[i];  //假設當前灰度i為閾值, 0~i 灰度的畫素(假設畫素值在此範圍的畫素叫做前景畫素) 所佔整幅影象的比例  
        u += i * histogram[i];  // 灰度i 之前的畫素(0~i)的平均灰度值： 前景畫素的平均灰度值  

        float t = avgValue * w - u;
        float variance = t * t / (w * (1 - w));
        if (variance > maxVariance)
        {
            maxVariance = variance;
            threshold = i;
        }
    }

    return threshold;
}