1## Copyright (C) 2008, 2009 Bill Denney
3## This file is part of Octave.
5## Octave is free software; you can redistribute it and/or modify it
6## under the terms of the GNU General Public License as published by
7## the Free Software Foundation; either version 3 of the License, or (at
8## your option) any later version.
10## Octave is distributed in the hope that it will be useful, but
11## WITHOUT ANY WARRANTY; without even the implied warranty of
12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13## General Public License for more details.
15## You should have received a copy of the GNU General Public License
16## along with Octave; see the file COPYING. If not, see
17## <http://www.gnu.org/licenses/>.
20## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b})
22## Compute the area of intersection of rectangles in @var{a} and
23## rectangles in @var{b}. Rectangles are defined as [x y width height]
24## where x and y are the minimum values of the two orthogonal
27## If @var{a} or @var{b} are matrices, then the output, @var{area}, is a
28## matrix where the i-th row corresponds to the i-th row of a and the j-th
29## column corresponds to the j-th row of b.
34## Author: Bill Denney <bill@denney.ws>
36function area = rectint (a, b)
40 elseif (ndims (a) != 2 || ndims (b) != 2)
41 error ("rectint: expecting arguments to be 2-d arrays");
42 elseif (columns (a) != 4)
43 error ("rectint: a must have 4 columns");
44 elseif (columns (b) != 4)
45 error ("rectint: b must have 4 columns");
46 elseif any ([a(:,3:4);b(:,3:4)](:) < 0)
47 error ("rectint: all widths and heights must be > 0");
50 ## This runs faster if the number of rows of a is greater than the
51 ## number of rows of b. Swap them and transpose to make it run
53 swapinputs = false ();
54 if (rows (a) > rows (b))
61 area = zeros (rows (a), rows (b));
62 r1 = [a(:,1:2) a(:,1:2)+a(:,3:4)];
63 r2 = [b(:,1:2) b(:,1:2)+b(:,3:4)];
64 for i = 1:columns (area)
65 ## Find the location of each point relative to the other points.
66 r1x1small = r1(:,1) < r2(i,1);
67 r1x1large = r1(:,1) > r2(i,3);
68 r1x1mid = (r1(:,1) >= r2(i,1)) & (r1(:,1) <= r2(i,3));
69 r1x2small = r1(:,3) < r2(i,1);
70 r1x2large = r1(:,3) > r2(i,3);
71 r1x2mid = (r1(:,3) >= r2(i,1)) & (r1(:,3) <= r2(i,3));
73 r1y1small = r1(:,2) < r2(i,2);
74 r1y1large = r1(:,2) > r2(i,4);
75 r1y1mid = (r1(:,2) >= r2(i,2)) & (r1(:,2) <= r2(i,4));
76 r1y2small = r1(:,4) < r2(i,2);
77 r1y2large = r1(:,4) > r2(i,4);
78 r1y2mid = (r1(:,4) >= r2(i,2)) & (r1(:,4) <= r2(i,4));
80 ## determine the width of the rectangle
81 ## r1 completely encloses r2
82 area(r1x1small & r1x2large,i) = r2(i,3) - r2(i,1);
83 ## the range goes from r2x min to r1x max
84 mask = r1x1small & r1x2mid;
85 area(mask,i) = r1(mask,3) - r2(i,1);
86 ## the range goes from r1x min to r2x max
87 mask = r1x1mid & r1x2large;
88 area(mask,i) = r2(i,3) - r1(mask,1);
89 ## the range goes from r1x min to r1x max
90 mask = r1x1mid & r1x2mid;
91 area(mask,i) = r1(mask,3) - r1(mask,1);
93 ## determine the height of the rectangle
94 ## r1 completely encloses r2
95 area(r1y1small & r1y2large,i) .*= r2(i,4) - r2(i,2);
96 ## the range goes from r2y min to r1y max
97 mask = r1y1small & r1y2mid;
98 area(mask,i) .*= r1(mask,4) - r2(i,2);
99 ## the range goes from r1y min to r2y max
100 mask = r1y1mid & r1y2large;
101 area(mask,i) .*= r2(i,4) - r1(mask,2);
102 ## the range goes from r1x min to r1x max
103 mask = r1y1mid & r1y2mid;
104 area(mask,i) .*= r1(mask,4) - r1(mask,2);
115## Exactly overlapping
116%!assert(rectint([0 0 1 1], [0 0 1 1]), 1)
117## rect2 completely enclosed by rect1
118%!assert(rectint([-1 -1 3 3], [0 0 1 1]), 1)
119## rect1 completely enclosed by rect2
120%!assert(rectint([0 0 1 1], [-1 -1 3 3]), 1)
121## rect1 right and top in rect2
122%!assert(rectint([-1 -1 1.5 1.5], [0 0 1 1]), 0.25)
123## rect2 right and top in rect1
124%!assert(rectint([0 0 1 1], [-1 -1 1.5 1.5]), 0.25)
125## no overlap - shared corner
126%!assert(rectint([0 0 1 1], [1 1 2 2]), 0)
127## no overlap - shared edge
128%!assert(rectint([0 0 1 1], [0 1 2 2]), 0)
129## Correct orientation of output
130%!assert(rectint([0 0 1 1;0.5 0.5 1 1;-1 -1 2 2], [1 1 2 2]), [0;0.25;0])
131%!assert(rectint([1 1 2 2], [0 0 1 1;0.5 0.5 1 1;-1 -1 2 2]), [0 0.25 0])
