8.4.3 Fitting a Rectangle

Fitting a rectangle requires four sets of points:

$\displaystyle P_i, i=1, \ldots, p, \quad Q_j, j= 1, \ldots, q, \quad
R_k, k= 1, \ldots, r, \quad S_l, l= 1, \ldots, s.

Since the sides of the rectangle are parallel and orthogonal we can proceed very similarly as before. The four sides will have the equations

a: & c_1 + n_1 x + n_2 y & =& 0 \\
b: &...
...4 - n_2 x + n_1 y & =& 0 \\
& n_1^2 + n_2^2 &=&1.

Inserting the sets of points we get the following constrained least squares problem:

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

subject to

$\displaystyle \left( \begin{array}{cccccc} 1 &0 &0 & 0 & x_{P_1} & y_{P_1} \\ \...
... = \left( \begin{array}{c} r_1\\ r_2\\ \vdots\\ r_{p+q+r+s} \end{array} \right)$   $\displaystyle \mbox{and $n_1^2 + n_2^2 = 1$.}$ (8.14)

Instead of explicitly giving the coordinates of the four sets of points, we will now enter the points with the mouse using the ginput function in MATLAB. [X,Y] = ginput(N) gets $ N$ points from the current axes and returns the $ x$- and $ y$-coordinates in the vectors $ X$ and $ Y$ of length $ N$. The points have to be entered clock- or counter clock wise in the same order as the sides of the rectangle: the next side is always orthogonal to the previous.

   % rectangle.m
   clf; hold on;
   axis([0 10 0 10])
   p=100; q=100; r=100; s=100;
   disp('enter points P_i belonging to side A')
   disp('by clicking the mouse in the graphical window.')
   disp('Finish the input by pressing the Return key')
   [Px,Py] = ginput(p); plot(Px,Py,'o')
   disp('enter points Q_i for side B ')
   [Qx,Qy] = ginput(q); plot(Qx,Qy,'x')
   disp('enter points R_i for side C ')
   [Rx,Ry] = ginput(r); plot(Rx,Ry,'*')
   disp('enter points S_i for side D ')
   [Sx,Sy] = ginput(s); plot(Sx,Sy,'+')
   zp = zeros(size(Px)); op =  ones(size(Px));
   zq = zeros(size(Qx)); oq =  ones(size(Qx));
   zr = zeros(size(Rx)); or =  ones(size(Rx));
   zs = zeros(size(Sx)); os =  ones(size(Sx));
   A = [ op zp zp zp Px  Py
         zq oq zq zq Qy -Qx
         zr zr or zr Rx  Ry
         zs zs zs os Sy -Sx]
   [c, n] = clsq(A,2)
   % compute the 4 corners of the rectangle
   B  = [n  [-n(2) n(1)]']
   X = -B* [c([1 3 3 1])'; c([2 2 4 4])']
   X = [X X(:,1)]
   plot(X(1,:), X(2,:))
   % compute the individual lines, if possible
   if all([sum(op)>1 sum(oq)>1 sum(or)>1 sum(os)>1]),
     [c1, n1] = clsq([op Px Py],2)
     [c2, n2] = clsq([oq Qx Qy],2)
     [c3, n3] = clsq([or Rx Ry],2)
     [c4, n4] = clsq([os Sx Sy],2)
     % and their intersection points
     aaa = -[n1(1) n1(2); n2(1) n2(2)]\[c1; c2];
     bbb = -[n2(1) n2(2); n3(1) n3(2)]\[c2; c3];
     ccc = -[n3(1) n3(2); n4(1) n4(2)]\[c3; c4];
     ddd = -[n4(1) n4(2); n1(1) n1(2)]\[c4; c1];
     plot([aaa(1) bbb(1) ccc(1) ddd(1) aaa(1)], ...
          [aaa(2) bbb(2) ccc(2) ddd(2) aaa(2)],':')
   hold off;

Abbildung 8.4: Fitting a Rectangle.
Image gvlsq4
The Program rectangle not only computes the rectangle but also fits the individual lines to the set of points. The result is shown as Figure 8.4. Some comments may be needed to understand some statements. To find the coordinates of a corner of the rectangle, we compute the point of intersection of the two lines of the corresponding sides. We have to solve the linear system

$\displaystyle \begin{array}{rcl}
n_1 x + n_2 y &=& -c_1\\
-n_2 x + n_1 y &=& -...
...\quad C =
n_1 & n_2 \\
-n_2 & n_1

Since $ C$ is orthogonal, we can simply multiply the right hand side by $ B = C^T$ to obtain the solution. By arranging the equations for the 4 corners so that the matrix $ C$ is always the system matrix, we can compute the coordinates of all 4 corners simultaneously with the compact statement
     X = -B* [c([1 3 3 1])'; c([2 2 4 4])'].

Peter Arbenz 2008-09-24