## 8.4.3 Fitting a Rectangle

Fitting a rectangle requires four sets of points:

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

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

subject to

 (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 points from the current axes and returns the - and -coordinates in the vectors  and  of length . 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])
axis('equal')
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)],':')
end
hold off;


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

Since is orthogonal, we can simply multiply the right hand side by to obtain the solution. By arranging the equations for the 4 corners so that the matrix 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