Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

sequence/parameters optimization #440

Open
zrazi opened this issue Apr 8, 2021 · 12 comments
Open

sequence/parameters optimization #440

zrazi opened this issue Apr 8, 2021 · 12 comments

Comments

@zrazi
Copy link

zrazi commented Apr 8, 2021

How would you optimize your sequence/parameters for T1/T2 mapping using IR for different T1/T2 series? What would you set up as your ground truth? What is the best way #to evaluate your model afterward quantitively?
@AlexMorello
Copy link

Same question!

@agahkarakuzu
Copy link
Member

For T1 mapping, I suggest these blog posts.

As for T2, e.g., for the mono_t2 model:

image

You can run single voxel curve simulation for the TEs shown on the protocol panel to fit data simulated for T2/M0 values you set:

image

It also reports Cramer-Rao lower bound. By changing TEs, T2/M0 and SNR, you can see how estimation changes. It is also important that your data acquisition complies with the signal representations the models are build on. For example, mono_t2 is based on a multi-echo spin echo acquisition.

You can expand this simulation to a sensitivity analysis to see the behaviour for a range of values. The simulations are available for inversion_recovery as well.

I highly recommend blog post on vfa_t1 for a more complete understanding.

@agahkarakuzu
Copy link
Member

What is the best way #to evaluate your model afterward quantitively?

For relaxometry measurements, having a ground truth is the most convenient. If you have a range of T1/T2 reference values, then you can look at percent error to see if your protocol is acceptable.

If not, one way is to check how data was fitted to the model per voxel (there's a sanity check section on the GUI). Although this gives an overall information about the appropriateness of the fit, not an actual measurement of accuracy.

If you think that your quantifications makes sense for whichever substance you are imaging, then another good practice is to test the stability of your measurements. I can offer some useful stats for such evaluation, if it is of concern.

On the other hand, these are broad pointers. A more complete answer often calls for details of the application, which sequence was used, what was the target etc. I hope that some of these answers are useful for you. Cheers!

@zrazi
Copy link
Author

zrazi commented Jan 19, 2022 via email

@agahkarakuzu
Copy link
Member

agahkarakuzu commented Jan 19, 2022 •
edited

As for stability, you could say that reproducibility over time. More common terms are reliability/consistency/agreement. As it is quite common in the literature, intraclass correlations would be a good start, this article is insightful for that. Again, it depends on your experimental design and what you'd like to achieve.

Have you seen the documentation about IR? It shows that for 9 TIs the input data is [128, 128, 1, 9]. In which format are your multiple inversion images? We can add a simple tool to merge and sort volumes.

@zrazi
Copy link
Author

zrazi commented Jan 19, 2022 via email

@agahkarakuzu
Copy link
Member

Can you send me a set of example files via email, I can quickly come up with something for you.

For mask, only volumetric dimensions should match (e.g., [10x1] TI, [128,128,2,10] IRData, [128,128,2] mask)

@zrazi
Copy link
Author

zrazi commented Jan 19, 2022 via email

@agahkarakuzu
Copy link
Member

I think you replied to the comment via email with an attachment, but they don't show up here. You'll need to send it to my address, I think you have it.

@zrazi
Copy link
Author

zrazi commented Jan 20, 2022 via email

@agahkarakuzu
Copy link
Member

agahkarakuzu commented Jan 21, 2022 •
edited

@zrazi

  1. Create a new matlab script with the following code and save it as dicom2ir.m.
  2. Make sure that qMRLab is added to your path
  3. Run the code, and select multiple DICOM files (as many TIs as you want) and confirm

The code will:

  • sort and merge files according to TI
  • create two masks for your phantom using simple morphometry, one binary and one indexed
  • fit data using merged volume, set protocol values and the binary mask
  • will save FitResults.mat and a csv (mean STD for T1 with the respective mask index).
  • You will find the outputs in the same directory where your input images are.

The mask detection is far from robust, yet still may perform well for this particular purpose.

function dicom2ir

[list,folder]=uigetfile('*','Select input DICOM images','MultiSelect','on');
lut = cell(2,length(list));
lut(1,:) = list;
for ii=1:length(list)
    in = dicominfo([folder filesep list{ii}]);
    lut(2,ii) = num2cell(in.InversionTime);
    im = dicomread(in);
    if ii==1
       sz = size(im);
       if ndims(sz) < 3
           ims = zeros(sz(1),sz(2),1,length(list)); 
       else
           ims = zeros(sz(1),sz(2),sz(3),length(list)); 
       end
    end
    
    ims(:,:,:,ii) = double(im);
end


[TI,idx] = sort([lut{2,:}]);
IRData = ims(:,:,:,idx);


Model = inversion_recovery;
Model.Prot.IRData.Mat = TI';
Model.Prot.TimingTable.Mat = in.RepetitionTime;
data = struct();
data.IRData = IRData; 

% A simple mask generator only for this phantom 
[a,b,~] = imfindcircles(mat2gray(ims(:,:,1,end)),[10 20],'Sensitivity', 0.98);
b = b - 2; % Erode radius by 2 pixels

% A simple sort logic assuming that all 12 circles are detected
% If not, will skip sorting
try
[a,idx] = sortrows(a,1);
a = round(a);
b = b(idx);
yall = repmat(sort(a(1:3,2)),[4,1]);
catch
yall = a(:,2);    
end
theta = (0:1000-1)*(2*pi/1000);
x = a(:,1) + b.*cos(theta);
y = yall + b.*sin(theta);

mask = zeros(sz(1),sz(2));
maskIdx = zeros(sz(1),sz(2));
for ii = 1:length(a)
    cur = poly2mask(x(ii,:),y(ii,:),sz(1),sz(2));
    mask = mask + cur ;
    maskIdx = maskIdx + cur*ii;
end

data.Mask= double(mask);
FitResults = FitData(data,Model,0);

FitResults.maskIdx = maskIdx;

% In case there are not 12 circles captured
idxs = unique(maskIdx(maskIdx(:)>0));
stat = zeros(length(idxs),3);
for ii=1:length(idxs)
    stat(ii,1) = ii; 
    stat(ii,2) = nanmean(FitResults.T1(maskIdx==ii));
    stat(ii,3) = nanstd(FitResults.T1(maskIdx==ii));
end

stat = sortrows(stat,2);
writematrix(stat,[folder filesep 'simplestat.csv']) 

save([folder filesep 'FitResults.mat'],'FitResults');
end

image

@agahkarakuzu
Copy link
Member

agahkarakuzu commented Feb 15, 2022 •
edited

@zrazi the script below works the same way, but for the mono_t2 model for T2 mapping:

function dicom2monot2

[list,folder]=uigetfile('*','Select input DICOM images','MultiSelect','on');
lut = cell(2,length(list));
lut(1,:) = list;
for ii=1:length(list)
    in = dicominfo([folder filesep list{ii}]);
    lut(2,ii) = num2cell(in.EchoTime);
    im = dicomread(in);
    if ii==1
       sz = size(im);
       if ndims(sz) < 3
           ims = zeros(sz(1),sz(2),1,length(list)); 
       else
           ims = zeros(sz(1),sz(2),sz(3),length(list)); 
       end
    end
    
    ims(:,:,:,ii) = double(im);
end


[TE,idx] = sort([lut{2,:}]);
SEdata = ims(:,:,:,idx);


Model = mono_t2;
Model.Prot.SEdata.Mat = TE';

data = struct();
data.SEdata = SEdata; 

% A simple mask generator only for this phantom 
[a,b,~] = imfindcircles(mat2gray(ims(:,:,1,end)),[10 20],'Sensitivity', 0.98);
b = b - 2; % Erode radius by 2 pixels

% A simple sort logic assuming that all 12 circles are detected
% If not, will skip sorting
try
[a,idx] = sortrows(a,1);
a = round(a);
b = b(idx);
yall = repmat(sort(a(1:3,2)),[4,1]);
catch
yall = a(:,2);    
end
theta = (0:1000-1)*(2*pi/1000);
x = a(:,1) + b.*cos(theta);
y = yall + b.*sin(theta);

mask = zeros(sz(1),sz(2));
maskIdx = zeros(sz(1),sz(2));
for ii = 1:length(a)
    cur = poly2mask(x(ii,:),y(ii,:),sz(1),sz(2));
    mask = mask + cur ;
    maskIdx = maskIdx + cur*ii;
end

data.Mask= double(mask);
FitResults = FitData(data,Model,0);

FitResults.maskIdx = maskIdx;

% In case there are not 12 circles captured
idxs = unique(maskIdx(maskIdx(:)>0));
stat = zeros(length(idxs),3);
for ii=1:length(idxs)
    stat(ii,1) = ii; 
    stat(ii,2) = nanmean(FitResults.T2(maskIdx==ii));
    stat(ii,3) = nanstd(FitResults.T2(maskIdx==ii));
end

stat = sortrows(stat,2);
writematrix(stat,[folder filesep 'simplestat.csv']) 

save([folder filesep 'FitResults.mat'],'FitResults');
end
  • Create a new matlab script with the following code and save it as dicom2monot2.m.
  • Make sure that qMRLab is added to your path
  • Run the code, and select multiple DICOM files (as many TEs as you want) and confirm
    image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@agahkarakuzu @zrazi @AlexMorello