
平台论坛博客文库

› 论坛 › 程序设计 › C/C++ › GUI编程 › 圆和椭圆求交算法

圆和椭圆求交算法 [复制链接]

醉卧水云间

富甲一方

论坛徽章:: 39

电梯直达

1楼 [收藏(0)] [报告]

发表于 2006-11-20 17:43 |只看该作者 |倒序浏览

写这个函数似乎不那么容易，需要的可以参考。

UINT CGeoFunc::GetCircleIntersectEllipse(CGeoCircle &Circle, CGeoEllipse &Elli, C3DPoint pArray[])
{
C3DVector OA(Elli.GetCenter(),Elli.GetCornerA());
C3DVector OB(Elli.GetCenter(),Elli.GetCornerB());
C3DVector OO(Elli.GetCenter(),Circle.GetCenter()); //椭圆心到圆心的向量
double a=OA.mod();
double b=OB.mod();
double r=Circle.GetRadius();
//检查椭圆是否是圆？如果是圆就简化处理
if(_IsEqual(a,b))
{
CGeoCircle EC;
EC.SetCenter(Elli.GetCenter());
EC.SetNormal(OA^OB);
EC.SetRadius(a);
return GetCircleIntersectCircle(Circle,EC,pArray);
}
C3DVector EV=OA^OB;
C3DVector EVXCV=Circle.GetNormal()^EV;
double modEVXCV=EVXCV.mod()/EV.mod()/Circle.GetNormal().mod();
if(_IsEqual(modEVXCV,0.0))
{//平行或共面
CGeoPlane CP(Circle);
if(_IsPointOnPlane(CP,Elli.GetCenter()))
{//共面，最难解的一个
double C2C=PointToPointDistance(Circle.GetCenter(),Elli.GetCenter());
if(C2C > Circle.GetRadius() + MAX(a,b))
{
return 0;
}
else
{//可能有交的情况,也是最复杂的情况
if(_IsEqual(Elli.GetCenter(),Circle.GetCenter()))
{
//特殊情况特别处理，这种情况下交点不需要求4次方程，也不能用4次方程来求
//1、圆心重合
if(_IsEqual( Circle.GetRadius() , a))
{
pArray[0]=Elli.GetCornerA();
pArray[1]=GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerA());
return 2;
}
else if(_IsEqual( Circle.GetRadius() , b))
{
pArray[0]=Elli.GetCornerB();
pArray[1]=GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerB());
return 2;
}
else if(Circle.GetRadius() < MIN(a,b))
{
return 0;
}
else if(Circle.GetRadius() > MAX(a,b))
{
return 0;
}
else
{
//*/
//2006/02/17重新计算得到x*x=a*a*(b*b-r*r)/(b*b-a*a);
C3DPoint Intersect[4];
double r=Circle.GetRadius();
double x=sqrt(a*a*(b*b-r*r)/(b*b-a*a));
double y=sqrt(r*r-x*x);
Intersect[0]=C3DPoint(x,y,0);
Intersect[1]=C3DPoint(-x,y,0);
Intersect[2]=C3DPoint(x,-y,0);
Intersect[3]=C3DPoint(-x,-y,0);
CUserCoordinate ucs;
ucs.SetSystem(Circle.GetCenter(),OA,OB);
CCSTranslater trans;
trans.SetCoordinate(ucs);
for(int i=0;i<4;i++)
{
pArray[i]=Circle.GetCenter()+
GetFormatVector(OA,Intersect[i].GetX()) +
GetFormatVector(OB,Intersect[i].GetY()) ;
C3DPoint compare=trans.ReverseTransPoint(Intersect[i]);
}
return 4;
}
}
//2、圆心在椭圆OA轴上
CGeoLine LOA(Elli.GetCenter(),Elli.GetCornerA());
LOA.SetType(GEO_LINE_INF);
if(_IsPointOnLine(LOA,Circle.GetCenter()))
{
double dist=OO.mod();
if(_IsEqual(dist,a+r))
{
//选哪个端点呢？假设是原A点
if((OA*OO)<0)
{
pArray[0]=Elli.GetCornerA();
}
else
{
pArray[0]=GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerA());
}
return 1;
}
else if(dist > a+r)
{
return 0;
}
else
{
//以圆心为原点，OA为x轴，OB为y轴建立坐标，则交点连接原点和x轴的夹角满足
//(b*b-a*a)*r*r*t^2-2*ex*r*b*b*t+ex*ex*b*b+a*a*r*r-a*a*b*b=0;
//其中t=cos(alhpa);
//计算前保证椭圆心的横坐标大于0，处于圆的右边
CGeoEllipse ElliCopy(Elli);
if((OO*OA) < 0)
ElliCopy.SetCornerA(GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerA()));
if((OO*OB) < 0)
ElliCopy.SetCornerB(GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerB()));
C3DVector OA(ElliCopy.GetCenter(),ElliCopy.GetCornerA());
C3DVector OB(ElliCopy.GetCenter(),ElliCopy.GetCornerB());
double A=(b*b-a*a)*r*r;
double B=-2*dist*r*b*b;
double C=dist*dist*b*b+a*a*r*r-a*a*b*b;
double delta=B*B-4*A*C;
if(delta < 0) return 0; //无解
double cosalpha1=(-B+sqrt(delta))/(2*A);
double cosalpha2=(-B-sqrt(delta))/(2*A);
bool value1=true;
bool value2=true;
if(cosalpha1 < -1 || cosalpha1 > 1)
value1=false;
if(cosalpha2 < -1 || cosalpha2 > 1)
value2=false;
double sinalpha1;
if(value1) sinalpha1=sqrt(1-cosalpha1*cosalpha1);
double sinalpha2;
if(value2) sinalpha2=sqrt(1-cosalpha2*cosalpha2);
C3DPoint Intersect[4]; //最多4个交
UINT valid=0;
if(value1)
{
Intersect[valid].SetX(r*cosalpha1);
Intersect[valid].SetY(r*sinalpha1);
valid++;
Intersect[valid].SetX(r*cosalpha1);
Intersect[valid].SetY(-r*sinalpha1);
valid++;
}
if(value2)
{
Intersect[valid].SetX(r*cosalpha2);
Intersect[valid].SetY(r*sinalpha2);
valid++;
Intersect[valid].SetX(r*cosalpha2);
Intersect[valid].SetY(-r*sinalpha2);
valid++;
}
//过滤相同的点
valid=FillterOffSamePoint(Intersect,valid);
//转换到世界坐标系
for(UINT i=0;i<valid;i++)
{
Intersect[i]=Circle.GetCenter()+
GetFormatVector(OA,Intersect[i].GetX()) +
GetFormatVector(OB,Intersect[i].GetY()) ;
}
//在世界坐标系检查
UINT onElli=0;
for(i=0;i<valid;i++)
{
if(_IsPointOnEllipse(Elli,Intersect[i]))
{
Intersect[onElli++]=Intersect[i];
}
}
//返回结果
for(i=0;i<onElli;i++)
{
pArray[i]=Intersect[i];
}
return onElli;
}
}
//3、圆心在椭圆OB轴上
CGeoLine LOB(Elli.GetCenter(),Elli.GetCornerB());
//LOB.SetType(GEO_LINE_INF);
if(_IsPointOnLine(LOB,Circle.GetCenter()))
{
double dist=OO.mod();
if(_IsEqual(dist,b+r))
{
//选哪个端点呢？假设是原A点
if((OB*OO)<0)
{
pArray[0]=Elli.GetCornerB();
}
else
{
pArray[0]=GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerB());
}
return 1;
}
else if(dist > b+r)
{
return 0;
}
else
{
//以圆心为原点，OB为x轴，OA为y轴建立坐标，则交点连接原点和x轴的夹角满足
//(a*a-b*b)*r*r*t^2-2*ex*r*a*a*t+ex*ex*a*a+b*b*r*r-a*a*b*b=0;
//其中t=cos(alhpa);
CGeoEllipse ElliCopy(Elli);
if((OO*OA) < 0)
ElliCopy.SetCornerA(GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerA()));
if((OO*OB) < 0)
ElliCopy.SetCornerB(GetMirrorPointByPoint(Elli.GetCenter(),Elli.GetCornerB()));
C3DVector OA(ElliCopy.GetCenter(),ElliCopy.GetCornerA());
C3DVector OB(ElliCopy.GetCenter(),ElliCopy.GetCornerB());
double A=(a*a-b*b)*r*r;
double B=-2*dist*r*a*a;
double C=dist*dist*a*a+b*b*r*r-a*a*b*b;
double delta=B*B-4*A*C;
if(delta < 0) return 0; //无解
double cosalpha1=(-B+sqrt(delta))/(2*A);
double cosalpha2=(-B-sqrt(delta))/(2*A);
bool value1=true;
bool value2=true;
if(cosalpha1 < -1 || cosalpha1 > 1)
value1=false;
if(cosalpha2 < -1 || cosalpha2 > 1)
value2=false;
double sinalpha1=sqrt(1-cosalpha1*cosalpha1);
double sinalpha2=sqrt(1-cosalpha2*cosalpha2);
C3DPoint Intersect[4]; //最多4个交
UINT valid=0;
if(value1)
{
Intersect[valid].SetX(r*cosalpha1);
Intersect[valid].SetY(r*sinalpha1);
valid++;
Intersect[valid].SetX(r*cosalpha1);
Intersect[valid].SetY(-r*sinalpha1);
valid++;
}
if(value2)
{
Intersect[valid].SetX(r*cosalpha2);
Intersect[valid].SetY(r*sinalpha2);
valid++;
Intersect[valid].SetX(r*cosalpha2);
Intersect[valid].SetY(-r*sinalpha2);
valid++;
}
//过滤相同的点
valid=FillterOffSamePoint(Intersect,valid);
//转换到世界坐标系
for(UINT i=0;i<valid;i++)
{
Intersect[i]=Circle.GetCenter()+
GetFormatVector(OB,Intersect[i].GetX()) +
GetFormatVector(OA,Intersect[i].GetY()) ;
}
//在世界坐标系检查
UINT onElli=0;
for(i=0;i<valid;i++)
{
if(_IsPointOnEllipse(Elli,Intersect[i]))
{
Intersect[onElli++]=Intersect[i];
}
}
//返回结果
for(i=0;i<onElli;i++)
{
pArray[i]=Intersect[i];
}
return onElli;
}
}
//一般情况cx!=0,cy!=0
//建立以椭圆中心为原点,A轴方向为X方向,B轴方向为Y方向的坐标
//有方程组...
// x*x/a/a+y*y/b/b=1; -----------------------1
// (x-cx)*(x-cx)+(y-cy)*(y-cy)=r*r; -----------------2
// 由1得 y=b*sqrt(1-(x*x)/(a*a)); -------------------3
// 由2得 y-cy=sqrt(r*r-(x-cx)*(x-cx)); ---------------4
// 3代入4 得 b*sqrt(1-(x*x)/(a*a))-cy=sqrt(r*r-(x-cx)*(x-cx)); ---------5
// 5两边平方得
// (b*b/a/a-1)*x^2+2*cx*x+r*r-cx*cx-cy*cy-b*b=-2*b*cy*sqrt(1-x*x/a/a); -----6
// 设
// K=b*b/a/a-1;
// L=2*cx;
// M=r*r-cx*cx-cy*cy-b*b;
// N=-2*b*cy;
// 则6可以写成
// K*x^2+L*x+M=N*sqrt(1-x*x/a/a); -----------7
// 7两边平方
// K*K*x^4 + 2*K*L*x^3 + (2*K*M+L*L+N*N/a/a)*x^2 + 2*L*M*x + M*M-N*N=0; -------8
// 设
// A=K*K;
// B=2*K*L;
// C=2*K*M+L*L+N*N/a/a;
// D=2*L*M;
// E=M*M-N*N;
// 8式可写成
// A*x^4+B*x^3+C*x^2+D*x+E=0; ------------9
// 直接解方程9就得到x
double r=Circle.GetRadius();
double cx,cy;
//圆中心的坐标
//利用OO向OA,OB投影求坐标
cx=OO*OA/OA.mod();
cy=OO*OB/OB.mod();
double K=b*b/a/a-1;
double L=2*cx;
double M=r*r-cx*cx-cy*cy-b*b;
double N=-2*b*cy;
double A=K*K;
double B=2*K*L;
double C=2*K*M+L*L+N*N/a/a;
double D=2*L*M;
double E=M*M-N*N;
double root[4]; //所有可能的根
int rootcount=CGeoFuncMathLib::SoluteFunction(A,B,C,D,E,root);
if(rootcount==0) return 0; //没有根
//检验所有可能的解共8个
C3DPoint Intersect[4];
int pc=0;
for(int i=0;i<rootcount;i++)
{
//在第一个方程求出y;
double y1=+b*sqrt(1-root[i]*root[i]/a/a);
double y2=-b*sqrt(1-root[i]*root[i]/a/a);
//代进第2方程检查误差
double c1=(root[i]-cx)*(root[i]-cx)+(y1-cy)*(y1-cy)-r*r;
double c2=(root[i]-cx)*(root[i]-cx)+(y2-cy)*(y2-cy)-r*r;
if(IsEqual(c1,0.0))
{
Intersect[pc].x=root[i];
Intersect[pc].y=y1;
Intersect[pc].z=0;
pc++;
}
else if(IsEqual(c2,0.0))
{
Intersect[pc].x=root[i];
Intersect[pc].y=y2;
Intersect[pc].z=0;
pc++;
}
else continue;
}
//现在这些点都是以椭圆圆心为原点的局部坐标，
//下面转换这些坐标到世界坐标系返回
//这个转换和采用正规坐标系转换相同
OA.Normalize();OB.Normalize();
for(i=0;i<pc;i++)
{
pArray[i]=
Elli.GetCenter() +
OA*Intersect[i].x +
OB*Intersect[i].y ;
}
return pc;
}
}
else
{//平行不共面
return 0;
}
}
else
{//不平行
CGeoPlane PC(Circle);
C3DVector Normal=OA^OB;
CGeoPlane PE(Normal,Elli.GetCenter());
CGeoLine L;
UINT nret=GetPlaneIntersectPlane(PC,PE,L);
//ASSERT(nret==1);
//C3DPoint P1,P2;
C3DPoint kArray[2];
nret=GetLineIntersectCircle(Circle,L,kArray);
UINT Counter=0;
for(UINT i=0;i<nret;i++)
{
if(_IsPointOnEllipse(Elli,kArray[i]))
{
pArray[Counter++]=kArray[i];
}
}
return Counter;
}
return 0;
}

复制代码

文库|博客

返回列表

Chinaunix › 论坛 › 程序设计 › C/C++ › GUI编程 › 圆和椭圆求交算法

圆和椭圆求交算法 [复制链接]

浏览过的版块