OpenCV:橢圓上點的計算方程

橢圓

        橢圓（Ellipse）是平面內到定點F1、F2的距離之和等於常數（大於|F1F2|）的動點P的軌跡，F1、F2稱為橢圓的兩個焦點。其數學表達式為：                                                                                                                |PF1|+|PF2|=2a（2a>|F1F2|）。         橢圓是圓錐曲線的一種，即圓錐與平面的截線。 橢圓在開普勒行星運行三定律中扮演了重要角色，即恆星是橢圓兩焦點中的一個，是數學科重點研究的一個項目。^[

標准方程：

        橢圓的標准方程共分兩種情況：               當焦點在x軸時，橢圓的標准方程是：x^2/a^2+y^2/b^2=1，(a>b>0)；               當焦點在y軸時，橢圓的標准方程是：y^2/a^2+x^2/b^2=1，(a>b>0)；                                                   其中a^2-c^2=b^2

參數方程：

      橢圓上點的參數方程為：

            y = a *sin(  alp ) 

            x=  a *cos( alp )  (a>b>0)；

 此時的角度alp不是中心點到橢圓上點的角度，而是橢圓的仿射圓上的點到圓心的角度，計算角度應考慮到壓縮。

壓縮方向：

        Height方向拉伸；

        計算變化後的beta；

        計算坐標：

            y = a *sin(  beta ) 

            x=  a *cos( beta )  (a>b>0)；

        Height方向壓縮；

              y = a *sin(  beta ) *（b/a）
              x=  a *cos( beta )                  (a>b>0)；

      計算距離。

橢圓上點的計算方程：

        對於 (a>b>0)；

              對應的圓的方程： R = a;

              圓上的點的坐標： x2 = R * sin(Beta)    y2 = R * cos(beta);

              不變性：  alp = beta                

             

      對應橢圓點的坐標：

              角度：  alp = beta

              角度：  alp = beta

計算橢圓上點的代碼：

    代碼是錯誤的，不能把點壓縮到橢圓上

		//調整橢圓邊緣到標准橢圓；在角度方向上進行拉伸
		//angleOfDip 為橢圓的偏斜角，弧度值！
		//增加邊界檢查
		template <class T1,class T2>
		float AdjustEllipseEdge(
			std::vector<std::pair< T1, T2 > >  &closeEdgeIn,
			std::vector<std::pair< T1, T2 > >  &closeEdgeOut,
			const cv::RotatedRect &ecf,
			const cv::Point2f &rfCentroidS,
			const double angleOfDipSrc,
			const int ww,
			const int hh)
		{
			assert(closeEdgeIn.size() == closeEdgeOut.size() );
			int w = ww -1;
			int h = hh -1;

			const cv::Point2f rfCentroid = ecf.center;
			//cv::Point2f rfCentroid(0,0);
			std::vector< double > angleListS;//為點橢圓角度，用於求取 橢圓點到中心的距離
			angleListS.resize( closeEdgeIn.size() );

			int vOrH = 0;//水平或者豎直？
			vOrH = ecf.size.width > ecf.size.height? 0:1;//若0，則為V;或者為1，水平

			double angleOfDip = 0;
			if (0 == vOrH )
			{//若為水平//width 的傾角
				angleOfDip  = angleOfDipSrc; 
			} 
			else
			{
				angleOfDip  = angleOfDipSrc - PI_1_2; 
			}

			double a = max(ecf.size.height/2.0,ecf.size.width /2.0);//長軸//固定後使用方程
			double b = min(ecf.size.height/2.0,ecf.size.width /2.0);

#ifdef SHOW_TEMP
			cv::Mat canvasSrc = cv::Mat::zeros(200,200,CV_8UC3);
			cv::bitwise_not(canvasSrc,canvasSrc);
			cv::ellipse(canvasSrc,ecf,cv::Scalar(0,0,255),1,8);
#endif

			//在此測試，cos計算的代碼
#ifdef SHOW_TEMP

			cv::RotatedRect ecT = RotatedRect(Point2f(100,100), Size2f(50,100), 30);
			std::vector<std::pair< cv::Point2f, double > >  PointCosTest(0);
			cvWish::polygon::GetElipseEdge(ecT, PointCosTest, (ecT.size.height + ecT.size.height)/5.0 );
			cv::ellipse(canvasSrc, ecT, cv::Scalar(0,0,255), 1, 8);
			
			for ( int i=0; i< PointCosTest.size(); ++i)
			{
				cv::circle( canvasSrc, PointCosTest[i].first, 1, cv::Scalar(255,0,0), 1, 8, 0 );
				double af = cvWish::cosCv(ecT.center,PointCosTest[i].first);//cosCv出現計算問題
				std::cout<< "Cos:" << af<< std::endl;
				std::cout<< "Angle:" << PointCosTest[i].second << std::endl;
				cv::imshow("PointCosTest",canvasSrc);
				cv::waitKey(1);
			}

#endif

			for ( int i=0; i<closeEdgeIn.size(); ++i )
			{
				closeEdgeIn[i].second  = cvWish::cosCv( rfCentroid, closeEdgeIn[i].first );

				angleListS[i]  = closeEdgeIn[i].second;
				angleListS[i] -= angleOfDip;//旋轉
				angleListS[i]  = angleListS[i]> PI_4_2 ? angleListS[i] - PI_4_2:angleListS[i];

				//探測距離
				double disPC    = cvWish::disCv(rfCentroid,closeEdgeIn[i].first);

				double alp =  angleListS[i];
				//alp =  alp *180/M_PI;
				double disShould = 
					sqrt( b*sin(alp ) *b*sin(alp ) + a*cos(alp) *a*cos(alp) );//公式無誤，角度出現問題？
				     //sqrt( b*cos(alp ) *b*cos(alp ) + a*sin(alp) *a*sin(alp) );//公式無誤，角度出現問題？
				//可能問題，方向角度出現往長軸極點的方向進行壓縮，導致生成距離變大。
				
				//double disShould = sqrt( 
				//	ecf.size.width*cos(angleListS[i]) *ecf.size.width*cos(angleListS[i]) /4
				//	+ ecf.size.height*sin(angleListS[i]) *ecf.size.height*sin(angleListS[i])/4 );
				std::cout<< alp << std::endl;
				std::cout<< cos(alp)  << std::endl;
				std::cout<<"disPc:" <<disPC << std::endl;
				std::cout<< "disShould:" << disShould << std::endl;

#ifdef SHOW_TEMP
				//cv::Mat canvasSrc(100,100,CV_8UC3);
				cv::circle(canvasSrc,closeEdgeIn[i].first,1,cv::Scalar(255,0,0),1,8,0);
				cv::imshow("edgeEvolution",canvasSrc);
				cv::waitKey(1);
#endif
				//調整點到橢圓上
				//adjustPoint2Elipse();

				//根據距離 往角度方向上拉伸點//角度其實產生了偏離//偏角使用圖片偏角
				cvWish::PullPoint2Out( closeEdgeIn[i].first, closeEdgeIn[i].second, ( disPC - disShould ) );	

				closeEdgeOut[i].first  = closeEdgeIn[i].first;
				////已確認大於0，此時確認不超邊界
				closeEdgeOut[i].first.x = min(closeEdgeOut[i].first.x,w);
				closeEdgeOut[i].first.y = min(closeEdgeOut[i].first.y,h);

				closeEdgeOut[i].second = closeEdgeIn[i].second;
#ifdef SHOW_TEMP
				cv::circle(canvasSrc,closeEdgeOut[i].first,1,cv::Scalar(0,255,0),1,8,0);
				cv::imshow("edgeEvolution",canvasSrc);
				cv::waitKey(1);
#endif
			}

			return 1.0;
		}

代碼修改：

    使用一個仿射變換

		//調整橢圓邊緣到標准橢圓；在角度方向上進行拉伸
		//angleOfDip 為橢圓的偏斜角，弧度值！
		//增加邊界檢查
		template <class T1,class T2>
		float AdjustEllipseEdge(
			std::vector<std::pair< T1, T2 > >  &closeEdgeIn,
			std::vector<std::pair< T1, T2 > >  &closeEdgeOut,
			const cv::RotatedRect &ecf,
			const cv::Point2f &rfCentroidS,
			const double angleOfDipSrc,
			const int ww,
			const int hh)
		{
			assert(closeEdgeIn.size() == closeEdgeOut.size() );
            int w = ww -1;
            int h = hh -1;

            const cv::Point2f rfCentroid = ecf.center;
            //cv::Point2f rfCentroid(0,0);
            std::vector< double > angleListS;//為點橢圓角度，用於求取 橢圓點到中心的距離
            angleListS.resize( closeEdgeIn.size() );

            int vOrH = 0;//水平或者豎直？
            vOrH = ecf.size.width > ecf.size.height? 0:1;//若0，則為V;或者為1，水平

            double angleOfDip = 0;
            if (0 == vOrH )
            {//若為水平//width 的傾角
                angleOfDip  = angleOfDipSrc; 
            } 
            else
            {
                angleOfDip  = angleOfDipSrc - PI_1_2; 
            }

            //double a = max(ecf.size.height/2.0,ecf.size.width /2.0);//長軸//固定後使用方程
            //double b = min(ecf.size.height/2.0,ecf.size.width /2.0);
            double b = ecf.size.height/2.0//長軸//固定後使用方程
            double a = ecf.size.width /2.0;
            double compressFactor = b /a ;//壓縮或者縮放因子

#ifdef SHOW_TEMP
            cv::Mat canvasSrc = cv::Mat::zeros(200,200,CV_8UC3);
            cv::bitwise_not(canvasSrc,canvasSrc);
            cv::ellipse(canvasSrc,ecf,cv::Scalar(0,0,255),1,8);
#endif

            //在此測試，cos計算的代碼
#ifdef SHOW_TEMP

            cv::RotatedRect ecT = RotatedRect(Point2f(100,100), Size2f(50,100), 30);
            std::vector<std::pair< cv::Point2f, double > >  PointCosTest(0);
            cvWish::polygon::GetElipseEdge(ecT, PointCosTest, (ecT.size.height + ecT.size.height)/5.0 );
            cv::ellipse(canvasSrc, ecT, cv::Scalar(0,0,255), 1, 8);
            
            for ( int i=0; i< PointCosTest.size(); ++i)
            {
                cv::circle( canvasSrc, PointCosTest[i].first, 1, cv::Scalar(255,0,0), 1, 8, 0 );
                double af = cvWish::cosCv(ecT.center,PointCosTest[i].first);//cosCv出現計算問題
                std::cout<< "Cos:" << af<< std::endl;
                std::cout<< "Angle:" << PointCosTest[i].second << std::endl;
                cv::imshow("PointCosTest",canvasSrc);
                cv::waitKey(1);
            }

#endif

            for ( int i=0; i<closeEdgeIn.size(); ++i )
            {
                closeEdgeIn[i].second  = cvWish::cosCv( rfCentroid, closeEdgeIn[i].first );

                //壓縮方向
                angleListS[i]  = closeEdgeIn[i].second;
                angleListS[i] -= angleOfDip;//旋轉
                angleListS[i]  = angleListS[i]> PI_4_2 ? angleListS[i] - PI_4_2:angleListS[i];

                //探測距離
                double disPC    = cvWish::disCv(rfCentroid,closeEdgeIn[i].first);

                //double alp =  angleListS[i];
                //alp =  alp *180/M_PI;
                //double disShould = sqrt( b*sin(alp ) *b*sin(alp ) + a*cos(alp) *a*cos(alp) );//公式無誤，角度出現問題？
                //可能問題，方向角度出現往長軸極點的方向進行壓縮，導致生成距離變大。

                //計算對應仿射圓的角度
                double xDeta =  closeEdgeIn[i].first.x - rfCentroid.x;
                double yDeta =  closeEdgeIn[i].first.y - rfCentroid.y;
                yDeta /= compressFactor;

                //計算角度
                double beta = cvWish::cosCv( rfCentroid, cv::Point2f( rfCentroid.x + xDeta, rfCentroid.y+ yDeta ) );
                double r = a;                                                                
                xDeta  =  r* cos(beta);
                yDeta  =  r* sin(beta);
                yDeta *= compressFactor;

                //直接計算距離
                double disShould = sqrt( xDeta*xDeta + yDeta*yDeta );//公式無誤，角度出現問題？

                std::cout<<"disPc:" <<disPC << std::endl;
                std::cout<< "disShould:" << disShould << std::endl;

#ifdef SHOW_TEMP
                //cv::Mat canvasSrc(100,100,CV_8UC3);
                cv::circle(canvasSrc,closeEdgeIn[i].first,1,cv::Scalar(255,0,0),1,8,0);
                cv::imshow("edgeEvolution",canvasSrc);
                cv::waitKey(1);
#endif
                //調整點到橢圓上
                //adjustPoint2Elipse();

                //根據距離 往角度方向上拉伸點//角度其實產生了偏離//偏角使用圖片偏角
                cvWish::PullPoint2Out( closeEdgeIn[i].first, closeEdgeIn[i].second, ( disPC - disShould ) );    

                closeEdgeOut[i].first  = closeEdgeIn[i].first;
                ////已確認大於0，此時確認不超邊界
                closeEdgeOut[i].first.x = min(closeEdgeOut[i].first.x,w);
                closeEdgeOut[i].first.y = min(closeEdgeOut[i].first.y,h);

                closeEdgeOut[i].second = closeEdgeIn[i].second;
#ifdef SHOW_TEMP
                cv::circle(canvasSrc,closeEdgeOut[i].first,1,cv::Scalar(0,255,0),1,8,0);
                cv::imshow("edgeEvolution",canvasSrc);
                cv::waitKey(1);
#endif
            }

            return 1.0;
        }

從一個橢圓上面獲取特定個數的點的函數：

//參數描述：橢圓；輸出的點集；欲獲取的點的個數
		int polygon::GetElipseEdge(
			const cv::RotatedRect &ecf,  
			std::vector<std::pair< cv::Point2f, double > >  &ellipseEdge,
			const int numPs,
			cv::Rect &roiRestrict,
			bool openEdgeRestrict )
		{
			if ( numPs == 0 )
			{
				return numPs;
			}
			else
			{
				ellipseEdge.resize( numPs );
			}

			//對橢圓進行劃分
			const double angleGap = PI_4_2/numPs;
			const double cx = ecf.center.x;
			const double cy = ecf.center.y;
			const float angleOfDip =  PI_1_2 + ecf.angle*3.1415926 /180.0;//為何偏移了 半個pi
			//const double angleOfDip =0- ecf.angle*3.1415926 /180.0;//

			double w = ecf.size.width /2.0;
			double h = ecf.size.height/2.0;
			for (int i=0 ;i< numPs;++i )
			{
				double as = i*angleGap ;

				double a = as ;
				a += angleOfDip;
				a = a>PI_4_2? a-PI_4_2:a;

				double y = (w) *sin( a );
				double x = (h) *cos( a );

				//旋轉
				float xDeta = x*cos( angleOfDip ) - y*sin( angleOfDip );
				float yDeta = x*sin( angleOfDip ) + y*cos( angleOfDip );

				cv::Point2f p( cx+xDeta, cy+yDeta);
				//ellipseEdge[i] = (std::pair< T1, T2 >)(std::make_pair( p,as ) );
				//ellipseEdge[i] = (std::pair< cv::Point2f, double >)(std::make_pair( p,as ) );//此處代碼只為運行於GCC修改，有問題，模板庫不能使用！！！wishchin!!!
				ellipseEdge[i].first.x = p.x;
				ellipseEdge[i].first.y = p.y;
				ellipseEdge[i].second  =  as;
			}

			
			if (openEdgeRestrict)
			{
				float x,y;
				float xS(roiRestrict.x), yS(roiRestrict.y), xE(roiRestrict.x+roiRestrict.width), yE(roiRestrict.y+roiRestrict.height );
				
				for (int i=0 ;i< numPs;++i )
				{
					x = ellipseEdge[i].first.x;
					y = ellipseEdge[i].first.y;

					x = (std::min)( (std::max)(x,xS),xE );
					y = (std::min)( (std::max)(y,yS),yE );

					//ellipseEdge[i].first = cv::Point2f(x,y);
					ellipseEdge[i].first.x = x;
					ellipseEdge[i].first.y = y;
				}
			} 
			else
			{
			}

			return 1;
		}

 結果顯示：

      原始結果：                                                                                    修改後結果：