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

16 views (last 30 days)

Commented: 雄大 on 4 Jul 2023

ご教授いただきたいです。

Hiroshi Iwamura on 25 Feb 2023

polyshape を使ってやってみました。

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 = 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')
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 でのデバッグ方法に慣れると、そういったことも理解しやすいかと思います。

もし、終了を内角の大きさで定義するとしたら、大変でしょうか。

Edited: 交感神経優位なあかべぇ on 25 Feb 2023

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)
areaImg = repmat(0xFF, [size(img, [1,2]), 3]);
color = uint8(colors(i,:) * 255);
tmpImg = reshape(areaImg, prod(size(img, [1,2])), 3);
areaImg = reshape(tmpImg, size(img));
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));
while true
areamap(row, col) = true;
if img(row - 1, col) || areamap(row - 1, col)
end
if img(row, col - 1) || areamap(row, col - 1)
end
if img(row + 1, col) || areamap(row + 1, col)
end
if img(row, col + 1) || areamap(row, col + 1)
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;
row = row + idxOffset(1);
col = col + idxOffset(2);
else
if ~isempty(row)
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
Image Processing Toolboxがありましたら、穴の塗りつぶしを行うimfill関数が使用できますので、もう少し簡潔に記述できます。
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)
areaImg = repmat(0xFF, [size(img, [1,2]), 3]);
color = uint8(colors(i,:) * 255);
tmpImg = reshape(areaImg, prod(size(img, [1,2])), 3);
areaImg = reshape(tmpImg, size(img));
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

Kenta on 25 Feb 2023
こんにちは、皆様の回答を限りなくシンプルにしたものを考えていました。あけべぇさんの結果（30708）と似た結果になっています。こういうのでもいいかもしれませんね(?)。
% グレースケールに変換
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

### Community Treasure Hunt

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

Start Hunting!