AABB aren't that difficult to deal with, but I'm having one heck of a time implementing oriented bounding regions (rectangular prisms) in 3-space. My current algorithm is an extension of my understanding of the separated axis theorem (SAT) method, but unfortunately the only discussion I've heard of it so far was related to 2-space. The algorithm I use is to first compute the 8 vertices ("points" in code) and 12 edges (point-point[k]) as vectors. I find an axis perpendicular to the edge using the cross product of the edge vector and the shortest principle axis (code below), and project each point of both colliders onto the result. This "projection" is then used to compute the minimum and maximum bounds (as per the SAT method). Last I check if the bounds of each collider overlap. If they do, I continue with SAT and use the axis as a direction vector and the magnitude of the lowest overlap as the magnitude of the separating vector. If any do not overlap, the test ends.

It all sounds good and well, but I've made at least one serious error in my code because its behavior is sporradic. I do not get false positives, but the separating vector is never the shortest and testing some collider a vs b yields true but b vs a yields false (so, for instance, this algorithm always returns false unless I do NOT swap a and b in the first for loop but can return true if I don't) and rarely sufficient for separation (it's slightly too short). I'm not sure where I'm going wrong, either.

The root of the problem:

bool Collider::collisionCheck(ColliderPtr other,Vec3* separation){
getPoints();//update points
other->getPoints();//also update OTHER's points
F64 overlap=10000;//10000 is an arbitrarily large number. No collision should ever occur at this range.
//for every edge on EACH collider, get the axis, compute the bounds, and compare. If no overlap, return FALSE.
Collider* col[2] = {this,other.get()};//simplifies code by reducing duplicates. "this" = col[0], "other" = col[1]
for(U8 n=0;n<2;n++){
U8 a=n,b=-1*(n-1);//flip between a,b = 1,0 and 0,1 to eliminate redundant code
for(U8 i=0;i<12;i++){
AxisEdge axis = col[a]->getAxisEdge(i,true);//get the axis of the i'th edge
//project bounds onto 1d axis
F64 max=-10000,omax=-10000,min=10000,omin=10000,projection;
for(U8 j=0;j<8;j++){
projection = axis.project(col[a]->points[j]);//project the j'th point onto the axis
if(projection
min = projection;
if(projection>max)
max = projection;
}
//project other bounds onto 1d axis
for(U8 k=0;k<8;k++){
projection = axis.project(col->points[k]);//project the k'th point onto the axis
if(projection
omin = projection;
if(projection>omax)
omax = projection;
}
//check and report overlap
if(omin>=min&&omin<=max){
if(separation!=NULL&&overlap>=max-omin){
overlap = max-omin;
(*separation)=axis.axis.normal()*overlap;
}
}else if(omax>=min&&omax<=max){
if(separation!=NULL&&overlap>=omax-min){
overlap = omax-min;
(*separation)=axis.axis.normal()*overlap;
}
}else
return false;//Separation = no collision
}
}
return true;
}

The principle axis method for computing an edge's axis. Note that edge is the aforementioned point[x]-point[y] vector.
Vec3 AxisEdge::perpendicular(void){
Vec3 principle;
F64 min = std::min(edge.x,std::min(edge.y,edge.z));
if(compared(min,edge.x))
principle = Vec3(1,0,0);
else if(compared(min,edge.y))
principle = Vec3(0,1,0);
else
principle = Vec3(0,0,1);
return (edge%principle).normal();
}

I can provide any of the functions mentioned in the code. Also note that Vec3 operator% is the cross product of two vectors, while compared is a simple function that returns if abs(a-b)<1E-14?true:false.

I really appreciate any help on the matter, or any links that might be informative. I hope I've been clear, and I'll gladly clarify any points which I've left ambiguous. Thank you for your time and assistance.

Edit: After considerable consideration and experimentation, I believe the problem is with the axis selection code itself. While "any" perpendicular vector works in 2D (since there is only one and its parallels/opposites), it seems like only a limited set of orthogonal vectors are useful for ascertaining the direction of separation. Thus, the collision detection portion works, but causes the separation vector to be wrong in some cases, which is why separation fails, does not follow the shortest route, or behaves unexpectedly. Unfortunately, in every discussion on the topic I find, they simplify the problem by explaining it in 2D, which isn't helpful for axis selection in 3D. At least, as I understand it.

I have tried grouping the different edges by face, and using the normal of that face as the orthogonal vector to the edge. On half the edges, this works well. On the other half, the problem persists and the separation does not follow the shortest path (sometimes taking even the longest path). Is there a way to intelligently chose these orthogonal vectors?

Edit: I've resolved this issue, but thank you for reading. If anyone has questions, I'll gladly post the code that resolved the situation, or explain the logic I used to arrive at the solution. Thanks again, even if it was just for reading.