/*
這個程式的基本思想是：對輸入的影象進行濾波去掉噪音，然後進行canny邊緣檢測，之後進行膨脹，然後尋找輪廓，對輪廓進行多邊形的逼近，檢測多邊形的點數是否是4而且各個角的的餘弦是否是小於某個值，程式中認為是0.3，然後就判斷該多邊形是四邊形，之後根據這四個點畫出該影象。

ps：我對程式中餘弦定理的使用 感覺公式用錯了
*/

#include "stdafx.h"   
  
#include "cv.h"   
#include "highgui.h"   
#include <stdio.h>   
#include <math.h>   
#include <string.h>   
  
  
int thresh = 50;   
IplImage* img = 0;   
IplImage* img0 = 0;   
CvMemStorage* storage = 0;   
CvPoint pt[4];   
const char* wndname = "Square Detection Demo";   
  
// helper function:   
// finds a cosine of angle between vectors   
// from pt0->pt1 and from pt0->pt2    
double angle( CvPoint* pt1, CvPoint* pt2, CvPoint* pt0 )   
{   
       double dx1 = pt1->x - pt0->x;   
       double dy1 = pt1->y - pt0->y;   
       double dx2 = pt2->x - pt0->x;   
       double dy2 = pt2->y - pt0->y;   
       //1e-10就是“aeb”的形式，表示a乘以10的b次方。   
       //其中b必須是整數，a可以是小數。   
       //?餘弦定理CosB=(a^2+c^2-b^2)/2ac？？所以這裡的計算似乎有問題   
       return (dx1*dx2 + dy1*dy2)/sqrt((dx1*dx1 + dy1*dy1)*(dx2*dx2 + dy2*dy2) + 1e-10);   
}   
  
// returns sequence of squares detected on the image.   
// the sequence is stored in the specified memory storage   
  
CvSeq* findSquares4( IplImage* img, CvMemStorage* storage )   
{   
       CvSeq* contours;   
       int i, c, l, N = 11;   
       CvSize sz = cvSize( img->width & -2, img->height & -2 );   
       IplImage* timg = cvCloneImage( img ); // make a copy of input image   
       IplImage* gray = cvCreateImage( sz, 8, 1 );    
       IplImage* pyr = cvCreateImage( cvSize(sz.width/2, sz.height/2), 8, 3 );   
       IplImage* tgray;   
       CvSeq* result;   
       double s, t;   
       // create empty sequence that will contain points -   
       // 4 points per square (the square's vertices)   
       CvSeq* squares = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvPoint), storage );   
       
       // select the maximum ROI in the image   
       // with the width and height divisible by 2   
       cvSetImageROI( timg, cvRect( 0, 0, sz.width, sz.height ));   
       
       // down-scale and upscale the image to filter out the noise   
       //使用gaussian金字塔分解對輸入影象向下取樣，首先對它輸入的影象用指定濾波器   
       //進行卷積，然後通過拒絕偶數的行與列來下采樣   
       cvPyrDown( timg, pyr, 7 );   
       //函式 cvPyrUp 使用Gaussian 金字塔分解對輸入影象向上取樣。首先通過在影象中插入 0 偶數行和偶數列，然後對得到的影象用指定的濾波器進行高斯卷積，其中濾波器乘以4做插值。所以輸出影象是輸入影象的 4 倍大小。   
       cvPyrUp( pyr, timg, 7 );   
       tgray = cvCreateImage( sz, 8, 1 );   
       
       // find squares in every color plane of the image   
       for( c = 0; c < 3; c++ )   
       {   
           // extract the c-th color plane   
           //函式 cvSetImageCOI 基於給定的值設定感興趣的通道。值 0 意味著所有的通道都被選定, 1 意味著第一個通道被選定等等。   
           cvSetImageCOI( timg, c+1 );   
           cvCopy( timg, tgray, 0 );   
           
           // try several threshold levels   
           for( l = 0; l < N; l++ )   
           {   
               // hack: use Canny instead of zero threshold level.   
               // Canny helps to catch squares with gradient shading      
               if( l == 0 )   
               {   
                   // apply Canny. Take the upper threshold from slider   
                   // and set the lower to 0 (which forces edges merging)    
                   cvCanny( tgray, gray,60, 180, 3 );   
                   // dilate canny output to remove potential   
                   // holes between edge segments    
                   //使用任意結構元素膨脹影象   
                   cvDilate( gray, gray, 0, 1 );   
               }   
               else   
               {   
                   // apply threshold if l!=0:   
                   //        tgray(x,y) = gray(x,y) < (l+1)*255/N ? 255 : 0   
                   //cvThreshold( tgray, gray, (l+1)*255/N, 255, CV_THRESH_BINARY );   
       cvThreshold( tgray, gray, 50, 255, CV_THRESH_BINARY );   
               }   
               
               // find contours and store them all as a list   
               cvFindContours( gray, storage, &contours, sizeof(CvContour),   
                   CV_RETR_LIST, CV_CHAIN_APPROX_SIMPLE, cvPoint(0,0) );   
               
               // test each contour   
               while( contours )   
               {   
                   // approximate contour with accuracy proportional   
                   // to the contour perimeter   
                   //用指定精度逼近多邊形曲線   
                   result = cvApproxPoly( contours, sizeof(CvContour), storage,   
                       CV_POLY_APPROX_DP, cvContourPerimeter(contours)*0.02, 0 );   
                   // square contours should have 4 vertices after approximation   
                   // relatively large area (to filter out noisy contours)   
                   // and be convex.   
                   // Note: absolute value of an area is used because   
                   // area may be positive or negative - in accordance with the   
                   // contour orientation   
                   //cvContourArea 計算整個輪廓或部分輪廓的面積   
                   //cvCheckContourConvexity測試輪廓的凸性                     
                   if( result->total == 4 &&   
                       fabs(cvContourArea(result,CV_WHOLE_SEQ)) > 1000 &&   
                       cvCheckContourConvexity(result) )   
                   {   
                       s = 0;   
                       
                       for( i = 0; i < 5; i++ )   
                       {   
                           // find minimum angle between joint   
                           // edges (maximum of cosine)   
                           if( i >= 2 )   
                           {   
                               t = fabs(angle(   
                               (CvPoint*)cvGetSeqElem( result, i ),   
                               (CvPoint*)cvGetSeqElem( result, i-2 ),   
                               (CvPoint*)cvGetSeqElem( result, i-1 )));   
                               s = s > t ? s : t;   
                           }   
                       }   
                       
                       // if cosines of all angles are small   
                       // (all angles are ~90 degree) then write quandrange   
                       // vertices to resultant sequence    
                       if( s < 0.3 )   
                           for( i = 0; i < 4; i++ )   
                               cvSeqPush( squares,   
                                   (CvPoint*)cvGetSeqElem( result, i ));   
                   }   
                   
                   // take the next contour   
                   contours = contours->h_next;   
               }   
           }   
       }   
       
       // release all the temporary images   
       cvReleaseImage( &gray );   
       cvReleaseImage( &pyr );   
       cvReleaseImage( &tgray );   
       cvReleaseImage( &timg );   
       
       return squares;   
}   
  
  
// the function draws all the squares in the image   
void drawSquares( IplImage* img, CvSeq* squares )   
{   
       CvSeqReader reader;   
       IplImage* cpy = cvCloneImage( img );   
       int i;   
       
       // initialize reader of the sequence   
       cvStartReadSeq( squares, &reader, 0 );   
       
       // read 4 sequence elements at a time (all vertices of a square)   
       for( i = 0; i < squares->total; i += 4 )   
       {   
           CvPoint* rect = pt;   
           int count = 4;   
           
           // read 4 vertices   
           memcpy( pt, reader.ptr, squares->elem_size );   
           CV_NEXT_SEQ_ELEM( squares->elem_size, reader );   
           memcpy( pt + 1, reader.ptr, squares->elem_size );   
           CV_NEXT_SEQ_ELEM( squares->elem_size, reader );   
           memcpy( pt + 2, reader.ptr, squares->elem_size );   
           CV_NEXT_SEQ_ELEM( squares->elem_size, reader );   
           memcpy( pt + 3, reader.ptr, squares->elem_size );   
           CV_NEXT_SEQ_ELEM( squares->elem_size, reader );   
           
           // draw the square as a closed polyline    
           cvPolyLine( cpy, &rect, &count, 1, 1, CV_RGB(0,255,0), 3, CV_AA, 0 );   
       }   
       
       // show the resultant image   
       cvShowImage( wndname, cpy );   
       cvReleaseImage( &cpy );   
}   
  
  
void on_trackbar( int a )   
{   
       if( img )   
           drawSquares( img, findSquares4( img, storage ) );   
}   
  
char* names[] = { "E:\\1.jpg", "E:\\2.jpg", "E:\\3.jpg",   
                     "E:\\4.jpg", "E:\\5.jpg", 0 };   
  
int main(int argc, char** argv)   
{   
       int i, c;   
       // create memory storage that will contain all the dynamic data   
       storage = cvCreateMemStorage(0);   
  
       for( i = 0; names[i] != 0; i++ )   
       {   
           // load i-th image   
           img0 = cvLoadImage( names[i], 1 );   
           if( !img0 )   
           {   
               printf("Couldn't load %s\n", names[i] );   
               continue;   
           }   
           img = cvCloneImage( img0 );   
           
           // create window and a trackbar (slider) with parent "image" and set callback   
           // (the slider regulates upper threshold, passed to Canny edge detector)    
           cvNamedWindow( wndname,0 );   
           cvCreateTrackbar( "canny thresh", wndname, &thresh, 1000, on_trackbar );   
           
           // force the image processing   
           on_trackbar(0);   
           // wait for key.   
           // Also the function cvWaitKey takes care of event processing   
           c = cvWaitKey(0);   
           // release both images   
           cvReleaseImage( &img );   
           cvReleaseImage( &img0 );   
           // clear memory storage - reset free space position   
           cvClearMemStorage( storage );   
           if( c == 27 )   
               break;   
       }   
       
       cvDestroyWindow( wndname );   
       
       return 0;   
}