Standard Deviational Ellipsoid, Average Direction of 3D Points

Question

I need help to calculate "average direction" from a set of 3D points, located around (0,0,0).

In other words I want the direction vector for biggest axis of the 3d ellipsoid covering the points.

I was able to make it work for 2D points. Based on this information: http://resources.esri.com/help/9.3/arcgisengine/java/gp_toolref/spatial_statistics_tools/how_directional_distribution_colon_standard_deviational_ellipse_spatial_statistics_works.htm

Using Standard Deviational Ellipse

My C++ code ended up like this:

void AvgDirU(Vec2 &dir, C Array<Vec2> &points) // // get average direction from array of points, points should be located around center Vec2(0,0)  , direction is not normalized
{
   float x2=0, y2=0, xy=0;
   for(int i=0; i<points.elms(); i++)
   {
    C Vec2 &p=points[i];
      x2+=Sqr(p.x);
      y2+=Sqr(p.y);
      xy+=p.x*p.y;
   }
   float a=x2-y2,
         c=2*xy,
         b=sqrt(sqr(a)+sqr(c));
   dir.x=a+b; // if we would normalize, then dir.x=cos
   dir.y=c  ; // if we would normalize, then dir.y=sin
}

I've tested this and it works great.

'dir' points in the average direction of an array of 2D points.

Note that instead of calculating tangent of the angle as in the link above, I prefer to specify Vec2 direction vector (cos and sin), this is much better because it avoids infinite values and division by zero.

However how to extend this to 3D Ellipsoid?

I need:

void AvgDirU(Vec3 &dir, C Array<Vec3> &points)
{
..
}

Many thanks!