22.txt

椭圆拟合的算法

void CVCTiSView::OnFitellips()
{
	// TODO: 在此添加命令处理程序代码
	if(HotPtArray.GetSize()>10000)
	{
		AfxMessageBox(L"数据量太大，请清除一些选择点!");
		return;
	}
	if(HotPtArray.GetSize()>2)
	{		
		double X1=0;
		double Y1=0;
		double X1cY1c=0;
		double X2c=0;
        double Y2c=0;
		
		double m_pointsX[10000],m_pointsY[10000];
		     for (int j=0;j<HotPtArray.GetSize();j++)
		      {	
			     m_pointsX[j] = (((CPolygon*)HotPtArray.GetAt(j))->Polygonpoint.x);
			     m_pointsY[j] = ((CPolygon*)HotPtArray.GetAt(j))->Polygonpoint.y;
		      }
		      for (int i=0;i<HotPtArray.GetSize();i++)
		      {
			      X1 = X1 + m_pointsX[i];
			      Y1 = Y1 + m_pointsY[i];			
		      }
		int num = HotPtArray.GetSize();
		double Xc,Yc;          //椭圆的圆心即质心；
		Xc = X1 / num ;
		Yc = Y1 / num ;
		for (int i=0;i<HotPtArray.GetSize();i++)
		{
		 X1cY1c = X1cY1c + (m_pointsX[i]-Xc)*(m_pointsY[i]-Yc);

		 X2c = X2c + (m_pointsX[i]-Xc)*(m_pointsX[i]-Xc);
		 Y2c = Y2c + (m_pointsY[i]-Yc)*(m_pointsY[i]-Yc);
		}
         
		 double longAngle,angleY,angle0;
		 angle0 = 0.5* atan2((2*X1cY1c),( X2c-Y2c));              //范围是正负90度
         double pie = 3.14159;          
		 if (angle0<0)
		 {
			 longAngle= angle0 +0.5* pie;
			 angleY= angle0 +pie;
		 }
		 else
		 {
			 longAngle= angle0+ pie;
			 angleY= angle0 + 1.5*pie;
		 }
		 double longaxis,shortaxis;	
		 distanceX =0.0,distanceY=0.0;
		 for (int i=0;i<HotPtArray.GetSize();i++)
		 {
             distanceX = distanceX + pow(((m_pointsY[i]-Yc)*cos(longAngle)-(m_pointsX[i]-Xc)*sin(longAngle)),2);
			 distanceY = distanceY + pow(((m_pointsY[i]-Yc)*cos(angleY)-(m_pointsX[i]-Xc)*sin(angleY)),2);
		 }
		 longaxis =1.2*sqrt (distanceX/num);
		 shortaxis =1.2*sqrt (distanceY/num);

		
		 			/*CString str;
			str.Format(_T("%d整型用d;%f其他类型都用f"),i,m_nAngle);
			AfxMessageBox(str);*/

		 if (longaxis<shortaxis)
		 {
			double mmm = longaxis;
		 longaxis = shortaxis;
		 shortaxis = mmm;
		// longAngle = angleY+0.5*pie;
		 }
		  longaxis =1.23*longaxis;
		 shortaxis =1.2*shortaxis;

		 EllipsCenter.x=Xc+2;
		 EllipsCenter.y=Yc;
		

		CFitEllips *pEllips=new CFitEllips(EllipsCenter,longaxis,shortaxis,longAngle);
		FitEllipsArray.Add(pEllips);
		
		 Invalidate();
		 m_nFittype=4;
		 FitEllipseErrorAnalyse( longaxis,shortaxis,Xc,Yc);

    }
}