Changes to CEST_EVAL, mainly for Bruker compatibility

BIG_BATCH2_cest-sources.m

@@ -22,6 +22,7 @@
 clear all; close all; clc
 
 %% LOAD CEST-DATA
+% for Philips or Bruker data, use 'LOAD_BATCH_cest-sources' script
 [M0_stack, Mz_stack, P] = LOAD('USER');
 % dimensions: M0_stack(x,y,z) or M0_stack(x,y,z,w) ;  Mz_stack(x,y,z,w) ;
 % make sure they are double; offets (deltaomega in ppm) are stored in P.SEQ.w
@@ -48,10 +49,6 @@
 
 [Z_corrExt] = NORM_ZSTACK(Mz_CORR,M0_stack,P,Segment);
 
-
-%% save
-save matlab.mat ;
-
 %% start imgui
 close(imgui); imgui
 
@@ -70,7 +67,7 @@
 p0 = [ 1    0.9   1.4   0       0.025  0.5   3.5     0.02   3    -3.5    0.1   25   -1     0.01  1.5  2.2 ];
 P.FIT.lower_limit_fit = lb; P.FIT.upper_limit_fit = ub; P.FIT.start_fit = p0;
 
-Segment= make_Segment(M0_stack, 'free', mean(M0_stack(M0_stack>0)).*[0.3]); % choose smalle ROI for testing
+% Segment= make_Segment(M0_stack, 'free', mean(M0_stack(M0_stack>0)).*[0.3]); % choose smalle ROI for testing
 close all
 
 tic
@@ -79,13 +76,19 @@
 
 [Zlab, Zref] = get_FIT_LABREF(popt,P,Segment);  % create Reference values Z_Ref=(Z_lab - Li)
 
-imgui
+close(imgui); imgui
+
+%% save
+save matlab.mat ;
 
 %% WASSR2/WASAB1 EVAL
+% Warning: will overwrite variables Mz_stack, P, Segment, etc.
+
+% for Philips or Bruker data, use 'LOAD_BATCH_cest-sources' script
 [M0_stack, Mz_stack, P] = LOAD('USER');
 
 %%
-Segment= make_Segment(M0_stack, 'free', mean(M0_stack(M0_stack>0)).*[0.20]); 
+% Segment= make_Segment(M0_stack, 'free', mean(M0_stack(M0_stack>0)).*[0.20]); 
 [Z_uncorr] = NORM_ZSTACK(Mz_stack,M0_stack,P,Segment);
 
 tic
@@ -97,10 +100,10 @@
 toc
 
 B1map=popt(:,:,1)/P.SEQ.B1;
- dB0_stack_ext=popt(:,:,2);
+dB0_stack_ext=popt(:,:,2);
 
- figure, subplot(1,2,1), imagesc(dB0_stack_ext(:,:,1),[-0.9 0.9]); title('\DeltaB0 map in ppm');
-         subplot(1,2,2), imagesc(B1map,[0.6 1.4]); title('relative B1 map');
+figure, subplot(1,2,1), imagesc(dB0_stack_ext(:,:,1),[-0.9 0.9]); title('\DeltaB0 map in ppm'); colorbar;
+        subplot(1,2,2), imagesc(B1map,[0.6 1.4]); title('relative B1 map'); colorbar;
 %%
 save B0_B1.mat ;
 
@@ -113,32 +116,51 @@
 Z_stack(:,:,:,:,n)=Z_corrExt; %% do this for all B1 measurements
 n=n+1;
 
+% for Bruker data (where inputs are sequence descriptions of B1 measurements):
+Z_stack = make_5D_B1_stack(protocol, directory_M0, M0_stack, Segment, 'example_B1seqDescr_1', 'example_B1seqDesc_2');
+
 %% B1 correction step2: run correction -rqures relative B1_map
 tic % etwa 100s
+
 B1_input=[0.3 0.6 0.8]; % give the nominal B1 values set at the scanner
+% or, for Bruker:
+% B1_input = get_protocol_B1_values(directory, protocol, 'example_B1seqDescr_1', 'example_B1seqDesc_2');
+
 B1_output=[0.3 0.4 0.6 0.7 0.8]; % choose which values you want to reconstruct, e.g. B1_output=[ 1 2 ], or B1_output=B1_input
+
 [Z_stack_corr] = Z_B1_correction(Z_stack,B1map,B1_input,B1_output,Segment,'linear');
 % [Z_stack_corr] = Z_B1_correction(Z_stack,B1_map,B1_input,B1_output,Segment);
 toc
 
-Z_corrExt=Z_stack_corr(:,:,:,:,2); % pick second reconstructed value (e.g. 2 for B1_output(2)=0.4�T) 
+Z_corrExt=Z_stack_corr(:,:,:,:,2); % pick second reconstructed value (e.g. 2 for B1_output(2)=0.4?T) 
 
-imgui
+% See WASABI spectrum (choose map 'Z_uncorr')
+close(imgui); imgui
 
 
 %% T1 mapping
 TI=[100 200 400 600 800 1000 1300 1600 2000 2500 3000 3500 4000 4500 5000 10000 15000];
 
-% TI=[20:20:100  200 400 600 800 1000 1200 1400 1600 1800 2000 2500 3000 3500 4000 4500 5000 10000 ];
+% for Bruker: make folder with T1 image DICOM files and change DICOM header
+ix_T1 = 7:12;   % protocol entry numbers of T1 images
+T1dirpath = make_T1_directory(directory, protocol, TI, ix_T1)
 
 % information about fit
 P_T1.FIT.options   = [1E-04, 1E-15, 1E-10, 1E-04, 1E-06];
 P_T1.FIT.nIter     = 100;
 P_T1.FIT.modelnum  = 031011;
-P_T1.SEQ.w = TI;
-% number of ROIS, mapflag (should complete T1map be calculated), P_T1 struct, Segment
-
-[T1info T1map popt_T1] = T1eval_levmar(1,2,P_T1);
+P_T1.SEQ.w = TI';
+
+% starting parameters (optional)
+%     T1          a      c
+lb = [0         -5000   0       ];
+ub = [50000      5000   5000    ];
+p0 = [1000      -2000   1000    ];
+P_T1.FIT.lower_limit_fit = lb; P_T1.FIT.upper_limit_fit = ub; P_T1.FIT.start_fit = p0;
+StartValues=p0;
+
+% mapflag (should complete T1map be calculated), number of ROIS, P_T1 struct, Segment
+[T1info T1map popt_T1] = T1eval_levmar(1,1,P_T1,Segment,StartValues);
 
 
 

LOAD_BATCH_cest-sources.m

@@ -27,10 +27,19 @@
 P.EVAL.lowerlim_slices=1;
 P.EVAL.upperlim_slices=size(M0_stack,3);
 clearvars -except  P Mz_stack M0_stack image_z x X ROI_def Segment dB0_stack_ext dB0_stack_int
-
-
-
-
-
-
 
+%% LOAD Bruker CEST-DATA
+% Given the patient directory, makes a 'protocol' structure containing
+% sequence names/descriptions (set by you at the scanner) and their
+% respective directory name (the sequence number, according to Bruker's
+% "E[n]" naming convention). Make sure sequence descriptions are
+% informative.
+
+directory = uigetdir;
+protocol = readprotocol(directory);
+Mz_name = 'example_sequence_description';
+M0_name = 'example_sequence_description';
+[directory_Mz, directory_M0] = getfolderpath(directory, protocol, Mz_name, M0_name);
+Mz_stack = load_Mz(directory_Mz);
+M0_stack = load_M0(directory_M0);
+P = wipread_modified(directory_Mz, directory_M0);   % writes all parameters into P structure

T1_evaluation/T1eval_levmar.m

@@ -22,14 +22,18 @@
 % read images from folder
 cd(uigetdir)
 listoffiles=dir('*.IMA');
-numfiles=size(listoffiles,1);
+% added CT 20161212
+if isempty(listoffiles)
+    listoffiles=dir('*.dcm');
+end
+numfiles=numel(listoffiles);
 
 for i=1:numfiles
     filect=listoffiles(i).name;
     dicom_info=dicominfo(filect);
     info.Instance(i,1)=dicom_info.InstanceNumber;
     info.Acquisition(i,1)=dicom_info.AcquisitionNumber;
-    image(:,:,info.Instance(i),info.Acquisition(i))=dicomread(listoffiles(i).name);
+    image(:,:,info.Instance(i),info.Acquisition(i))=double(dicomread(listoffiles(i).name));
 end
 clear i filect listoffiles ans
 fclose('all');
@@ -83,7 +87,7 @@
         selected_slice=1;
     end
 
-    ROIS = ROItool(squeeze(image(:,:,selected_slice,end)),ROInumber,'ellipse');
+    ROIS = ROItool(squeeze(image(:,:,selected_slice,end)),ROInumber,'free');
 
      % preview plot
     if mapflag

T1_evaluation/T1eval_levmar.m~

@@ -0,0 +1,178 @@
+function [T1info T1map popt] = T1eval_levmar(mapflag,ROInumber,P,Segment,StartValues)
+% T1eval_levmar(mapflag,ROInumber,TI,Segment,StartValues)
+% mapflag: if 1 fits full given image (inside Segment)
+% ROInumber: number of ROIs fitted
+% P : P struct
+
+
+if nargin < 5
+    StartValues = false;
+    sprintf('No StartValues were given')
+else
+    sprintf('Benutze folgende Startwerte:')
+    StartValues
+end
+
+try
+    TI=P.SEQ.w;
+catch
+    error('No inversion times given!');
+end
+
+% read images from folder
+cd(uigetdir)
+listoffiles=dir('*.IMA');
+% added CT 20161212
+if isempty(listoffiles)
+    listoffiles=dir('*.dcm');
+end
+numfiles=numel(listoffiles);
+
+for i=1:numfiles
+    filect=listoffiles(i).name;
+    dicom_info=dicominfo(filect);
+    info.Instance(i,1)=dicom_info.InstanceNumber;
+    info.Acquisition(i,1)=dicom_info.AcquisitionNumber;
+    image(:,:,info.Instance(i),info.Acquisition(i))=double(dicomread(listoffiles(i).name));
+end
+clear i filect listoffiles ans
+fclose('all');
+
+if (nargin < 4 && mapflag==1)
+    sprintf('No Segment was given to function T1eval')
+    Segment=make_Segment(image,'free');
+    assignin('base','Segment',Segment);
+end
+
+% check if number of images is equal to number of inversion times
+if not(dicom_info.AcquisitionNumber==numel(TI))
+    error('Number of inversion times is different to number of files');
+end
+
+
+
+
+if mapflag
+
+    [popt, P] = FIT_3D(image,P,Segment,StartValues);
+    T1map=popt(:,:,:,1);
+
+    figure,
+    for ii=1:size(image,3)
+        subplot(1,size(image,3),ii);
+        imagesc(T1map(:,:,ii),[0 5000]);
+        axis image
+        set(gca,'xtick',[])
+        set(gca,'ytick',[])
+        colormap jet;
+        if (ii == size(image,3))
+            colorbar
+        end
+    end
+
+else
+    T1map=1;
+    popt=0;
+end;
+
+% ROI eval
+if (ROInumber~=0)
+    % define ROI_def for all ROIs
+    if (size(image,3)>1)
+        prompt = {'Select slice for ROI analysis:'};
+        dlg_title = 'Input';
+        answer = inputdlg(prompt,dlg_title,1,{'1'});    
+        selected_slice=abs(str2num(answer{1}));
+    else
+        selected_slice=1;
+    end
+
+    ROIS = ROItool(squeeze(image(:,:,selected_slice,end)),ROInumber,'free');
+
+     % preview plot
+    if mapflag
+        if (ndims(T1map)==2)
+            im = T1map;
+        else
+            im = T1map(:,:,selected_slice);
+        end
+    else
+        im=squeeze(image(:,:,selected_slice,end));
+    end
+
+    im=double(im);
+    im(isnan(im))=0;
+
+    figure;
+    subplot(1,2,1);
+    imagesc(im, [0.1*mean(im(im~=0)) 2*mean(im(im~=0))]), axis image;
+    title('ROI definition')
+
+    for ii=1:ROInumber
+        % get ROI_def for current ROI
+        ROI_def=ROIS{ii}.ROI_def;
+        % fit every pixel in ROI
+        [ROI_popt, P] = FIT_3D(image(:,:,selected_slice,:),P,ROI_def,StartValues);
+        T1mapROI=squeeze(ROI_popt(:,:,:,1));
+
+        % write values (mean,std,...) for current ROI
+        T1info{ii}.mean = mean(T1mapROI(ROI_def==1));
+        T1info{ii}.std = std(T1mapROI(ROI_def==1));
+        T1info{ii}.img = T1mapROI;
+        T1info{ii}.ROI_def = ROI_def;
+
+        % calculate mean values of fit-result for preview plot
+        for kk=1:size(ROI_popt,4);
+            temp=ROI_popt(:,:,1,kk);
+            ROI_popt_mean(kk)=mean(temp(ROI_def==1));
+        end
+
+        % calculate mean signal and std for all TI in current ROI
+        for jj=1:numel(TI)
+            buffer=squeeze(image(:,:,selected_slice,jj));
+            ROI_mean(jj)= mean(buffer(ROI_def == 1));
+            ROI_std(jj)= std(double(buffer(ROI_def == 1)));
+        end;
+
+        % calculate y_fit values for current ROI
+        [f] = fitmodelfunc_ANA(ROI_popt_mean,P);
+
+        % save all ROI results into S-struct
+        S{ii}.y = ROI_mean;
+        S{ii}.y_std = ROI_std;
+        S{ii}.popt = ROI_popt_mean;
+        S{ii}.T1 = T1info{ii}.mean;
+        S{ii}.dT1 = T1info{ii}.std;
+        S{ii}.y_fit = f(TI);
+
+        TI_inter=[min(TI):1:max(TI)];
+
+        subplot(1,2,1);
+        hold on;
+        contour(T1info{ii}.ROI_def,1,'m-','LineWidth',2);
+        [yy xx]=find(T1info{ii}.ROI_def);
+        t_h=text(fix(max(xx)),fix(min(yy)),'# ');
+        set(t_h,'String',sprintf('ROI # %d',ii),'BackgroundColor',[1 1 1])
+        leg1{ii}=sprintf('# ROI %d',ii);
+
+        subplot(1,2,2)
+        hold on
+        leg2{ii}=sprintf('ROI%i: T1= %4.0f +- %4.0f ms',ii,S{ii}.T1,S{ii}.dT1);
+        %leg{numel(leg)+1}=sprintf('ROI%i: data',ii);
+        plot(TI_inter,f(TI_inter),'Color',cl(ii,ROInumber));
+        hold on;
+        h2(ii)=plot(TI,ROI_mean,':o','Color',cl(ii,ROInumber));
+
+    end;
+%    subplot(1,2,1)
+%    legend(leg1,'location','SouthEast'); 
+   subplot(1,2,2)
+   legend(h2,leg2,'location','East'); 
+   box on 
+
+
+else
+    T1info = 1;
+end
+
+

bruker/get_protocol_B1_values.m

@@ -0,0 +1,20 @@
+function B1_vector = get_protocol_B1_values(directory, protocol, varargin)
+% ** function B1_vector = get_protocol_B1_values(directory, protocol, sequenceDescription1, ... sequenceDescriptionN)
+%
+% Returns B1 values for the specified sequence(s).
+% Input is one or more sequence description strings.
+%     e.g. 'protocol(3:6).SequenceDescription' for protocol entries 3-6.
+% Output is a vector of B1 field strength (in microTesla) of each sequence.
+%
+% CT 20170111
+
+if nargin<3
+    error('Must specify one or more input strings')
+end
+
+B1_vector = nan(1, nargin-2);
+for i=1:nargin-2
+    seqdirec = getfolderpath(directory, protocol, varargin{i});
+    methodpath = [seqdirec, '/method'];
+    B1_vector(i) = mineExecFile(methodpath, '##$PVM_MagTransPower');
+end