Sphere - AABOX Edge detection

Question

I am trying to implement a particle-AABOX edge collision, the below images represent, two timesteps(dt) where the spherical particle is accelerated by gravity.

I have the particle center C and its radius r, along with AABOX coordinates as B_min and Bmax.

I have to check if the particle collides with any edges, and then bounce back, if a collision happens.

If I use Arvo's AABB vs Sphere collision algorithm, then I can find if there is an intersection within the radius, but I can't get the hitpoint and hitdistance.

And after getting the hitpoint and hitdistance, how to reflect it back? (for every timestep, I will have the current position C, and also the previous position P).

I tried to calculate it by implementing,

BOX_N[0] = make_vector(-1, 0, 0);// x-min;
BOX_P[0] = make_point(bb_min.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N[1] = make_vector(1, 0, 0); // x-max
BOX_P[1] = make_point(bb_max.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N[2] = make_vector(0, -1, 0);// y-min;
BOX_P[2] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N[3] = make_vector(0, 1, 0); // y-max
BOX_P[3] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_max.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);

BOX_N[4] = make_vector(0, 0, -1);// z-min;
BOX_P[4] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z);

BOX_N[5] = make_vector(0, 0, 1); // z-max
BOX_P[5] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_max.z);



__device__ int box_collision(const Point& previous_position, const Point& current_position, const Vector& direction, const float& radius, Point& hitPoint, Vector& normal, float& hit_distance) {
    Point boxPoint;
    int index;

         if (current_position.x < (bb_min.x + radius)) { index = 1; normal = BOX_N[0]; boxPoint = BOX_P[0]; boxPoint.x += radius;}
    else if (current_position.x > (bb_max.x - radius)) { index = 1; normal = BOX_N[1]; boxPoint = BOX_P[1]; boxPoint.x -= radius;}

    else if (current_position.y < (bb_min.y + radius)) { index = 2; normal = BOX_N[2]; boxPoint = BOX_P[2]; boxPoint.y += radius;}
    else if (current_position.y > (bb_max.y - radius)) { index = 2; normal = BOX_N[3]; boxPoint = BOX_P[3]; boxPoint.y -= radius;}

    else if (current_position.z < (bb_min.z + radius)) { index = 3; normal = BOX_N[4]; boxPoint = BOX_P[4]; boxPoint.z += radius;}
    else if (current_position.z > (bb_max.z - radius)) { index = 3; normal = BOX_N[5]; boxPoint = BOX_P[5]; boxPoint.z -= radius;}
    else return 0;

    auto denom = vdot(direction, normal);
    if (denom < 1e-6) return 0;

    hit_distance = vdot(normal, boxPoint - previous_position) / denom;
    if (hit_distance < 0) return 0;

    hitPoint = previous_position + hit_distance * direction;
    
    return index;
}

inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius) {
    Point hitPoint;
    Vector normal;
    float hit_distance;

    Vector direction = vnormalize(current_position - previous_position);

    int collision = box_collision(previous_position, current_position, direction, radius, hitPoint, normal, hit_distance);

    if (!collision) return true;

    Vector damping{1, 1, 1};
    if (collision == 1) damping.x = bounce_factor;
    else if (collision == 2) damping.y = bounce_factor;
    else if (collision == 3) damping.z = bounce_factor;
        
    float relection_length = vlength(hitPoint - current_position);
    
    Vector R = vnormalize(direction - 2.0f * vdot(normal, direction) * normal) * damping;
            
    current_position = hitPoint + R * relection_length;
    previous_position = hitPoint - R * hit_distance;

    return false;
}

by calling the 2nd function to check for collision,

inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius)

In the above code, I used the ray-plane interesection, and then use reflection of the direction to get the bounce direction, but after 500-550 timesteps, the solution deviates and becomes numerically unstable(possibly they settle at the bottom and collides continously). Is there something, I am missing here? Or am I doing something wrong?

Note: I am using simiulating 10k particles, so the solution deviation compounds by a large amount.

Answer 1

You're in a very simple case here. You can simply shrink the disc by its radius to make it a point, and shrink the box by the same radius on each side to make it a smaller box. Now you can test for collisions of this point against the smaller box, and extract from those the equivalent collisions of the disc versus the original box.

Instead of doing a sphere versus box check, now you can perform a raycast against the box (or against each side as a line segment), which is significantly simpler.

When you get a hit with the raycast, you can determine the point and normal of the collision. That's the yellow dot and arrow where the red displacement vector crosses the inner yellow box.

The normal will be the same for both the raycast hit and the disc hit (the blue arrow points in the same direction as the yellow one).

To convert the collision point from the raycast to the disc (the red point of contact in the diagram above), subtract the (unit) normal times the disc radius.

Since it looks like you're using Verlet integration, what you want to do next is reflect the start and end positions across the line of collision - the yellow line segment the ray collided with.

You can compute this by expressing your begin/end positions as offsets from the yellow point where the ray hits the shrunken box. Take the dot product of that offset vector with your normal vector, and subtract twice the result times the normal vector from the offset, and add it back to the collision point. In pseudo-code it would look a bit like...

Vector2 ReflectAcrossCollision(Vector2 position, Vector2 collisionPoint, Vector2 collisionNormal) {
    Vector2 offset = position - collisionPoint;
    Vector2 dot = offset.x * collisionNormal.x + offset.y * collisionNormal.y;
    offset -= 2 * dot * collisionNormal;
    return collisionPoint + offset;
}

Note that after doing this reflection in this example, our new reflected path also crosses the box. So you'll want to repeat your collision checks with the new start and end points (disregarding places where the ray "enters" the box, and counting only places where it exits)