找出加盖矩形物体的方向，长度和半径 [英] find out the orientation, length and radius of capped rectangular object
问题描述
我有一张如图1所示的图像。我试图适应这个带有上限矩形的二进制图像（图2）以确定：
I have a image as shown as fig.1. I am trying to fit this binary image with a capped rectangular (fig.2) to figure out:
- 方向（长轴和水平轴之间的角度）
- 长度（l）和半径（ R）对象。最好的方法是什么？
感谢您的帮助。
我非常天真的想法是使用最小二乘法来找出这些信息，但我发现了没有上限矩形的等式。在matlab中有一个叫做矩形的函数可以完美地创建上限矩形，但它似乎只是为了绘图目的。
My very naive idea is using least square fit to find out these information however I found out there is no equation for capped rectangle. In matlab there is a function called rectangle can create the capped rectangle perfectly however it seems just for the plot purpose.
推荐答案
我解决了这2种不同的方式，并在下面的每种方法有说明。每种方法的复杂程度各不相同，因此您需要为您的应用决定最佳交易。
I solved this 2 different ways and have notes on each approach below. Each method varies in complexity so you will need to decide the best trade for your application.
第一种方法：最小二乘优化：
这里我通过Matlab的fminunc（）函数使用无约束优化。看一下Matlab的帮助，看看在优化之前你可以设置的选项。我做了一些相当简单的选择，只是为了让这种方法适合你。
First Approach: Least-Squares-Optimization: Here I used unconstrained optimization through Matlab's fminunc() function. Take a look at Matlab's help to see the options you can set prior to optimization. I made some fairly simple choices just to get this approach working for you.
总之，我设置了一个带帽的矩形模型作为参数L的函数， W，和theta。如果你愿意，你可以包括R，但我个人认为你不需要;通过检查每个边缘的半半圆的连续性，我认为通过检查模型几何结构让R = W就足够了。这也将优化参数的数量减少了一个。
In summary, I setup a model of your capped rectangle as a function of the parameters, L, W, and theta. You can include R if you wish but personally I don't think you need that; by examining continuity with the half-semi-circles at each edge, I think it may be sufficient to let R = W, by inspection of your model geometry. This also reduces the number of optimization parameters by one.
我使用布尔图层制作了加盖矩形的模型，请参阅下面的cappedRectangle（）函数。因此，我需要一个函数来计算模型相对于L，W和θ的有限差分梯度。如果你没有为fminunc（）提供这些渐变，它会尝试估计这些，但我发现Matlab的估计对这个应用程序不起作用，所以我提供了自己作为错误函数的一部分，由fminunc调用（）（见下文）。
I made a model of your capped rectangle using boolean layers, see the cappedRectangle() function below. As a result, I needed a function to calculate finite difference gradients of the model with respect to L, W, and theta. If you don't provide these gradients to fminunc(), it will attempt to estimate these but I found that Matlab's estimates didn't work well for this application, so I provided my own as part of the error function that gets called by fminunc() (see below).
我最初没有您的数据，所以我只需右键点击上面的图片并下载：'aRIhm.png'
I didn't initially have your data so I simply right-clicked on your image above and downloaded: 'aRIhm.png'
为了读取你的数据，我做了这个（创建变量 cdata ）：
To read your data I did this (creates the variable cdata):
image = importdata('aRIhm.png');
vars = fieldnames(image);
for i = 1:length(vars)
assignin('base', vars{i}, image.(vars{i}));
end
然后我转换为双重类型并通过标准化清理数据。注意：此预处理对于使优化正常工作非常重要，并且可能因为我没有您的原始数据而需要（如上所述，我从网页上下载了此图片中的图片）：
Then I converted to double type and "cleaned-up" the data by normalizing. Note: this pre-processing was important to get the optimization to work properly, and may have been needed since I didn't have your raw data (as mentioned I downloaded your image from the webpage for this question):
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）的函数。
Note #1: you might consider adding the center of the capped rectangle, (xc,yc), as optimization parameters. These extra degrees of freedom will make a difference in the overall fitting results (see comment on final error function values below). I didn't set that up here but you can follow the approach I used for L, W, and theta, to add that functionality with the finite difference gradients. You will also need to setup the capped rectangle model as a function of (xc,yc).
编辑：出于好奇我在加盖上添加了优化矩形中心，请查看底部的结果。
注意＃2：表示封顶末端的连续性如果您愿意，可以在R，W，theta的示例之后包含R作为显式优化
参数。您甚至可能希望将每个端点的R1和R2作为变量？
Note #2: for "continuity" at the ends of the capped rectangle, let R = W. If you like, you can later include R as an explicit optimization parameter following the examples for L, W, theta. You might even want to have say R1 and R2 at each endpoint as variables?
下面是我用来简单说明示例优化的任意起始值。 我不知道您的申请中有多少信息，但一般来说，您应该尝试提供最佳的初步估算。
Below are arbitrary starting values that I used to simply illustrate an example optimization. I don't know how much information you have in your application but in general, you should try to provide the best initial estimates that you can.
L = 25;
W = L;
theta = 90;
params0 = [L W theta];
请注意，根据您的初步估算，您会得到不同的结果。
Note that you will get different results based on your initial estimates.
接下来显示起始估计值（稍后定义cappedRectangle（）函数）：
Next display the starting estimate (the cappedRectangle() function is defined later):
capRect0 = reshape(cappedRectangle(params0,nX,nY),nX,nY);
figure(2); imshow(capRect0);
为错误指标定义一个匿名函数（下面列出了errorFunc（））：
Define an anonymous function for the error metric (errorFunc() is listed below):
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的帮助）：
Below is the output from fminunc (for more details on each column see Matlab's help):
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
你可以看到最后的错误度量值相对于起始值没有减少那么多，这向我表明模型函数可能没有足够的degr自由地真正适应好的数据，所以考虑添加额外的优化参数，例如图像中心，如前所述。
You can see that the final error metric values have not decreased that much relative to the starting value, this indicates to me that the model function probably doesn't have enough degrees of freedom to really "fit" the data that well, so consider adding extra optimization parameters, e.g., image center, as discussed earlier.
编辑：已添加优化覆盖矩形中心，查看底部的结果。
现在打印结果（使用2011 Macbook Pro）：
Now print the results (using a 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的值更大。见底部的更新结果。
通过将数据与最终估计值进行比较，您可以看到即使是相对简单的模型也能很好地开始与数据相似。
You can see by comparing the data to the final estimate that even a relatively simple model starts to resemble the data fairly well.
您可以进一步计算估算的误差条，方法是设置自己的蒙特卡罗模拟，以检查作为噪声和其他降级因素的函数的准确度（使用已知输入生成模拟数据）。
You can go further and calculate error bars for the estimates by setting up your own Monte-Carlo simulations to check accuracy as a function of noise and other degrading factors (with known inputs that you can generate to produce simulated data).
下面是我用于加盖矩形的模型函数（注意：我进行图像旋转的方式是粗略的数字而不是数字ry对于有限差分而言非常强大，但它又快又脏，让你前进）：
Below is the model function I used for the capped rectangle (note: the way I did image rotation is kind of sketchy numerically and not very robust for finite-differences but its quick and dirty and gets you going):
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调用的错误函数：
And then you will also need the error function called by 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
以及计算有限差分梯度的函数：
As well as the function to calculate the finite difference gradients:
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（）。总的来说，我认为这可以通过一些调整和检查来确保它按照您的期望做到最好。所以这种方法就是这样称呼它（它肯定比第一种方法简单得多！）：
As others have pointed out you can also use Matlab's regionprops(). Overall I think this could work the best with some tuning and checking to insure that its doing what you expect. So the approach would be to call it like this (it certainly is a lot simpler than the first approach!):
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
这提供了足够的信息来输入加盖矩形的模型。乍一看，这似乎是要走的路，但看起来你的想法是另一种方法（也许是上面的第一种方法）。
This gives enough information to feed into a model of a capped rectangle. At first glance this seems like the way to go, but it seems like you have your mind set on another approach (maybe the first approach above).
无论如何，下面是结果图像（红色）覆盖在您可以看到的数据之上，看起来非常好：
Anyway, below is an image of the results (in red) overlaid on top of your data which you can see looks quite good:
编辑：我无法帮助自己，我怀疑通过将图像中心作为优化参数，可以获得更好的结果，所以我继续前进并做了它只是为了检查。果然，在最小二乘估计中使用的初始估计值相同，结果如下：
I couldn't help myself, I suspected that by including the image center as an optimization parameter, much better results could be obtained, so I went ahead and did it just to check. Sure enough, with the same starting estimates used earlier in the Least-Squares Estimation, here are the results:
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
与之前的价值相比，我们可以看到新的价值当包括图像中心时，误差平方误差值相当小，确认了我们之前所怀疑的（所以没什么大惊喜）。
By comparison with the earlier values we can see that the new least-square error values are quite a bit smaller when including the image center, confirming what we suspected earlier (so no big surprise).
加盖矩形的更新估计值因此参数是：
The updated estimates for the capped rectangle parameters are thus:
Convergence time (sec) = 96.0418
L0 = 25; L estimate = 89.0784
W0 = 25; W estimate = 80.4379
theta0 = 90; theta estimate = 31.614
相对于图像阵列中心，我们得到：
And relative to the image array center we get:
xc = -22.9107
yc = 35.9257
优化需要更长的时间，但视觉检查结果会有所改善：
The optimization takes longer but the results are improved as seen by visual inspection:
如果性能有问题，您可能需要考虑编写自己的优化器或首先尝试调整Matlab的优化参数，也许也使用不同的算法选项;请参阅上面的优化选项。
If performance is an issue you may want to consider writing your own optimizer or first try tuning Matlab's optimization parameters, perhaps using different algorithm options as well; see the optimization options above.
以下是更新型号的代码：
Here is the code for the updated model:
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（）之前我调整了参数列表：
and then prior to calling fminunc() I adjusted the parameter list:
L = 25;
W = L;
theta = 90;
% set image center to zero as intial guess:
xc = 0;
yc = 0;
params0 = [L W theta xc yc];
享受。
这篇关于找出加盖矩形物体的方向，长度和半径的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
