%	To perform a quantitation of snr, sfnr, stability and drift
%	includes weisskoff plot  MRM 36:643 (1996)
%
%	rev 0	3/3/00		original from noiseave and imgroi
%	rev 1	3/29/02		fix a few header things
%	rev 2	9/4/02		add weissnoise plot
%	rev 3	1/28/03		add phase drift plot 
%				.freq image is scaled 10x
%	rev 4	4/1/03		for fbirn

%	acq: 35	slice 64x64 TR 3.0 TE 30 (3T) 40 (1.5T) 200 frames
%	grecons -s 18 -e 18 Pxxxxx

more off;
clear roi;
clear roir;

% some defaults

NPIX = 64;
%NPIX = input('npix (cr = 64) = ');	% matrix
if(isempty(NPIX))
  NPIX = 64;
end

I1 = 3403;		% first image
I2 = 3600;		% last image

TR = 3.0;		% rep time, s

if(NPIX == 128)
  R = 30;			% ROI width
else
  R = 15;			% probably 64x64
end
npo2 = NPIX/2;
ro2 = fix(R/2);
X1 = npo2 - ro2;
X2 = X1 + R - 1;
Y1 = X1;
Y2 = X2;
r1 = 1;
r2 = R;

%  set up input and ROI mask

fkern = input('image file path = ','s');
mask = ones(R);
npx = sum(sum(mask));
img = zeros(NPIX);
mag = zeros(NPIX);
str = sprintf('gimme image numbers (cr = [%d %d]) = ', I1, I2);
foo = input(str);
if(~isempty(foo))
  i1 = foo(1); i2 = foo(2);
else
  i1 = I1;
  i2 = I2;
end

N = i2 - i1 + 1;                % num time frames
M = r2 - r1 + 1;                % num ROI's
roir = zeros(N,M);

%TR = input('gimme TR, s = ');

% check for phase images and make base complex image

fname = sprintf('%s.%03dp',fkern, i1);
fid = fopen(fname, 'r');
if(fid > 0) 
  numwin = 4;
  [buf, n] = fread(fid, 'short');
  fclose(fid);
  img(:) = buf;
  phase1 = img*.001;
  base = exp(i*phase1);
  TE = .001*input('gimme TE (ms) = ');
else
  numwin = 3;
end

%  begin loop through images

Iodd = zeros(NPIX*NPIX,1);
Ieven = zeros(NPIX*NPIX,1);
Syy = zeros(NPIX*NPIX,1);
Syt = zeros(NPIX*NPIX,1);
St = 0;
Stt = 0;
S0 = 0;
for j = i1:i2
    fname = sprintf('%s.%03d',fkern, j);
    fprintf('open file  %s\n', fname);
    fid = fopen(fname, 'r');
    I = fread(fid, 'short');
    fclose(fid);
    if(mod(j,2)==1)
      Iodd = Iodd + I;
    else
      Ieven = Ieven + I;
    end
    Syt = Syt + I*j;
    Syy = Syy + I.*I;
    S0 = S0 + 1;
    St = St + j;
    Stt = Stt + j*j;
    img(:) = I;
    sub = img(X1:X2,Y1:Y2);
    roi(S0) = sum(sum(sub))/npx;
    for r = r1:r2                 % each roi size
      ro2 = fix(r/2);
      x1 = npo2 - ro2;
      x2 = x1 + r - 1;
      sub = img(x1:x2,x1:x2);
      roir(j-i1+1, r) = mean(sub(:));
    end;
    if(numwin == 4)		% do the phase
      fname = sprintf('%s.%03dp', fkern, j);
      %fprintf('read %s ...\n', fname);
      fid = fopen(fname, 'r');
      [buf, n] = fread(fid, 'short');
      fclose(fid);
      img(:) = buf;
      phase = img*.001;
      img1 = exp(i*phase);
      z = img1./base; 
      phi = atan2(imag(z), real(z));
      freq = phi/(2*pi*TE);
      sub = freq(X1:X2,Y1:Y2);
      roip(j-i1+1) = mean(sub(:));
    end
end;

% write out last freq drift image

if (numwin==4)
  fname = sprintf('%s.freq', fkern); 
  fout = fopen(fname, 'w');
  fwrite(fout, 10*freq(:), 'short');
  fprintf('\nwrite file %s\n', fname);
  fclose(fout);
end

%  write out diff image

Isub = Iodd - Ieven;
img(:) = Isub;
sub = img(X1:X2,Y1:Y2);
varI = var(sub(:));
fname = sprintf('%s.nave', fkern); 
fout = fopen(fname, 'w');
fwrite(fout, Isub, 'short');
fprintf('\nwrite file %s\n', fname);
fclose(fout);
        
%  write out ave image

Sy = Iodd + Ieven;
Iave = Sy/N;
img(:) = Iave;
sub = img(X1:X2,Y1:Y2);
meanI = mean(sub(:));
fname = sprintf('%s.ave', fkern);
fout = fopen(fname, 'w');
fwrite(fout, Iave, 'short');
fprintf('write file %s\n', fname);
fclose(fout);

% find trend line at + b

D = (Stt*S0 - St*St);
a = (Syt*S0 - St*Sy)/D;
b = (Stt*Sy - St*Syt)/D;

% make sd image

Var = Syy + a.*a*Stt +b.*b*S0 + 2*a.*b*St - 2*a.*Syt - 2*b.*Sy;
Isd = sqrt(Var/(N-1));
fname = sprintf('%s.sd', fkern);
fout = fopen(fname, 'w');
fwrite(fout, 10*Isd, 'short');
fprintf('write file %s\n', fname);
fclose(fout);

% make sfnr image

sfnr = Iave./(Isd + eps);
img(:) = sfnr;
sub = img(X1:X2,Y1:Y2);
sfnrI = mean(sub(:));
fname = sprintf('%s.sfnr', fkern);
fout = fopen(fname, 'w');
fwrite(fout, 10*sfnr, 'short');
fprintf('write file %s\n', fname);
fclose(fout);

snr = meanI/sqrt(varI/N);
fprintf('\nmean, SNR, SFNR = %5.1f  %5.1f  %5.1f\n', meanI, snr, sfnrI);

%  Do fluctation analysis

x=(1:N);
p=polyfit(x,roi,2);
yfit = polyval(p, x);
y = roi - yfit;

subplot(numwin,1,1)
plot(x,roi,x,yfit);
xlabel('frame num');
ylabel('Raw signal');
grid
m=mean(roi);
sd=std(y);
drift = (yfit(N)-yfit(1))/m;
title(sprintf('%s   percent fluct (trend removed), drift= %5.2f %5.2f', fkern, 100*sd/m, 100*drift));

fprintf('std, percent fluc, drift = %5.2f  %6.2f %6.2f \n', sd, 100*sd/m, 100*drift);

z = fft(y);
fs = 1/TR;
nf = N/2+1;
f = 0.5*(1:nf)*fs/nf;
subplot(numwin,1,2);plot(f, abs(z(1:nf)));grid
ylabel('spectrum');
xlabel('frequency, Hz');
ax = axis;

fid = fopen(fkern, 'r');
if(fid ~= -1)
  hdr = fread(fid, 39984/2, 'short');
  R1 = hdr(207)*65536 + hdr(208);	% mps
  TG = hdr(211)*65536 + hdr(212);
  u1 = hdr(214);
  if(u1<0) u1 = u1 + 65536; end
  freq = .1*(hdr(213)*65536 + u1);  
  fprintf('R1, TG, freq = %d  %d  %9.0f\n', R1, TG, freq);
  text(ax(2)*.1, ax(4)*.8, sprintf('mean, SNR, SFNR = %5.1f  %5.1f  %5.1f    R1, TG, freq= %d  %d  %9.0f', meanI, snr, sfnrI, R1, TG, freq));
else
  text(ax(2)*.2, ax(4)*.8, sprintf('mean, SNR, SFNR = %5.1f  %5.1f  %5.1f', meanI, snr, sfnrI));
end

%  now do analysis for each roi size

t = (1:N);
for r = r1:r2
  y = roir(:, r)';
  yfit = polyval(polyfit(t, y, 2), t);  % 2nd order trend
  F(r) = std(y - yfit)/mean(yfit);
end
rr = (r1:r2);
F = 100*F;              % percent
fcalc = F(1)./rr;
rdc = F(1)/F(r2);	% decorrelation distance

subplot(numwin,1,3);
loglog(rr, F, '-x', rr, fcalc, '--');
grid
xlabel('ROI full width, pixels');
ylabel('Relative std, %');
axis([r1 r2 .01 1]);
text(6, 0.5, 'solid: meas   dashed: calc');
text(6, 0.25, sprintf('rdc = %3.1f pixels',rdc));

if(numwin==4)
  subplot(numwin,1,4);
  plot(x,roip)
  xlabel('frame num');
  ylabel('freq drift, Hz');
  grid
end

more on;