Given an original axis-aligned bounding box and a transformation matrix, you can compute a new axis-aligned bounding box from the extents of the original box, transformed by the matrix, like so:
// Keep a scratchpad array we can re-use for listing all 8 corners of a box.
Vector3[] corners = new Vector3[8];
// Fill our scratchpad array with the corners of an axis-aligned bounding box.
void PopulateCorners(AABB bounds) {
corners[0] = bounds.Min;
corners[1] = new Vector3(bounds.Min.x, bounds.Min.y, bounds.Max.z);
corners[2] = new Vector3(bounds.Min.x, bounds.Max.y, bounds.Min.z);
corners[3] = new Vector3(bounds.Max.x, bounds.Min.y, bounds.Min.z);
corners[4] = new Vector3(bounds.Min.x, bounds.Max.y, bounds.Max.z);
corners[5] = new Vector3(bounds.Max.x, bounds.Min.y, bounds.Max.z);
corners[6] = new Vector3(bounds.Max.x, bounds.Max.y, bounds.Min.z);
corners[7] = bounds.Max;
}
// Compute a new axis-aligned bounding box that will contain whatever the original
// bounds did, after an affine transformation. (Note this is a lossy operation)
AABB GetTransformedBounds(AABB original, Matrix transformation) {
PopulateCorners(original);
Vector3 min = new Vector3(1, 1, 1) * float.PositiveInfinity;
Vector3 max = new Vector3(1, 1, 1) * float.NegativeInfinity;
// Transform all of the corners, and keep track of the greatest and least
// values we see on each coordinate axis.
for(int i = 0; i < 8; i++) {
Vector3 transformed = transformation * corners[i];
min = Vector3.ComponentwiseMin(min, transformed);
max = Vector3.ComponentwiseMax(max, transformed);
}
return new AABB(min, max);
}
Note that you'll want to save the object's original "neutral orientation" bounding box separately. Anytime you want to compute a set of rotated bounds for the object, always work from this original, not from a previous rotated result. The reason is that an AABB might not be a tight fit - it might leave some room at the corners for example. If we rotate that box, then take the axis-aligned bounds of the rotated box, we'll tend to accumulate extra padding around the object, and our bounds will become less and less tight. Working always from the original ensures this error can't accumulate, and we get the best fit we can under the circumstances.
Of course, you can be more accurate still by always working from the source object/mesh/primitives. Transform all the points of the object (or its convex hull for simplicity) and then compute a new min & max from scratch off of those values. That will get the absolute tightest axis-aligned bounding box, like the one shown in the gif of the rotating knot in the article you link. But it's also often more work than just transforming the AABB.