找出加蓋矩形物體的方向，長度和半徑

Question

我有一個圖像，如圖1所示。 我試圖用封端的矩形（圖2），以適應該二進制圖像弄清楚：

1. 取向
2. 的長度（L（長軸和水平軸之間的角度） ）和物體的半徑（R）。什麼是最好的方式來做到這一點？ 感謝您的幫助。

我非常天真的想法是使用最小二乘法擬合來找出這些信息，但是我發現沒有方程式來限制矩形。在matlab中有一個稱爲矩形的函數，可以完美地創建頂點矩形，但它似乎只是爲了繪圖目的。

Answer 1

我解決了這2種不同的方式，並在下面的每種方法上有筆記。每種方法的複雜程度都不相同，因此您需要爲您的應用程序決定最佳交易。

第一種方法：最小二乘法優化： 這裏我使用Matlab的fminunc（）函數進行無約束優化。看看Matlab的幫助，看看你可以在優化之前設置的選項。我做了一些相當簡單的選擇，只是爲了讓這種方法適合你。

總之，我建立了一個封閉的矩形模型作爲參數L，W和θ的函數。如果你願意，你可以包括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值更大。請參見底部的更新結果。

通過比較數據和最終估計，即使是相對簡單的模型開始很好地模擬數據，您也可以看到。

您可以進一步計算誤差值，通過設置您自己的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 

通過與我們可以看到前面的值進行比較，新的最小均方誤差值是頗有幾分包括圖像中心在內，確認我們之前懷疑的是什麼（所以沒有什麼大驚喜）。

的最新估計數爲上限矩形的參數是這樣的：

xc = -22.9107 
yc = 35.9257 

優化需要較長的時間，但所看到的結果改進：

Convergence time (sec) = 96.0418 
    L0 = 25;  L estimate = 89.0784 
    W0 = 25;  W estimate = 80.4379 
theta0 = 90; theta estimate = 31.614 

相對於圖像陣列中心，我們得到

而且通過目視檢查：

如果性能是一個問題，您可能需要考慮編寫自己的優化程序，或者先嚐試調整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]; 

享受。

Answer 2

首先，我必須說，我沒有你所有問題的答案，但我可以幫助你定位。

我建議在二值圖像上使用主成分分析。關於PCA的好教程由Jon Shlens給出。在他的教程的圖2中有一個例子，它可以用於。在他的論文的第5部分中，您可以看到一些指導如何計算主要組成部分。利用奇異值分解，它更容易，如6.1節所示。

要使用PCA，您必須獲得要爲其計算主要組件的測量值。在你的情況下，每個白色像素是一個由像素位置(x, y)'表示的測量值。你將有N二維向量，給你的測量。因此，您的測量2xN矩陣X由級聯向量表示。

當你建立這個矩陣時，按照6.1節的規定進行。奇異值代表不同組件的「強度」。因此，最大奇異值代表橢圓的長軸。第二大奇異值（它應該只有兩個）表示另一個（垂直）軸。

請記住，如果橢圓是一個圓，你的奇異值應該是相等的，但與你的離散圖像表示，你不會得到一個完美的圓。