最近在項目進行中遇到要提取離散點邊界的問題，像我這樣的對于matlab不是特別熟練的朋友一開始肯定摸不著頭腦，到底選用哪種算法可以有效地提取到所有已知點的輪廓線呢。本人經(jīng)過大量的文獻搜索及代碼實驗找到了幾個效果比較好的輪廓提取代碼，在這里做個總結(jié)，并且希望能夠?qū)τ龅酵瑯訂栴}的朋友有所啟發(fā)。

關(guān)于離散點邊界提取的三種方法：

1.Convhull 離散點集獲得邊界

2.Alpha Shape算法檢測邊緣點

3.Delaunay 三角剖分算法

前兩種方法在之前的博客中已經(jīng)做了總結(jié)這里就不展開了，現(xiàn)在主要介紹第三種算法。

該算法的總體思路如下：

1、利用 delaunay 函數(shù)，對所有數(shù)據(jù)點進行 Delaunay 三角剖分處理，delaunay 函數(shù)的返回值是一個 N * 3 的矩陣，其中 N 為剖分出的三角形個數(shù)，3 為每個三角形的三個端點的序號。

2、根據(jù) triangles 矩陣，提取出所有 delaunay 三角剖分時所連接的邊，依次掃描 triangles 矩陣的每一行，將 delaunay 三角剖分時所連接的邊添加到一個新的矩陣中，最后構(gòu)成一個 M * 2 的矩陣，其中 M 是一共所連接的邊的條數(shù)。

3、顯然，最小凸多邊形上的邊應(yīng)該僅在以上矩陣中出現(xiàn)一次，因此，將以上矩陣中那些出現(xiàn)次數(shù)超過一次的邊全部去掉，最后保留的便是最小凸多邊形的邊。

4、根據(jù)最小凸多邊形的邊，很容易得到構(gòu)成最小凸多邊形的結(jié)點的順序，從而解決問題。

輸入?yún)?shù) points 是一個 2 * P 矩陣， P 為數(shù)據(jù)點的個數(shù)，第一行是這些數(shù)據(jù)點對應(yīng)的 x 坐標，第二行是對應(yīng)的 y 坐標；輸出參數(shù) polygon 是一個 2 * Q 矩陣， Q 為凸多邊形的頂點個數(shù)（首尾相連），第一行是這些頂點對應(yīng)的 x 坐標，第二行是對應(yīng)的 y 坐標。代碼實現(xiàn)如下：

function polygon = minimal_convex_polygon(points)
 % 進行 delaunay 三角剖分，將所有連接了的邊保存在矩陣 lines 中
 triangles = sort(delaunay(points(1, :), points(2, :)), 2);
 lines = zeros(size(triangles, 1) * 3, 2);
 for i = 1:size(triangles, 1)
 lines(3 * i - 2,:) = [triangles(i, 1), triangles(i, 2)];
 lines(3 * i - 1,:) = [triangles(i, 1), triangles(i, 3)];
 lines(3 * i,:) = [triangles(i, 2), triangles(i, 3)];
 end
 % 去掉 lines 中出現(xiàn)次數(shù)超過一次的邊
 [~, IA] = unique(lines, 'rows');
 lines = setdiff(lines(IA, :), lines(setdiff(1:size(lines, 1), IA), :), 'rows');
 % 跟蹤 lines 中的數(shù)據(jù)點，將凸多邊形的頂點編號保存在 seqs 中
 seqs = zeros(size(lines, 1) + 1,1);
 seqs(1:2) = lines(1, :);
 lines(1, :) = [];
 for i = 3:size(seqs)
 pos = find(lines == seqs(i - 1));
 row = rem(pos - 1, size(lines, 1)) + 1;
 col = ceil(pos / size(lines, 1));
 seqs(i) = lines(row, 3 - col);
 lines(row, :) = [];
 end
 % 根據(jù) seqs ， 得到凸多邊形頂點坐標
 polygon = points(:, seqs);
end

定義了實現(xiàn)函數(shù)，下面進行調(diào)用：

plot(Pp(1,:),Pp(2,:), '*r', 'LineWidth', 4);  % Pp第一行為x坐標，第二行為y坐標
polygon = minimal_convex_polygon(Pp);
hold on;
plot(polygon(1, :), polygon(2, :), 'LineWidth', 2);

效果圖片我還不會添加進來，有興趣的朋友可以試一試。

以上就是本文的全部內(nèi)容，希望對大家的學(xué)習(xí)有所幫助，也希望大家多多支持腳本之家。

最新評論