用OPENCV视觉解数独

用OPENCV视觉解数独

 

2010-06-29

 

看到增强视觉网站上介绍老外用视觉解SUDOKU(http://www.cvchina.info/2011/05/29/video-sudoku-solver/)，觉得应该不难，于是用OPENCV和训练好的数字分类器，也试着做一个，纯属娱乐

基本思路如下：

一，定位网格，
 1，寻找图像中的矩形，使用OPENCV中经典找矩形（因为速度问题使用C版本）代码就可以，将条件降低（CANNY的对比度和角度是否直角）以便找到更多的矩形。
 2，由于不一定能找到所有的矩形，需要根据找到的矩形对网格进行推算，推算出的矩形被称为“虚拟矩形”
 3，最后对网格中的虚拟矩形根据上下的矩形进行位置修正

二，数字分析
 网格中存在两种可能，有数据和没有数据。我们通过对比度和轮廓大小过滤掉没有数据的网格。
 对剩下的网格进行数字分析，使用SVM手写数字分类器，得到该网格的数据

 

三，求解
 剩下的问题就简单了，对得到的SUDOKU题目进行求解就可以啦

 

 

 

效果吗，对网格和数字比较清晰的图像效果不错，从模糊的图片就差一些了，主要是数字和边框不清晰，导致矩形寻找和数字判断错误，可以根据图像的特征进行优化，但如果能做到各种图片通吃就困难多了。

手写数字的图片（蓝色为找到的网格，绿色为推算出的虚拟网格）：

 

识别结果：

 

印刷数字的图片（蓝色为找到的网格，绿色为推算出的虚拟网格）：

 

失败的例子，由于网格不清晰，导致无法找到足够的矩形网格：

 

 

主要代码如下（不包括SUDOKU求解部分），

 

 

// Sudoku detect

#include "stdafx.h"

#ifdef _CH_
#pragma package <opencv>
#endif

#ifndef _EiC
#include "opencv2/opencv.hpp"

#include <stdio.h>
#include <math.h>
#include <string.h>
#endif

#include <vector>
#include <list>

using namespace cv;

#define SHOW_IMG 0

const double MIN_CONTOUR_PROPORTION = 0.02;
const int CANNY_LOW_THRESH = 0;
const int CANNY_HIGH_THRESH = 30;
//int thresh = 50;

const double MIN_ANGLE = 0.1;
//double MIN_ANGLE = 0.05;

IplImage* img = 0;
IplImage* img0 = 0;
CvMemStorage* storage = 0;
CvPoint pt[4];
const char* wndname = "Square Detection Demo";


class FilterRect
{
public:
	Rect rc;
	bool useFlag;

	FilterRect()
	{
		rc = Rect(0, 0, 0, 0);
		useFlag = false;
	}
};

vector<Rect> allRects;
vector<FilterRect> filterRects;

const int SUDOKU_GRID_NUM = 9;
const unsigned int DIFF_THRESHOLD = 16;

class RtTrainData
{
public:
	float data[64];
	int result;
};

class ValidRect
{
public:
	Rect rc;
	int row;
	int col;
	bool blankFlag;
	bool virtualFlag;
	Rect contourRc;

	double contrast;
	double areaProportion;

	ValidRect()
	{
		rc = Rect(0, 0, 0, 0);
		contourRc = Rect(0, 0, 0, 0);
		row = 0;
		col = 0;
		blankFlag = true;
		virtualFlag = true;

		contrast = 0.0;
		areaProportion = 0.0;
	}
};

ValidRect validRects[81];

typedef struct rectWithArea
{
	CvPoint pt[4];
	double area;

}RectWithArea;


bool compareRects(const Rect& rc1, const Rect& rc2, bool& replaceFirst)
{
	int dx = rc1.x - rc2.x;
	int dy = rc1.y - rc2.y;
	int dw = rc1.width - rc2.width;
	int dh = rc1.height - rc2.height;

	//if((dx*dx + dy*dy + dw*dw + dh*dh) < DIFF_THRESHOLD)
	if((dx*dx + dy*dy + dw*dw + dh*dh) < DIFF_THRESHOLD)
	{
		if(rc1.area() > rc2.area())
			replaceFirst = true;
		else
			replaceFirst = false;
		return true;
	}
	else
		return false;
}

double getBkInfo(Mat& src, float& highProportion, Scalar& contrast)
{
	Mat temp, high, low;
	Rect rc(src.cols / 10, src.rows / 10, src.cols * 0.9, src.rows * 0.9);
	temp = src(rc);

	double thresh = threshold(temp, high, 0, 0, CV_THRESH_OTSU); //求出分割阈值
	threshold(temp, high, thresh, 0, CV_THRESH_TOZERO);   //分割高亮部分
	threshold(temp, low, thresh, 0, CV_THRESH_TOZERO_INV);   //分割低暗部分

	int highCount = cv::countNonZero(high);
	highProportion = (float)highCount / src.total(); 
	//highProportion = (float)highCount / (src.rows * src.cols);

	Scalar lowValue = mean(low);
	Scalar highValue = mean(high);
	contrast = (highValue - lowValue) / lowValue;

	double ret = abs(contrast.val[0]);
	return ret;
}

// 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;
    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 )
{
	const int MAX_CONTOUR_SIZE = img->width * img->height / 80;
	const int MIN_CONTOUR_SIZE = MAX_CONTOUR_SIZE / 16;

	int contours_num = 0;
    CvSeq* contours;
    int i, c, l, N = 1; //11;
    CvSize sz = cvSize( img->width & -2, img->height & -2 );
	
    IplImage* timg = cvCloneImage( img ); // make a copy of input image
    //IplImage* timg2 = 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(RectWithArea), 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 ));
    
	//Only care RED
    //cvSetImageCOI( timg2, 3 );

	//cvCopy(timg2, gray, NULL);

    // down-scale and upscale the image to filter out the noise
    cvPyrDown( timg, pyr, 7 );
    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( 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, CANNY_LOW_THRESH, CANNY_HIGH_THRESH, 5 );
                // 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 );
            }

 			//imshow("gray", gray);
			//waitKey(0);
           
            // 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 )
            {
				contours_num++;

                // 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
                if( result->total == 4 &&
                    fabs(cvContourArea(result,CV_WHOLE_SEQ)) > MIN_CONTOUR_SIZE && fabs(cvContourArea(result,CV_WHOLE_SEQ)) < MAX_CONTOUR_SIZE &&
                    cvCheckContourConvexity(result) )
                {
                    s = 0;
                    
                    for( i = 0; i < 5; i++ )
                    {
                        // find minimum angle between joint 检查3个角就足够了（对四边形来说如果三个角为直角，则第四个也为直角）
                        // 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 < MIN_ANGLE )
					{
 						RectWithArea rwa;
						for(int i=0; i<4; i++)
							rwa.pt[i] = *(CvPoint*)cvGetSeqElem( result, i );
						rwa.area = fabs(cvContourArea(result,CV_WHOLE_SEQ));
                        cvSeqPush( squares, &rwa);

						//Add to all rects
						Rect rect(rwa.pt[0], rwa.pt[2]);
						allRects.push_back(rect);
					}
                }
                
                // 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

/* Sort 2d points in top-to-bottom left-to-right order */
static int cmp_func( const void* _a, const void* _b, void* userdata )
{
    RectWithArea* a = (RectWithArea*)_a;
    RectWithArea* b = (RectWithArea*)_b;
    
    return a->area < b->area ? -1 : a->area > b->area ? 1 : 0;
}

void drawSquares( IplImage* img, CvSeq* squares )
{
	cvSeqSort( squares, cmp_func, 0);

    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 ++ )
    {
		RectWithArea rwa;
        int count = 4;
        
        // draw the square as a closed polyline 
        memcpy(&rwa, reader.ptr, squares->elem_size);

		CvPoint* rect = pt;
		pt[0] = rwa.pt[0];
		pt[1] = rwa.pt[1];
		pt[2] = rwa.pt[2];
		pt[3] = rwa.pt[3];

        cvPolyLine( cpy, &rect, &count, 1, 1, CV_RGB(0,255,0), 3, CV_AA, 0 );
        CV_NEXT_SEQ_ELEM( squares->elem_size, reader );

		//if(i>10)
		//	break;
   }
    
    // show the resultant image
    cvShowImage( wndname, cpy );
    cvReleaseImage( &cpy );
}


void filterDupRect()
{
	for(unsigned int i=0; i<allRects.size(); i++)
	{
		bool insertFlag = true;
		for(unsigned int j=0; j<filterRects.size(); j++)
		{
			bool replaceFlag = false;
			if(compareRects(filterRects[j].rc, allRects[i], replaceFlag))
			{
				if(replaceFlag)
					filterRects[j].rc = allRects[i]; //Using smaller one

				insertFlag = false; //Have similiar, not save
				break;
			}
		}

		if(insertFlag)
		{
			FilterRect fr;
			fr.rc = allRects[i];
			filterRects.push_back(fr);
		}
	}
}

int findTopLeftRect(Mat& src)
{
	int min_dist = src.size().width + src.size().height;
	int index = -1;

	for(unsigned int j = 0; j < filterRects.size(); j ++ )
	{
		int dist = filterRects[j].rc.x + filterRects[j].rc.y;
		if(dist < min_dist)
		{
			min_dist = dist;
			index = j;
		}
	}

	return index;
}

void updateValidRect(int col, int row, Rect rc, bool virtualFlag)
{
	int index = col + row * 9;
	validRects[index].col = col;
	validRects[index].row = row;
	validRects[index].rc = rc;
	//validRects[index].processFlag = true;
	validRects[index].virtualFlag = virtualFlag;
}

int findNearest(Rect& rc, int offsetX, int offsetY)
{
	int min_dist = 1000;
	int index = -1;

	for(unsigned int j = 0; j < filterRects.size(); j ++ )
	{
		if(filterRects[j].useFlag == true)
			continue;

		int dx = abs(filterRects[j].rc.x -offsetX);
		int dy = abs(filterRects[j].rc.y -offsetY);
		int dw = abs(filterRects[j].rc.width - rc.width);
		int dh = abs(filterRects[j].rc.height - rc.height);

		int dist = dx+dy+dw+dh;
		if( dist < min_dist)
		{
			index = j;
			min_dist = dist;
		}
	}

	int maxDiff = rc.width / 4 + rc.height / 4;
	if(min_dist < maxDiff)
		return index;
	else
		return -1;
}


bool findValidRects(Mat& src)
{
	int maxW = src.size().width;
	int maxH = src.size().height;

	int dx = 10, dy = 10;
	for(int row = 0; row < 9; row++)
	{
		for(int col = 0; col < 9; col++)
		{
			if((row == 0) && (col == 0))
			{
				int index = findTopLeftRect(src);
				if(index>0)
				{
					updateValidRect(col, row, filterRects[index].rc, false);
					filterRects[index].useFlag = true;
				}
				else
					return false;
			}
			else if(row == 0)
			{
				Rect leftRc = validRects[col -1].rc; //Align to left
				int offsetX = leftRc.x + leftRc.width + dx;
				int offsetY = leftRc.y;
				int index = findNearest(leftRc, offsetX, offsetY);
				if(index>0) //Find rect
				{
					updateValidRect(col, row, filterRects[index].rc, false);
					filterRects[index].useFlag = true;
				}
				else if( (offsetX + leftRc.width) > maxW) //over the scale
				{
					return false;
				}
				else
				{	
					leftRc.x = offsetX;
					leftRc.y = offsetY;
					updateValidRect(col, row, leftRc, true);
				}
			}
			else if(col == 0)
			{
				Rect upRc = validRects[9*(row-1) + col ].rc; //Align to top
				int offsetX = upRc.x;
				int offsetY = upRc.y + upRc.height + dy;
				int index = findNearest(upRc, offsetX, offsetY);
				if(index>0)
				{
					updateValidRect(col, row, filterRects[index].rc, false);
					filterRects[index].useFlag = true;
				}
				else if( (offsetY + upRc.height) > maxH) 
				{
					return false;
				}
				else
				{
					upRc.x = offsetX;
					upRc.y = offsetY;
					updateValidRect(col, row, upRc, true);
				}
			}
			else
			{
				Rect leftRc = validRects[9*row + col -1].rc; //Align to left and top
				Rect upRc = validRects[9*(row-1) + col ].rc;
				int offsetX = upRc.x;
				int offsetY = leftRc.y;
				int index = findNearest(leftRc, offsetX, offsetY); //Widht and height can use left or top
				if(index>0)
				{
					updateValidRect(col, row, filterRects[index].rc, false);
					filterRects[index].useFlag = true;
				}
				else if( (offsetX + leftRc.width) > maxW) 
				{
					return false;
				}
				else if( (offsetY + upRc.height) > maxH) 
				{
					return false;
				}
				else
				{					
					leftRc.x = offsetX;
					leftRc.y = offsetY;
					leftRc.width =  upRc.width;
					updateValidRect(col, row, leftRc, true);
				}
			}
		}
	}
}

int getVirtualRectsCount()
{
	int count = 0;
	for(int i = 0; i < 81; i++)
	{
		if(validRects[i].virtualFlag == true)
			count++;
	}
	return count;
}

bool modifyValidRect()
{
	for(int row = 0; row < 9; row++)
	{
		for(int col = 0; col < 9; col++)
		{
			if(validRects[9*row+col].virtualFlag == false)
				continue;

			ValidRect leftRc, upRc, rightRc, bottomRc;
			bool left, up, right, bottom;
			left = false;
			up = false;
			right = false;
			bottom = false;

			if(row != 0)
			{
				 upRc = validRects[9*(row-1) + col ];
				 up = true;
			}
			if(row != 8)
			{
				bottomRc = validRects[9*(row + 1) + col];
				bottom = true;
			}
			if(col != 0)
			{
				leftRc = validRects[9*row + col -1];
				left = true;
			}
			if(col != 8)
			{
				rightRc = validRects[9*row + col + 1]; 
				right = true;
			}

			if(left && (leftRc.virtualFlag == false))
			{
				validRects[9*row+col].rc.y = leftRc.rc.y;
				validRects[9*row+col].rc.height = leftRc.rc.height;
			}
			else if(right && (rightRc.virtualFlag == false))
			{
				validRects[9*row+col].rc.y = rightRc.rc.y;
				validRects[9*row+col].rc.height = rightRc.rc.height;
			}

			if(up && (upRc.virtualFlag == false))
			{
				validRects[9*row+col].rc.x = upRc.rc.x;
				validRects[9*row+col].rc.width = upRc.rc.width;
			}
			else if(bottom && (bottomRc.virtualFlag == false))
			{
				validRects[9*row+col].rc.x = bottomRc.rc.x;
				validRects[9*row+col].rc.width = bottomRc.rc.width;
			}
		}
	}
	return true;
}

//int findNearest(Mat& src, int offsetX, int offsetY)
//{
//	int max_dist = src.size().width * src.size().width + src.size().height * src.size().height;
//	int index = -1;
//
//	for(unsigned int j = 0; j < stardRect.size(); j ++ )
//	{
//		int dx = stardRect[j].x -offsetX;
//		int dy = stardRect[j].y -offsetY;
//
//		int dist = dx*dx + dy*dy;
//		if( dist < max_dist)
//		{
//			index = j;
//			max_dist = dist;
//		}
//	}
//
//	if(max_dist < 1800)
//		return index;
//	else
//		return -1;
//}
//

void GetROI(Mat& src, Mat& dst)
{
	int left, right, top, bottom;
	left = src.cols;
	right = 0;
	top = src.rows;
	bottom = 0;

	for(int i=0; i<src.rows; i++)
	{
		for(int j=0; j<src.cols; j++)
		{
			if(src.at<uchar>(i, j) > 0)
			{
				if(j<left) left = j;
				if(j>right) right = j;
				if(i<top) top = i;
				if(i>bottom) bottom = i;
			}
		}
	}

	cv::Point center;
	center.x = (left + right) / 2;
	center.y = (top + bottom) / 2;

	int width = right - left;
	int height = bottom - top;
	int len = (width < height) ? height : width;

	dst = Mat::zeros(len, len, CV_8UC1);

	cv::Rect dstRect((len - width)/2, (len - height)/2, width, height);
	cv::Rect srcRect(left, top, width, height);
	Mat dstROI = dst(dstRect);
	Mat srcROI = src(srcRect);
	srcROI.copyTo(dstROI);
}


//bool findValidRects(Mat& src, double stardWidth, double stardHeight)
//{
//	int minX, minY, maxX = 0, maxY = 0;
//	minX = src.cols;
//	minY = src.rows;
//
//	for(unsigned int j = 0; j < stardRect.size(); j ++ )
//	{
//		if(stardRect[j].x < minX) minX = stardRect[j].x;
//		if(stardRect[j].x > maxX) maxX = stardRect[j].x;
//		if(stardRect[j].y < minY) minY = stardRect[j].y;
//		if(stardRect[j].y > maxY) maxY = stardRect[j].y;
//	}
//
//	int cols = (maxX - minX) / stardWidth + 0.5;
//	int rows = (maxY - minY) / stardHeight + 0.5;
//
//
//	//Not enough valid rect
//	if( (cols != 9) || (rows != 9) )
//		return false;
//
//	int realWidth = (maxX - minX) / 8;
//	int realHight = (maxY - minY) / 8;
//
//	for(int row = 0; row<rows; row++)
//	{
//		for(int col = 0; col<cols; col++)
//		{
//			int offsetX = col * realWidth + minX;
//			int offsetY = row * realHight + minY;
//
//			int index = findNearest(src, offsetX, offsetY);
//			int vIndex = col + row *9;
//
//			if(index >= 0)
//			{
//				validRects[vIndex].rc = stardRect[index];
//				validRects[vIndex].flag = true;
//			}
//			else
//			{
//				validRects[vIndex].rc = Rect( offsetX, offsetY, stardWidth, stardHeight);
//				validRects[vIndex].flag = false;
//			}
//			validRects[vIndex].col = col;
//			validRects[vIndex].row = row;
//
//			//Only care about 90% rect(ignore edge noise)
//			//validRects[vIndex].rc.x = validRects[vIndex].rc.x + validRects[vIndex].rc.width / 10;
//			//validRects[vIndex].rc.y = validRects[vIndex].rc.y + validRects[vIndex].rc.height / 10;
//			//validRects[vIndex].rc.width = validRects[vIndex].rc.width / 10 * 9;
//			//validRects[vIndex].rc.height = validRects[vIndex].rc.height / 10 * 9;
//
//		}
//	}
//
//	return true;
//}
//
//bool findStardardRect(Mat& dst)
//{
//	int allCount = rects.size();
//	 
//	if(allCount < 81)
//		return false;
//
//    // read 4 sequence elements at a time (all vertices of a square)
//	int step = 20;
//	//int step = 8;
//	int start1 = 0, start2 = step / 2;
//
//	Mat widths = Mat::zeros(2, dst.cols / step + 1, CV_16UC1);
//	Mat heights = Mat::zeros(2, dst.rows / step + 1, CV_16UC1);
//
//	int col1, col2, row1, row2;
//
//    for(unsigned int i = 0; i < rects.size(); i ++ )
//    {
//		col1 = rects[i].width / step;
//		col2 = (rects[i].width - start2) / step;
//
//		row1 = rects[i].height / step;
//		row2 = (rects[i].height - start2) / step;
//
//		widths.at<short>(0, col1) =  widths.at<short>(0, col1) + 1;
//		widths.at<short>(1, col2) =  widths.at<short>(1, col2) + 1;
//		heights.at<short>(0, row1) =  heights.at<short>(0, row1) + 1;
//		heights.at<short>(1, row2) =  heights.at<short>(1, row2) + 1;
//
//		//cv::rectangle(dst, rects[i], CV_RGB(0,255,0));
//	}
//
//	double maxCol, maxRow;
//
//	Point pt1, pt2;
//
//	cv::minMaxLoc(widths, 0, &maxCol, 0, &pt1, Mat());
//	cv::minMaxLoc(heights, 0, &maxRow, 0, &pt2, Mat());
//
//	double stardWidth = pt1.x * step + step / 2 + pt1.y * start2;
//	double stardHeight = pt2.x * step + step / 2 + pt2.y * start2;
//
//    for(unsigned int i = 0; i < rects.size(); i ++ )
//    {
//		if( (abs(rects[i].width - stardWidth) < step / 2) && (abs(rects[i].height - stardHeight) < step / 2 ) ) //Sharp check
//		{
//			bool flag = true;
//			for(unsigned int j = 0; j < stardRect.size(); j ++ )
//			{
//				if( (abs(rects[i].x - stardRect[j].x) < (stardWidth - step / 2)) && (abs(rects[i].y - stardRect[j].y) < (stardHeight - step / 2)) )
//				{
//					flag = false;
//					break;
//				}
//			}
//
//			if(flag)
//				stardRect.push_back(rects[i]);
//		}
//	}
//
//	return findValidRects(dst, stardWidth, stardHeight);
//}
//
void drawAllRect(Mat& dst)
{
	//for(unsigned int i = 0; i < stardRect.size(); i ++ )
 //   {
	//	cv::rectangle(dst, stardRect[i], CV_RGB(0,255,0));
	//}
	for(unsigned int i = 0; i < allRects.size(); i ++ )
    {
			cv::rectangle(dst, allRects[i], CV_RGB(0,255,128), 5);
	}
}

void drawFilterRect(Mat& dst)
{
	//for(unsigned int i = 0; i < stardRect.size(); i ++ )
 //   {
	//	cv::rectangle(dst, stardRect[i], CV_RGB(0,255,0));
	//}
	for(unsigned int i = 0; i < filterRects.size(); i ++ )
    {
			cv::rectangle(dst, filterRects[i].rc, CV_RGB(255,255,0), 5);
	}
}

void drawValidRect(Mat& dst)
{
	//for(unsigned int i = 0; i < stardRect.size(); i ++ )
 //   {
	//	cv::rectangle(dst, stardRect[i], CV_RGB(0,255,0));
	//}
	for(unsigned int i = 0; i < 81; i ++ )
    {
		if(validRects[i].virtualFlag)
			cv::rectangle(dst, validRects[i].rc, CV_RGB(0,255,0), 5);
		else
			cv::rectangle(dst, validRects[i].rc, CV_RGB(0,0,255), 5);
	}
}

void newFindSquares4(Mat& src)
{
    storage = cvCreateMemStorage(0);
	IplImage* img = &IplImage(src);
	findSquares4(img, storage);
}

int getMaxContour(vector<vector<Point>> &contours, double& maxArea)
{
	int index = -1;
	maxArea = -1;

	for(unsigned int i=0; i<contours.size(); i++)
	{
		double area = cv::contourArea(Mat(contours[i]));
		if(area>maxArea)
		{
			index = i;
			maxArea = area;
		}
	}

	return index;
}


void removeBlank(Mat& src)
{
	Mat gray, tempGray;
	cvtColor(src, gray, CV_RGB2GRAY);

	cv::pyrDown(gray, tempGray); //, 7 );
	cv::pyrUp(tempGray, gray ); //, 7 );

	//cv::GaussianBlur(gray, tempGray, cvSize(3,3), 3);
	//cv::GaussianBlur(tempGray, gray, cvSize(3,3), 3);
	//medianBlur(gray, tempGray, 1);   
	//medianBlur(tempGray, gray, 1);  

	//dilate(dst, temp, Mat(), Point(-1,-1), 3);   
	//dilate(temp, dst, Mat(), Point(-1,-1), 3);   

	//erode(gray, tempGray, Mat(), Point(-1,-1), 1);   
	//dilate(tempGray, gray, Mat(), Point(-1,-1), 1);   

	//cv::subtract(dst, back, dst, Mat());
	//equalizeHist(dst, dst);

#if(0)
		Mat gray2;
		resize(gray, gray2, cvSize(gray.size().width / 2, gray.size().height / 2));
		imshow("gray", gray2);
		if(waitKey(0) == 27)
			return;
#endif

	for(unsigned int i = 0; i < 81; i ++ )
    {
		validRects[i].blankFlag = true;

		//Mat dst, back, temp;
		Mat dst = gray(validRects[i].rc);

		float proportion;
		Scalar contrast;
		validRects[i].contrast = getBkInfo(dst, proportion, contrast);

		if(validRects[i].contrast<2) // || proportion>0.95)
		{
			//stringstream ss;
			cout << "No number at " << i << ", contrast " << validRects[i].contrast << endl;
			continue;
		}

		//Make inner ellipse to reduce edge noise
		Point2f centerPt;
		centerPt.x = dst.size().width / 2;
		centerPt.y = dst.size().height / 2;
		cv::RotatedRect rrc(centerPt, cvSize(dst.size().width, dst.size().height),  0);

		Mat mask = Mat::zeros(dst.size(), CV_8UC1);
		cv::ellipse(mask, rrc, Scalar(255,255,255), -1);

		Mat temp, temp2;
		temp2 = Mat::zeros(dst.size(), CV_8UC1);


		double thresh = threshold(dst, temp, 0, 0, CV_THRESH_OTSU); //求出分割阈值
		threshold(dst, temp, thresh, 255, CV_THRESH_BINARY_INV);   //分割

		temp.copyTo(temp2, mask & 1);
		temp2.copyTo(temp);

		//if(i==60)
		//{
		//	imshow("BIN", temp);
		//		if(waitKey(0) == 27)
		//			break;
		//}

		//cv::adaptiveThreshold(dst, temp, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY_INV, 3, 5);
		vector<vector<Point>> contours;   
		vector<Vec4i> hierarchy;   

#if(0)
			imshow("ellipse", temp);
			if(waitKey(0) == 27)
				break;
#endif

		//Mat mColor = Mat::zeros(temp.size(), CV_8UC3);
		findContours(temp, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE); //CV
		//findContours(temp, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE); //CV
		Scalar color(0, 255, 255 );   

		double maxArea = -1;
		int maxContourIndex = getMaxContour(contours, maxArea);
		validRects[i].areaProportion = maxArea / validRects[i].rc.area();

		if(maxContourIndex < 0)
		{
			cout << "No number at " << i << ", no contour find " << endl;
			continue;
		}
		else if(validRects[i].areaProportion < MIN_CONTOUR_PROPORTION) //(validRects[i].rc.area() / 50) )
		{
			cout << "No number at " << i << ", contour area too small, proportion is " << validRects[i].areaProportion << endl;
			continue;
		}
		else
		{
			validRects[i].blankFlag = false;
			validRects[i].contourRc = boundingRect(Mat(contours[maxContourIndex]));
			//if(validRects[i].contrast > 20)
			//{
			//	imshow("ellipse", temp);
			//	if(waitKey(0) == 27)
			//		break;
			//}
		}
	}
}


void Predict(Mat& src)
{
	char result[81];
	memset(result, '-', sizeof(result));

	//CvRTrees forest;
	//forest.load( "D:\\SourceCode\\2011\\CV\\CvChessMove\\SvmClassifierForNum\\new_rtrees_10000.xml" );
	int featureLen = 64;

	CvSVM svm;
	svm.load( "D:\\SourceCode\\2011\\CV\\CvChessMove\\SvmClassifierForNum\\SVM_DATA_60000.xml" );

	RtTrainData rtd;
	memset(rtd.data, 0, sizeof(rtd.data));

	Mat m = Mat::zeros(1, featureLen, CV_32FC1);

	Mat gray, tempGray;
	cvtColor(src, gray, CV_RGB2GRAY);

#if(0)
		imshow("gray", gray);
		if(waitKey(0) == 27)
			return;
#endif

	for(unsigned int i = 0; i < 81; i ++ )
    {
		Mat dst = gray(validRects[i].rc);
		if(validRects[i].blankFlag == false)
		{
			Mat target, temp;

			double thresh = threshold(dst, target, 0, 0, CV_THRESH_OTSU); //求出分割阈值
			threshold(dst, target, thresh, 255, CV_THRESH_BINARY_INV);   //分割
			GetROI(target(validRects[i].contourRc), temp);

#if(0)
			imshow("NUM", temp);
			if(waitKey(0) == 27)
				break;
#endif

			Mat temp2 = Mat::zeros(8, 8, CV_8UC1);
			resize(temp, temp2, temp2.size());
		
			for(int i = 0; i<8; i++)
			{
				for(int j = 0; j<8; j++)
				{
					rtd.data[j + i*8] = temp2.at<uchar>(i, j);
					//注意m.data为UCHAR*类型，因此m的值不能使用下面的表达式
					//m.data[j + i*8] = temp3.at<uchar>(i, j);
				}
			}

			memcpy(m.data, rtd.data, featureLen * sizeof(float));

			////GetRectFeatures(rtd);
			//memcpy(m.data, rtd.data, featureLen * sizeof(float));
			//result[i] = forest.predict(m); 
			normalize(m, m);
			float svm_ret = svm.predict(m); 
			char svm_char = (char)svm_ret;

			cout << svm_char << " at " << i << ", contour proportion " << validRects[i].areaProportion << ", contrast " << validRects[i].contrast << endl;
			result[i] = svm.predict(m); 
		}
	}

	cout << endl << "SUDOKU RESULT: "<< endl;
	for(int i=0; i<9; i++)
	{
		for(int j=0; j<9; j++)
			cout << result[9*i+j] << '\t';
		cout << endl;
	}
}

int main(int argc, char** argv)
{
	string fileName = "../../res/sudoku/svg_sudoku7.jpg";
	Mat src = imread(fileName, 1);
	Mat dst; // = Mat::zeros(cvSize(src.size().width / 2, src.size().height / 2), CV_8UC3);
	//Mat org = imread(fileName, 1);

	//flip(src, src, -1);

	//Mat src = Mat::zeros(768, 1024, CV_8UC3);
	//resize(org, src, src.size());

	newFindSquares4(src);
	filterDupRect();

	//drawAllRect(src);
	//drawFilterRect(src);
	//resize(src, dst, cvSize(src.size().width / 2, src.size().height / 2));
	//imshow("all", src);
	//waitKey(0);

	if(findValidRects(src))
	{
		int vc = getVirtualRectsCount();
		cout << "Virtual rects count is " << vc << endl;

		if(vc < 30)
		{
			modifyValidRect();
			removeBlank(src);
			Predict(src);
		}
		else
			cout << "Too much virtual rects to continue predict" << endl;
	}
	else
		cout << "Find rects error" << endl;

	drawValidRect(src);

	resize(src, dst, cvSize(src.size().width / 2, src.size().height / 2));
	imshow("IMG", dst);
	waitKey(0);

	return 0;
}