I have used least squares methods in the past, but they are axis dependent and give only approximizations. (Which means the result improves if you use the resulting line as the reference coordinate system for the next iteration).
Not sure if PCA is the same thing i did, but guess so.
So what i came up with is using average weighted angles, utilizing periodical nature of angles to get average lines or crosses. (I use degree 2 for lines to find curvature directions on meshes, and degree 4 to build crossfields from that.)
It calculates exact solution in one step and should work well for you too:
std::vector<vec2> points;
points.push_back(vec2(-3,1));
points.push_back(vec2(-1,2));
points.push_back(vec2(1,-2));
points.push_back(vec2(3,-1));
points.push_back(vec2(5,-2));
points.push_back(vec2(7,0));
vec2 com(0,0);
for (int i=0; i<points.size(); i++)
{
com += points[i];
}
com /= float(points.size());
//for (int i=0; i<points.size(); i++) // some code i've used to proof axis independence
//{
// quat r; r.FromAxisAndAngle (vec(0,0,1), -20.0f/180.0f*PI);
// points[i] = r.Rotate(points[i] - com) + com;
//}
float sinA = 0;
float cosA = 0;
float const degree = 2; // use 1 for vectors, 2 for lines, 4 for crosses etc...
for (int i=0; i<points.size(); i++)
{
vec2 d = points[i] - com;
float weight = d.Dot(d); // can use anything here, but i use squared distance as the weight, which has some backing from physics laws; you could multiply this also with point mass for axample.
float angle = atan2(d[0], d[1]);
angle *= degree;
sinA += sin(angle) * weight;
cosA += cos(angle) * weight;
}
float averageAngle = atan2 (sinA, cosA) / degree;
vec2 averageLine (sin(averageAngle), cos(averageAngle));
// visualize result
ImGui::Text("averageAngle: %f", averageAngle/PI*180.0f);
RenderVector (com, averageLine, 1,1,1);
RenderVector (com, -averageLine, 1,1,1);
for (int i=0; i<points.size(); i++) RenderPoint (points[i], 1,1,1);
I guess this is a common method but i don't know how people call it.