8.4 Fitting Lines, Rectangles and Squares in the Plane

Fitting a line to a set of points in such a way that the sum of squares of the distances of the given points to the line is minimized, is known to be related to the computation of the main axes of an inertia tensor. Often this fact is used to fit a line and a plane to given points in the 3d space by solving an eigenvalue problem for a $3\times3$ matrix.

In this section we will develop a different algorithm, based on the singular value decomposition, that will allow us to fit lines, rectangles and squares to measured points in the plane and that will also be useful for some problems in 3d space.

Let us first consider the problem of fitting a straight line to a set of given points $P_1$, $P_2$, ..., $P_m$ in the plane. We denote their coordinates with  $(x_{P_1},y_{P_1})$, $(x_{P_2},y_{P_2})$, ..., $(x_{P_m},y_{P_m})$. Sometimes it is useful to define the vectors of all $x$- and $y$-coordinates. We will use ${\bf x}_{P}$ for the vector $(x_{P_1}, x_{P_2}, \ldots, x_{P_m})$ and similarly  ${\bf y}_{P}$ for the $y$ coordinates.

The problem we want to solve is not linear regression. Linear regression means to fit the linear model $y = ax+b$ to the given points, i.e. to determine the two parameters $a$ and $b$ such that the sum of squares of the residual is minimized:

$$\sum_{i=1}^m r_i^2 = \min,$$   where$$\quad r_i = y_{P_i} - a x_{P_i} -b.$$

This simple linear least squares problem is solved in MATLAB by the following statements (assuming that x = ${\bf x}_{P}$ and y = ${\bf y}_{P}$):

   p = [ x ones(size(x))]\y;
   a = p(1); b = p(2);

In the case of the linear regression the sum of squares of the differences of the $y$ coordinates of the points ${P_i}$ to the fitted linear function is minimized. What we would like to minimize now, however, is the sum of squares of the distances of the points from the fitted straight line.

In the plane we can represent a straight line uniquely by the equations

$$c + n_1 x + n_2 y = 0,\quad n_1^2 + n_2^2 = 1.$$ (8.8)

The unit vector $(n_1 , n_2 )$ is orthogonal to the line. A point is on the line if its coordinates $(x,y)$ satisfy the first equation. On the other hand if $P= (x_P,y_P)$ is some point not on the line and we compute

$$r = c + n_1 x_P + n_2 y_P$$

then $\vert r\vert$ is its distance from the line. Therefore if we want to determine the line for which the sum of squares of the distances to given points is minimal, we have to solve the constrained least squares problem

$$\vert\vert{ \bf r} \vert\vert = \sum_{i=1}^m r_i^2 = \min$$

subject to

$$\left( \begin{array}{ccc} 1 & x_{P_1} & y_{P_1} \\ 1 & x_{P_2} & ... ... \right) = \left( \begin{array}{c} r_1\\ r_2\\ \vdots\\ r_m \end{array} \right) \mbox{and $n_1^2 + n_2^2 = 1$.}$$ (8.9)

Let A be the matrix of the linear system (8.9), x denote the vector of unknowns  $(c, n_1, n_2)^T$ and r the right hand side. Since orthogonal transformations  ${\bf y}={\bf Q}^T\ {\bf r}$ leave the norm invariant ( $\vert\vert{\bf y}\vert\vert _2 = \vert\vert{\bf r}\vert\vert _2$ for an orthogonal matrix $\bf Q$), we can proceed as follows to solve problem (8.9).

First we compute the $QR$ decomposition of A and reduce our problem to solving a small system:

$${\bf A} = {\bf Q R} \quad \Longrightarrow \quad {\bf Q}^T {\bf A ... ... \left( \begin{array}{c} c \\ n_1\\ n_2 \end{array} \right) = {\bf Q}^T {\bf r}$$ (8.10)

Since the nonlinear constraint only involves two unknowns we now have to solve

$$\left( \begin{array}{cc} r_{22} & r_{23} \\ 0 & r_{33} \\ \end{ar... ...\end{array} \right) \approx \left( \begin{array}{c} 0 \\ 0 \end{array} \right), \quad n_1^2 + n_2^2 = 1.$$ (8.11)

Problem (8.11) is of the form $\vert\vert{\bf B x} \vert\vert _2 = \min$, subject to $\vert\vert{\bf x}\vert\vert _2=1$. The value of the minimum is the smallest singular value of B, and the solution is given by the corresponding singular vector. Thus we can determine $n_1$ and $n_2$ by a singular value decomposition of a 2-by-2 matrix. Inserting the values into the first component of (8.10) and setting it to zero, we then can compute $c$. As a slight generalization we denote the dimension of the normal vector $\bf n$ by dim. Then, the MATLAB function to solve problem (8.9) is given by the folowing Algorithm cslq.

    function [c,n] = clsq(A,dim);
    % solves the constrained least squares Problem
    % A (c n)' ~ 0 subject to norm(n,2)=1
    % length(n) = dim
    % [c,n] = clsq(A,dim)
    [m,p] = size(A);
    if p < dim+1, error ('not enough unknowns'); end;
    if m < dim, error ('not enough equations'); end;
    m = min (m, p);
    R = triu (qr (A));
    [U,S,V] = svd(R(p-dim+1:m,p-dim+1:p));
    n = V(:,dim);
    c = -R(1:p-dim,1:p-dim)\R(1:p-dim,p-dim+1:p)*n;

Let us test the function clsq with the following main program:

   % mainline.m
   Px = [1:10]'
   Py = [ 0.2 1.0 2.6 3.6 4.9 5.3 6.5 7.8 8.0 9.0]'
   A = [ones(size(Px)) Px Py]
   [c, n] = clsq(A,2)

The line computed by the program mainline has the equation $0.4162 -0.7057 x + 0.7086y = 0$. We would now like to plot the points and the fitted line. For this we need the function plotline,

   function plotline(x,y,s,c,n,t)
   % plots the set of points (x,y) using the symbol s
   % and plots the straight line c+n1*x+n2*y=0 using
   % the line type defined by t
   plot(x,y,s)
   xrange = [min(x) max(x)];
   yrange = [min(y) max(y)];
   if n(1)==0, % c+n2*y=0  => y = -c/n(2)
      x1=xrange(1); y1 = -c/n(2);
      x2=xrange(2); y2 = y1
   elseif n(2) == 0, % c+n1*x=0  => x = -c/n(1)
      y1=yrange(1); x1 = -c/n(1);
      y2=yrange(2); x2 = x1;
   elseif xrange(2)-xrange(1)> yrange(2)-yrange(1),
      x1=xrange(1); y1 = -(c+n(1)*x1)/n(2);
      x2=xrange(2); y2 = -(c+n(1)*x2)/n(2);
   else
      y1=yrange(1); x1 = -(c+n(2)*y1)/n(1);
      y2=yrange(2); x2 = -(c+n(2)*y2)/n(1);
   end
   plot([x1, x2], [y1,y2],t)

The picture is generated by adding the commands

   clf; hold on;
   axis([-1, 11 -1, 11])
   plotline(Px,Py,'o',c,n,'-')
   hold off;