octave_packages/m/geometry/rectint.m

   1 ## Copyright (C) 2008-2012 Bill Denney
   2 ##
   3 ## This file is part of Octave.
   4 ##
   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.
   9 ##
  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.
  14 ##
  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/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b})
  21 ##
  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
  25 ## dimensions.
  26 ##
  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.
  30 ##
  31 ## @seealso{polyarea}
  32 ## @end deftypefn
  33
  34 ## Author: Bill Denney <bill@denney.ws>
  35
  36 function area = rectint (a, b)
  37
  38   if (nargin != 2)
  39     print_usage ();
  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");
  48   endif
  49
  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
  52   ## faster.
  53   swapinputs = false ();
  54   if (rows (a) > rows (b))
  55     tmp = a;
  56     a = b;
  57     b = tmp;
  58     swapinputs = true ();
  59   endif
  60
  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));
  72
  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));
  79
  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);
  92
  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);
 105
 106   endfor
 107
 108   if swapinputs
 109     area = area';
 110   endif
 111
 112 endfunction
 113
 114 ## Tests
 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])