Find Gaps and Overlaps in Rectangular Prisms

I have rectangular prisms at fixed locations defined below, and I'm trying to see if they fit in a box without any gaps or overlaps. The perfect no-gaps, no-overlaps solution looks like the picture of the attached Jenga tower.

% Fit the following rectangular prisms in a box:
prisms.p1.x = [75 200];     % Min/max of x range for prism p1
prisms.p1.y = [4 100];      % Min/max of y range for prism p1
prisms.p1.z = [1e3 1e4];    % Min/max of z range for prism p1

prisms.p2.x = [20 80];
prisms.p2.y = [0.1 2.2];
prisms.p2.z = [0 1e3];

prisms.p3.x = [80 115];
prisms.p3.y = [0.1 4];
prisms.p3.z = [0 1e3];

prisms.p4.x = [20 180];
prisms.p4.y = [2.2 8];
prisms.p4.z = [0 1e3];

prisms.p5.x = [200 1000];
prisms.p5.y = [0 1];
prisms.p5.z = [0 1e4];

% The box that the prisms must fit in, given the values above, is:
box.x = [20 1000];  % Min/max of x range for the box of prisms
box.y = [0 100];    % Min/max of x range for the box of prisms
box.z = [0 1e4];    % Min/max of x range for the box of prisms

Above is a struct called prisms. It has 5 fields; each field describes one prism (p1,p2,p3,p4,p5) using x,y,z coordinates. I want to identify if the rectangular prisms have any overlaps or gaps between them, and if so, report where those gaps and/or overlaps occur. The ranges are all aligned with the axes, and are always >=0.

The size of the box that the prisms must fit in is defined to be the min and max values over all prisms for a particular axis. If two prisms touch at only one point in their range, I don't count that as an "overlap". I don't want to move the prisms, but instead want to report where two prisms overlap, or where there are gaps in the box.

I've found all the overlapping regions by testing pairs of rectangular prisms for overlap, and it seems to work:

%% Get all combinations of classes (combinations since order does not matter)
nFields = length(fieldnames(prisms));
combos = nchoosek(1:nFields,2);
nCombos = nchoosek(nFields,2);

% Index the struct by index instead of field name
fns = fieldnames(prisms);

%% Check for overlap between each pair of classes
for ic = 1:nCombos
    x11 = prisms.(fns{combos(ic,1)}).x(1);
    x12 = prisms.(fns{combos(ic,1)}).x(2);
    y11 = prisms.(fns{combos(ic,1)}).y(1);
    y12 = prisms.(fns{combos(ic,1)}).y(2);
    z11 = prisms.(fns{combos(ic,1)}).z(1);
    z12 = prisms.(fns{combos(ic,1)}).z(2);

    x21 = prisms.(fns{combos(ic,2)}).x(1);
    x22 = prisms.(fns{combos(ic,2)}).x(2);
    y21 = prisms.(fns{combos(ic,2)}).y(1);
    y22 = prisms.(fns{combos(ic,2)}).y(2);
    z21 = prisms.(fns{combos(ic,2)}).z(1);
    z22 = prisms.(fns{combos(ic,2)}).z(2);

    if( (x12>x21 && y12>y21 && z12>z21) || (x21>x12 && y21>y12 && z21>z12) )
        % There is either overlap, or the edge of classes align
        
        % Check for overlap in x
        if( ((x12>=x21) && (x11<=x22))  ||  ((x11<=x22) && (x21<=x21)) )
            if((x21==x12) || (x11==x22))
                continue;
            else
                fprintf('Overlap in x between %s and %s\n', fns{combos(ic,1)}, fns{combos(ic,2)});
                fprintf('  Limits for %s x: [%g %g]\n', fns{combos(ic,1)}, x11, x12);
                fprintf('  Limits for %s x: [%g %g]\n', fns{combos(ic,2)}, x21, x22);
            end
        end

        % Check for overlap in y
        if( ((y12>=y21) && (y11<=y22))  ||  ((y11<=y22) && (y21<=y21)) )
            if((y21==y12) || (y11==y22))
                continue;
            else
                fprintf('Overlap in y between %s and %s\n', fns{combos(ic,1)}, fns{combos(ic,2)});
                fprintf('  Limits for %s y: [%g %g]\n', fns{combos(ic,1)}, y11, y12);
                fprintf('  Limits for %s y: [%g %g]\n', fns{combos(ic,2)}, y21, y22);
            end
        end

        % Check for overlap in z
        if( ((z12>=z21) && (z11<=z22))  ||  ((z11<=z22) && (z21<=z21)) )
            if((z21==z12) || (z11==z22))
                continue;
            else
                fprintf('Overlap in z between %s and %s\n', fns{combos(ic,1)}, fns{combos(ic,2)});
                fprintf('  Limits for %s z: [%g %g]\n', fns{combos(ic,1)}, z11, z12);
                fprintf('  Limits for %s z: [%g %g]\n', fns{combos(ic,2)}, z21, z22);
            end
        end
    end
end

I just can't get the last part which will identify where the gaps happen. I think I can modify the above code to push this function across the finish line, but I need some help.

[Catalytic]

How are you defining a gap? Is it when one of the prisms has no contact with any of the other prisms?

[Alex]

I define a gap as any location in the box that isn't filled by a prism. If two prisms touch at their edges, like how p1 and p5's x ranges touch, I don't count that as an overlap.

prisms.p1.x = [75 200];
prisms.p5.x = [200 1000];

% No gap or overlap here since p1 ends at 200 and p5 starts at 200. 

[Alex]

However, if two prisms don't touch, then there is a gap. A prism doesn't have to be touched on every axis to form a gap.

prisms.p1.x = [75 200];
prisms.p5.x = [200 1000];
% No gap here since p1 ends at 200 and p5 starts at 200. 

prisms.p1.x = [75 199];
prisms.p5.x = [200 1000];
% There is a gap here since the values between 199 and 200 aren't filled
% in. 

Answer 1

[Matt J]

If you can get R2025a with the PDE Toolbox, you gain access to tools that are very well-suited to this problem. In particular with subtract(),

https://www.mathworks.com/help/pde/ug/fegeometry.subtract.html

and union(),

https://www.mathworks.com/help/pde/ug/fegeometry.union.html

you can subtract the union of cuboids from a solid box, thus exposing the cavities.

[Alex]

Oh wow, that's interesting. I was hoping to port this into Python and going down this route might make that very difficult. I'll try this out regardless, as it looks like a great tool.

[Matt J]

I don't think it should be hard, e.g.,

%MATLAB
function model = buildBricksModel(boxDims, bricks)
% buildBricksModel(boxDims, bricks)
% boxDims = [L, W, H] of the main box
% bricks = struct array with fields: x, y, z, w, h, d
% where (x, y, z) is the corner, and w/h/d are dimensions

model = createpde();
box = multicuboid(boxDims(1), boxDims(2), boxDims(3));

% Convert bricks to geometry objects
brickSolids = [];

for i = 1:length(bricks)
    b = bricks(i);
    % Create and position the brick
    brick = multicuboid(b.w, b.h, b.d);
    brick = translate(brick, [b.x + b.w/2, b.y + b.h/2, b.z + b.d/2]);
    brickSolids = [brickSolids, brick];
end

% Union all bricks
if ~isempty(brickSolids)
    brickUnion = brickSolids(1);
    for i = 2:length(brickSolids)
        brickUnion = union(brickUnion, brickSolids(i));
    end
    gm = subtract(box, brickUnion);
else
    gm = box;
end

% Assign geometry
model.Geometry = gm;

% Optional: generate mesh here
% generateMesh(model);

end

#PYTHON

import matlab.engine

# Start MATLAB engine
eng = matlab.engine.start_matlab()

# Define box dimensions (length, width, height)
box_dims = [100, 100, 100]

# Define list of bricks
brick_data = [
    {"x": 10, "y": 10, "z": 10, "w": 30, "h": 30, "d": 30},
    {"x": 50, "y": 20, "z": 20, "w": 20, "h": 40, "d": 20},
    {"x": 10, "y": 60, "z": 10, "w": 40, "h": 30, "d": 30}
]

# Convert to MATLAB struct array
brick_structs = []
for b in brick_data:
    brick_structs.append(eng.struct(**{
        "x": b["x"], "y": b["y"], "z": b["z"],
        "w": b["w"], "h": b["h"], "d": b["d"]
    }))

brick_matlab_array = matlab.struct([brick_structs])  # wrap as MATLAB struct array

# Call MATLAB function
eng.cd('path/to/your/matlab/scripts')  # Update to correct path
model = eng.buildBricksModel(box_dims, brick_matlab_array, nargout=1)

# Optional: visualize in MATLAB
eng.eval("pdegplot(model, 'FaceAlpha', 0.5, 'FaceLabels', 'on')", nargout=0)

Answer 2

[Matt J]

Does a "gap" not occur whenever the two prisms being compared neither touch, nor overlap? If so, why not as below?

% Fit the following rectangular prisms in a box:
prism(1).x = [75 200];     % Min/max of x range for prism 1
prism(1).y = [4 100];      % Min/max of y range for prism 1
prism(1).z = [1e3 1e4];    % Min/max of z range for prism 1

prism(2).x = [20 80];
prism(2).y = [0.1 2.2];
prism(2).z = [0 1e3];

prism(3).x = [80 115];
prism(3).y = [0.1 4];
prism(3).z = [0 1e3];

prism(4).x = [20 180];
prism(4).y = [2.2 8];
prism(4).z = [0 1e3];

prism(5).x = [200 1000];
prism(5).y = [0 1];
prism(5).z = [0 1e4];

ij=nchoosek(1:5,2);

clear contactType
for k=height(ij):-1:1
    i=ij(k,1);
    j=ij(k,2);
    contactType(k,1)=comparePrisms(prism(i),prism(j));
end

T=table(ij,contactType)

T = 10x2 table
      ij       contactType 
    ______    _____________

    1    2    "unconnected"
    1    3    "touching"   
    1    4    "touching"   
    1    5    "unconnected"
    2    3    "touching"   
    2    4    "touching"   
    2    5    "unconnected"
    3    4    "overlapping"
    3    5    "unconnected"
    4    5    "unconnected"

function contactType=comparePrisms(pa,pb)

   f=@(a,b) any(b(1)<=a & a<=b(2));
   f=@(a,b) f(a,b) | f(b,a);
   Q(1)=f(pa.x,pb.x);
   Q(2)=f(pa.y,pb.y);
   Q(3)=f(pa.z,pb.z);

   f=@(q) numel(unique(q));
   N(1)=f([pa.x,pb.x]);
   N(2)=f([pa.y,pb.y]);
   N(3)=f([pa.z,pb.z]);

   if all(Q) && ~any(N==3)
       contactType="overlapping"; %prisms intersect with positive volume
   elseif all(Q) && any(N==3)
       contactType="touching";    %prisms intersect, but with zero volume
   else
       contactType="unconnected"; %prisms do not intersect
   end

end

[Catalytic]

If a gap is supposed to mean there is an isolated brick, then you can use @Matt J's code to form a graph of the bricks and analyze that -

A=false(numel(prisms));
for k=height(ij):-1:1
    i=ij(k,1);
    j=ij(k,2);
    contactType(k,1)=comparePrisms(prism(i),prism(j));
    
    A(i,j)=contactType(k)=="touching"; %adjacency matrix
end

G=graph(A,'upper');

isolatedBricks=find(degree(G)==0); %Gap?

[Matt J]

If a gap is supposed to mean there is an isolated brick...

Possibly, a gap means that the prisms don't form a single connected tree. That would take something more general than identifying prisms that are isolated. You would need to decompose the grap into connected components:

%test if all bricks are connected
isConnected = max(conncomp(G))==1 & ~any(T.contactType=="overlapping") 

[Alex]

For this problem, a gap occurs whenever the face of any prism doesn't completely touch another prism, or the box the prisms are contained in. If one prism's face touches multiple faces of other prisms or the box, then there are still no gaps. This must be true for all the prisms.

The problem is similar to taking a few bricks of different size and seeing if they fit into a box with no wiggle room, but the bricks are fixed in place, and the box is made around the bricks! Like a Jenga tower (see attached .png for a perfect no-gaps, no-overlap picture).

[Matt J]

a gap occurs whenever the face of any prism doesn't completely touch another prism, or the box the prisms are contained in.

That sounds equivalent to saying that the prisms must completely fill the volume of the box. If so, that should be a simple thing to test, just by adding up the volumes of all the prisms. What are you trying to do beyond that?

[Alex]

Ah, yes, that is better wording. It should be easy to test if the box is completely filled, but I need to report the specific locations where there are gaps in the box, or overlaps between the prisms.

[Matt J]

How are those "locations" to be specified? A cavity in the box won't consist of a single point. It could have some arbitrary non-convex shape.

Similarly, in what way are the overlaps to be specified? The intersection of two overlapping prisms will form another prism. Do you just want the vertices of that prism?

[Alex]

Exactly - the cavity might not be a rectangular prism. The overlaps are specified by the vertices of the prisms that make the overlap (it should suffice to send the user the dimensions of the prisms (p1,p2,p3,p4,p5) that make the overlap). The specification for a gap is a little harder to define. The goal for gaps is to give the user an idea of where the gaps are, and then they would manually adjust the dimensions of the prisms to fill the cavity.

[Alex]

I tried your code with a different set of prisms (see attached diagram)

% Fit the following rectangular prisms in a box:
prism(1).x = [0 3];     % Min/max of x range for prism 1
prism(1).y = [0 6];     % Min/max of y range for prism 1
prism(1).z = [0 1];     % Min/max of z range for prism 1

prism(2).x = [0 5];
prism(2).y = [6 8];
prism(2).z = [0 1];

prisms(3).x = [5 10];
prisms(3).y = [6 8];
prisms(3).z = [0 1];

prisms(4).x = [3 10];
prisms(4).y = [3 6];
prisms(4).z = [0 1];

prisms(5).x = [3 10];
prisms(5).y = [0 3];
prisms(5).z = [0 1];

And it produced the following output:

T =

  10×2 table

      ij       contactType 
    ______    _____________

    1    2    "touching"   
    1    3    "unconnected"
    1    4    "touching"   
    1    5    "touching"   
    2    3    "touching"   
    2    4    "unconnected"
    2    5    "touching"   
    3    4    "unconnected"
    3    5    "touching"   
    4    5    "touching"   

The ideal solution here would be for essentially no output at all, as all of the prisms fit perfectly into a box.

[Alex]

Here's a better diagram of the box of prisms, attached.

[Matt J]

it should suffice to send the user the dimensions of the prisms (p1,p2,p3,p4,p5) that make the overlap

If so, I think the table T generated by my code gives you that.

and then they would manually adjust the dimensions of the prisms to fill the cavity.

I wonder if it might be better to undertake this using a 3D image, assuming it's possible to voxelize all the prisms. Then you can just form a logical map of the shape they form, together with the box that encloses them. The complement of this can be segmented into separate holes using bwconncomp or regionprops.

[Alex]

It's possible - I was also thinking about using a 3D grid of finely spaced points and just seeing if any of them fall in either multiple prisms, or no prisms. I figured I'd go for gold and try something more mathematically-based first.

[Alex]

For the record, the initial code that I posted doesn't work with this new set of prisms; it reports too many overlaps.