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

二値化した画像に対して、ある領域に内接する最大の楕円を自動で設定したいです。

現状、面積を求めたいのですが画像が欠けており囲まれていないため、面積がうまく算出できません。そこでその領域に最大の楕円を設定できれば近似的に面積を算出できると考えています。マニュアルで楕円を設定することは可能ですがそれが最大のものかはわかりません。なので、自動的に最大のものを設定したいです。

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

Answer 1

[Hiroshi Iwamura]

色々なやり方がありそうで、コンペとかやっても面白そうな題材ですね。

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

埋めたい領域の中心辺りをクリックし、面積が表示されるまで待ちます。

何かキーを押すと終了します。

条件によって、 tunable parameters を調整する必要があるとは思います。

参考：MATLAB で計算幾何学 ～ 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 = 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')

[John D'Errico]

I removed the spam flag. The site bots sometimes get overzealous.

[Hiroshi Iwamura]

Thank you for your quick response!

[雄大]

回答ありがとうございます。

お手数をおかけしますがもしよろしければこのシステムがどのような原理になっているのか教えていただけないでしょうか。

[Hiroshi Iwamura]

polyshape 関数を使って、その領域内に元の白点が含まれているかどうかを判定しながら頂点座標を更新し、面積が増えなくなるか凸度が大きくなったら停止しています。

pgon.Vertices = pax; が polyshape の座標更新、

p.Shape.Vertices = pax; は表示の更新用です。

TFin = isinterior(pgon,bx,by); が当たり判定関数です。

リンク先に polyshape の基本的な説明は書いてありますが、上のコードをライブエディターにコピペし、TFin = ～辺りにブレイクポイントを置いて実行すると、マウスオーバーで各変数も全て見れるので動きが理解できるかと思います。

[雄大]

お時間空いてしまいましたが、こちらの方法で変数を色々と変えながら実行してみたのですが、

・白で繋がっていない場所から、飛び出してしまうことがあります。

　そこで、凸度を調整することで例えば、鈍角にすることができればなくなるのではないかと考えています。調整の仕方を教えていただけないでしょうか。

もしくは凸型の多角形しか生成しないようにするにはどのような変更を加えたらよろしいでしょうか。

・頂点座標を更新する1番初めは、とても小さい正多角形を設置しているという解釈であってますでしょうか。

[Hiroshi Iwamura]

形状の自由度（と実行速度）とのトレードオフになるので調整が難しそうですが、以下の辺りをお試しください。

・N = 16;　→　20　多角形の角数を増やす

・convMax = 1.5; % max convexity　→　1.0　最大凸度を減らす

・minStep = 2;　→　1　最小ステップサイズを減らす

・if convexity > convMax で、最大凸度超えた場合にブレイクするのではなくステップサイズを小さくしてやり直す

break ;　を　

-----

pax = paxPre;

p.Shape.Vertices = pax;

dr = dr / 2;

dr = max(dr,minStep);

-----

　に変更

> 頂点座標を更新する1番初めは、とても小さい正多角形を設置しているという解釈であってますでしょうか。

そうですね

while (~w) の

dr = minStep;

で最初の半径が決まります。

[雄大]

ありがとうございます。

・convMax=1.5や1.0というのは度数法で表すことはできますでしょうか。

・半径は中心座標から、minStep=1であれば1ずつ大きくしていくという解釈であってますか？つまり2週目では半径2、3週目では半径3...でしょうか

[Hiroshi Iwamura]

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

MATLAB コード ファイルのデバッグ

[雄大]

何から何までありがとうございます。

面積比と言いますと、新しく更新したものとその一つ前のものとの比較が1.５や設定した値になるということでしょうか。

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

Answer 2

[交感神経優位なあかべぇ]

楕円の作成ではないですが、画像内のそれぞれの囲まれてる黒の部分を塗りつぶしする方法でなら、面積が求められるかなぁと思い、試してみました。

画像内の黄色のエリアの面積は、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

[交感神経優位なあかべぇ]

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

Answer 3

[Kenta]

こんにちは、皆様の回答を限りなくシンプルにしたものを考えていました。あけべぇさんの結果（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)