matlab - 找出加盖矩形物体的方向、长度和半径

Question

我有一个如图 1 所示的图像。 我试图用一个带帽的矩形(图 2) 来拟合这个二值图像弄清楚:

方向(长轴和水平轴之间的角度)

物体的长度 (l) 和半径 (R)。最好的方法是什么？
谢谢您的帮助。

我非常天真的想法是使用最小二乘法来找出这些信息，但是我发现没有带帽矩形的方程。在 matlab 中有一个名为 rectangle 的函数可以完美地创建加盖矩形，但它似乎只是为了绘图目的。

Answer 1

我解决了这两种不同的方法，并在下面对每种方法都有注释。每种方法的复杂性各不相同，因此您需要为您的应用程序决定最佳交易。

第一种方法:最小二乘优化:
这里我通过 Matlab 的 fminunc() 函数使用了无约束优化。查看 Matlab 的帮助以查看您可以在优化之前设置的选项。我做了一些相当简单的选择，只是为了让这种方法适合你。

总之，我设置了一个带帽矩形的模型，作为参数 L、W 和 theta 的函数。如果您愿意，您可以包含 R，但我个人认为您不需要它；通过检查每个边缘的半圆的连续性，我认为通过检查模型几何形状，让 R = W 可能就足够了。这也将优化参数的数量减少了一个。

我使用 bool 图层制作了一个带帽矩形的模型，请参阅下面的 chedRectangle() 函数。因此，我需要一个函数来计算模型关于 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

然后我转换为 double 类型并通过规范化“清理”数据。注意:此预处理对于使优化正常工作很重要，并且可能需要，因为我没有您的原始数据(如前所述，我从该问题的网页下载了您的图像):

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 = 
           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]; 

享受。