计算旋转矩形中的最大矩形答案

【问题标题】：Calculate largest inscribed rectangle in a rotated rectangle计算旋转矩形中的最大矩形
【发布时间】：2011-08-12 22:52:07
【问题描述】：

我正在尝试找到计算可包含在旋转矩形内的最大（面积）矩形的最佳方法。

一些图片应该有助于（我希望）可视化我的意思：

输入矩形的宽度和高度是给定的，旋转的角度也是如此。输出矩形未旋转或倾斜。

我正在走漫长的路线，我什至不确定它是否能处理极端情况（不是双关语）。我确信有一个优雅的解决方案。有什么建议吗？

编辑：输出矩形点不一定要接触输入矩形的边缘。（感谢 E 先生）

【问题讨论】：

“最大的矩形”是指面积最大的那个吗？
@Sven 是的，就是这个意思。我会做一个编辑...谢谢。
@George Profenza 唯一的选择是写三千字......
这不是编程问题吗？
@zaf 看这里的图片：i.imgur.com/22yAQ.jpg，可能稍微旋转了一些。你怎么能把这样一个矩形放在这个里面？

标签： algorithm math language-agnostic geometry

【解决方案1】：

我只是来这里寻找相同的答案。在想到涉及如此多的数学时不寒而栗，我想我会求助于一个半受过教育的猜测。稍微涂鸦一下，我得出了一个（直观且可能不完全准确）的结论，即最大的矩形与外部生成的矩形成比例，并且它的两个对角位于外部矩形的对角线与最长边的交点旋转的矩形。对于正方形，任何对角线和边都可以……我想我对此很满意，现在将开始清除我生锈的三角技巧（可悲，我知道）。

小更新...设法进行了一些三角计算。这适用于图像的高度大于宽度的情况。

更新。得到了整个工作。这是一些js代码。它连接到一个更大的程序，并且大多数变量都在函数范围之外，并且直接在函数内部进行修改。我知道这不好，但我在孤立的情况下使用它，不会与其他脚本混淆：redacted

我冒昧地清理了代码并将其提取到一个函数中：

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
    };
}

我在gamma-计算中遇到了一些问题，并对其进行了修改以考虑原始框在哪个方向上最长。

-- 马格努斯霍夫

【讨论】：

漂亮的图形。我会考虑这个想法。如果您设法生成代码，请在此处发布！
我现在正在解决同样的问题。尝试为一些基于服务器的图像旋转和裁剪构建一个所见即所得的前端。我也做了一些计算。把它们贴在这里。作为图像......我还没有编码任何东西。
我最终使用了这个。谢谢！在此过程中，我重写了您的代码。我将其发布为编辑，因为我认为它更好，但请随时恢复或进一步编辑:)
哇，我差点忘了这件事。谢谢你的重写。
这个功能太棒了！我只是在一个黑客马拉松项目中使用它，如果没有它，我会迷路的。谢谢你们俩！ :)

【解决方案2】：

尽量不打破传统，将问题的解决方案作为图片：)

编辑： 第三个方程是错误的。正确的是：

3.w * cos(α) * X + w * sin(α) * Y - w * w * sin(α) * cos(α) - w * h = 0

要求解线性方程组，您可以使用Cramer rule 或Gauss method。

【讨论】：

如何将 P、Q、R、S 代入方程 1、2、3 和 4？请给出一个替换为 4 个方程之一的样本。谢谢。
P 应放在第一个方程中（即 (A, B) 线方程）。并且因为 P(x1, y1) 在那条线上，所以 x1 和 y1 应该是等式 w * cos(a) * x1 + w * sin(a) * y1 -w * w * sin(a) * cos(a) = 0 成立。
@Mihran Hovsepyan 谢谢。我会研究一下，看看我能不能摸到它。
Mihran - 我已经更新了我的答案，并提供了一个解决您问题的研究论文的链接。请在下面查看我的更新答案。
对不起@Jason Moore 你的意思是什么？我没有问题。

【解决方案3】：

首先，我们处理角度为零或 pi/2 的倍数的简单情况。那么最大的矩形就和原来的矩形一样了。

一般来说，内矩形在外矩形的边界上会有 3 个点。如果不是，则可以移动它，使一个顶点位于底部，一个顶点位于左侧。然后，您可以放大内部矩形，直到剩下的两个顶点之一碰到边界。

我们称外部矩形的边为 R1 和 R2。不失一般性，我们可以假设 R1

H cos a + W sin a <= R1
H sin a + W cos a <= R2

由于边界上至少有 3 个点，因此这些不等式中至少有一个实际上必须是等式。让我们使用第一个。不难看出：

W = (R1 - H cos a) / sin a

所以面积是

A = H W = H (R1 - H cos a) / sin a

我们可以采用 wrt 的导数。 H 并要求它等于 0：

dA/dH = ((R1 - H cos a) - H cos a) / sin a

求解 H 并使用上面的 W 表达式，我们发现：

H = R1 / (2 cos a)
W = R1 / (2 sin a)

在第二个不等式中代入，经过一些操作，

R1 (tan a + 1/tan a) / 2 <= R2

左边的因子总是至少为 1。如果不等式得到满足，那么我们就有了解决方案。如果不满足，则解决方案是将两个不等式都满足为等式的解决方案。换句话说：它是接触外部矩形的所有四个边的矩形。这是一个有 2 个未知数的线性系统，很容易求解：

H = (R2 cos a - R1 sin a) / cos 2a
W = (R1 cos a - R2 sin a) / cos 2a

根据原始坐标，我们得到：

x1 = x4 = W sin a cos a
y1 = y2 = R2 sin a - W sin^2 a 
x2 = x3 = x1 + H
y3 = y4 = y2 + W

【讨论】：

我会尽量找时间检查您的解决方案。你能看到一种快速获取目标内部矩形的 xy 位置（如果有多个位置就可以了）的方法吗？
确实，这似乎是正确区分这两种情况的唯一解决方案 1) R2 的长度足以获得关于 R1 的最优解（并且最优矩形不接触第四边） 2)最佳矩形接触所有 4 个边。案例 1) 有一个有趣的特性：面积最大的矩形在短边的中点与外矩形相接触。
我尝试了这个解决方案（我的问题在这里发布：stackoverflow.com/questions/16702966/…），但无法重现您的结果 - 您能否更新您的答案以包含完整的伪代码函数列表？
例如“外部矩形”是什么意思？ R1 和 R2 是原始矩形的尺寸吗？还是限定旋转矩形的较大矩形？
@aaronsnoswell 查看问题中的第二张图片。外部矩形是红色的。还要注意条件R1 <= R2。如果不是这种情况，您必须进行相应的调整。

【解决方案4】：

编辑：我在下面的 Mathematica 答案是错误的 - 我解决的问题与我认为您真正要问的问题略有不同。

为了解决您真正提出的问题，我将使用以下算法：

On the Maximum Empty Rectangle Problem

使用此算法，表示形成旋转矩形边界的有限数量的点（可能是 100 左右，并确保包括角） - 这些将是论文中描述的集合 S。

.

为了后代的缘故，我把我原来的帖子留在了下面：

面积最大的内部矩形始终是矩形的中下角（图表上靠近 alpha 的角）等于外部矩形宽度的一半的矩形。

我有点作弊并使用 Mathematica 为我解代数：

由此可以看出，内部矩形的最大面积等于 1/4 宽度^2 * 角的余割乘以角的正割。

现在我需要弄清楚这个最佳条件下底角的 x 值是多少。在我的面积公式上使用 mathematica 中的 Solve 函数，我得到以下信息：

这表明底角的x坐标等于宽度的一半。

现在只是为了确保，我将根据经验测试我们的答案。通过下面的结果，您可以看到我所有测试的最高区域（肯定不是详尽的，但您明白了）是底角的 x 值 = 外部矩形宽度的一半时。

【讨论】：

我从未使用过 Mathematica。您能否再扩展一些，以便我了解发生了什么？
Jason，请考虑参与mathematica标签。
@belisarius 当然，为什么不呢？ :-)
@Jason Moore 你能发布代码（不仅仅是代码的图像）吗

【解决方案5】：

@Andri 在我测试的width > height 的图像中无法正常工作。所以，我通过这种方式修复和优化了他的代码（只有两个三角函数）：

calculateLargestRect = function(angle, origWidth, origHeight) {
    var w0, h0;
    if (origWidth <= origHeight) {
        w0 = origWidth;
        h0 = origHeight;
    }
    else {
        w0 = origHeight;
        h0 = origWidth;
    }
    // Angle normalization in range [-PI..PI)
    var ang = angle - Math.floor((angle + Math.PI) / (2*Math.PI)) * 2*Math.PI; 
    ang = Math.abs(ang);      
    if (ang > Math.PI / 2)
        ang = Math.PI - ang;
    var sina = Math.sin(ang);
    var cosa = Math.cos(ang);
    var sinAcosA = sina * cosa;
    var w1 = w0 * cosa + h0 * sina;
    var h1 = w0 * sina + h0 * cosa;
    var c = h0 * sinAcosA / (2 * h0 * sinAcosA + w0);
    var x = w1 * c;
    var y = h1 * c;
    var w, h;
    if (origWidth <= origHeight) {
        w = w1 - 2 * x;
        h = h1 - 2 * y;
    }
    else {
        w = h1 - 2 * y;
        h = w1 - 2 * x;
    }
    return {
        w: w,
        h: h
    }
}

更新

我还决定发布以下用于比例矩形计算的函数：

calculateLargestProportionalRect = function(angle, origWidth, origHeight) {
    var w0, h0;
    if (origWidth <= origHeight) {
        w0 = origWidth;
        h0 = origHeight;
    }
    else {
        w0 = origHeight;
        h0 = origWidth;
    }
    // Angle normalization in range [-PI..PI)
    var ang = angle - Math.floor((angle + Math.PI) / (2*Math.PI)) * 2*Math.PI; 
    ang = Math.abs(ang);      
    if (ang > Math.PI / 2)
        ang = Math.PI - ang;
    var c = w0 / (h0 * Math.sin(ang) + w0 * Math.cos(ang));
    var w, h;
    if (origWidth <= origHeight) {
        w = w0 * c;
        h = h0 * c;
    }
    else {
        w = h0 * c;
        h = w0 * c;
    }
    return {
        w: w,
        h: h
    }
}

【讨论】：

感谢您的修复。我的答案是由 Magnus Hoff 在某个时候编辑的，我还没有测试过新版本。我知道旧的（丑陋的）版本可以工作，因为我已经使用它没有问题了 ~2 年了。
这种方法可以用来计算旋转矩形的边界框并进行一些调整吗？在我的项目中，我需要在旋转矩形时同时计算最大的矩形和 bbox，如果这能同时返回就太好了！
不适用于矩形（宽度和高度不相等）:(
已修复和清理...解决方案一点也不明显，如果没有您的实施，我将无法到达那里，所以谢谢！

【解决方案6】：

Coproc 在另一个线程 (https://stackoverflow.com/a/16778797) 上以简单有效的方式解决了这个问题。另外，他在那里给出了很好的解释和python代码。

下面是他的解决方案的我的 Matlab 实现：

function [ CI, T ] = rotateAndCrop( I, ang )
%ROTATEANDCROP Rotate an image 'I' by 'ang' degrees, and crop its biggest
% inner rectangle.

[h,w,~] = size(I);
ang = deg2rad(ang);

% Affine rotation
R = [cos(ang) -sin(ang) 0; sin(ang) cos(ang) 0; 0 0 1];
T = affine2d(R);
B = imwarp(I,T);

% Largest rectangle
% solution from https://stackoverflow.com/a/16778797

wb = w >= h;
sl = w*wb + h*~wb;
ss = h*wb + w*~wb;

cosa = abs(cos(ang));
sina = abs(sin(ang));

if ss <= 2*sina*cosa*sl
    x = .5*min([w h]);
    wh = wb*[x/sina x/cosa] + ~wb*[x/cosa x/sina];
else
    cos2a = (cosa^2) - (sina^2);
    wh = [(w*cosa - h*sina)/cos2a (h*cosa - w*sina)/cos2a]; 
end

hw = flip(wh);

% Top-left corner
tl = round(max(size(B)/2 - hw/2,1));

% Bottom-right corner
br = tl + round(hw);

% Cropped image
CI = B(tl(1):br(1),tl(2):br(2),:);

【讨论】：

【解决方案7】：

很抱歉这里没有给出推导，但我几天前在 Mathematica 中解决了这个问题，并提出了以下程序，非 Mathematica 的人应该能够阅读。如有疑问，请咨询http://reference.wolfram.com/mathematica/guide/Mathematica.html

下面的过程返回一个最大面积的矩形的宽度和高度，该矩形适合另一个宽度为 w 和高度为 h 的矩形已被 alpha 旋转。

CropRotatedDimensionsForMaxArea[{w_, h_}, alpha_] := 
With[
  {phi = Abs@Mod[alpha, Pi, -Pi/2]},
  Which[
   w == h, {w,h} Csc[phi + Pi/4]/Sqrt[2],
   w > h, 
     If[ Cos[2 phi]^2 < 1 - (h/w)^2, 
       h/2 {Csc[phi], Sec[phi]}, 
       Sec[2 phi] {w Cos[phi] - h Sin[phi], h Cos[phi] - w Sin[phi]}],
   w < h, 
     If[ Cos[2 phi]^2 < 1 - (w/h)^2, 
       w/2 {Sec[phi], Csc[phi]}, 
       Sec[2 phi] {w Cos[phi] - h Sin[phi], h Cos[phi] - w Sin[phi]}]
  ]
]

【讨论】：

感谢您的回答，欢迎堆栈溢出！
什么是Abs@Mod 和Csc 和Sec？

【解决方案8】：

这是最简单的方法... :)

Step 1
//Before Rotation

int originalWidth = 640;
int originalHeight = 480;

Step 2
//After Rotation
int newWidth = 701;  //int newWidth = 654;  //int newWidth = 513;
int newHeight = 564; //int newHeight = 757; //int newHeight = 664;

Step 3
//Difference in height and width
int widthDiff ;
int heightDiff;
int ASPECT_RATIO = originalWidth/originalHeight; //Double check the Aspect Ratio

if (newHeight > newWidth) {

    int ratioDiff = newHeight - newWidth;
    if (newWidth < Constant.camWidth) {
        widthDiff = (int) Math.floor(newWidth / ASPECT_RATIO);
        heightDiff = (int) Math.floor((originalHeight - (newHeight - originalHeight)) / ASPECT_RATIO);
    }
    else {
        widthDiff = (int) Math.floor((originalWidth - (newWidth - originalWidth) - ratioDiff) / ASPECT_RATIO);
        heightDiff = originalHeight - (newHeight - originalHeight);
    }

} else {
    widthDiff = originalWidth - (originalWidth);
    heightDiff = originalHeight - (newHeight - originalHeight);
}

Step 4
//Calculation
int targetRectanleWidth = originalWidth - widthDiff;
int targetRectanleHeight = originalHeight - heightDiff;

Step 5
int centerPointX = newWidth/2;
int centerPointY = newHeight/2;

Step 6
int x1 = centerPointX - (targetRectanleWidth / 2); 
int y1 = centerPointY - (targetRectanleHeight / 2);
int x2 = centerPointX + (targetRectanleWidth / 2);
int y2 = centerPointY + (targetRectanleHeight / 2);

Step 7
x1 = (x1 < 0 ? 0 : x1);
y1 = (y1 < 0 ? 0 : y1);

【讨论】：

widthDiff = originalWidth - (originalWidth); 永远为 0

【解决方案9】：

这只是Jeffrey Sax's solution above的插图，供我以后参考。

参考上图，解决方法是：

（我用了tan(t) + cot(t) = 2/sin(2t)这个身份）

【讨论】：