diff --git a/BIG_BATCH2_cest-sources.m b/BIG_BATCH2_cest-sources.m index 7b0a1ea..63dd7e3 100644 --- a/BIG_BATCH2_cest-sources.m +++ b/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); diff --git a/LOAD_BATCH_cest-sources.m b/LOAD_BATCH_cest-sources.m index 9fa433b..71ad90b 100644 --- a/LOAD_BATCH_cest-sources.m +++ b/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 \ No newline at end of file diff --git a/T1_evaluation/T1eval_levmar.m b/T1_evaluation/T1eval_levmar.m old mode 100644 new mode 100755 index 47d2d5d..748a673 --- a/T1_evaluation/T1eval_levmar.m +++ b/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 diff --git a/T1_evaluation/T1eval_levmar.m~ b/T1_evaluation/T1eval_levmar.m~ new file mode 100644 index 0000000..748a673 --- /dev/null +++ b/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 + + diff --git a/bruker/get_protocol_B1_values.m b/bruker/get_protocol_B1_values.m new file mode 100755 index 0000000..7099d3c --- /dev/null +++ b/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 \ No newline at end of file diff --git a/bruker/getfolderpath.m b/bruker/getfolderpath.m new file mode 100755 index 0000000..06c63b3 --- /dev/null +++ b/bruker/getfolderpath.m @@ -0,0 +1,28 @@ +function varargout = getfolderpath(directory, protocol, varargin) +% ** function [path1, ..., pathN] = getfolderpath(directory, protocol, SequenceDescription1, ..., SequenceDescriptionN) +% +% Retrieves path of a sequence's directory. A sequence is identified by +% its 'SequenceDescription' string, stored in protocol structure 'protocol' +% created by function 'readprotocol'. +% Input is one or more 'SequenceDescription' strings. +% Output is one or more 'path' strings. +% +% CT 20170113 + +if nargin<3 + error('Must specify one or more input strings') +end + +nseq = nargin-2; +if nargout~=nseq + error('Number of outputs must be equal to number of input strings.') +end + +for n=1:nseq + ix = find(strcmp({protocol.SequenceDescription}, varargin{n})); + if ~isempty(ix) + varargout{n} = fullfile(directory, protocol(ix).FolderNum); + else + error(['No folder number is associated with sequence description ''', varargin{n}, '''.']); + end +end \ No newline at end of file diff --git a/bruker/load_M0.m b/bruker/load_M0.m new file mode 100755 index 0000000..a1dfe19 --- /dev/null +++ b/bruker/load_M0.m @@ -0,0 +1,73 @@ +function M0_stack = load_M0(directory_M0) +% ** function M0_stack = load_M0(directory_M0) + +cd(directory_M0); + +% read acqp +file_acqp = textread('method','%s','delimiter','=','whitespace',''); +size_file_acqp = size(file_acqp); +for cont_acqp = 1:size_file_acqp(1) + if strcmp(file_acqp(cont_acqp),'##$PVM_NRepetitions')==1 + num_rip = round(str2num(char(file_acqp(cont_acqp+1)))); + end +end + +% read reco +file_reco = textread('pdata/1/reco','%s','delimiter','\n','whitespace',''); +size_reco = size(file_reco); +for cont_reco = 1:size_reco(1) + if strcmp(file_reco(cont_reco),'##$RECO_size=( 2 )') + size_xy = str2num(char(file_reco(cont_reco+1))); + size_x_res = size_xy(1); + size_y_res = size_xy(2); + break + end +end + +file_meta=textread('pdata/1/reco','%s','delimiter','=','whitespace',''); +size_meta=size(file_meta); + +for cont_meta=1:size_meta(1) + if strcmp(file_meta(cont_meta),'##$RECO_wordtype') + reco_word=(char(file_meta(cont_meta+1))); + end +end + +if strcmp(reco_word,'_32BIT_SGN_INT') + reco_bit='uint32'; +elseif strcmp(reco_word,'_16BIT_SGN_INT') + reco_bit='uint16'; +end + +file_reco_slope=textread('pdata/1/reco','%s','delimiter','=','whitespace',''); + +for cont_reco=1:size(file_reco_slope,1) + if strcmp(file_reco_slope(cont_reco),'##$RECO_map_slope') + slope=str2num(char(file_reco_slope(cont_reco+2))); % correction factor + break + end +end + +size_zx=size_x_res; +size_zy=size_y_res; +size_zz=num_rip; + +M0_stack = zeros(size_zx,size_zy,1,1); + +id_read = waitbar(0,'Reading 2dseq M0...'); +fid=fopen('pdata/1/2dseq'); + +image_M0 = zeros(size_zx,size_zy,1,size_zz); +for cont_z=1:size_zz + waitbar(cont_z/size_zz); + for cont_i=1:size_zx + for cont_j=1:size_zy + image_M0(cont_i,cont_j,1,cont_z)=fread(fid,1,reco_bit); + end + end +end +fclose(fid); +close(id_read); + +image_M0 = image_M0/slope(1); % scale image by correction factor +M0_stack = mean(image_M0(:,:,1,:),4); \ No newline at end of file diff --git a/bruker/load_Mz.m b/bruker/load_Mz.m new file mode 100755 index 0000000..f83e03d --- /dev/null +++ b/bruker/load_Mz.m @@ -0,0 +1,116 @@ +function Mz_stack = load_Mz(directory_Mz) +% ** function Mz_stack = load_Mz(directory_Mz) + +cd(directory_Mz); + +num_ppm=0; +cont_ppm=1; +x=[]; +file_acqp=textread('acqp','%s','delimiter','=','whitespace',''); +size_file_acqp=size(file_acqp); +for cont_acqp=1:size_file_acqp(1) + if strcmp(file_acqp(cont_acqp),'##$ACQ_institution')==1 + institution=(char(file_acqp(cont_acqp+2))); + end + if strcmp(file_acqp(cont_acqp),'##$ACQ_station')==1 + scanner=(char(file_acqp(cont_acqp+2))); + end + if strcmp(file_acqp(cont_acqp),'##$BF1')==1 + B0_field=round(str2num(char(file_acqp(cont_acqp+1)))); + B0_field_nonround=str2num(char(file_acqp(cont_acqp+1))); % CT 20161205 + end + if strcmp(file_acqp(cont_acqp),'##$ACQ_O2_list_size')==1 + num_ppm=round(str2num(char(file_acqp(cont_acqp+1)))); + end + if strcmp(file_acqp(cont_acqp),'##$ACQ_O2_list')==1 + cont_acqp=cont_acqp+1; + + offset=(str2num(char(file_acqp(cont_acqp+1)))/B0_field_nonround); + x=horzcat(x,offset); + punti = size(x); + cont_acqp=cont_acqp+1; + + while punti(2) < num_ppm + offset=(str2num(char(file_acqp(cont_acqp+1)))/B0_field_nonround); + x=horzcat(x,offset); + punti = size(x); + cont_acqp=cont_acqp+1; + end + end + if strcmp(file_acqp(cont_acqp),'##$ACQ_coils')==1 + coil_info=char(file_acqp(cont_acqp+2)); + coil_info2=textscan(coil_info,'%s','delimiter',',', 'whitespace','()'); + coil_1H=coil_info2{1}{1} + end +end +size_zz=num_ppm; + + +file_imnd=textread('method','%s','delimiter','=','whitespace',''); +size_file_imnd=size(file_imnd); +for s=1:size_file_imnd(1) + + if strcmp(file_imnd(s),'##$PVM_MagTransPower')==1 + power=str2num(char(file_imnd(s+1))); + end + if strcmp(file_imnd(s),'##$PVM_MagTransModuleTime')==1 + durata=str2num(char(file_imnd(s+1))); + end +end + +% leggi nel file reco +file_reco=textread('pdata/1/reco','%s','delimiter','\n','whitespace',''); +size_reco=size(file_reco); +for cont_reco=1:size_reco(1) + if strcmp(file_reco(cont_reco),'##$RECO_size=( 2 )') + size_xy=str2num(char(file_reco(cont_reco+1))); + size_zx=size_xy(1); + size_zy=size_xy(2); + break + end +end + +file_meta=textread('pdata/1/reco','%s','delimiter','=','whitespace',''); +size_meta=size(file_meta); + +for cont_meta=1:size_meta(1) + if strcmp(file_meta(cont_meta),'##$RECO_wordtype') + reco_word=(char(file_meta(cont_meta+1))); + end +end + +if strcmp(reco_word,'_32BIT_SGN_INT') + reco_bit='uint32'; +elseif strcmp(reco_word,'_16BIT_SGN_INT') + reco_bit='uint16'; +end + +file_reco_slope=textread('pdata/1/reco','%s','delimiter','=','whitespace',''); + +% leggi la scala applicata ai valori in questa acquisizione +for cont_reco=1:size(file_reco_slope,1) + if strcmp(file_reco_slope(cont_reco),'##$RECO_map_slope') + slope=str2num(char(file_reco_slope(cont_reco+2))); + break + end +end + +id_read = waitbar(0,'Reading Mz stack...'); +image_Mz=zeros(size_zx,size_zy,1,size_zz); + +fid=fopen('pdata/1/2dseq'); +for cont_z=1:size_zz + waitbar(cont_z/size_zz); + for cont_i=1:size_zx + for cont_j=1:size_zy + image_Mz(cont_i,cont_j,1,cont_z)=fread(fid,1,reco_bit); + end + end +end +fclose(fid); +close(id_read); + +image_Mz = image_Mz/slope(1); % scale image by correction factor + +[~,IX] = sort(x'); +Mz_stack = double(image_Mz(:,:,:,IX)); \ No newline at end of file diff --git a/bruker/make_5D_B1_stack.m b/bruker/make_5D_B1_stack.m new file mode 100755 index 0000000..fe1f21c --- /dev/null +++ b/bruker/make_5D_B1_stack.m @@ -0,0 +1,25 @@ +function Z_stack = make_5D_B1_stack(protocol, directory_M0, M0_stack, Segment, varargin) +% ** function Z_stack = make_5D_B1_stack(protocol, directory_M0, M0_stack, Segment, seqDescr1, ..., seqDescrN) +% +% Creates 5D stack of z-spectra for each saturation field power (B1). +% Input is sequence descriptions of Mz images at different B1 values. +% Output 'Z_stack' dimensions are (x,y,z,w,B1). +% +% CT 20170111 + +nsat = nargin-4; +directory = regexprep(directory_M0, '/\d{1,2}$', ''); + +for i = 1:nsat + fprintf('Processing z-spectrum for ''%s'' ...\n', varargin{i}); + directory_Mz = getfolderpath(directory, protocol, varargin{i}); + Mz_stack = load_Mz(directory_Mz); + P = wipread_modified(directory_Mz, directory_M0); + + % calculate B0-corrected z-spectra + [dB0_stack_int,~] = MINFIND_SPLINE_3D(Mz_stack, Segment, P); + Mz_CORR = B0_CORRECTION(Mz_stack, dB0_stack_int, P, Segment); + Z_corrExt = NORM_ZSTACK(Mz_CORR, M0_stack, P, Segment); + + Z_stack(:,:,:,:,i) = Z_corrExt; +end \ No newline at end of file diff --git a/bruker/make_T1_directory.m b/bruker/make_T1_directory.m new file mode 100644 index 0000000..e0c554a --- /dev/null +++ b/bruker/make_T1_directory.m @@ -0,0 +1,30 @@ +function T1dirpath = make_T1_directory(directory, protocol, TI, ix_T1) +% ** function T1dirpath = make_T1_directory(TI, directory, protocol, ix_T1) +% +% Make folder containing DICOM files of TI images for T1 mapping. Changes +% DICOM header for compatibility with fitting function T1eval_levmar. +% Inputs: 'TI' is a vector of inversion recovery times. +% 'ix_T1' is a vector of protocol entry numbers of T1 images. +% Output is the full path to the new directory of T1 images. +% +% CT 20170115 + +T1dirpath = fullfile(directory, 'T1dicoms'); +mkdir(T1dirpath); +cd(T1dirpath); + +if length(ix_T1)==length(TI) + for i=1:length(TI) + pathT1 = fullfile(directory, 'T1dicoms', sprintf('TI%04d.dcm', TI(i))); + copyfile(fullfile(getfolderpath(directory,protocol,protocol(ix_T1(i)).SequenceDescription),... + 'pdata/1/dicom/MRIm1.dcm'), pathT1); + dmat = dicomread(pathT1); + dhead = dicominfo(pathT1); + dhead.AcquisitionNumber = i; + dhead.InstanceNumber = 1; + dhead.ImagesInAcquisition = length(TI); + dicomwrite(dmat, sprintf('TI%04d.dcm', TI(i)), dhead, 'CreateMode', 'copy'); + end +else + error('Number of TI images does not correspond to number of specified inversion times.') +end \ No newline at end of file diff --git a/bruker/mineExecFile.m b/bruker/mineExecFile.m new file mode 100644 index 0000000..f2525c1 --- /dev/null +++ b/bruker/mineExecFile.m @@ -0,0 +1,43 @@ +function paramval = mineExecFile(filepath, paramname) +% ** function paramval = mineExecFile(filepath, paramname) +% +% Extracts parameter from a Bruker executable file, e.g. 'method', 'acqp' +% or 'reco' file in the scan directory. +% e.g. B1value = mineExecFile(methodpath, '##$PVM_MagTransPower') +% See Paravision manual for list/description of (most) parameters. +% +% CT 20170113 + +FID = fopen(filepath,'r'); +A = textscan(FID,'%s','delimiter','=','whitespace',''); +fclose(FID); + +pos_param = find(strcmp(A{1}, paramname)) + 1; +if ~isempty(pos_param) + paramval = A{1}{pos_param}; + + % check for bad value (e.g. '( 1 )') + if regexp(paramval, '^\( \d.* \)$') + pos_param = pos_param+1; + paramval = A{1}{pos_param}; + end + % check for values (e.g. row vector) continuing on multiple cells + while isempty(regexp(A{1}{pos_param+1}, '^(#|\$)', 'once')) + pos_param = pos_param+1; + paramval = [paramval A{1}{pos_param}]; + end + % check for parameter list + if regexp(paramval, '^\(\d.*\)$') + ix = 1; % index of parameter in the list (not in manual--check in Paravision GUI!) + paramlist = regexp(paramval, '([^\(\)\s]+?)(?:,|$)', 'match'); + paramlist = regexprep(paramlist,',',''); + paramval = paramlist{ix}; + end + + paramdouble = str2double(cellstr(strsplit(paramval))); + if ~isnan(paramdouble) + paramval = paramdouble; + end +else + error(['No value found for parameter name ''',paramname,'''']); +end \ No newline at end of file diff --git a/bruker/readprotocol.m b/bruker/readprotocol.m new file mode 100755 index 0000000..d06ce1c --- /dev/null +++ b/bruker/readprotocol.m @@ -0,0 +1,41 @@ +function P = readprotocol(directory) +% ** function P = readprotocol(directory) +% +% Given patient directory, creates structure 'P' with fields: +% 'FolderNum': name of folder (in Bruker's directory structure, a number) +% 'SequenceDescription': extracted from DICOM header file as either the +% 'SeriesDescription' (sequence name set by the user at the scanner) +% or the 'ProtocolName' parameter. +% +% CT 20170113 + +P = dir(directory); +[P.FolderNum] = P.name; +P = rmfield(P, {'name', 'date', 'bytes', 'isdir', 'datenum'}); + +nonnumericmask = cellfun('isempty', regexp({P.FolderNum},'^\d+$')); +P(nonnumericmask) = []; % sequence folder names must be numbers +for i=1:length(P) + nzeros=0; + ERR=1; + while ERR + nzeros = nzeros+1; + try + dicomhdr = dicominfo(fullfile(directory,P(i).FolderNum,'pdata/1/dicom',sprintf('MRIm%0*d.dcm',nzeros,1))); + ERR=0; + catch + end + end + if isfield(dicomhdr, 'SeriesDescription') + P(i).SequenceDescription = dicomhdr.SeriesDescription; + elseif isfield(dicomhdr, 'ProtocolName') + P(i).SequenceDescription = dicomhdr.ProtocolName; + else + P(i).SequenceDescription = ''; + warning(sprintf('Could not extract sequence description of folder ''%s''.', P(i).FolderNum)); + end +end + +% sort entries by increasing folder number (to preserve sequence order) +[~,ixsort] = sort(str2double({P.FolderNum})); +P = P(ixsort); \ No newline at end of file diff --git a/bruker/wipread_modified.m b/bruker/wipread_modified.m new file mode 100755 index 0000000..b399e06 --- /dev/null +++ b/bruker/wipread_modified.m @@ -0,0 +1,53 @@ +function P = wipread_modified(directory_Mz, directory_M0) +% ** function P = wipread_modified(directory_Mz, directory_M0) +% +% Writes parameters of Mz and M0 images to structure 'P'. +% Modified for compatibility with Bruker/Paravision data format. +% +% CT 20161205 + +try + methodpath = [directory_Mz, '/method']; + acqppath = [directory_Mz, '/acqp']; + acqppathM0 = [directory_M0, '/acqp']; + recopath = [directory_Mz, '/pdata/1/reco']; + dicompath = [directory_Mz, '/pdata/1/dicom/MRIm01.dcm']; + dicomhdr = dicominfo(dicompath); +catch ME + ME +end + +% cell strings for the string parameters +pulse_shape = {'Gauss' 'Sinc' 'Rect' 'Spinlock' 'Adia Fullpass' 'Adia Spinlock' 'AdiaInvRec' 'Fermi'}; +spoiling = {'none' 'constant' 'alternating' 'varying'}; +sampling = {'regular' 'alternating' 'List' 'SingleOffset'}; + +% write all parameters into P struct +try + P.SEQ.imageflipangle = dicomhdr.FlipAngle; + P.SEQ.averages = dicomhdr.NumberOfAverages; + P.SEQ.nominalB0 = dicomhdr.MagneticFieldStrength; + P.SEQ.FS = round(P.SEQ.nominalB0); + P.SEQ.FREQ = dicomhdr.ImagingFrequency; + P.SEQ.measurements = dicomhdr.ImagesInAcquisition; % number of repetitions + P.SEQ.TR = dicomhdr.RepetitionTime; + P.SEQ.TE = dicomhdr.EchoTime; + P.SEQ.B1 = mineExecFile(methodpath, '##$PVM_MagTransPower'); + P.SEQ.n = mineExecFile(methodpath, '##$PVM_MagTransPulsNumb'); + P.SEQ.tp = mineExecFile(methodpath, '##$PVM_MagTransPulse1')/1000; % pulse duration in seconds + P.SEQ.DC = mineExecFile(acqppath, '##$HDDUTY'); % homodecoupling duty cycle? (CHECK) + P.SEQ.recovertime = mineExecFile(acqppath, '##$ACQ_recov_time'); + P.SEQ.recovertimeM0 = mineExecFile(acqppathM0, '##$ACQ_recov_time'); + P.SEQ.sampling = sampling{1}; % let's assume regular sampling (CHECK) + P.SEQ.w = mineExecFile(acqppath, '##$ACQ_O2_list')'/P.SEQ.FREQ; + P.SEQ.Offset = max(P.SEQ.w); + P.SEQ.stack_dim = [mineExecFile(recopath, '##$RECO_size'), 1, P.SEQ.measurements]; + P.EVAL.N_asym = fix((P.SEQ.measurements-1)/2); + P.EVAL.w_fit = (min(P.SEQ.w):0.01:max(P.SEQ.w))'; + P.EVAL.w_interp = P.SEQ.w; + P.EVAL.lowerlim_slices = 1; + P.EVAL.upperlim_slices = 1; % we're working on single-slice data +catch ME + warning(sprintf('Could not write all parameters into P structure.\n')); + ME +end \ No newline at end of file diff --git a/levmar_fit/fitmodelfunc_ANA.m b/levmar_fit/fitmodelfunc_ANA.m index a6f5cea..66b96f1 100644 --- a/levmar_fit/fitmodelfunc_ANA.m +++ b/levmar_fit/fitmodelfunc_ANA.m @@ -259,7 +259,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -elseif (strcmp(modelname,'lorentzfit6pool_BW')) % fόgt einen Lorentz ein, der ein plateau der breite BW hat +elseif (strcmp(modelname,'lorentzfit6pool_BW')) % f?gt einen Lorentz ein, der ein plateau der breite BW hat fZi=@(k) k*0+A(1) ; f1=@(k) L0_BW(A(1:4),A(20),k); f2=@(k) A(5).*A(6).^2/4./ (A(6).^2/4+(k-A(7)).^2) ; @@ -282,7 +282,7 @@ elseif ( strcmp(modelname,'WASABIFIT') ) freq=P.SEQ.FREQ; - t_p=P.SEQ.tp/(1E6); + t_p=P.SEQ.tp; B1=A(1); offset=A(2); c=A(3); @@ -294,7 +294,7 @@ elseif ( strcmp(modelname,'WASABIFIT_2') ) freq=P.SEQ.FREQ; - t_p=P.SEQ.tp/(1E6); + t_p=P.SEQ.tp; B1=A(1); offset=A(2); c=A(3); diff --git a/levmar_fit/matlabfit.m b/levmar_fit/matlabfit.m old mode 100644 new mode 100755 index eaf8819..ff11761 --- a/levmar_fit/matlabfit.m +++ b/levmar_fit/matlabfit.m @@ -6,13 +6,13 @@ %to be split up in seperate variables, that is a bit stupid to programm to %be general working, i implement it for the optimization toolbox firt. % -% Algorithm — Algorithm used for the fitting procedure: -% Trust-Region — This is the default algorithm and must be used if you specify Lower or Upper coefficient constraints. -% Levenberg-Marquardt — If the trust-region algorithm does not produce a reasonable fit, and you do not have coefficient constraints, try the Levenberg-Marquardt algorithm. +% Algorithm ? Algorithm used for the fitting procedure: +% Trust-Region ? This is the default algorithm and must be used if you specify Lower or Upper coefficient constraints. +% Levenberg-Marquardt ? If the trust-region algorithm does not produce a reasonable fit, and you do not have coefficient constraints, try the Levenberg-Marquardt algorithm. % Finite Differencing Parameters % -% DiffMinChange — Minimum change in coefficients for finite difference Jacobians. The default value is 10-8. -% DiffMaxChange — Maximum change in coefficients for finite difference Jacobians. The default value is 0.1. +% DiffMinChange ? Minimum change in coefficients for finite difference Jacobians. The default value is 10-8. +% DiffMaxChange ? Maximum change in coefficients for finite difference Jacobians. The default value is 0.1. % Note that DiffMinChange and DiffMaxChange apply to: % % Any nonlinear custom equation, that is, a nonlinear equation that you write @@ -21,17 +21,17 @@ % % Fit Convergence Criteria % -% MaxFunEvals — Maximum number of function (model) evaluations allowed. The default value is 600. -% MaxIter — Maximum number of fit iterations allowed. The default value is 400. -% TolFun — Termination tolerance used on stopping conditions involving the function (model) value. The default value is 10-6. -% TolX — Termination tolerance used on stopping conditions involving the coefficients. The default value is 10-6. +% MaxFunEvals ? Maximum number of function (model) evaluations allowed. The default value is 600. +% MaxIter ? Maximum number of fit iterations allowed. The default value is 400. +% TolFun ? Termination tolerance used on stopping conditions involving the function (model) value. The default value is 10-6. +% TolX ? Termination tolerance used on stopping conditions involving the coefficients. The default value is 10-6. % Coefficient Parameters % -% Unknowns — Symbols for the unknown coefficients to be fitted. -% StartPoint — The coefficient starting values. The default values depend on the model. For rational, Weibull, and custom models, default values are randomly selected within the range [0,1]. For all other nonlinear library models, the starting values depend on the data set and are calculated heuristically. See optimized starting points below. -% Lower — Lower bounds on the fitted coefficients. The tool only uses the bounds with the trust region fitting algorithm. The default lower bounds for most library models are -Inf, which indicates that the coefficients are unconstrained. However, a few models have finite default lower bounds. For example, Gaussians have the width parameter constrained so that it cannot be less than 0. See default constraints below. -% Upper — Upper bounds on the fitted coefficients. The tool only uses the bounds with the trust region fitting algorithm. The default upper bounds for all library models are Inf, which indicates that the coefficients are unconstrained. -% For more information about these fit options, see the lsqcurvefit function in the Optimization Toolbox™ documentation. +% Unknowns ? Symbols for the unknown coefficients to be fitted. +% StartPoint ? The coefficient starting values. The default values depend on the model. For rational, Weibull, and custom models, default values are randomly selected within the range [0,1]. For all other nonlinear library models, the starting values depend on the data set and are calculated heuristically. See optimized starting points below. +% Lower ? Lower bounds on the fitted coefficients. The tool only uses the bounds with the trust region fitting algorithm. The default lower bounds for most library models are -Inf, which indicates that the coefficients are unconstrained. However, a few models have finite default lower bounds. For example, Gaussians have the width parameter constrained so that it cannot be less than 0. See default constraints below. +% Upper ? Upper bounds on the fitted coefficients. The tool only uses the bounds with the trust region fitting algorithm. The default upper bounds for all library models are Inf, which indicates that the coefficients are unconstrained. +% For more information about these fit options, see the lsqcurvefit function in the Optimization Toolbox? documentation. % fitable_func =@(x,w,Z) feval(fitfunc,x, w, P); @@ -69,7 +69,11 @@ [w, M] = prepareCurveData(w, M ); -fitable_func =@(x,w,Z) feval(fitfunc,x, w, P)'; +if strcmp(fitfunc,'T1recovery') || strcmp(fitfunc,'T1recovery_biex') + fitable_func =@(x,w,Z) feval(fitfunc,x, w, P); +else + fitable_func =@(x,w,Z) feval(fitfunc,x, w, P)'; +end % fit-options matlab_options = optimset('TolFun',options(4),'TolX',options(3), 'MaxIter',nIter,'Display','off'); @@ -80,12 +84,13 @@ popt=0; try -[popt resnorm RES,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN] = lsqcurvefit(fitable_func,p0,w,M,lb,ub,matlab_options); -ret=EXITFLAG; -info=OUTPUT; -catch - warning('lsqrfit retrned with an error'); - end; + [popt resnorm RES,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN] = lsqcurvefit(fitable_func,p0,w,M,lb,ub,matlab_options); + ret=EXITFLAG; + info=OUTPUT; +catch ME + warning('lsqrfit returned with an error'); + ME +end; try [ci, varb, corrb, varinf] = nlparci(popt,RES,JACOBIAN,0.95); diff --git a/ultravist_CEST/get_FIT_LABREF_modified.m b/ultravist_CEST/get_FIT_LABREF_modified.m new file mode 100755 index 0000000..00e8f0b --- /dev/null +++ b/ultravist_CEST/get_FIT_LABREF_modified.m @@ -0,0 +1,62 @@ +function [Zlab, Zref] = get_FIT_LABREF_modified(popt, P, Segment, mode, x_inter) +% function [Zlab, Zref] = get_FIT_LABREF_modified(popt, P, Segment, mode, x_inter) +% +% Version of get_FIT_LABREF modified for Ultravist analysis. +% +% CT 20170111 + + +if nargin < 5 + x_inter=P.SEQ.w; +end + +switch mode + case 'invivo' + Zref_fields = {'Background', 'Amide', 'NOE', 'MT', 'Amine', 'onlyWater'}; + case 'ultravist' + Zref_fields = {'Background', 'ppm4p2', 'ppm5p6', 'onlyWater'}; + otherwise + error('Input ''mode'' can only be ''invivo'' or ''ultravist''') +end + +sizes=size(popt); +Zlab=zeros(sizes(1),sizes(2),sizes(3),numel(x_inter)); +for nfield = 1:length(Zref_fields) + Zref.(Zref_fields{nfield}) = zeros(sizes(1),sizes(2),sizes(3),numel(x_inter)); +end + +if ndims(Segment) == 2 + for ii=1:sizes(3) + Segment(:,:,ii)=Segment(:,:,1); + end +end +clear ii + +for kk=1:sizes(3)%z + for ii=1:sizes(1) %x + h=waitbar(ii/sizes(1)); + for jj=1:sizes(2) %y + + if Segment(ii,jj,kk)==1 + + [f, fZi, f1, f2, f3, f4, f5, f6, g, g2, g3, g4, g5]= fitmodelfunc_ANA(popt(ii,jj,kk,:),P); + Zlab(ii,jj,kk,:) = f(x_inter); + Zref.(Zref_fields{1})(ii,jj,kk,:) = g(x_inter); + Zref.(Zref_fields{2})(ii,jj,kk,:) = g2(x_inter); + Zref.(Zref_fields{3})(ii,jj,kk,:) = g3(x_inter); + if strcmp(mode,'invivo') + Zref.(Zref_fields{4})(ii,jj,kk,:) = g4(x_inter); + Zref.(Zref_fields{5})(ii,jj,kk,:) = g5(x_inter); + end + Zref.(Zref_fields{end})(ii,jj,kk,:) = fZi(x_inter)-f1(x_inter); + else + Zlab(ii,jj,kk,:)=NaN(numel(x_inter),1); + for nfield = 1:length(Zref_fields) + Zref.(Zref_fields{nfield})(ii,jj,kk,:) = NaN(numel(x_inter),1); + end + end; + + end; + end; + end; +close(h); diff --git a/ultravist_CEST/lorentzianfit_main.m b/ultravist_CEST/lorentzianfit_main.m new file mode 100755 index 0000000..35bf5fb --- /dev/null +++ b/ultravist_CEST/lorentzianfit_main.m @@ -0,0 +1,46 @@ +function [Zlab, Zref, Pnew, popt] = lorentzianfit_main(Z_corrExt, P, Segment, mode) +% ** function [Zlab, Zref, Pnew, popt] = lorentzianfit_main(Z_corrExt, P, Segment, mode) +% 'mode' can be: +% 'invivo' (default) uses standard 5-pool Lorentzian model (water, amide, +% NOE, MT, amine). +% 'ultravist' uses 3 pools corresponding to Ultravist z-spectrum peaks +% (water, @4.2ppm, @5.6ppm). +% +% CT 20170111 + +if nargin<4 + mode='invivo'; +end + +Pnew = P; +Pnew.FIT.options = [1E-04, 1E-15, 1E-10, 1E-04, 1E-06]; +Pnew.FIT.nIter = 50; +Pnew.FIT.extopt = 1; % change parameters explicitly + +switch mode + case 'invivo' + Pnew.FIT.modelnum = 5; % number of Lorentzian pools (possible 1-5) + %Lorentzian line LI defined by amplitude Ai, width Gi [ppm] and offset dwi[ppm]: Li=Ai.*Gi^2/4./ (Gi^2/4+(dw-dwi).^2) ; + %1=water; 2=amide; 3=NOE; 4=MT; 5=amine + %const.Zi A1 G1 dw1 A2 G2 dw2 A3 G3 dw3 A4 G4 dw4 A5 G5 dw5 + lb = [ 0.5 0.02 0.3 -1 0.0 0.4 +3 0.0 1 -4.5 0.0 10 -4 0.0 1 1 ]; + ub = [ 1 1 10 +1 0.2 3 +4 0.4 5 -2 1 100 -1 0.2 3.5 2.5 ]; + 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 ]; + case 'ultravist' + Pnew.FIT.modelnum = 3; + %Lorentzian line LI defined by amplitude Ai, width Gi [ppm] and offset dwi[ppm]: Li=Ai.*Gi^2/4./ (Gi^2/4+(dw-dwi).^2) ; + %1=water; 2=@4.2ppm; 3=@5.6ppm + %const.Zi A1 G1 dw1 A2 G2 dw2 A3 G3 dw3 + lb = [ 0.5 0.02 0.3 -1 0.0 1 3.8 0.0 1 5 ]; + ub = [ 1 1 10 +1 0.6 5 4.8 0.6 5 6 ]; + p0 = [ 1 0.9 1.4 0 0.1 3 4.2 0.1 3 5.6 ]; + otherwise + error('Input ''mode'' can only be ''invivo'' or ''ultravist''') +end + +Pnew.FIT.lower_limit_fit = lb; +Pnew.FIT.upper_limit_fit = ub; +Pnew.FIT.start_fit = p0; + +[popt, Pnew] = FIT_3D(Z_corrExt, Pnew, Segment); % perform the fit pixelwise +[Zlab, Zref] = get_FIT_LABREF_modified(popt, Pnew, Segment, mode); % create Reference values Z_Ref=(Z_lab - Li) \ No newline at end of file diff --git a/ultravist_CEST/main_ultravist.m b/ultravist_CEST/main_ultravist.m new file mode 100644 index 0000000..b4bfcdd --- /dev/null +++ b/ultravist_CEST/main_ultravist.m @@ -0,0 +1,58 @@ +% ======================= main_ultravist ======================= % +% Script for pH mapping with Ultravist phantom (work in progress) +% CT 20170113 + +%% calculate log10ratio +CESTratiomap = log10ratio(popt, P, Segment, [4.2 5.6]); + +%% -------------------- CEST-pH calibration -------------------- +%% Load phantom scans and fit z-spectra +LOAD_BATCH_cest-sources +Segment = make_Segment(M0_stack, 'ellipse', mean(M0_stack(M0_stack>0)).*[0.3]); +Z_uncorr = NORM_ZSTACK(Mz_stack,M0_stack,P,Segment); +[~, ~, P, popt] = lorentzianfit_main(Z_uncorr, P, Segment, 'ultravist'); + +%% log10 ratio of MTR_Rex +delta = [4.2 5.6]; +CESTratio = log10ratio(popt, P, Segment, delta); +figure; imagesc(CESTratio); colorbar; +title('log_{10}[MTR_{4.2ppm} / MTR_{5.6ppm}]') + +%% extract mean log10ratio from phantom ROIs +W = evalin('base','whos'); +doesSegmentExist = ismember('Segment_ph',[W(:).name]); +if ~doesSegmentExist + Segment_ph = nan(size(CESTratio,1), size(CESTratio,2), 7); + for i=1:7 + Segment_ph(:,:,i) = make_Segment(CESTratio, 'ellipse'); + end +end + +CESTratio_ROIs = nan(1,7); +for i=1:7 + roivals = CESTratio .* Segment_ph(:,:,i); + CESTratio_ROIs(i) = nanmean(roivals(:)); +end +clear roivals + +%% pH-CEST curve +ph = [6.1 6.3 6.5 6.7 6.9 7.1 7.3]; +% phantom# 7 6 5 4 3 2 1 +CESTratio_ROIs = flip(CESTratio_ROIs); +figure; +plot(ph, CESTratio_ROIs, '+', 'markersize', 8) +set(gca,'xtick',ph); +xlabel('pH') +ylabel('log_{10}[((M0-Mz)/Mz)_{4.2ppm} / ((M0-Mz)/Mz)_{5.6ppm}]') +title(sprintf('pH-CEST calibration at B1=%duT', P.SEQ.B1)) +hold on +lsline + +%% linear regression fit +% calibration model (pH ~ CESTratio) +calibmdl = fitlm(CESTratio_ROIs,ph); + +%% ----------------------- pH map ------------------------ +pH_pred = predict_pH(calibmdl,CESTratiomap); +figure; imagesc(pH_pred); colorbar; +title('pH map') \ No newline at end of file diff --git a/ultravist_CEST/predict_pH.m b/ultravist_CEST/predict_pH.m new file mode 100755 index 0000000..329024b --- /dev/null +++ b/ultravist_CEST/predict_pH.m @@ -0,0 +1,12 @@ +function y = predict_pH(mdl,x) +% ** function y = predict_pH(mdl,x) +% +% 'mdl' is a fit object specifying the linear model (CESTratio ~ pH). +% 'x' is the CESTratio map. +% 'y' is the predicted pH map. + +imsize = size(x); +x = x(:); +y = predict(mdl,x); +y(y<5 | y>8) = NaN; +y = reshape(y,imsize); \ No newline at end of file