I just came here looking for the same answer. After shuddering at the thought of so much math involved, I thought I would resort to a semi-educated guess. Doodling a bit I came to the (intuitive and probably not entirely exact) conclusion that the largest rectangle is proportional to the outer resulting rectangle, and its two opposing corners lie at the intersection of the diagonals of the outer rectangle with the longest side of the rotated rectangle. For squares, any of the diagonals and sides would do... I guess I am happy enough with this and will now start brushing the cobwebs off my rusty trig skills (pathetic, I know).

Minor update... Managed to do some trig calculations. This is for the case when the Height of the image is larger than the Width.

Update. Got the whole thing working. Here is some js code. It is connected to a larger program, and most variables are outside the scope of the functions, and are modified directly from within the functions. I know this is not good, but I am using this in an isolated situation, where there will be no confusion with other scripts: redacted

I took the liberty of cleaning the code and extracting it to a function:

function getCropCoordinates(angleInRadians, imageDimensions) {
    var ang = angleInRadians;
    var img = imageDimensions;

    var quadrant = Math.floor(ang / (Math.PI / 2)) & 3;
    var sign_alpha = (quadrant & 1) === 0 ? ang : Math.PI - ang;
    var alpha = (sign_alpha % Math.PI + Math.PI) % Math.PI;

    var bb = {
        w: img.w * Math.cos(alpha) + img.h * Math.sin(alpha),
        h: img.w * Math.sin(alpha) + img.h * Math.cos(alpha)
    };

    var gamma = img.w < img.h ? Math.atan2(bb.w, bb.h) : Math.atan2(bb.h, bb.w);

    var delta = Math.PI - alpha - gamma;

    var length = img.w < img.h ? img.h : img.w;
    var d = length * Math.cos(alpha);
    var a = d * Math.sin(alpha) / Math.sin(delta);

    var y = a * Math.cos(gamma);
    var x = y * Math.tan(gamma);

    return {
        x: x,
        y: y,
        w: bb.w - 2 * x,
        h: bb.h - 2 * y
    };
}

I encountered some problems with the gamma-calculation, and modified it to take into account in which direction the original box is the longest.

-- Magnus Hoff

Solution

Given bounding box dimensions bx by by and t being the anticlockwise rotation of rectangle sized x by y:

x = (1/(cos(t)^2-sin(t)^2)) * (  bx * cos(t) - by * sin(t))
y = (1/(cos(t)^2-sin(t)^2)) * (- bx * sin(t) + by * cos(t))

Derivation

Why is this?

First, consider that the length bx is cut in two pieces, a and b, by the corner of the rectangle. Use trigonometry to express bx in terms of x, y, and theta:

bx = b          + a
bx = x * cos(t) + y * sin(t)            [1]

and similarly for by:

by = c          + d
by = x * sin(t) + y * cos(t)            [2]

1 and 2 can be expressed in matrix form as:

[ bx ] = [ cos(t)  sin(t) ] * [ x ]     [3]
[ by ]   [ sin(t)  cos(t) ]   [ y ]

Note that the matrix is nearly a rotation matrix (but not quite - it's off by a minus sign.)

Left-divide the matrix on both sides, giving:

[ x ] = inverse ( [ cos(t)  sin(t) ]    * [ bx ]                        [4]
[ y ]             [ sin(t)  cos(t) ] )    [ by ]

The matrix inverse is easy to evaluate for a 2x2 matrix and expands to:

[ x ] = (1/(cos(t)^2-sin(t)^2)) * [ cos(t) -sin(t) ] * [ bx ]           [5]
[ y ]                             [-sin(t)  cos(t) ]   [ by ]

[5] gives the two formulas:

x = (1/(cos(t)^2-sin(t)^2)) * (  bx * cos(t) - by * sin(t))             [6]
y = (1/(cos(t)^2-sin(t)^2)) * (- bx * sin(t) + by * cos(t))

Easy as pie!