我有一个如图 1 所示的图像。
我试图用一个带顶盖的矩形（图 2）来拟合这个二进制图像
来弄清楚：
- 方向（长轴和水平轴之间的角度）
- 物体的长度 (l) 和半径 (R)。最好的方法是什么？谢谢您的帮助。
我非常天真的想法是使用最小二乘拟合来找出这些信息，但是我发现没有带盖矩形的方程。在 matlab 中有一个名为 rectangle 的函数可以完美地创建有盖的矩形，但它似乎只是为了绘图目的。
问问题
5611 次
2 回答
16
我解决了这两种不同的方法，并对下面的每种方法都有注释。每种方法的复杂性各不相同，因此您需要为您的应用程序确定最佳交易方式。
第一种方法：最小二乘优化： 这里我通过 Matlab 的 fminunc() 函数使用无约束优化。查看 Matlab 的帮助，了解您可以在优化之前设置的选项。我做了一些相当简单的选择，只是为了让这种方法适合你。
总之，我根据参数 L、W 和 theta 设置了一个带盖矩形的模型。如果您愿意，您可以包含 R，但我个人认为您不需要它；通过检查每个边缘的半圆的连续性，我认为通过检查模型几何形状让 R = W 可能就足够了。这也将优化参数的数量减少了一个。
我使用布尔图层制作了一个带帽矩形的模型，请参见下面的 cappedRectangle() 函数。因此，我需要一个函数来计算模型关于 L、W 和 theta 的有限差分梯度。如果您不向 fminunc() 提供这些梯度，它将尝试估计这些，但我发现 Matlab 的估计不适用于此应用程序，因此我提供了自己的估计作为 fminunc 调用的错误函数的一部分（） （见下文）。
我最初没有你的数据，所以我只是右键单击上面的图像并下载：'aRIhm.png'
为了读取您的数据，我这样做了（创建变量cdata）：
image = importdata('aRIhm.png');
vars = fieldnames(image);
for i = 1:length(vars)
assignin('base', vars{i}, image.(vars{i}));
end
然后我转换为双精度类型并通过规范化“清理”数据。注意：此预处理对于使优化正常工作很重要，并且可能需要，因为我没有您的原始数据（如前所述，我从网页下载了您的图像以解决此问题）：
data = im2double(cdata);
data = data / max(data(:));
figure(1); imshow(data); % looks the same as your image above
现在获取图像尺寸：
nY = size(data,1);
nX = size(data,2);
注意 #1：您可能会考虑添加有盖矩形的中心 (xc,yc) 作为优化参数。这些额外的自由度将对整体拟合结果产生影响（请参阅下面对最终误差函数值的评论）。我没有在这里设置，但您可以按照我用于 L、W 和 theta 的方法来添加具有有限差分梯度的功能。您还需要将封顶矩形模型设置为 (xc,yc) 的函数。
编辑：出于好奇，我在有盖矩形中心添加了优化，请参阅底部的结果。
注意 #2：对于有盖矩形末端的“连续性”，令 R = W。如果您愿意，您可以稍后将 R 作为显式优化参数包含在 L、W、theta 的示例之后。您甚至可能希望将每个端点的 R1 和 R2 作为变量？
下面是我用来简单说明示例优化的任意起始值。我不知道你的申请中有多少信息，但总的来说，你应该尽量提供最好的初步估计。
L = 25;
W = L;
theta = 90;
params0 = [L W theta];
请注意，根据您的初始估计，您将获得不同的结果。
接下来显示起始估计值（cappedRectangle() 函数稍后定义）：
capRect0 = reshape(cappedRectangle(params0,nX,nY),nX,nY);
figure(2); imshow(capRect0);
为错误度量定义一个匿名函数（errorFunc() 如下所示）：
f = @(x)errorFunc(x,data);
% Define several optimization parameters for fminunc():
options = optimoptions(@fminunc,'GradObj','on','TolX',1e-3, 'Display','iter');
% Call the optimizer:
tic
[x,fval,exitflag,output] = fminunc(f,params0,options);
time = toc;
disp(['convergence time (sec) = ',num2str(time)]);
% Results:
disp(['L0 = ',num2str(L),'; ', 'L estimate = ', num2str(x(1))]);
disp(['W0 = ',num2str(W),'; ', 'W estimate = ', num2str(x(2))]);
disp(['theta0 = ',num2str(theta),'; ', 'theta estimate = ', num2str(x(3))]);
capRectEstimate = reshape(cappedRectangle(x,nX,nY),nX,nY);
figure(3); imshow(capRectEstimate);
以下是 fminunc 的输出（有关每列的更多详细信息，请参阅 Matlab 的帮助）：
Iteration f(x) step optimality CG-iterations
0 0.911579 0.00465
1 0.860624 10 0.00457 1
2 0.767783 20 0.00408 1
3 0.614608 40 0.00185 1
.... and so on ...
15 0.532118 0.00488281 0.000962 0
16 0.532118 0.0012207 0.000962 0
17 0.532118 0.000305176 0.000962 0
您可以看到最终的误差度量值相对于起始值没有减少太多，这表明模型函数可能没有足够的自由度来真正“拟合”数据，所以考虑添加如前所述，额外的优化参数，例如图像中心。
编辑：在有盖矩形中心上添加了优化，请参阅底部的结果。
现在打印结果（使用 2011 Macbook Pro）：
Convergence time (sec) = 16.1053
L0 = 25; L estimate = 58.5773
W0 = 25; W estimate = 104.0663
theta0 = 90; theta estimate = 36.9024
并显示结果：
编辑：上面拟合结果的夸大“厚度”是因为模型试图拟合数据，同时保持其中心固定，导致 W 值更大。请参阅底部的更新结果。
通过将数据与最终估计值进行比较，您可以看到，即使是相对简单的模型也开始与数据非常相似。
您可以通过设置您自己的 Monte-Carlo 模拟来检查精度作为噪声和其他退化因素的函数（使用您可以生成以生成模拟数据的已知输入）来进一步计算估计的误差线。
下面是我用于加盖矩形的模型函数（注意：我进行图像旋转的方式在数值上有点粗略，对于有限差分来说不是很健壮，但它又快又脏，可以让你继续前进）：
function result = cappedRectangle(params, nX, nY)
[x,y] = meshgrid(-(nX-1)/2:(nX-1)/2,-(nY-1)/2:(nY-1)/2);
L = params(1);
W = params(2);
theta = params(3); % units are degrees
R = W;
% Define r1 and r2 for the displaced rounded edges:
x1 = x - L;
x2 = x + L;
r1 = sqrt(x1.^2+y.^2);
r2 = sqrt(x2.^2+y.^2);
% Capped Rectangle prior to rotation (theta = 0):
temp = double( (abs(x) <= L) & (abs(y) <= W) | (r1 <= R) | (r2 <= R) );
cappedRectangleRotated = im2double(imrotate(mat2gray(temp), theta, 'bilinear', 'crop'));
result = cappedRectangleRotated(:);
return
然后你还需要 fminunc 调用的错误函数：
function [error, df_dx] = errorFunc(params,data)
nY = size(data,1);
nX = size(data,2);
% Anonymous function for the model:
model = @(params)cappedRectangle(params,nX,nY);
% Least-squares error (analogous to chi^2 in the literature):
f = @(x)sum( (data(:) - model(x) ).^2 ) / sum(data(:).^2);
% Scalar error:
error = f(params);
[df_dx] = finiteDiffGrad(f,params);
return
以及计算有限差分梯度的函数：
function [df_dx] = finiteDiffGrad(fun,x)
N = length(x);
x = reshape(x,N,1);
% Pick a small delta, dx should be experimented with:
dx = norm(x(:))/10;
% define an array of dx values;
h_array = dx*eye(N);
df_dx = zeros(size(x));
f = @(x) feval(fun,x);
% Finite Diff Approx (use "centered difference" error is O(h^2);)
for j = 1:N
hj = h_array(j,:)';
df_dx(j) = ( f(x+hj) - f(x-hj) )/(2*dx);
end
return
第二种方法：使用 regionprops()
正如其他人指出的那样，您也可以使用 Matlab 的 regionprops()。总的来说，我认为这可以通过一些调整和检查来确保它达到你的预期效果。所以方法是这样称呼它（它肯定比第一种方法简单得多！）：
data = im2double(cdata);
data = round(data / max(data(:)));
s = regionprops(data, 'Orientation', 'MajorAxisLength', ...
'MinorAxisLength', 'Eccentricity', 'Centroid');
然后是结构结果s：
>> s
s =
Centroid: [345.5309 389.6189]
MajorAxisLength: 365.1276
MinorAxisLength: 174.0136
Eccentricity: 0.8791
Orientation: 30.9354
这提供了足够的信息来输入有盖矩形的模型。乍一看，这似乎是要走的路，但似乎您已经决定采用另一种方法（可能是上面的第一种方法）。
无论如何，下面是覆盖在您的数据之上的结果图像（红色），您可以看到它看起来相当不错：
编辑： 我忍不住，我怀疑通过将图像中心作为优化参数，可以获得更好的结果，所以我继续做只是为了检查。果然，使用前面最小二乘估计中使用的相同起始估计，结果如下：
Iteration f(x) step optimality CG-iterations
0 0.911579 0.00465
1 0.859323 10 0.00471 2
2 0.742788 20 0.00502 2
3 0.530433 40 0.00541 2
... and so on ...
28 0.0858947 0.0195312 0.000279 0
29 0.0858947 0.0390625 0.000279 1
30 0.0858947 0.00976562 0.000279 0
31 0.0858947 0.00244141 0.000279 0
32 0.0858947 0.000610352 0.000279 0
通过与之前的值相比，我们可以看到新的最小二乘误差值在包括图像中心时要小得多，这证实了我们之前的猜测（所以没什么大惊小怪的）。
因此，带盖矩形参数的更新估计值为：
Convergence time (sec) = 96.0418
L0 = 25; L estimate = 89.0784
W0 = 25; W estimate = 80.4379
theta0 = 90; theta estimate = 31.614
并且相对于我们得到的图像数组中心：
xc = -22.9107
yc = 35.9257
优化需要更长的时间，但通过目测可以看出结果得到了改善：
如果性能是一个问题，您可能需要考虑编写自己的优化器或首先尝试调整 Matlab 的优化参数，也可能使用不同的算法选项；请参阅上面的优化选项。
这是更新模型的代码：
function result = cappedRectangle(params, nX, nY)
[X,Y] = meshgrid(-(nX-1)/2:(nX-1)/2,-(nY-1)/2:(nY-1)/2);
% Extract params to make code more readable:
L = params(1);
W = params(2);
theta = params(3); % units are degrees
xc = params(4); % new param: image center in x
yc = params(5); % new param: image center in y
% Shift coordinates to the image center:
x = X-xc;
y = Y-yc;
% Define R = W as a constraint:
R = W;
% Define r1 and r2 for the rounded edges:
x1 = x - L;
x2 = x + L;
r1 = sqrt(x1.^2+y.^2);
r2 = sqrt(x2.^2+y.^2);
temp = double( (abs(x) <= L) & (abs(y) <= W) | (r1 <= R) | (r2 <= R) );
cappedRectangleRotated = im2double(imrotate(mat2gray(temp), theta, 'bilinear', 'crop'));
result = cappedRectangleRotated(:);
然后在调用 fminunc() 之前我调整了参数列表：
L = 25;
W = L;
theta = 90;
% set image center to zero as intial guess:
xc = 0;
yc = 0;
params0 = [L W theta xc yc];
享受。
于 2013-11-27T01:23:55.657 回答
1
首先，我不得不说，我没有你所有问题的答案，但我可以帮助你进行定向。
我建议使用二值图像的主成分分析。Jon Shlens给出了关于 PCA 的一个很好的教程。在他的教程的图 2 中有一个示例，它可以用于什么。在他论文的第 5 节中，您可以看到一些关于如何计算主成分的说明。如第 6.1 节所示，使用奇异值分解要容易得多。
要使用 PCA，您必须获得要计算主成分的测量值。在您的情况下，每个白色像素都是由像素位置表示的测量值(x, y)'。您将有N二维向量来提供您的测量值。因此，您的测量2xN矩阵X由级联向量表示。
当你建立了这个矩阵后，按照第 6.1 节中的说明进行操作。奇异值代表不同组件的“强度”。因此，最大的奇异值代表椭圆的长轴。第二大奇异值（它应该只有两个）代表另一个（垂直）轴。
请记住，如果椭圆是圆形，则您的奇异值应该相等，但是使用离散图像表示，您将不会得到完美的圆形。
于 2013-11-22T22:45:40.137 回答