%mrsnr -- Written June, 2001 by Michael Worden
%This is a script to algorithmically calculate signal-to-noise ratios on sets of 
%MRI images.  The hope was that this method would be less sensitive to variability
%resulting from placement of ROIs by hand onto the images and would prove useful
%for comparing SNR across different sets of images and possibly across different sites.
%USAGE: Set pathname to be a directory containing MR images.  The default behavior
%is to strip off a standard GE header.  If your data have no header or have a header
%with a different size, change the headersize variable accordingly.  The script outputs
%a file in the directory with the images named 'snr' containing statistics for each
%slice individually and calculated on the volume as a whole.  Additionally, two figures
%are generated.  Each figure is a collection of images graphically depicting which
%voxels were considered background, which were considered to be the imaged object, and 
%which had intermediate values and were not included.  One set of images reflects the 
%slice by slice values and the other set of images corresponds to the thresholds
%calculated from the entire volume at once (GT appended to filename).  If these images 
%don't look reasonable, you may want to try changing the 'edgetol' parameter.

pathname = '/d1/3T/mike/211/s2/';
edgetol = .50;     %Edge tolerance to reduce the number of voxels with partial voluming
headersize = 8432; %set to size of header or 0 for no header.  GE is 8432.

cmap = [0 0 0;1 0 0; 0 0 1;0 1 0]; %create color map for visualization later
outfile=[pathname 'snr'] 
outfid=fopen(outfile,'w');
fprintf(outfid,'fname bgthreshold bgmean bgsd objthreshold objmean objsd snr\n');


f1 = figure('Position',[100 100 600 800]);
f2 = figure('Position',[200,100 600 800]);

oldpath = pwd;
cd(pathname);
flist = dir('e*s*i*');   %change this if you want to read in files differently named
nfiles = size(flist,1);

dat = []; alldat = [];
for ifile = 1:nfiles
  fname =flist(ifile).name;
  fid=fopen(fname);
  fseek(fid,headersize,'bof');                            %lop off the header, if any
  [dat,count]=fread(fid,'uint16');
  fclose(fid);
  alldat = [alldat; dat];

  [n,xout]=hist(dat,max(dat));
  outbound = round(mean(dat) + (1*std(dat)));             %boundary for outliers

  %figure out what value occurs the smallest number of times (thresh)
  %note that we stop looking before the very end because you could 
  %get a stray small value there if there is an outlier
  [c,thresh]=min(n(1:outbound));
  phthresh = thresh + round(thresh * edgetol);
  bgthresh = thresh - round(thresh * edgetol);
  datsort=sort(dat);
  bg = datsort(datsort < bgthresh);
  ph = datsort(datsort > phthresh);
  bgmean = mean(bg); bgsd = std(bg);
  phmean = mean(ph); phsd = std(ph);
  snr = phmean/bgsd;

  % output
  fprintf(1,'File = %s\n',fname);
  fprintf(1,'Threshold= %g\n',thresh);
  fprintf(1,'bgmean= %g bgstd= %g phmean= %g phsd= %g snr= %g\n',bgmean,bgsd,phmean,phsd,snr);
  fprintf(outfid,'%s %g %g %g %g %g %g %g\n',fname,bgthresh,bgmean,bgsd,phthresh,phmean,phsd,snr);

  %just some visualizaton for the thresholding
  figure(f1);
  subplot(6,4,ifile);           %assumes there are 24 images
  imgsize = sqrt(size(dat,1));  %assumes a square image
  img=reshape(dat,imgsize,imgsize);
  phimg=img > phthresh;
  phimg=phimg*2;
  bgimg=img < bgthresh;
  timg = bgimg+phimg;
  image(timg+2);
  colormap(cmap); axis equal; axis tight; axis off; title(fname);

end

%This repeated code really ought to be a function call
%We're going to do exactly as above but with the whole volume at once
[n,xout]=hist(alldat,max(alldat));
outbound = round(mean(alldat) + (1*std(alldat)));              %boundary for outliers

[c,thresh]=min(n(1:outbound));
phthresh = thresh + round(thresh * edgetol);
bgthresh = thresh - round(thresh * edgetol);
datsort=sort(alldat);
bg = datsort(datsort < bgthresh);
ph = datsort(datsort > phthresh);
bgmean = mean(bg); bgsd = std(bg);
phmean = mean(ph); phsd = std(ph);
snr = phmean/bgsd;

fprintf(1,'Volume\n');
fprintf(1,'Threshold= %g\n',thresh);
fprintf(1,'bgmean= %g bgstd= %g phmean= %g phsd= %g snr= %g\n',bgmean,bgsd,phmean,phsd,snr);
fprintf(outfid,'VOLUME %g %g %g %g %g %g %g\n',bgthresh,bgmean,bgsd,phthresh,phmean,phsd,snr);

fclose(outfid);
%print; % print the figure

%It would be a good idea to look at each slice thresholded with the volume threshold value
figure(f2);
for ifile = 1:nfiles
  fname =flist(ifile).name;
  fid=fopen(fname);
  fseek(fid,8432,'bof');
  [dat,count]=fread(fid,'uint16');
  fclose(fid);

  subplot(6,4,ifile);
  imgsize = sqrt(size(dat,1));  %assumes a square image
  img=reshape(dat,imgsize,imgsize);
  phimg=img > phthresh;
  phimg=phimg*2;
  bgimg=img < bgthresh;
  timg = bgimg+phimg;
  image(timg+2);
  colormap(cmap); axis equal; axis tight; axis off; title([fname ' GT']);
end

cd(oldpath);  %change back to the directory where we started