实例介绍
【实例简介】根据图像上的椭圆计算出拟合椭圆的5个参数：长轴、短轴、圆心x坐标、圆心y坐标、倾角
【实例截图】
【核心代码】
// C_EDlg.cpp : implementation file
//
#include "stdafx.h"
#include "C_E.h"
#include "C_EDlg.h"
#include "math.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define pi 3.14159265358979
struct FXY{double x,y;};
int all_gauss( int n, double *a, double *b);
/////////////////////////////////////////////////////////////////////////////
// CC_EDlg dialog
CC_EDlg::CC_EDlg(CWnd* pParent /*=NULL*/)
: CDialog(CC_EDlg::IDD, pParent)
{
//{{AFX_DATA_INIT(CC_EDlg)
m_xc = 0.0f;
m_a = 0.0f;
m_b = 0.0f;
m_yc = 0.0f;
m_bt = 0.0f;
//}}AFX_DATA_INIT
m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}
void CC_EDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
//{{AFX_DATA_MAP(CC_EDlg)
DDX_Control(pDX, IDC_LIST1, m_lst);
DDX_Text(pDX, IDC_EDIT1, m_xc);
DDX_Text(pDX, IDC_EDIT3, m_a);
DDX_Text(pDX, IDC_EDIT4, m_b);
DDX_Text(pDX, IDC_EDIT2, m_yc);
DDX_Text(pDX, IDC_EDIT5, m_bt);
//}}AFX_DATA_MAP
}
BEGIN_MESSAGE_MAP(CC_EDlg, CDialog)
//{{AFX_MSG_MAP(CC_EDlg)
ON_WM_PAINT()
ON_WM_QUERYDRAGICON()
ON_BN_CLICKED(IDC_BUTTON1, OnButton1)
ON_BN_CLICKED(IDC_BUTTON2, OnButton2)
ON_EN_CHANGE(IDC_EDIT1, OnChangeEdit1)
ON_EN_CHANGE(IDC_EDIT2, OnChangeEdit2)
ON_EN_CHANGE(IDC_EDIT3, OnChangeEdit3)
ON_EN_CHANGE(IDC_EDIT4, OnChangeEdit4)
ON_EN_CHANGE(IDC_EDIT5, OnChangeEdit5)
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CC_EDlg message handlers
BOOL CC_EDlg::OnInitDialog()
{
CDialog::OnInitDialog();
SetIcon(m_hIcon, TRUE); // Set big icon
SetIcon(m_hIcon, FALSE); // Set small icon
// TODO: Add extra initialization here
m_xc=m_yc=m_bt=0.;
m_a=m_b=100.;UpdateData(false);
return TRUE; // return TRUE unless you set the focus to a control
}
// If you add a minimize button to your dialog, you will need the code below
// to draw the icon. For MFC applications using the document/view model,
// this is automatically done for you by the framework.
void CC_EDlg::OnPaint()
{
if (IsIconic())
{
CPaintDC dc(this); // device context for painting
SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);
// Center icon in client rectangle
int cxIcon = GetSystemMetrics(SM_CXICON);
int cyIcon = GetSystemMetrics(SM_CYICON);
CRect rect;
GetClientRect(&rect);
int x = (rect.Width() - cxIcon 1) / 2;
int y = (rect.Height() - cyIcon 1) / 2;
// Draw the icon
dc.DrawIcon(x, y, m_hIcon);
}
else
{
CDialog::OnPaint();
}
}
HCURSOR CC_EDlg::OnQueryDragIcon()
{
return (HCURSOR) m_hIcon;
}
void CC_EDlg::OnButton1()
{
// TODO: Add your control notification handler code here
int i,j,k,ret;
double r,xc,yc,x,y,z;
double af;
char st[80];
FXY xy[360];
r=m_a;xc=m_xc,yc=m_yc;
m_lst.ResetContent();
sprintf(st,"r=%9.4f xc=%9.4f yc=%9.4f",r,xc,yc);m_lst.AddString(st);
m_lst.AddString("");
for(k=0;k<360;k )
{
af=k*pi/180.;
xy[k].x=xc r*cos(af);
xy[k].y=yc r*sin(af);
sprintf(st,"k=%3d %9.4f %9.4f",k,xy[k].x,xy[k].y);m_lst.AddString(st);
}
m_lst.AddString("");
//-----------------------
double *XX=(double *)new double[12];
FillMemory(XX,sizeof(double)*12,0);
for(k=0;k<360;k )
{
x=xy[k].x;y=xy[k].y;z=x*x y*y;
XX[0] =x*x;XX[4] =y*y;XX[8] =1;
XX[1] =x*y;XX[2] =x;XX[5] =y;
XX[9] =x*z;XX[10] =y*z;XX[11] =z;
}
XX[3]=XX[1];XX[6]=XX[2];XX[7]=XX[5];
ret=all_gauss(3,XX,&XX[9]);
if(ret==0)return;
//x0=XX[9]/2.;y0=XX[10]/2.; c=XX[11]=>r0*r0=c x0*x0 y0*y0
sprintf(st,"%9.4f %9.4f %9.4f",XX[9],XX[10],XX[11]);m_lst.AddString(st);
xc=XX[9]/2.;yc=XX[10]/2.;
// r=sqrt(XX[11]-x*x-y*y);
r=sqrt(XX[9]*XX[9] XX[10]*XX[10] 4*XX[11])/2.;
sprintf(st,"xc=%9.4f yc=%9.4f r0=%9.4f",xc,yc,r);m_lst.AddString(st);
delete [] XX;
}
int all_gauss( int n, double *a, double *b)
{
int l,k,i,j,is,p,q;//js 用于记忆列交换信息的空间
double d,t;
int *js = (int *) new int[n];
l = 1 ;
for(k=0; k<n-1; k )
{
d = 0.0;
for(i=k; i<n; i )for(j=k; j<n; j )
{
t = fabs(*(a i*n j));
if(t>d){ d = t; js[k] = j; is = i;}
} //选主元
if(d 1.0==1.0){return 0;} //奇异标志
else
{
if(js[k]!=k)
for(i=0; i<n; i )
{ //列交换
p = k;
q = js[k];
t = *(a i*n p);
*(a i*n p)=*(a i*n q);
*(a i*n q)=t;
}
if(is!=k)
{
for(j=k; j<n; j )//行交换
{
t = *(a k*n j);
*(a k*n j)=*(a is*n j);
*(a is*n j) = t;
}
t = b[k];
b[k] = b[is];
b[is] = t;
}
}
d = *(a k*n k);
for(j=k 1; j<n; j )
{
*(a k*n j) /= d ;
}
b[k]/=d;
for(i=k 1; i<n; i )
{ //消元
for(j=k 1; j<n; j )
{
p = j;
*(a i*n p)-= (*(a i*n k)) * (*(a k*n j));
}
b[i] -= (*(a i*n k))*b[k];
}
}
d =*(a (n-1)*n n-1);
if(fabs(d) 1.0==1.0){delete [] js;return 0;}
b[n-1]/=d;
for(i=n-2; i>=0; i--)
{
t = 0.0;
for(j=i 1; j<n; j ) t = *(a i*n j)*b[j];
b[i] -= t;
}
js[n-1] = n-1;
for(k=n-1; k>=0; k--)if(js[k]!=k)
{
t = b[k];
b[k] = b[js[k]];
b[js[k]] = t; //恢复解向量
}
delete [] js;
return 1;
} // all_gauss
void CC_EDlg::OnButton2()
{
// TODO: Add your control notification handler code here
double xc,yc,a,b;
double x,y,r,af,bt;
double A,B,C,D,E,F,G;
int i,j,k,h2,w2,i0,j0;
RECT rct2;
char st[80];
FXY p9,xy[360];
long hw;
xc=m_xc;yc=m_yc;
a=m_a;b=m_b;
bt=m_bt*pi/180.;
m_lst.ResetContent();
sprintf(st,"xc=%9.4f yc=%9.4f bt=%9.4f %9.4f°",xc,yc,bt,bt*180/pi);m_lst.AddString(st);
sprintf(st,"a=%9.4f b=%9.4f",a,b);m_lst.AddString(st);
for(k=0;k<360;k ){af=k*pi/180.;xy[k].x=a*cos(af);xy[k].y=b*sin(af);}
for(k=0;k<360;k )
{
x=xy[k].x;y=xy[k].y;r=sqrt(x*x y*y);
af=asin(y/r);if(x<0.)af=pi-af;
xy[k].x=r*cos(af-bt) xc;xy[k].y=r*sin(af-bt) yc;
}
for(k=0;k<360;k ){sprintf(st,"%3d %9.4f %9.4f",k,xy[k].x,xy[k].y);m_lst.AddString(st);}
//拟合
double U,V,W,X,Y,Z,z,xx,yy,Xy,Q,C2,C4;
int ret;
double *XX=(double *)new double[30];
FillMemory(XX,sizeof(double)*30,0);
for(k=0;k<360;k )
{
xx=xy[k].x*xy[k].x;
Xy=xy[k].x*xy[k].y;
yy=xy[k].y*xy[k].y;
x=xy[k].x;
y=xy[k].y;
z=-xx;
XX[0] =yy*yy;
XX[6] =Xy*Xy;
XX[12] =xx;
XX[18] =yy;
XX[24] =1;
XX[1] =yy*Xy;
XX[2] =yy*x;
XX[3] =yy*y;
XX[4] =yy;
XX[7] =Xy*x;
XX[8] =Xy*y;
XX[9] =Xy;
XX[13] =Xy;
XX[14] =x;
XX[19] =y;
XX[25] =yy*z;
XX[26] =Xy*z;
XX[27] =x*z;
XX[28] =y*z;
XX[29] =z;
}
//对称
XX[5]=XX[1];
XX[10]=XX[2];
XX[11]=XX[7];
XX[15]=XX[3];
XX[16]=XX[8];
XX[17]=XX[13];
XX[20]=XX[4];
XX[21]=XX[9];
XX[22]=XX[14];
XX[23]=XX[19];
ret=all_gauss(5,XX,&XX[25]);
if(ret==0){SetWindowText("all_gauss fail!");return;}
sprintf(st,"%9.4f %9.4f %9.4f %9.4f %9.4f",XX[25],XX[26],XX[27],XX[28],XX[29]);m_lst.AddString(st);
// X=-W*yc-2*xc
// Y=-2*V*yc-W*xc
// 2*Y/W=-4*V/W*yc-2*xc
//2*Y/W-X=(W-4*V/W)*yc
//(W-4*V/W)*yc=2*Y/W-X
//yc=(2*Y/W-X)/(W-4*V/W)
yc=(2*XX[28]/XX[26]-XX[27])/(XX[26]-4*XX[25]/XX[26]);
xc=(-XX[26]*yc-XX[27])/2.;
sprintf(st,"xc=%9.4f yc=%9.4f",xc,yc);m_lst.AddString(st);
Q=(1.-XX[25])/XX[26];
C=Q sqrt(Q*Q 1.);bt=atan(C);
sprintf(st,"bt=%9.4f°",bt*180/pi);m_lst.AddString(st);
U=1./(xc*xc XX[25]*yc*yc XX[26]*xc*yc-XX[29]);
sprintf(st,"U=%15.12f C=%9.4f",U,C);m_lst.AddString(st);
V=XX[25]*U;
W=XX[26]*U;
X=XX[27]*U;
Y=XX[28]*U;
Z=XX[29]*U-1;
// sprintf(st,"%15.12f %15.12f %15.12f",U,V,W);m_lst.AddString(st);
// sprintf(st,"%15.12f %15.12f %15.12f",X,Y,Z);m_lst.AddString(st);
C2=C*C;
C4=C2*C2;
A=sqrt((U-C2*V)/(1.-C4));
B=sqrt((V-C2*U)/(1.-C4));
sprintf(st,"A=%9.6f B=%9.6f",A,B);m_lst.AddString(st);
a=cos(bt)/A;
b=cos(bt)/B;
sprintf(st,"a=%9.4f b=%9.4f",a,b);m_lst.AddString(st);
delete [] XX;
}
void CC_EDlg::OnChangeEdit1() {UpdateData(true);}
void CC_EDlg::OnChangeEdit2() {UpdateData(true);}
void CC_EDlg::OnChangeEdit3() {UpdateData(true);}
void CC_EDlg::OnChangeEdit4() {UpdateData(true);}
void CC_EDlg::OnChangeEdit5() {UpdateData(true);}
标签： 椭圆拟合算法
好例子网口号：伸出你的我的手 — 分享！
小贴士
感谢您为本站写下的评论，您的评论对其它用户来说具有重要的参考价值，所以请认真填写。
- 类似“顶”、“沙发”之类没有营养的文字，对勤劳贡献的楼主来说是令人沮丧的反馈信息。
- 相信您也不想看到一排文字/表情墙，所以请不要反馈意义不大的重复字符，也请尽量不要纯表情的回复。
- 提问之前请再仔细看一遍楼主的说明，或许是您遗漏了。
- 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦，又没人会搭理，于人于己都无利。
关于好例子网
本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享，对搜索内容的合法性不具有预见性、识别性、控制性，仅供学习研究，请务必在下载后24小时内给予删除，不得用于其他任何用途，否则后果自负。基于互联网的特殊性，平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查；无论平台是否已进行审查，用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场，基于网友分享，根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定，若资源存在侵权或相关问题请联系本站客服人员，点此联系我们。关于更多版权及免责申明参见 版权及免责申明
网友评论
我要评论