几何矩求解椭圆

勒让德惯性椭圆求解

1.matlab利用二阶矩求解椭圆长轴、短轴、离心率、长轴与x轴夹角

xbar=stats(k).Centroid(1);%区域的重心坐标
ybar = stats(k).Centroid(2);


 x = list(:,1) - xbar;
 y = -(list(:,2) - ybar); % This is negative for the
 % orientation calculation (measured in the
 % counter-clockwise direction).

 N = length(x);


 % Calculate normalized second central moments for the region. 1/12 is
  % the normalized second central moment of a pixel with unit length.
  uxx = sum(x.^2)/N + 1/12;%X的二阶矩，即方差
  uyy = sum(y.^2)/N + 1/12;%Y的二阶矩，即方差
  uxy = sum(x.*y)/N;%X,Y的混合二阶矩，协方差

其中，a,b为长轴和短轴，lanmuda为特征值

% Calculate major axis length, minor axis length, and eccentricity.
   common = sqrt((uxx - uyy)^2 + 4*uxy^2);
    stats(k).MajorAxisLength = 2*sqrt(2)*sqrt(uxx + uyy + common);
    stats(k).MinorAxisLength = 2*sqrt(2)*sqrt(uxx + uyy - common);
    stats(k).Eccentricity = 2*sqrt((stats(k).MajorAxisLength/2)^2 - ...
        (stats(k).MinorAxisLength/2)^2) / ...
         stats(k).MajorAxisLength;

% Calculate orientation
 if (uyy > uxx)
       num = uyy - uxx + sqrt((uyy - uxx)^2 + 4*uxy^2);
       den = 2*uxy;
 else
      num = 2*uxy;
      den = uxx - uyy + sqrt((uxx - uyy)^2 + 4*uxy^2);
 end
 if (num == 0) && (den == 0)
      stats(k).Orientation = 0;
 else
      stats(k).Orientation = (180/pi) * atan(num/den);
end

C++实现：

 void get_LegendreEllipse(vector<Point>&contour, double&major_axis, double&minor_axis)
 {
	 Point2f c;
	 get_centroid(contour, c);
	 int length = contour.size();
	 //求x和y的二阶矩，即方差
	 double uxx = 0, uyy = 0, uxy = 0;
	 for (int i = 0; i < length; i++)
	 {
		 uxx += (contour[i].x - c.x)*(contour[i].x - c.x)/length;
		 uyy += (c.y - contour[i].y)*(c.y - contour[i].y)/length;
		 uxy += (contour[i].x - c.x)*(c.y - contour[i].y)/length;
	 }
	 uxx = uxx+1/12;
	 uyy = uxy+1/12;
	 uxy = uyy;

	 double common = sqrt((uxx - uyy)*(uxx - uyy) + 4 * uxy*uxy);
	 major_axis = 2 * sqrt(2)*sqrt(uxx + uyy + common);
	 minor_axis = 2 * sqrt(2)*sqrt(uxx + uyy - common);

     //离心率
	 double eccentricity = 2 * sqrt((major_axis*major_axis) / 4 - (minor_axis*minor_axis) / 4) / major_axis;
	 cout << eccentricity << endl;

	 //计算与x轴夹角
	 double orientation = 0;
	 double num = 0;
	 double den = 0;
	 if (uyy > uxx)
	 {
		 num = uyy - uxx + sqrt((uyy - uxx)*(uyy - uxx) + 4 * uxy*uxy);
		 den = 2 * uxy;
	 }
	 else {
		 num = 2 * uxy;
		 den = uxx - uyy + sqrt((uxx - uyy)*(uxx - uyy) + 4 * uxy*uxy);
	 }
	 if (num == 0 && den == 0)
	 {
		 orientation = 0;
	 }
	 else {
		 orientation = (180 / PI)*atan(num / den);
	 }
	 //cout << orientation << endl;
 }