| [S, MG, Nmax]=seg_fuzzy(I,smooth,rgbImage,Nm,Nth,weighted) |
function [S, MG, Nmax]=seg_fuzzy(I,smooth,rgbImage,Nm,Nth,weighted)
%
% SEG_FUZZY Fuzzy thresholding segmentation of 2D or 3D data
% V. 3.0
%
% INPUT
% I input image (2D o 3D data)
% smooth (=1) gaussian smoothing of data (default)
% rgbImage (=1) if the image is rgb (default 0)
% (for rgb size MxNx3)
% Nm: normalization
% =1 absolute max
% =2 relative max (default)
% Nth: Maximum of sets to use
% (default is the maximum posible)
% =0 means the maximum posible
% weighted (=0) no histogram ponderation (default)
%
% OUTPUT
% S: Output value
% MG: Local Agregation (based on median)
% Nmax Number of maxima in histogram
% (Number of fuzzy sets)
%
% Implementation done up to Nth maxima in histogram output sets.
% If histogram shows more than Nth maxima reduced to the top Nth.
% In this case output Nmax has two values:
%
% Nmax=[Real_Nmax, 5]
%
% Algorithm proposed in:
%
% Santiago Aja-Fernández, Gonzalo Vegas-Sánchez-Ferrero,
% Miguel A. Martín Fernández, Soft thresholding for medical
% image segmentation, EMBC'2010, Buenos Aires, Sept. 2010.
%
% EXAMPLE
%
% [S,MG,Nmax]=seg_fuzzy(I);
%
% DEFAULT
%
% seg_fuzzy(I)=seg_fuzzy(I,0,1,2,0,0)
%
%
% Santiago Aja-Fernandez (V1.0), Ariel Hernan Curiale (V3.0)
% LPI V3.0
% www.lpi.tel.uva.es/~santi
% LPI Valladolid, Spain
% Original: 06/05/2012,
if(~exist('rgbImage','var'))
rgbImage = 0;
end
if(~exist('smooth','var'))
smooth = 1;
end
if(~exist('Nm','var'))
Nm = 2;
end
if(~exist('Nth','var'))
Nth = 0;
end
if(~exist('weighted','var'))
weighted = 0;
end
dim = ndims(I);
if(rgbImage == 1)
dim = dim -1;
end
if(rgbImage)
I=mean(double(I),dim+1);
else
I = double(I);
end
%
%Smooth original image for better thresholding
if smooth==1
if(dim == 2)
I2=I;
I2b=[I2(1,:); I2(1,:);I2(1,:);I2(1,:);I2(1,:);
I2;
I2(end,:); I2(end,:);I2(end,:);I2(end,:);I2(end,:)];
I2b=[I2b(:,1) I2b(:,1) I2b(:,1) I2b(:,1) I2b(:,1) ...
I2b ...
I2b(:,end) I2b(:,end) I2b(:,end) I2b(:,end) I2b(:,end)];
window = fspecial('gaussian', 11, 1.5);
I3=filter2(window,I2b);
I3=I3(6:end-5,6:end-5);
I=I3;
elseif(dim == 3)
% Gaussian smooth
siz=[5 5 5];% 11 in each dimension
sig = [1.5 1.5 1.5]; % 1.5 in each dimension
siz = (siz-1)/2;
[x,y,z] = ndgrid(-siz(1):siz(1),-siz(2):siz(2),-siz(3):siz(3));
h = exp(-(x.*x/2/sig(1)^2 + y.*y/2/sig(2)^2 + z.*z/2/sig(3)^2));
window = h/sum(h(:));
I = imfilter(I,window,'replicate');
else
error('Unsupported Image Dimension !');
end
end
%Normalization-------------
if Nm==1
Mw=max(I(:));
else
S1=std(I(:));
M1=mean(I(:));
Mw=max(I((I>(M1-2.*S1))&(I<(M1+2.*S1))));
end
I=I./Mw;
%Histogram of the Image
[h1,x]=hist(I(:),100);
h2=conv(ones([7,1])/7,h1);
h=h2(4:(end-3));
h=(h./sum(h))./(x(2)-x(1));
Dr=diff(h);
% Find the points when the function (dDr) is 0.
M=[Dr(1:end-1).*Dr(2:end)];
Maxt=find(double(M<0).*(Dr(1:end-1)>0))+1;
Nmax=length(Maxt); %(Number of fuzzy sets)
%MAX search
%Modeling the Histogram as a Gaussian mixture
if Nmax==1
error('Only one maximum find. Segmentation is no possible');
elseif Nth == 0 || Nth > Nmax
Nth = Nmax;
elseif Nmax > Nth
Nmax=[Nmax,Nth];
[A,B]=sort(Dr(Maxt));
Maxt=Maxt(B);
Maxt=Maxt(1:Nth);
end
total_sum = sum(h(Maxt(1:Nth)));%weighted
alpha_i = 1/Nth; %uniform
sigma_i = 0.1;
X0 = [];
for iNmax=1:Nth
%[X0,alpha_i,mu_i,sigma_i]
if(weighted == 1)
X0 = [X0, h(Maxt(iNmax))/total_sum,x(Maxt(iNmax)),sigma_i];
else
X0 = [X0, alpha_i,x(Maxt(iNmax)),sigma_i];
end
end
X=fminsearch(@(y)minimizationFunction(y,h,x,Nth),X0);
for iNmax=1:Nth
id = 3*(iNmax-1)+1;
Ci=normpdf(x,X(id+1),X(id+2));
Mci=max(Ci);
Ci=Ci./Mci;
%Membership of data (image) to each of the classes
Aci=Mci.*normpdf(I,X(id+1),X(id+2));
if(dim == 2)
TM(:,:,iNmax)=Aci;
else
TM(:,:,:,iNmax)=Aci;
end
end
[Vx Mx]=max(TM,[],dim+1);%AGREGRATION by MAX
%--------------------------------------
if(dim == 2)
id=0;
for i=-1:1:1
M_tmp = shiftmat(Mx,i,1);
for j=-1:1:1
id = id+1;
Mx2(:,:,id) = shiftmat(M_tmp,j,2);
end
end
else
id=0;
for i=-1:1:1
M_tmp = shiftmat(Mx,i,1);
for j=-1:1:1
M_tmp2 = shiftmat(M_tmp,j,2);
for k=-1:1:1
id = id+1;
Mx2(:,:,:,id) = shiftmat(M_tmp2,k,3);
end
end
end
end
MG=median(Mx2,dim+1);
S=Mx;
|
|