２次元グラフの面積取得方法

２次元の座標データをもっていて

点に囲まれた座標の面積を求めたいです．

polyareaを使ってみましたが，形状が変わってしまいます．　（sample.png)

左図：もとめたい形状

右図：plot した画像　→　polyareaはこの形状の面積を求めている？

左図の形状を保った状態で，面積を取得する方法はありますか？

座標の分解能は保持したいです．

[Hernia Baby]

面積というのは総面積ですか？それとも各領域ごとの面積ですか？

[H.O]

総面積です。 欲を言えば各領域ごとの面積です。

[Atsushi Ueno]

凸包 - MATLAB convhull - MathWorks 日本

横やり失礼します。

点群の凸包やその面積の演算はMATLABの関数で出来ますが、部分的な領域をどう認識させるかが問題ですね。

予めどの塊がどういう順番で来るかが把握できていれば楽ですが、そうでなければ認識させる必要があります。また点群の塊が凸形状でない場合や、点群の密度が荒い場合は、正確な面積が出せないおそれがあります。

n = 100;
for k = 0:9 % ランダム値でサンプル点群データを生成
    x(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
    y(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
end
[k,av] = convhull(x,y); % 凸包を計算
av % 面積

av = 4.0018e+05

plot(x,y,'*')
hold on
plot(x(k),y(k))

[Atsushi Ueno]

2 次元および 3 次元の点で構成される多角形および多面体 - MATLAB - MathWorks 日本

上記なら上手く計算出来そうです。

n = 100;
for k = 0:9 % ランダム値でサンプル点群データを生成
    x(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
    y(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
end
shp = alphaShape(x',y'); % alphaShape オブジェクトを作成
A = area(shp) % 面積を計算

A = 7.7024e+04

plot(shp)

[Atsushi Ueno]

領域境界のタイプ - MATLAB & Simulink - MathWorks 日本

convhull関数はboundary関数の係数sが0の時であると知りました。同係数sを0⇒1に変化させていくと、2次元点群をラバーバンドで括った形状が徐々に締め付けられていくのが判ります。

n = 100;
for k = 0:9 % ランダム値でサンプル点群データを生成
    x(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
    y(k*n+1:(k+1)*n) = rand(100,1)*n + randi(10,1,1)*n;
end
plot(x,y,'*')
hold on
h = plot(x(k),y(k));
for s = 0:0.01:1
    [k,av] = boundary(x',y',s); % 境界を計算%[k,av] = convhull(x,y); % 凸包を計算
    %av % 面積
    h.XData = x(k);
    h.YData = y(k);
    drawnow
end

[H.O]

ありがとうございます．

実はアルファシェイプは最初に試したのですが，

各領域の面積や領域間との距離が小さいためか，アルファシェイプだと形状の変形が大きかったので，

ほかの方法を探していました．

申し訳ありません．

Answer 1

[Hiroshi Iwamura]

外枠いらなければ、最後のループも不要ですね。

I = imread("https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1345389/sample_c.png");

BW = imbinarize(I(:,:,2));

se = strel('disk',2);
BW2 = imclose(BW,se);
BW3 = imopen(BW2,se);
montage({I,BW,BW2,BW3})

stats = regionprops(BW3,'basic');
L = bwlabel(BW3);

T = struct2table(stats);


imshow(L,[],Colormap=jet)

idx = (T.Area > 100);
text(T.BoundingBox(idx,1),T.Centroid(idx,2),num2str(T.Area(idx)),'Color','white','FontSize',10);

fprintf('Total Area = %d\n',sum(T.Area(T.Area>10)))

Total Area = 4289

[Atsushi Ueno]

的確な回答ですが、問題は要の regionprops 関数等に ImageProcessingToolbox が必要な点です。

と思ったのですが、以前の@H.Oさんの質問でImageProcessingToolboxを使う回答にOKのメッセージがあるので、良いみたいですね。

[Hiroshi Iwamura]

なるほど

Computer Vision Toolbox もお持ちのようですね。

Answer 2

[Hiroshi Iwamura]

なかなか正確に求めるのは難しく調整が必要になりますが、Image Processing Toolbox をお持ちであればモフォロジーを使う手はあります。Simulink (Computer Vision Toolbox) でやった方が色々と調整が簡単かもしれません。

I = imread("sample_c.png");

BW = imbinarize(I(:,:,2));

se = strel('disk',2);
BW2 = imclose(BW,se);
BW3 = imopen(BW2,se);

montage({I,BW,BW2,BW3})

stats = regionprops(BW3,'basic');
L = bwlabel(BW3);

T = struct2table(stats);
T = sortrows(T,'Area','descend');

% imshow(I)
imshow(L,[],Colormap=jet)
hold on
n = 1;
while T.Area(n) > 100
    rectangle('Position',T.BoundingBox(n,:),EdgeColor=[1 0.2 0])
    text(T.BoundingBox(n,1),T.Centroid(n,2),num2str(T.Area(n)),'Color','white','FontSize',10)
    n = n + 1;
end
hold off

fprintf('Total Area = %d\n',sum(T.Area(T.Area>10)))

Total Area = 4289

[H.O]

ありがとうございます。 必要なtoolboxも運良く所有していました。 挙げて頂いたスクリプトで実行、面積取得(pixel)できました。