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

Question

我的图像如图1所示. 我试图用带帽的矩形拟合这个二进制图像(图2) 弄明白:

1. 方向(长轴和水平轴之间的角度)
2. 物体的长度(l)和半径(R).最好的方法是什么？谢谢您的帮助.

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

Answer 1

我解决了这两种不同的方式,并对下面的每种方法都有说明.每种方法的复杂程度各不相同,因此您需要确定应用程序的最佳交易.

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

总之,我根据参数L,W和theta设置了加盖矩形的模型.如果你愿意,你可以包括R,但我个人认为你不需要; 通过检查每个边缘的半半圆的连续性,我认为通过检查模型几何结构让R = W就足够了.这也将优化参数的数量减少了一个.

我使用布尔图层制作了加盖矩形的模型,请参阅下面的cappedRectangle()函数.结果,我需要一个函数来计算模型相对于L,W和θ的有限差分梯度.如果你没有为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的值更大.请参见底部的更新结果.

通过将数据与最终估计进行比较,您可以看到即使是相对简单的模型也能很好地开始与数据相似.

您可以进一步计算估算的误差条,方法是设置自己的蒙特卡罗模拟,以检查作为噪声和其他降级因素(您可以生成用于生成模拟数据的已知输入)的函数的准确度.

下面是我用于加盖矩形的模型函数(注意:我进行图像旋转的方式在数值上是粗略的,对于有限差分而言不是很稳健,但它的快速和肮脏并让你前进):

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'); 

然后是struct结果:

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

请享用.