% Segmentation of two-tissue images
% Probabilities are not weighted as usual, that is, by the voxel count
% for each class allover the image, but for the voxel count IN THE
% NEIGHBORHOOD of each voxel. 
% Jose A. Rodriguez Serrano 2002
% This time with 4 gaussians...
% Matlab code 6.1.0.450 (R12.1)

fprintf ('\n');
disp('*****************************************');
disp('*   Segmentation of two-tissue images   *');
disp('* with local weighting of probabilities *');
disp('*      This time with 4 gaussians       *');
disp('*    Jose A. Rodriguez-Serrano 2002     *');
disp('*****************************************');
fprintf ('\n');

% Initialization
slice=13;
if exist('slice')==0
  slice=0
end  
paramsfilename = '/home/josear/scratch/gaussian4.dat';
directory = '/home/josear/images/verhoeven';
outputfileprefix = '/home/josear/results/segmentation/muscle_verhoeven';
filename = sprintf('%s%d%s', directory, slice,'.w')
outputfilename = sprintf('%s%d%s', outputfileprefix, slice,'.w');

WhatDoISegment=2 % 1=fat 2=muscle
DoIWriteTheResult=1;
DoIApplyAThreshold=1;
DoINormalize=0;
NormalizationDepth=255;  

lo=800;
hi=1170;
window=3; % Must be an even number!
miniside=(window-1)/2;  

  % Read parameters from input file
  fid = fopen(paramsfilename, 'rt');
  params=fscanf(fid, '%f');
  fclose (fid);
  mu1=params(1);
  mu2=params (2);
  mu3=params(3);
  mu4=params(4);
  sigma1=params(5);
  sigma2=params(6);
  sigma3=params(7);
  sigma4=params(8);
  p1=params(9);
  p2=params(10);
  p3=params(11);
  p4=params(12);
  disp ('Initialised with');
  initialization=[mu1 mu2 mu3 mu4 sigma1 sigma2 sigma3 sigma4 p1 p2 ...
        p3 p4]
  clear params initialization

  % Open the image
  fid=fopen(filename, 'r');
  A=fread(fid, 'short');
  fclose (fid);
  clear fid;
  side=sqrt(size(A));
  side=side(1);
  A=reshape(A,[side, side]);
  P=zeros(side); % Matrix which will store the probabilities    
  
% Classify: Compute the probabilities of belonging to each tissue
for i=1:side,
  for j=1:side,    
    x=A(i,j);      
    gauss1=p1/(sqrt(2*pi)*sigma1)*exp(-(x-mu1)*(x-mu1)/(2*sigma1*sigma1));
    gauss2=p2/(sqrt(2*pi)*sigma2)*exp(-(x-mu2)*(x-mu2)/(2*sigma2*sigma2));
    gauss3=p3/(sqrt(2*pi)*sigma3)*exp(-(x-mu3)*(x-mu3)/(2*sigma3*sigma1));
    gauss4=p4/(sqrt(2*pi)*sigma4)*exp(-(x-mu4)*(x-mu4)/(2*sigma4*sigma2));    
    totalgauss=gauss1+gauss2+gauss3+gauss4;
    if totalgauss ~= 0
      if WhatDoISegment == 1, 
        P(i,j)=(gauss1+gauss2)/totalgauss; 
      else 
        P(i,j)=(gauss3+gauss4)/totalgauss;
      end  
    end      
  end
end  
clear x gauss1 gauss2 gauss3 gauss4 totalgauss WhatDoISegment
disp('Probabilities computed');

% Weight these probabilities with the voxel count for each class BUT
% ONLY in the 3x3 neighborhood
Q=zeros(side);
for i=1+miniside:side-miniside,
  for j=1+miniside:side-miniside,
    if (A(i,j)>lo) & (A(i,j)<hi)
      weight=0;
      % Count fat and muscle voxels in the neighborhood
      tissue1=0;
      tissue2=0;    
      for k=-miniside:miniside
        for l=-miniside:miniside
          distance = sqrt(l*l+k*k);
          if distance~=0
            if P(i+k,j+l)>0.5
              tissue1=tissue1+1/distance;
            else  
              tissue2=tissue2+1/distance;
            end %if                
          end %if 
        end                   
      end  
      % Weight!
      tissue1=tissue1/8;
      tissue2=tissue2/8;     
      total=P(i,j)*tissue1+(1-P(i,j))*tissue2;
      if total~=0
        Q(i,j)=P(i,j)*tissue1/total;        
      end  
    else 
      Q(i,j)=0; % Mask effect!! (It also makes the process faster)       
    end %if (A(i,j)..)  
  end  
end  
clear tissue1 tissue2 total A 
disp('Classification done');

% Threshold
if DoIApplyAThreshold==1  
  for i=1:side,
    for j=1:side,
      if Q(i,j)>0.5 
        Q(i,j)=1;
      else  
        Q(i,j)=0;
      end  
    end  
  end  
  disp ('Threshold done');
else  
  disp ('Threshold not applied'); 
end

% Normalize matrix Q to some new integer (!) range
if DoINormalize == 1
  Q=round(Q*NormalizationDepth);  
end  

% Write the result of the threshold in an image, if desired
if DoIWriteTheResult == 1
  fid=fopen(outputfilename, 'w');
  fwrite (fid, Q, 'short');
  fclose (fid);
  fprintf ('Result written in %s\n', outputfilename);
else
  disp ('Result not written');
end  
clear paramsfilename filename outputfilename directory ...
    outputfileprefix DoIWriteTheResult DoINormalize DoIApplyThreshold ...
    NormalizationDepth lo hi window miniside

imshow(mat2gray(Q));