f

#### Best least-squares line fit

```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;
}
}

```
 0
ukram (3)
11/18/2003 8:34:45 PM
comp.graphics.algorithms 6674 articles. 1 followers.

1 Replies
743 Views

Similar Articles

[PageSpeed] 52

```"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

```
 0
eberly (231)
11/18/2003 9:41:23 PM