--- a/NEWS	Wed Apr 22 08:41:50 2015 -0700
+++ b/NEWS	Wed Apr 22 20:54:14 2015 +0100
@@ -179,7 +179,7 @@
 
  ** The following functions now support N-dimensional arrays:
 
-      fliplr   flipud   rot90
+      fliplr   flipud   rot90   rectint
 
  ** The new warning ID "Octave:data-file-in-path" replaces the three
     previous separate warning IDs "Octave:fopen-file-in-path",

--- a/scripts/geometry/rectint.m	Wed Apr 22 08:41:50 2015 -0700
+++ b/scripts/geometry/rectint.m	Wed Apr 22 20:54:14 2015 +0100
@@ -1,4 +1,4 @@
-## Copyright (C) 2008-2015 Bill Denney
+## Copyright (C) 2015 Carnë Draug
 ##
 ## This file is part of Octave.
 ##
@@ -18,100 +18,64 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b})
+## Compute area or volume of intersection of rectangles or ND boxes.
 ##
 ## Compute the area of intersection of rectangles in @var{a} and
-## rectangles in @var{b}.  Rectangles are defined as [x y width height]
-## where x and y are the minimum values of the two orthogonal
-## dimensions.
+## rectangles in @var{b}.  N dimensional boxes are supported in which
+## case the volume, or hypervolume is computed according to the number
+## of dimensions.
 ##
-## If @var{a} or @var{b} are matrices, then the output, @var{area}, is a
-## matrix where the i-th row corresponds to the i-th row of a and the j-th
-## column corresponds to the j-th row of b.
+## 2 dimensional rectangles are defined as @code{[xpos ypos width height]}
+## where xpos and ypos are the position of the bottom left corner.
+## Higher dimensions are supported where the coordinates for the minimum
+## value of each dimension follow the length of the box in that dimension,
+## e.g., @code{[xpos ypos zpos kpos @dots{} width height depth k_length @dots{}]}.
+##
+## Each row of @var{a} and @var{b} define a rectangle, and if both define
+## multiple rectangles, then the output, @var{area}, is a matrix where
+## the i-th row corresponds to the i-th row of a and the j-th column
+## corresponds to the j-th row of b.
 ##
 ## @seealso{polyarea}
 ## @end deftypefn
 
-## Author: Bill Denney <bill@denney.ws>
+## Author: 2015 Carnë Draug <carandraug@octave.org>
 
-function area = rectint (a, b)
+function dists = rectint (a, b)
 
   if (nargin != 2)
     print_usage ();
-  elseif (ndims (a) != 2 || ndims (b) != 2)
-    error ("rectint: expecting arguments to be 2-d arrays");
-  elseif (columns (a) != 4)
-    error ("rectint: A must have 4 columns");
-  elseif (columns (b) != 4)
-    error ("rectint: B must have 4 columns");
-  elseif (any ([a(:,3:4);b(:,3:4)](:) < 0))
-    error ("rectint: all widths and heights must be > 0");
-  endif
-
-  ## This runs faster if the number of rows of a is greater than the
-  ## number of rows of b.  Swap them and transpose to make it run
-  ## faster.
-  swapinputs = false ();
-  if (rows (a) > rows (b))
-    tmp = a;
-    a = b;
-    b = tmp;
-    swapinputs = true ();
+  elseif (columns (a) != columns (b))
+    error ("rectint: A and B must have same number of columns");
+  elseif (ndims (a) > 2)
+    error ("rectint: A and B must be 2-d arrays");
+  elseif (mod (columns (a), 2))
+    error ("rectint: number of columns of A and B must be a multiple of two");
   endif
 
-  area = zeros (rows (a), rows (b));
-  r1 = [a(:,1:2) a(:,1:2)+a(:,3:4)];
-  r2 = [b(:,1:2) b(:,1:2)+b(:,3:4)];
-  for i = 1:columns (area)
-    ## Find the location of each point relative to the other points.
-    r1x1small = r1(:,1) < r2(i,1);
-    r1x1large = r1(:,1) > r2(i,3);
-    r1x1mid = (r1(:,1) >= r2(i,1)) & (r1(:,1) <= r2(i,3));
-    r1x2small = r1(:,3) < r2(i,1);
-    r1x2large = r1(:,3) > r2(i,3);
-    r1x2mid = (r1(:,3) >= r2(i,1)) & (r1(:,3) <= r2(i,3));
+  nd = columns (a) / 2;
+  na = rows (a);
+  nb = rows (b);
 
-    r1y1small = r1(:,2) < r2(i,2);
-    r1y1large = r1(:,2) > r2(i,4);
-    r1y1mid = (r1(:,2) >= r2(i,2)) & (r1(:,2) <= r2(i,4));
-    r1y2small = r1(:,4) < r2(i,2);
-    r1y2large = r1(:,4) > r2(i,4);
-    r1y2mid = (r1(:,4) >= r2(i,2)) & (r1(:,4) <= r2(i,4));
+  a_start = a(:,1:nd);
+  b_start = b(:,1:nd);
+
+  a_end = a_start + a(:,nd+1:end);
+  b_end = b_start + b(:,nd+1:end);
 
-    ## determine the width of the rectangle
-    ## r1 completely encloses r2
-    area(r1x1small & r1x2large,i) = r2(i,3) - r2(i,1);
-    ## the range goes from r2x min to r1x max
-    mask = r1x1small & r1x2mid;
-    area(mask,i) = r1(mask,3) - r2(i,1);
-    ## the range goes from r1x min to r2x max
-    mask = r1x1mid & r1x2large;
-    area(mask,i) = r2(i,3) - r1(mask,1);
-    ## the range goes from r1x min to r1x max
-    mask = r1x1mid & r1x2mid;
-    area(mask,i) = r1(mask,3) - r1(mask,1);
+  a_start = reshape (a_start, [na 1 nd]);
+  b_start = reshape (b_start, [1 nb nd]);
 
-    ## determine the height of the rectangle
-    ## r1 completely encloses r2
-    area(r1y1small & r1y2large,i) .*= r2(i,4) - r2(i,2);
-    ## the range goes from r2y min to r1y max
-    mask = r1y1small & r1y2mid;
-    area(mask,i) .*= r1(mask,4) - r2(i,2);
-    ## the range goes from r1y min to r2y max
-    mask = r1y1mid & r1y2large;
-    area(mask,i) .*= r2(i,4) - r1(mask,2);
-    ## the range goes from r1x min to r1x max
-    mask = r1y1mid & r1y2mid;
-    area(mask,i) .*= r1(mask,4) - r1(mask,2);
+  a_end   = reshape (a_end,   [na 1 nd]);
+  b_end   = reshape (b_end,   [1 nb nd]);
 
-  endfor
-
-  if (swapinputs)
-    area = area';
-  endif
+  ## We get a 3d matrix where each dimension is in the 3rd dimension
+  dists = bsxfun (@min , a_end, b_end) - bsxfun (@max, a_start, b_start);
+  dists(dists < 0) = 0;
+  dists = prod (dists, 3);
 
 endfunction
 
-
 ## Exactly overlapping
 %!assert (rectint ([0 0 1 1], [0 0 1 1]), 1)
 ## rect2 completely enclosed by rect1
@@ -130,3 +94,23 @@
 %!assert (rectint ([0 0 1 1;0.5 0.5 1 1;-1 -1 2 2], [1 1 2 2]), [0;0.25;0])
 %!assert (rectint ([1 1 2 2], [0 0 1 1;0.5 0.5 1 1;-1 -1 2 2]), [0 0.25 0])
 
+## bug #44904
+%!assert (rectint ([0 0 5 5], [6 6 5 5]), 0)
+%!assert (rectint ([0 0 5 5], [0 6 5 5]), 0)
+%!assert (rectint ([0 0 5 5], [6 0 5 5]), 0)
+%!assert (rectint ([0 0 0 5 5 5], [0 0 6 5 5 5]), 0)
+
+## Test volumes
+%!shared r1, r2, r3, r4, r5
+%! r1 = [  5   3 0  7   5 2];
+%! r2 = [  2   5 0  4   2 2];
+%! r3 = [ 10   7 0 10   3 2];
+%! r4 = [ 10  -5 0  5   7 2];
+%! r5 = [-10   0 0 40  11 2];
+
+%!assert (rectint (r5, r1), 70)
+%!assert (rectint (r5, r4), 20)
+%!assert (rectint (r5, [r1; r2; r3; r4]), [70 16 60 20])
+
+## Test multiple volumes in both A and B
+%!assert (rectint ([r2; r5], [r1; r3; r4]), [4 0 0; 70 60 20])