qualicheck_7T_20071219_ER.m

%	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
I1 = 1;		% first image
I2 = 500;	% last image
numwin = 4;

%Datenpfad festlegen
%[fdummy,fpath] = uigetfile('*.IMA');
fpath = '/scr/mrincoming7t/QA/QA_180104_32CH.QA/S2_EPI_time_series_500Rep_NEW_position/';
fdummy = 'QA_180104_32CH.MR.TA_QA.0002.0001.2018.01.04.13.58.41.578125.19504220.IMA';
proband=strrep(fdummy,'.', ' ');
fnummern=sscanf(proband,'%*s %*s %*s %d %d %d %d %d %d %d %d %d %d IMA');
proband=sscanf(proband,'%s %*d %*d %*d.IMA');

fname=sprintf('%s%s',fpath,fdummy);
header=dicominfo(fname);

NPIX = header.AcquisitionMatrix(1);  % No. of Pix Read
TR = header.RepetitionTime/1000;		% rep time, s

if(NPIX == 128)     % H? ??
  R = 30;			% ROI width
else
  R = 15;			% probably 64x64
end
r1 = 1;
r2 = R;
%  set up input and ROI mask
 mask = ones(R);
 npx = sum(sum(mask));
 img = zeros(header.AcquisitionMatrix(4),header.AcquisitionMatrix(1),1);
    
% Mosaic Dimensionen ********Nur wenn phase H/F?********
Mosaic_rows = header.Height / header.AcquisitionMatrix(4);
Mosaic_columns = header.Width / header.AcquisitionMatrix(1);



% Anzahl Rep.
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

% Schichtanzahl
slices = 20;
str = sprintf('gimme No. of slices (cr = %d) = ', slices);
foo = input(str);
if(~isempty(foo))
  slices = foo(1);
end


% Schichtauswahl f???r eval
slice_eval_start = 1; slice_eval_stop = slices;
str = sprintf('gimme slice No. for eval. (cr =[%d %d]) = ', slice_eval_start, slice_eval_stop);
foo = input(str);
if(~isempty(foo))
  slice_eval_start = foo(1); slice_eval_stop = foo(2);
end

N = i2 - i1 + 1;                % num time frames
M = r2 - r1 + 1;                % num ROI's


%
% Auswertung f???r alle alsgew???hlten Schichten
%
for slice_eval = slice_eval_start:slice_eval_stop
    
    %Einzelbildkoordinaten im Mosaic
    imgX = fix((slice_eval -1) / double(Mosaic_columns)) * header.AcquisitionMatrix(4) + 1;
    imgY = mod(slice_eval - 1,Mosaic_columns) * header.AcquisitionMatrix(1) +1;
        
% init ROIs    
    roir = zeros(N,M);

%  begin loop through images

Iodd  = zeros(header.AcquisitionMatrix(4),header.AcquisitionMatrix(1),1); % Nur f???r ein Einzelbild im Mosaic
Ieven = zeros(header.AcquisitionMatrix(4),header.AcquisitionMatrix(1),1);
Syy   = zeros(header.AcquisitionMatrix(4),header.AcquisitionMatrix(1),1); 
Syt   = zeros(header.AcquisitionMatrix(4),header.AcquisitionMatrix(1),1);
St = 0;
Stt = 0;
S0 = 0;

% Phantomposition im ersten Bild finden
%    fname = sprintf('%s%s -%04d-%04d-%04d.dcm',fpath,proband,fnummern(1),i1,(i1-1)*slices+1);
fname = deblank(ls(sprintf('%s%s.*.%04d.%04d.%04d.%02d.%02d.%02d.%02d.%02d.*',fpath,proband,fnummern(1),i1,fnummern(3),fnummern(4),fnummern(5),fnummern(6),fnummern(7),fnummern(8))));
dummy = double(dicomread(fname));
    FitI = dummy(imgX:imgX+header.AcquisitionMatrix(4)-1,imgY:imgY+header.AcquisitionMatrix(1)-1)/4096; % Nur Einzelbild skaliert auf 0-1
    threshold = graythresh(FitI);
    BW = im2bw(FitI,threshold);
    dim = size(BW);
%    col = 86;
    col = 80;
    row = min(find(BW(:,col)));
    connectivity = 8;
    num_points   = 180;
    contour = bwtraceboundary(BW, [row, col], 'N', connectivity, num_points);
    x = contour(:,2);
    y = contour(:,1);

    % solve for parameters a, b, and c in the least-squares sense by
    % using the backslash operator
    abc=[x y ones(length(x),1)]\[-(x.^2+y.^2)];
    a = abc(1); b = abc(2); c = abc(3);

    % calculate the location of the center and the radius
    xc = -a/2;
    yc = -b/2;
    radius  =  sqrt((xc^2+yc^2)-c);
    
    xc=round(xc);
    yc=round(yc);
      
    % ROI Koordinaten
     X1 = xc - fix(R/2);
     X2 = X1 + R -1; 
     Y1 = yc - fix(R/2); 
     Y2 = Y1 + R - 1; 
    
for j = i1:i2
    %fname = sprintf('%s%s -%04d-%04d-%04d.dcm',fpath,proband,fnummern(1),j,(j-1)*slices+1);
    fname = deblank(ls(sprintf('%s%s.*.%04d.%04d.%04d.%02d.%02d.%02d.%02d.%02d.*',fpath,proband,fnummern(1),j,fnummern(3),fnummern(4),fnummern(5),fnummern(6),fnummern(7),fnummern(8))));
    %fname = sprintf('%s%s-%04d-%04d-%05d.dcm',fpath,proband,fnummern(1),j,j);
    dummy = double(dicomread(fname));
    I = dummy(imgX:imgX+header.AcquisitionMatrix(4)-1,imgY:imgY+header.AcquisitionMatrix(1)-1); % Nur Einzelbild
    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;
    sub = I(X1:X2,Y1:Y2);
    roi(S0) = sum(sum(sub))/npx;
    for r = r1:r2                 % each roi size
      x1 = xc - fix(r/2);   
      x2 = x1 + r - 1; % Reihe ende unten
      y1 = yc - fix(r/2); 
      y2 = y1 + r - 1; % Spalte ende rechts
      sub = I(x1:x2,y1:y2);
      roir(j-i1+1, r) = mean(sub(:));
    end;
end;

%  write out diff image

Isub = Iodd - Ieven;
sub = Isub(X1:X2,Y1:Y2);
varI = var(sub(:));
%fname=sprintf('%s%s-%04d-%02d-%s.nave','/tmp/',proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate);
fname=sprintf('%s%s-%04d-%02d-%s.nave','/scr/erhard1/home_overflow/data/7t/stability/tmp/',proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate);
%fname = sprintf('%s%s.MR.QA_8CHANNEL.%d.%d.%04d.%02d.%02d.%d.%d.%d.*',fpath,proband,fnummern(1),slice_eval,fnummern(3),fnummern(4),fnummern(5),fnummern(6),fnummern(7),fnummern(8))))

%fout = fopen(fname, 'w');
%fwrite(fout, Isub, 'short');
%fprintf('\nwrite file %s\n', fname);
%fclose(fout);
        
%  write out ave image
dpath = '/scr/erhard1/home_overflow/data/7t/stability/tmp/'
Sy = Iodd + Ieven; % MosaicBild
Iave = Sy/N;
sub = Iave(X1:X2,Y1:Y2);
meanI = mean(sub(:));
fname = sprintf('%s%s-%04d-%02d-%s.ave',dpath,proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate)
fout = fopen(fname, 'w');
fwrite(fout, Iave, 'short');
fprintf('write file %s\n', fname);
fclose(fout);

% find trend line at + b Einzelbild

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%s-%04d-%02d-%s.sd',dpath,proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate);
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 = sfnr(X1:X2,Y1:Y2);
sfnrI = mean(sub(:));
fname = sprintf('%s%s-%04d-%02d-%s.sfnr',dpath,proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate);
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);

%
figure (slice_eval);

%  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-%d-slice%02d-%s   percent fluct (trend removed), drift= %5.2f %5.2f', proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate, 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;

text(ax(2)*.2, ax(4)*.8, sprintf('mean, SNR, SFNR = %5.1f  %5.1f  %5.1f', meanI, snr, sfnrI));


%  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

% write log file
logname = sprintf('%s%s-%04d-%s.log',dpath,proband,fnummern(1),header.PerformedProcedureStepStartDate);
logout = fopen(logname, 'a');
fprintf(logout,'Schicht %04d %6.2f %6.2f %5.1f  %5.1f  %5.1f %3.1f \n',slice_eval, 100*sd/m, 100*drift, meanI, snr, sfnrI,rdc );
fclose(logout);

% plot
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));

subplot(numwin,1,4);
imshow(FitI);
hold on;
plot(contour(:,2),contour(:,1),'g','LineWidth',1);
% display the calculated center
plot(xc,yc,'yx','LineWidth',2);
   
% plot the entire circle
theta = 0:0.01:2*pi;
% use parametric representation of the circle to obtain coordinates
% of points on the circle
Xfit = radius*cos(theta) + xc;
Yfit = radius*sin(theta) + yc;

plot(Xfit, Yfit,'LineWidth',2);
        
% plot the ROI
plot(xc-fix(R/2):xc+fix(R/2)-1,yc-fix(R/2):yc-fix(R/2),'y','Linewidth',2);
plot(xc+fix(R/2)-1:xc+fix(R/2)-1,yc-fix(R/2):yc+fix(R/2)-1,'y','Linewidth',2);
plot(xc-fix(R/2):xc+fix(R/2)-1,yc+fix(R/2)-1:yc+fix(R/2)-1,'y','Linewidth',2);
plot(xc-fix(R/2):xc-fix(R/2),yc-fix(R/2):yc+fix(R/2)-1,'y','Linewidth',2); 
hold off;

fname=sprintf('%s%s-%04d-%02d-%s.jpg',dpath,proband,fnummern(1),slice_eval,header.PerformedProcedureStepStartDate);
print ('-djpeg',fname);
end % slice_eval 
more on;