Code covered by the BSD License  

Highlights from
Soft thresholding for image segmentation

image thumbnail
from Soft thresholding for image segmentation by SANTIAGO AJA-FERNANDEZ
Image segmentation based on histogram soft thresholding

[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;

Contact us