I copied the best least squares line fit algorithm from Graphics Gems V book, 
but now i'm having problems with vertical lines, although the book claims that 
the algorithm can handle vertical lines.

For example using 0,0 and 0,2 as test points gives a=0 and b=0 and c=0 for the 
line ax + by = c, although it should be something like a=1, b=0, c=0.

Here is my code:

void best_fit(void) 
{
	if(n_points<2)
	{
		a=1; // dummy values
		b=1;
		c=0; 
		return; 
	}

	double sx=0;
	double sy=0;
	double sxx=0;
	double syy=0;
	double sxy=0;

	for(point *p=points.first(); p; p=p->next)
	{
		sx += p->x;
		sy += p->y;
	}
	double xmid = sx/n_points;
	double ymid = sy/n_points;

	for(p=points.first(); p; p=p->next)
	{
		double x1 = p->x-xmid;
		double y1 = p->y-ymid;
		sxx += x1*x1;
		syy += y1*y1;
		sxy += x1*y1;
	}
	double a1 = sxx-syy;
	double b1 = sxy;
	a = 2*b1;
	b = - (a1 + sqrt(a1*a1 + 4*b1*b1));
	c = - a*xmid - b*ymid;

	if(a==0.0 && b==0.0) // try to solve the vertical line problem
	{
		a=1; c=-xmid;
	}
}		

"Markku Poysti" <ukram@puosu.dna.fi_nospam> wrote in message
news:3fba8264$0$6509$1b6aedd2@news.songnet.fi...
> I copied the best least squares line fit algorithm from Graphics Gems V
book,
> but now i'm having problems with vertical lines, although the book claims
that
> the algorithm can handle vertical lines.
<snip>

The derivation has problems.  For example, consider equation (3)
which has x*sin(theta)+y*cos(theta)+p = 0.  Figure 1 shows what
p and theta mean.  Now choose p = 0 and theta = 45 degrees.
The figure indicates a line y = x.  The equation says x + y = 0.
Does this give you warm fuzzies about the rest of the construction?

Let t = ((sxx+syy)-sqrt((sxx-syy)*(sxx-syy)+4*sxy*sxy))/2.
The least squares line contains the average (xmid,ymid) of the
points (using your notation).  The direction is chosen to be
either (-sxy, t-syy) or (t-sxx, -sxy), whichever has larger length.
(The GG5 construction essentially chose the wrong direction
that, as you noticed, was (0,0).)

Grundgy details in
http://www.magic-software.com/Documentation/LeastSquaresFitting.pdf

--
Dave Eberly
eberly@magic-software.com
http://www.magic-software.com
http://www.wild-magic.com