領域に内接する最大の楕円を自動設定する方法

16 views (last 30 days)
雄大
雄大 on 20 Feb 2023
Commented: 雄大 on 4 Jul 2023
二値化した画像に対して、ある領域に内接する最大の楕円を自動で設定したいです。
現状、面積を求めたいのですが画像が欠けており囲まれていないため、面積がうまく算出できません。そこでその領域に最大の楕円を設定できれば近似的に面積を算出できると考えています。マニュアルで楕円を設定することは可能ですがそれが最大のものかはわかりません。なので、自動的に最大のものを設定したいです。
ご教授いただきたいです。

Accepted Answer

Hiroshi Iwamura
Hiroshi Iwamura on 25 Feb 2023
色々なやり方がありそうで、コンペとかやっても面白そうな題材ですね。
polyshape を使ってやってみました。
埋めたい領域の中心辺りをクリックし、面積が表示されるまで待ちます。
何かキーを押すと終了します。
条件によって、 tunable parameters を調整する必要があるとは思います。
close all
clear
tEn = true; % text on the area
% tunable parameters
N = 16; % number of vertices
convMax = 1.5; % max convexity
maxStep = 10;
minStep = 2;
areaNum = 1; areas = {}; % for keep each area
BW = imread('https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1300565/%E4%BA%8C%E5%80%A4%E5%8C%96%E5%89%8D%E3%81%AE%E7%94%BB%E5%83%8F.png');
BW = imbinarize(BW);
BW = BW(:,:,1);
imSize = size(BW);
% draw border box
BW(1,:) = 1; BW(end,:) = 1;
BW(:,1) = 1; BW(:,end) = 1;
f = figure;
ax = axes;
imshow(BW)
[by,bx] = find(BW==1);
f.WindowButtonDownFcn = 'pos = ax.CurrentPoint(1,1:2);';
w = waitforbuttonpress; % wait for the mouse click or key press
hold on
while (~w)
pax = repmat(pos,N+1,1); % start position
dr = minStep; % increase step size of radious
for k=0:N
th = k * 2*pi/N;
[dx,dy] = pol2cart(th,dr);
pax(k+1,:) = pax(k+1,:) + [dx,dy];
end
pgon = polyshape(pax(:,1),pax(:,2),Simplify=false);
TFin = isinterior(pgon,bx,by);
if sum(TFin) == 0
f.Pointer = 'watch'; % change cursor shape
p = plot(pgon,FaceAlpha= 0.75);
paxPre = pax;
sPre = 0;
while (pgon.area > sPre)
sPre = pgon.area;
for k=0:N
th = k * 2*pi/N;
[dx,dy] = pol2cart(th,dr);
pax(k+1,:) = pax(k+1,:) + [dx,dy];
pgon.Vertices = pax;
p.Shape.Vertices = pax;
TFin = isinterior(pgon,bx,by);
if sum(TFin) > 0
pax = paxPre;
p.Shape.Vertices = pax;
dr = dr / 2;
dr = max(dr,minStep);
else
dr = dr * 2;
dr = min(dr,maxStep);
end
paxPre = pax;
end
drawnow
polyout = convhull(pgon);
convexity =polyout.area / pgon.area;
if convexity > convMax
break;
end
end
% renew border line
F = getframe;
BW = imbinarize(F.cdata(:,:,[2 3]),1e-6); % except for red(text)
BW = BW(:,:,1);
[by,bx] = find(BW==1);
areaText = sprintf('Area = %d',round(pgon.area));
areas{areaNum} = num2str(round(pgon.area));
areaNum = areaNum + 1;
col = [0.9 0 0];
if tEn
text(min(pax(:,1)),pos(2),areaText,Color=col,fontsize=12)
end
f.Pointer = 'arrow'; % restore cursor shape
end
w = waitforbuttonpress;
end
f.WindowButtonDownFcn = '';
hold off
l = legend(areas,fontsize=14);
title(l,'Areas')
  9 Comments
Hiroshi Iwamura
Hiroshi Iwamura on 2 Jul 2023
Edited: Hiroshi Iwamura on 2 Jul 2023
> convMax
-----
polyout = convhull(pgon);
convexity = polyout.area / pgon.area;
if convexity > convMax
break;
end
-----
 ここで定義している通り、凸包との面積比です。
> 半径は中心座標から、minStep=1であれば~
-----
if sum(TFin) > 0
pax = paxPre;
p.Shape.Vertices = pax;
dr = dr / 2;
dr = max(dr,minStep);
else
dr = dr * 2;
dr = min(dr,maxStep);
end
-----
dr の更新は /2 か *2 です。
minStep は初期ステップサイズかつ最小ステップサイズです。
整数である必要はありません。
MATLAB でのデバッグ方法に慣れると、そういったことも理解しやすいかと思います。
雄大
雄大 on 4 Jul 2023
何から何までありがとうございます。
面積比と言いますと、新しく更新したものとその一つ前のものとの比較が1.5や設定した値になるということでしょうか。
もし、終了を内角の大きさで定義するとしたら、大変でしょうか。

Sign in to comment.

More Answers (2)

交感神経優位なあかべぇ
Edited: 交感神経優位なあかべぇ on 25 Feb 2023
楕円の作成ではないですが、画像内のそれぞれの囲まれてる黒の部分を塗りつぶしする方法でなら、面積が求められるかなぁと思い、試してみました。
画像内の黄色のエリアの面積は、areaの2番目(30708)になります。(areaの1番目(93228)は,青いエリアになります)
img = imread('https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1300565/%E4%BA%8C%E5%80%A4%E5%8C%96%E5%89%8D%E3%81%AE%E7%94%BB%E5%83%8F.png');
allAreamap = GetAreamap(img); % 画像の黒い部分に塗りつぶしを行った範囲のそれぞれの2値の画像を取得
area = arrayfun(@(i) nnz(allAreamap(:,:,i)), 1 : size(allAreamap, 3));
area = sort(area, 'descend') % 面積の大きいもの順に表示
area = 1×28
93228 30708 44 40 36 36 28 28 28 20 20 20 20 16 16 12 8 8 4 4 4 4 4 4 4 4 4 4
% 塗りつぶし画像の描画
colors = lines;
imshow(img);
hold on;
for i = 1 : size(allAreamap, 3)
areaData = allAreamap(:,:,i);
areaImg = repmat(0xFF, [size(img, [1,2]), 3]);
color = uint8(colors(i,:) * 255);
tmpImg = reshape(areaImg, prod(size(img, [1,2])), 3);
tmpImg(areaData,:) = repmat(color, nnz(areaData), 1);
areaImg = reshape(tmpImg, size(img));
image(areaImg, 'AlphaData', areaData);
end
function allAreamap = GetAreamap(img)
tmpImg = img(:,:,1) > 50 | img(:,:,2) > 50 | img(:,:,3) > 50;
img = true(size(img, [1,2]) + 2);
img(2:end-1, 2:end-1) = tmpImg;
nonMap = ~img;
allAreamap = false([size(img), 0]);
setIdx = 0;
while true
[row,col] = find(nonMap, 1);
if ~isempty(row)
area = FillInArea(row, col, img);
nonMap = xor(nonMap, area);
setIdx = setIdx + 1;
allAreamap(:,:,setIdx) = area;
else
break;
end
end
allAreamap = allAreamap(2:end-1,2:end-1,:);
end
function areamap = FillInArea(row, col, img)
areamap = false(size(img));
spread = false([size(img), 4]);
while true
areamap(row, col) = true;
spread(row, col, :) = true(4,1);
if img(row - 1, col) || areamap(row - 1, col)
spread(row, col, 1) = false;
end
if img(row, col - 1) || areamap(row, col - 1)
spread(row, col, 2) = false;
end
if img(row + 1, col) || areamap(row + 1, col)
spread(row, col, 3) = false;
end
if img(row, col + 1) || areamap(row, col + 1)
spread(row, col, 4) = false;
end
spread(row - 1, col, 3) = false;
spread(row, col - 1, 4) = false;
spread(row + 1, col, 1) = false;
spread(row, col + 1, 2) = false;
if any(spread(row, col, :))
idxOffset = GetNextIdx(spread(row, col, :));
row = row + idxOffset(1);
col = col + idxOffset(2);
else
spreadTruemap = any(spread, 3);
[row, col] = find(spreadTruemap, 1);
if ~isempty(row)
idxOffset = GetNextIdx(spread(row, col, :));
row = row + idxOffset(1);
col = col + idxOffset(2);
else
break;
end
end
end
end
function idxOffset = GetNextIdx(logicalArray)
idx = find(logicalArray, 1);
switch idx
case 1
idxOffset = [-1, 0];
case 2
idxOffset = [0,-1];
case 3
idxOffset = [1,0];
case 4
idxOffset = [0,1];
end
end
  1 Comment
交感神経優位なあかべぇ
Image Processing Toolboxがありましたら、穴の塗りつぶしを行うimfill関数が使用できますので、もう少し簡潔に記述できます。
img = imread('https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1300565/%E4%BA%8C%E5%80%A4%E5%8C%96%E5%89%8D%E3%81%AE%E7%94%BB%E5%83%8F.png');
allAreamap = GetAreamap(img); % 画像の黒い部分に塗りつぶしを行った範囲のそれぞれの2値の画像を取得
area = arrayfun(@(i) nnz(allAreamap(:,:,i)), 1 : size(allAreamap, 3));
area = sort(area, 'descend') % 面積の大きいもの順に表示
area = 1×28
93228 30708 44 40 36 36 28 28 28 20 20 20 20 16 16 12 8 8 4 4 4 4 4 4 4 4 4 4
% 塗りつぶし画像の描画
colors = lines;
imshow(img);
hold on;
for i = 1 : size(allAreamap, 3)
areaData = allAreamap(:,:,i);
areaImg = repmat(0xFF, [size(img, [1,2]), 3]);
color = uint8(colors(i,:) * 255);
tmpImg = reshape(areaImg, prod(size(img, [1,2])), 3);
tmpImg(areaData,:) = repmat(color, nnz(areaData), 1);
areaImg = reshape(tmpImg, size(img));
image(areaImg, 'AlphaData', areaData);
end
function allAreamap = GetAreamap(img)
img = img(:,:,1) > 50 | img(:,:,2) > 50 | img(:,:,3) > 50; % 2値画像の作成
allAreamap = false([size(img), 0]); % 塗りつぶしたエリアを表す画像の宣言(3次元目のそれぞれの塗りつぶし画像をコピー)
while true
[row,col] = find(~img, 1); % 穴があるピクセルを検索
if ~isempty(row)
area = imfill(img, [row, col]);% 穴が存在していたら塗りつぶしの実施
allAreamap(:,:,end + 1) = xor(area, img); % 塗りつぶしたエリアだけを抽出しコピー
img = area;
else
break; % 穴が全て埋まったら、終了
end
end
end

Sign in to comment.


Kenta
Kenta on 25 Feb 2023
こんにちは、皆様の回答を限りなくシンプルにしたものを考えていました。あけべぇさんの結果(30708)と似た結果になっています。こういうのでもいいかもしれませんね(?)。
img = imread('https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1300565/%E4%BA%8C%E5%80%A4%E5%8C%96%E5%89%8D%E3%81%AE%E7%94%BB%E5%83%8F.png');
% グレースケールに変換
gray = rgb2gray(img);
% 穴埋めをして差分を求める
Ifill = imfill(gray,'holes');
% 変化のあった箇所を可視化
diff = -single(gray)+single(Ifill);
% 5%を閾値として、塗った面積を求める
length(find(diff>prctile(diff(diff>0),5)))
ans = 31184
% 可視化
figure;imagesc(diff)

Categories

Find more on Chemistry in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!