if your cubes are axis aligned then you may construct tree interval Trees one for each axis.
An intervall tree is a data structure that let you query for an intervall overllapping or a point hit of a set of intervall.
If you can ensure that no cube will ever have an edge less than ε then you can easely see if a cube touches other cubes by querying your trees.
To understand how this works lets scale the problem down to "1D cubes" i.e. segments
It the testing segments is AB then take the points A-ε,A+ε,B-ε,B+ε check the intervall tree for the hit. For every segment that matches there are a few of possible relation with your testing segment:
Testing Segment A- A+ B- B+
Case - - - - impossible : no match
Case X - - - touch: the segments are disjoint but very near
Case - X - - overlap
Case X X - - overlap: segments touches for more that 2 epsilon
Case - - X - overlap
Case X - X - impossible: no holes in segments
Case - X X - overlap: the segments are nearly the same
Case X X X - overlap: the segments are nearly the same
Case - - - X touch: disjointed segments but very near
Case X - - X impossible: no holes in segments
Case - X - X impossible: no holes in segments
Case X X - X impossible: no holes in segments
Case - - X X overlap: segments touches for more that 2 epsilon
Case X - X X impossible: no holes in segments
Case - X X X overlap: the segments are nearly the same
Case X X X X overlap: the segments are nearly the same
As you can see a segment touches our testing one if A- matches but A+ does not OR B+ matches but B- does not; in any other case (excluding the impossible ones) there is an overlap.
In interval trees other data can be stored alongside the start and the end of the intervall. This possibility can be used to give names to the segments so you can wrap the query function so given a couple of points it gives you a list of (segment_name,"touch") or (segment_name,"overlap").
Now we need a way to use 1D results to get 2D results, a way that possibly scale up to 3D and over.
We can say that a rectangle touches another if there is a touch in an axis and an overlap in the other (edge touch) or a touch in both axes (corner touch).
To do this we need to store our intervals into two separate trees taking care of naming intervals correctly.
Let say that we are looking if the rectangle (Ax,Ay)(Bx,By)
touches some of the stored rectangle, then we have to query the X-tree for Ax,Bx
and the Y-tree for Ay,By
.
Lets suppose that the result is [("alpha","touch"),("beta","touch")]
by X-tree and [("alpha","overlap")]
for Y-tree whe can say that our rectangle touches "alpha" ("edge" touches) but don't touches "beta" nor other rectangles.
The test can be easely managed using a dictionary that keep the rectangle name as key and stores the X,Y values for example: {"alpha":[0,1],"beta":[0,-1]}
means that "alpha" touches(0) in X and overlap(1) in Y while "beta" touches(0) in X but is disjoint in Y(-1).
This reasoning can be extended further adding other dimensions: there is a touch if at least one dimension touches and there are no disjoint dimensions:
3D (0:touch, 1:overlap)
X Y Z
0 0 0 corner touch
0 0 1 edge touch (over one of the four Z edge)
0 1 0 edge touch (over one of the four Y edge)
0 1 1 surface touch (over one of the two X surfaces)
1 0 0 edge touch (over one of the four X edge)
1 0 1 surface touch (over one of the two Y surfaces)
1 1 0 surface touch (over one of the two Z surfaces)
1 1 1 eerrr no touch... the cubes are nested (remeber, at least one dimension touch)
[Appendix]
To speed up the interval query match you can test only A- and B+; when you find a match with A- then you get the upper end of the segment that matches and measure the distance with A-: if the distance is greather then 2ε this means that the overlap is too large to be considered a touch so you don't have to continue.
The same has to be done with B+ but this time you have to check the value of B+ with the lower end of the matching segments, again to see if the overlap you find with the extended segment is too large to be considered a simple touch.
A balanced interval tree takes ~log(n) iterations to check if a point hits against n intervals, this means that if you keep (by using red black trees for example) your trees balanced you can check against a huge set in reasonable time