[FieldTrip] Incorporating White Matter conductivity anisotropy into FEM simbio
Marios Antonakakis
antonakakismar at gmail.com
Tue Nov 6 14:21:56 CET 2018
Hi Ravi,
assuming that you work under Linux and you have already the DTI data
registered to MRI, you need then to run the below script.
*%%1, calculate conductivity tensor for every voxel*
% mri_segmented_final, V1, ... L1 variables have ft MRI structure
%skin,skull compacta, skull spongiosa, CSF, GM, WM
conductivity = [0.4300, 0.0024, 0.0086, 1.7900, 0.3300,
0.1400];
cfg =[];
cfg.compartments = mri_segmented_final.anatomylabel; % anatomical labels,
skin: 0, skull spongiosa 2 and show on...
cfg.conductivity = conductivity;
[condcell, s, fail] =
sb_calcTensorCond_tuch(cfg,mri_segmented_final,V1,V2,V3,L1,L2,L3);
*%% 2. assign the conductivity tensor with the hex mesh *
[mesh.nodes,mesh.elements,mesh.labels]= read_vista_mesh('foo_mesh.v);
mask = mri_segmented_final;
mask.anatomy((mask.anatomy~=5) & (mask.anatomy~=6))=0;
mask.anatomy((mask.anatomy==5) | (mask.anatomy==6))=1;
tensors =
sb_assiTensorCond(mask,mesh.nodes,mesh.elements,mesh.labels,condcell);
write_vista_mesh('foo_mesh_aniso.v',mesh.nodes,mesh.elements,mesh.labels,tensors');
Do not hesitate to ask further questions.
**Note that write_vista_mesh and read_vista_mesh are private ft functions.*
Best regards,
Marios
ᐧ
Στις Τρί, 6 Νοε 2018 στις 12:10 μ.μ., ο/η Carsten Wolters <
carsten.wolters at uni-muenster.de> έγραψε:
> Hi Ravi,
>
> Marios (in CC) promised to send you the short Matlab-script that we use to
> transform the diffusion tensors
> to conductivity tensors using Tuch's effective medium approach.
>
> For this approach, please check e.g. the subsection
> "Calibrated Finite Element Head Model and Forward Solution"
> in
> https://link.springer.com/article/10.1007/s10548-017-0568-9
>
> BR
> Carsten
>
> Am 05.11.18 um 18:26 schrieb Ravi Mill:
>
> Hi Carsten and Johannes
>
>
> Many thanks for responding, and for developing these great tools!
>
>
> I'm in the process of acquiring a large diffusion MR dataset from which I
> can hopefully create an 'averaged' atlas. From your responses I think I
> have a sense of how to integrate the conductivity tensors derived from this
> atlas with the Fieldtrip FEM pipeline.
>
>
> But I was wondering if you had any advice on how to compute
> these conductivity tensors in the first place? From the paper that Carsten
> sent, It seems like the FDT program within FSL is what I need to compute
> diffusion tensors from the raw diffusion images (steps 1-6 from the FDT
> user guide
> https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/UserGuide#Processing_pipeline)?
> Seemingly, these diffusion tensors need to then be converted to
> conductivity tensors - any advice on how to do this (or if you could point
> me to some example code) would be greatly appreciated.
>
>
> Thanks
>
> Ravi
>
>
> FDT/UserGuide - FslWiki - University of Oxford
> <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/UserGuide#Processing_pipeline>
> fsl.fmrib.ox.ac.uk
> Eddy Current Correction. Eddy currents in the gradient coils induce
> (approximate) stretches and shears in the diffusion weighted images. These
> distortions are different for different gradient directions.
>
>
> Thanks again
> Ravi
>
>
> ------------------------------
> *From:* fieldtrip <fieldtrip-bounces at science.ru.nl>
> <fieldtrip-bounces at science.ru.nl> on behalf of Johannes Vorwerk
> <j.vorw01 at gmail.com> <j.vorw01 at gmail.com>
> *Sent:* Monday, October 29, 2018 10:14:16 AM
> *To:* FieldTrip discussion list
> *Subject:* Re: [FieldTrip] Incorporating White Matter conductivity
> anisotropy into FEM simbio
>
> Dear Ravi,
>
> as Carsten already said, calculating FEM with anisotropic conductivities
> is not directly supported by the FieldTrip-SimBio implementation. However,
> if you are willing to invest a bit of time it is possible to work around
> this.
>
> The „only“ thing that needs to be changed is the calculation of the FEM
> stiffness matrix, which is performed by the routine
> „calc_stiffness_matrix_val“ in the function sb_calc_stiff (usually called
> from ft_prepare_headmodel). The problem is that FieldTrip does not support
> anisotropic conductivities, so that you would have to call
> calc_stiffness_matrix_val directly. You can see the correct call in
> sb_calc_stiff. For anisotropic conductivities you have to replace the input
> „cond“ by a #elements x 6 matrix containing your anisotropic conductivities
> in the format "xx yy zz xy yz zx“. If you now follow the normal
> FieldTrip-SimBio workflow using the resulting stiffness matrix, you will
> get results for anisotropic conductivities.
>
> Best,
> Johannes
>
> Am 29.10.2018 um 12:31 schrieb Carsten Wolters <
> carsten.wolters at uni-muenster.de>:
>
> Dear Ravi,
>
> 1) You can use the pure SimBio-code from
> https://www.mrt.uni-jena.de/simbio/index.php/Main_Page
> <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.mrt.uni-jena.de%2Fsimbio%2Findex.php%2FMain_Page&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863059067&sdata=cFED89%2BXMUfro9URWv5pWzGC5SnzTqEHT%2FmXd%2F2eB8Q%3D&reserved=0>
> to treat WM anisotropy.
> While it would in principle also be possible to use anisotropic
> conductivities with FieldTrip-SimBio,
> this is currently not implemented using ft_prepare_headmodel. Johannes (in
> CC), who implemented
> Fieldtrip-SimBio, answered a same question by Junjie Wu in March 2018:
> "Depending on your matlab skills and your available time, I could help you
> to give it a
> try though. It should be possible with using some direct function calls
> instead of the high-level fieldtrip-functions."
>
> 2) We recommend
>
> http://www.sci.utah.edu/~wolters/PaperWolters/2012/RuthottoEtAl_PhysMedBiol_2012.pdf
> <https://na01.safelinks.protection.outlook.com/?url=http:%2F%2Fwww.sci.utah.edu%2F~wolters%2FPaperWolters%2F2012%2FRuthottoEtAl_PhysMedBiol_2012.pdf&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863069076&sdata=tYM5Oabmkx3Gr4R0I1Wqauly3lVENtp1H%2Fnebvcl0pw%3D&reserved=0>
> on individual data. I could imagine that an atlas does a reasonable job
> w.r.t. the main
> bigger fiber tracts such as corpus callosum or pyramidal tracts, but that
> the finer details
> in the cortices are individual. We always measure T1, T2 and DTI from each
> subject
> and I personally do not have experience with such a group-level anisotropy
> compared
> to the individual one. Might be interesting to hear from others what they
> think!?
>
> BR
> Carsten
>
>
>
> Am 25.10.18 um 23:05 schrieb Ravi Mill:
>
> Dear Fieldtrippers
>
> I have applied the FEM simbio head modeling pipeline implemented
> in Fieldtrip to my EEG data. My understanding is that this pipeline
> assumes isotropic conductivities for 5 head compartments (as specified by
> cfg.conductivity in ft_prepare_headmodel). After reading some papers
> (e.g. Vorwerk et al 2014 https://doi.org/10.1016/j.neuroimage.2014.06.040
> <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1016%2Fj.neuroimage.2014.06.040&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863069076&sdata=GbsiChTzmUZMEBsH1k44etRdVJUcUQ6SskLJtnzkbSk%3D&reserved=0>),
> it seems like incorporating white matter conductivity anisotropy has a
> relatively small albeit significant effect on the source solution. I am
> interested in comparing FEM results when treating white matter as
> anisotropic. My questions are as follows:
>
>
> 1. Is there a way to implement the FEM simbio head model whilst
> treating WM as anisotropic within Fieldtrip? If so, how would one do this
> (or are there any resources available that demonstrate this)?
> 2. From previous papers and some simbio documentation (
> https://www.mrt.uni-jena.de/simbio/index.php/SIMBIO/Releasenotes/Examples
> <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.mrt.uni-jena.de%2Fsimbio%2Findex.php%2FSIMBIO%2FReleasenotes%2FExamples&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863079077&sdata=LLPhLC8V85SmMAJbnV%2F15%2Bcs366ek6FQbx81p25mwqA%3D&reserved=0>)
> it seems like diffusion MRI data is required to calculate the WM
> conductivity for each individual subject. I only have T1 and T2 scans for
> my subjects. So would it be possible to use WM anisotropic information
> obtained from some kind of diffusion MRI group average/atlas instead
> (accepting some loss in subject-level precision)? If so, does such a group
> average/atlas exist?
>
>
> Any help would be greatly appreciated!
>
> Thanks
> Ravi
>
>
>
> _______________________________________________
> fieldtrip mailing listhttps://mailman.science.ru.nl/mailman/listinfo/fieldtrip <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmailman.science.ru.nl%2Fmailman%2Flistinfo%2Ffieldtrip&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863089093&sdata=M3y4f%2BGg7MpJXOr%2FF8DC9cTM0dUb2wS9hRyyEffxuRw%3D&reserved=0>https://doi.org/10.1371/journal.pcbi.1002202 <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1371%2Fjournal.pcbi.1002202&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863089093&sdata=1nMk28PGgmZtyX1M5Rl2wDkhBMv2vSXknpWYm0TGX9w%3D&reserved=0>
>
>
>
> --
> Prof. Dr.rer.nat. Carsten H. Wolters
> University of Münster
> Institute for Biomagnetism and Biosignalanalysis
> Malmedyweg 15
> 48149 Münster, Germany
>
> Phone:
> +49 (0)251 83 56904
> +49 (0)251 83 56865 (secr.)
>
> Fax:
> +49 (0)251 83 56874
>
> Email: carsten.wolters at uni-muenster.de
> Web: https://campus.uni-muenster.de/biomag/das-institut/mitarbeiter/carsten-wolters/ <https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fcampus.uni-muenster.de%2Fbiomag%2Fdas-institut%2Fmitarbeiter%2Fcarsten-wolters%2F&data=02%7C01%7Crdm146%40newark.rutgers.edu%7Cfd39adb746104deddfa808d63da91b81%7Cb92d2b234d35447093ff69aca6632ffe%7C1%7C0%7C636764193863099102&sdata=Xuoiy%2FV4Vfti4QM2bSxaitn%2FxnJeG3vG8UaWnd5XMXI%3D&reserved=0>
>
>
>
> _______________________________________________
> fieldtrip mailing listhttps://mailman.science.ru.nl/mailman/listinfo/fieldtriphttps://doi.org/10.1371/journal.pcbi.1002202
>
>
>
> --
> Prof. Dr.rer.nat. Carsten H. Wolters
> University of Münster
> Institute for Biomagnetism and Biosignalanalysis
> Malmedyweg 15
> 48149 Münster, Germany
>
> Phone:
> +49 (0)251 83 56904
> +49 (0)251 83 56865 (secr.)
>
> Fax:
> +49 (0)251 83 56874
>
> Email: carsten.wolters at uni-muenster.de
> Web: https://campus.uni-muenster.de/biomag/das-institut/mitarbeiter/carsten-wolters/
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20181106/1be638d5/attachment.html>
-------------- next part --------------
function condtensor = sb_assiTensorCond(mask,nodes,elem,labels,condcell)
%%
condtensor = zeros(9,size(elem,1));
count = 0;
for i = 1 : size(elem)
pos = nodes(elem(i,7),:);
pos = round(pos);
if (~(mask.anatomy(pos(1),pos(2),pos(3)) == 0))
if labels(i) ~= 5 && labels(i) ~= 6
labels(i)
disp('not 5 or 6')
end
count = count+1;
for j = 1 : 3
for k = 1 : 3
condtensor((j-1)*3 + k,i) = condcell{pos(1),pos(2),pos(3)}(j,k);
end
end
end
end
count
end
-------------- next part --------------
function [condtensor, s, fail] = sb_calcTensorCond_tuch(cfg,mask,V1,V2,V3,L1,L2,L3)
check1 = isequal(V1.dim,V2.dim,V3.dim);
check2 = isequal(V1.dim,size(V1.anatomy),size(V2.anatomy),size(V3.anatomy));
check3 = isequal(L1.dim,L2.dim,L3.dim);
check4 = isequal(L1.dim,size(L1.anatomy),size(L2.anatomy),size(L3.anatomy));
check5 = isequal(V1.dim(1:3),L1.dim,mask.dim);
check6 = ~(isempty(cfg)||isempty(cfg.conductivity));
check7 = ~(length(cfg.conductivity)<6);
if (check1 && check2 && check3 && check4 && check5 && check6 && check7)
fail = 0;
failnan = 0;
condtensor = cell(mask.dim);
N = zeros(1,3);
vol = zeros(1,3);
for i = 1 : V1.dim(1)
for j = 1 : V1.dim(2)
for k = 1 : V1.dim(3)
if (((mask.anatomy(i,j,k) == 5)|(mask.anatomy(i,j,k) == 6))&(L1.anatomy(i,j,k)>10^-4)&(L2.anatomy(i,j,k)>10^-4)&(L3.anatomy(i,j,k)>10^-5))
S = zeros(3);
S(:,1) = V1.anatomy(i,j,k,:);
S(:,2) = V2.anatomy(i,j,k,:);
S(:,3) = V3.anatomy(i,j,k,:);
if(norm(S'*S-diag([1,1,1]),2)>10e-7)
S(:,1) = S(:,1) / norm(S(:,1));
S(:,2) = S(:,2) - (S(:,1)'*S(:,2))*S(:,1);
S(:,2) = S(:,2) / norm(S(:,2));
S(:,3) = S(:,3) - (S(:,1)'*S(:,3))*S(:,1) - (S(:,2)'*S(:,3))*S(:,2);
S(:,3) = S(:,3) / norm(S(:,3));
failnan = failnan + 1;
end
if(sum(sum(isnan(S),1),2)>0)
condtensor{i,j,k} = cfg.conductivity(mask.anatomy(i,j,k))*diag([1,1,1]);
else
D = diag([L1.anatomy(i,j,k),L2.anatomy(i,j,k),L3.anatomy(i,j,k)]);
T = S * D * S';
%T = hier skalieren
condtensor{i,j,k} = T;
if (mask.anatomy(i,j,k) == 5)
N(1) = N(1) + 1;
vol(1) = vol(1) + L1.anatomy(i,j,k)*L2.anatomy(i,j,k)*L3.anatomy(i,j,k);
elseif (mask.anatomy(i,j,k) == 6)
N(2) = N(2) + 1;
vol(2) = vol(2) + L1.anatomy(i,j,k)*L2.anatomy(i,j,k)*L3.anatomy(i,j,k);
% elseif (mask.anatomy(i,j,k) == 3)
% N(3) = N(3) + 1;
% vol(3) = vol(3) + L1.anatomy(i,j,k)*L2.anatomy(i,j,k)*L3.anatomy(i,j,k);
end
end
elseif(((mask.anatomy(i,j,k) == 5)|(mask.anatomy(i,j,k) == 6)))
condtensor{i,j,k} = cfg.conductivity(mask.anatomy(i,j,k))*diag([1,1,1]);
if(((mask.anatomy(i,j,k) == 5)|(mask.anatomy(i,j,k) == 6))&(~((L1.anatomy(i,j,k)>10^-4)&(L2.anatomy(i,j,k)>10^-4)&(L3.anatomy(i,j,k)>10^-5))))
fail = fail + 1;
end
else
condtensor{i,j,k} = zeros(3);
end
end
end
end
d(1) = vol(1) / N(1);
d(1) = d(1)^(1/3);
d(2) = vol(2) / N(2);
d(2) = d(2)^(1/3);
% d(3) = vol(3) / N(3);
% d(3) = d(3)^(1/3);
s = d(1)*cfg.conductivity(5)+d(2)*cfg.conductivity(6);
s = s / (d(1)^2 + d(2)^2);
fprintf('s*d = %.6f\n',s*(d(1)+d(2)))
fprintf('s = %.6f\n',s)
% failper = fail / (N(1)+N(2));
mx = -100000;
mn = 100000;
for i = 1 : V1.dim(1)
for j = 1 : V1.dim(2)
for k = 1 : V1.dim(3)
if (((mask.anatomy(i,j,k) == 5)||(mask.anatomy(i,j,k) == 6))&(L2.anatomy(i,j,k)>10^-4)&(L3.anatomy(i,j,k)>10^-5))
condtensor{i,j,k} = s * condtensor{i,j,k};
%keep outliers wma bigger from the highest cond around cond of wm
if max(condtensor{i,j,k}(:)) > cfg.conductivity(4)
condtensor{i,j,k}(1:1+size(condtensor{i,j,k},1):end) = cfg.conductivity(6);
end
if mx < max(max(condtensor{i,j,k}))
mx = max(max(condtensor{i,j,k}));
fprintf('mx %d, %d, %d\n',i,j,k)
end
if mn > min(min(condtensor{i,j,k}))
mn = min(min(condtensor{i,j,k}));
fprintf('mn %d, %d, %d\n',i,j,k)
end
end
end
end
end
mn
mx
end
end
More information about the fieldtrip
mailing list