orzSphericalSepFilter.m

function [sep dif]= orzSphericalSepFilter(img, r,wi,wo)
% 分離度フィルタの関数
% input
%    img  : 入力画像(double型)
%    flt1 : フィルタ１(領域１)
%    flt2 : フィルタ２(領域２)
% output
%    sep  : 出力(分離度マップ)
% 
% ---example
% 円形分離度フィルタの場合
%    r  : 円の半径
%    wi : 内円の厚み
%    wo : 外円の厚み
[x,y,z] = meshgrid(-r-wo:r+wo);
val = x.^2 + y.^2 + z.^2;
flt1 = (r-wi)^2 <= val & val <= r^2;
flt2 = r^2 < val & val <= (r+wo)^2;
% sep = SepFilter(img,flt1,flt2);


n1 = sum(flt1(:));
n2 = sum(flt2(:));
n = n1 + n2;
if ~isequal(size(flt1),size(flt2))
    [h1,w1] = size(flt1); [h2,w2] = size(flt2);
    flt = zeros(max([h1,w1],[h2,w2]));
    [h,w] = size(flt); h0 = ceil(h/2); w0 = ceil(w/2);
    Y = h0-floor((h1-1)/2) : h0+ceil((h1-1)/2); X = w0-floor((w1-1)/2) : w0+ceil((w1-1)/2);
    flt(Y,X) = flt(Y,X) + flt1;
    Y = h0-floor((h2-1)/2) : h0+ceil((h2-1)/2); X = w0-floor((w2-1)/2) : w0+ceil((w2-1)/2);
    flt(Y,X) = flt(Y,X) + flt2;
    flt = flt / n;
else
    flt = (flt1 + flt2) / n;
end
flt1 = flt1 / n1;
flt2 = flt2 / n2;
M = imfilter(img, flt, 'replicate');
MM = imfilter(img.^2, flt, 'replicate');
M1 = imfilter(img, flt1, 'replicate');
M2 = imfilter(img, flt2, 'replicate');

SB = n1*(M1-M).^2 + n2*(M2-M).^2;
ST = n*(MM - M.^2);
dif = (M2-M1);
sep = SB ./ ST;
sep(isnan(sep)) = 0;
sep=sep.*sign(dif);