#pragma warning(disable : 4996)
#pragma comment( lib, "cxcore.lib" )
#pragma comment( lib, "highgui.lib" )
#pragma comment( lib, "cv.lib" )
#pragma comment( lib, "ml.lib" )
#pragma comment( lib, "cvaux.lib" )
#pragma comment( lib, "cvcam.lib" )

#include <cxcore.h>
#include <highgui.h>
#include <cv.h>
#include<iostream.h>
#include <math.h>
#include "cv.h"
#include "highgui.h"
#include <stdio.h>
#include <stdlib.h>


IplImage* inpaint_mask = 0;
IplImage* img0 = 0, *img = 0, *inpainted = 0;
CvPoint prev_pt = {-1,-1};
CvPoint points[100];
int number=0;
//a^2+b^2=1
struct Line
{
	double a;
	double b;
	double c;
};
Line LineFitting(CvPoint* p,int n);
void PrintMat(CvMat* A);
void on_mouse( int event, int x, int y, int flags, void* );
double MaxX(CvPoint *p,int n);
double MinX(CvPoint *p,int n);
void PrintLine(IplImage* img,Line l,double xfrom,double xto);

void main()
{

    cvNamedWindow( "image", 1 );
	img= cvCreateImage(cvSize(700,500), IPL_DEPTH_8U, 3);



//    cvZero(img);

    cvShowImage( "image", img );
    cvSetMouseCallback( "image", on_mouse, 0 );

    for(;;)
    {
        int c = cvWaitKey(0);

        if( (char)c == 27 )
        break;

        if( (char)c == 'r' )
        {
            cvZero(img);
			number=0;
            cvShowImage( "image", img );
        }

    }


}
void on_mouse( int event, int x, int y, int flags, void* )
{
	double xmax,xmin;
	Line l;
    if( !img )
    return;
	if( event == CV_EVENT_LBUTTONDOWN )
	{
		cvCircle(img,cvPoint(x,y),1,cvScalarAll(255));
		points[number]=cvPoint(x,y);
		number++;
		cvShowImage( "image", img );
		cout<<x<<endl;
	}
	if( event == CV_EVENT_RBUTTONUP )
	{

		xmax=MaxX(points,number);
		xmin=MinX(points,number);
		cout<<"xmax="<<xmax<<endl;
		cout<<"xmin="<<xmin<<endl;

		l=LineFitting(points,number);
		cout<<"a="<<l.a<<endl;
		cout<<"b="<<l.b<<endl;
		cout<<"c="<<l.c<<endl;

		PrintLine(img,l,xmin,xmax);
//		cvCircle(img,cvPoint(x,y),50,cvScalarAll(255));
		cvShowImage( "image", img );
	}
    
    
}
double distance(Line l,CvPoint p)
{
	double d;
	d=l.a*p.x+l.b*p.y+l.c;
	if(d>0)
		return d/sqrt(l.a*l.a+l.b*l.b);
	else
		return -d/sqrt(l.a*l.a+l.b*l.b);
}
double distance(Line l,CvPoint* p,int n)
{
	double sum=0;
	for(int i=0;i<n;i++)
	{
		sum=sum+distance(l,p[i]);
	}
	return sum;
}
double MaxX(CvPoint *p,int n)
{
	int max;
	max=p[0].x;
	for(int i=1;i<n;i++)
	{
		if(max<p[i].x)
			max=p[i].x;
	}
	return max;
}
double MinX(CvPoint *p,int n)
{
	int min;
	min=p[0].x;
	for(int i=1;i<n;i++)
	{
		if(min>p[i].x)
			min=p[i].x;
	}
	return min;
}
Line LineFitting(CvPoint *p,int n)
{
	CvMat Ma, Mb, Mc;
	Line l1,l2;
	double a[4],b[4],c[2];
	double avgx2=0;
	double avgx=0;
	double avgxy=0;
	double avgy=0;
	double avgy2=0;
	double t;
	
	for(int i=0;i<n;i++)
	{
		avgx=avgx+((p+i)->x);
		avgx2=avgx2+((p+i)->x)*((p+i)->x);
		avgxy=avgxy+((p+i)->x)*((p+i)->y);
		avgy=avgy+((p+i)->y);
		avgy2=avgy2+((p+i)->y)*((p+i)->y);
	}

	avgx=avgx/n;
	avgx2=avgx2/n;
	avgxy=avgxy/n;
	avgy=avgy/n;
	avgy2=avgy2/n;
	
	
    
	a[0]=avgx2-avgx*avgx;a[1]=avgxy-avgx*avgy;
	a[2]=avgxy-avgx*avgy;a[3]=avgy2-avgy*avgy;

    // 從緩存給矩陣賦值
    cvInitMatHeader( &Ma, 2, 2, CV_64FC1, a, CV_AUTOSTEP);
    cvInitMatHeader( &Mb, 2, 2, CV_64FC1, b, CV_AUTOSTEP);
    cvInitMatHeader( &Mc, 2, 1, CV_64FC1, c, CV_AUTOSTEP);
	cout<<"矩陣A："<<endl;

	PrintMat(&Ma);
	
	
	cvEigenVV(&Ma,&Mb,&Mc);
	

	l1.a=cvGetReal2D(&Mb,0,0);
	l1.b=cvGetReal2D(&Mb,0,1);

	cout<<"特征向量矩陣："<<endl;
	PrintMat(&Mb);
	
	t=sqrt(l1.a*l1.a+l1.b*l1.b);
	l1.a=l1.a/t;
	l1.b=l1.a/t;
	l1.c=-l1.a*avgx-l1.b*avgy;


	l2.a=cvGetReal2D(&Mb,1,0);
	l2.b=cvGetReal2D(&Mb,1,1);
	t=sqrt(l2.a*l2.a+l2.b*l2.b);
	l2.a=l2.a/t;
	l2.b=l2.b/t;
	l2.c=-l2.a*avgx-l2.b*avgy;

	cout<<"l1 p "<<distance(l1,p,n)<<endl;
	cout<<"l2 p "<<distance(l2,p,n)<<endl;

	if(distance(l1,p,n)<distance(l2,p,n))
	{
		cout<<"應(yīng)該取l1"<<endl;
		return l1;
	}
	else
	{
		cout<<"應(yīng)該取l2"<<endl;
		return l2;
	}
	
}


void PrintLine(IplImage* img,Line l,double xfrom,double xto)
{
	CvPoint pFrom,pTo;
	double x,y;

	x=xfrom;
	y=(-l.c-l.a*x)/l.b;

	pFrom=cvPoint(x,y);

	x=xto;
	y=(-l.c-l.a*x)/l.b;
	pTo=cvPoint(x,y);

	cvLine(img,pFrom,pTo,cvScalarAll(255),1);	
}
void PrintMat(CvMat* A)
{
    int i,j;
    //printf("\nMatrix = :");
    for(i=0;i<A->rows;i++)
    {
        printf("\n");
        
        switch( CV_MAT_DEPTH(A->type) )
        {
        case CV_32F:
        case CV_64F:
            for(j=0;j<A->cols;j++)
                printf("%9.3f ", (float) cvGetReal2D( A, i, j ));
            break;
        case CV_8U:
        case CV_16U:
            for(j=0;j<A->cols;j++)
                printf("%6d",(int)cvGetReal2D( A, i, j ));
            break;
        default:
            break;
        }
    }
    printf("\n");
}