Skip to content
This repository has been archived by the owner on May 18, 2021. It is now read-only.

Première version corrigé SWFM - fonctionnel #171

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Première version corrigé SWFM - fonctionnel #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
2fb3312
Première publication de l'algorithme SWFM (prototype initial):
alexandrebiomed Nov 9, 2020
e6f941e
Première version corrigée 1.2
alexandrebiomed Nov 13, 2020
28feb9a
Algo SWFM 1.3
alexandrebiomed Nov 14, 2020
9802b44
added a handedness variable to switch between left- and right-handed …
evaalonsoortiz Sep 29, 2020
676712c
Algo SWFM 1.3
alexandrebiomed Nov 14, 2020
11bffbd
Merge branch 'ab/SWFM' of https://github.com/shimming-toolbox/shimmin…
alexandrebiomed Nov 17, 2020
0beebb4
Trying to remove .ssh from git history
alexandrebiomed Nov 17, 2020
0735122
SWFM sortie de la somme et correction bruit
alexandrebiomed Nov 18, 2020
ff059e8
Algo SWFM fonctionnel avec fct atan 4 quad modifiée
alexandrebiomed Nov 24, 2020
4c2ff1f
Algo SWFM fonctionnel avec fct atan 4 quad modifiée et sup fichiers i…
alexandrebiomed Nov 24, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 41 additions & 0 deletions +imutils/+b0/B0mapping_sim_SWFM_Formulecorrige.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
% generate a "Shepp-Logan" phantom (ie; shape) in 2D
SheppLogan2d_vol = NumericalModel('Shepp-Logan2d',256);
figure; imagesc(SheppLogan2d_vol.starting_volume)
% generate a deltaB0 (B0 inhomogeneity) distribution in 2D (it will be in
% units of Hz)
SheppLogan2d_vol.generate_deltaB0('2d_linearIP', [5 0]);
% plot it to see how it looks
figure; imagesc(SheppLogan2d_vol.deltaB0)

% simulate MRI magnitude and phase data for your Shepp-Logan phantom in the
% presence of the deltaB0 field you generated
% indice = 1 ;
deltaTE=0.0001;
TE=[0.00100 0.00115 0.00120 0.00125 0.00130 0.00135 0.00140]
% for t=0:1
% TE(indice) = 0.001 + t*0.0005;
% indice = indice + 1 ;
% end
% echo times
SheppLogan2d_vol.simulate_measurement(15, TE, 100);

% get magnitude and phase data
magn = SheppLogan2d_vol.getMagnitude;
phase = SheppLogan2d_vol.getPhase;

% save the magnitude and phase data as a nifti file
SheppLogan2d_vol.save('Phase', 'SheppLogan2d_simulated_phase.nii');
SheppLogan2d_vol.save('Magnitude', 'SheppLogan2d_simulated_magnitude.nii');

% generate a complex data set from your magnitude and phase data
compl_vol = magn.*exp(1i*phase);

% calculate the deltaB0 map from the magnitude and phase data using the
% SWFM method (it will be in units of Hz)
[SWFM_delf] = +imutils.b0.SWFM_Formulecorrige(compl_vol(:,:,:,:), TE);
% plot it, it should look the same as the deltaB0 you simulated!
figure('Name','Champ magn�tique inhomog�ne'); imagesc(SWFM_delf)

% save your deltaB0 map as a nifti file
nii_vol = make_nii(SWFM_delf);
save_nii(nii_vol, ['dualechoB0_SheppLogan2d' '.nii'])
34 changes: 34 additions & 0 deletions +imutils/+b0/B0mapping_sim_dual_echo.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
% generate a "Shepp-Logan" phantom (ie; shape) in 2D
SheppLogan2d_vol = NumericalModel('Shepp-Logan2d',256);
figure; imagesc(SheppLogan2d_vol.starting_volume)
% generate a deltaB0 (B0 inhomogeneity) distribution in 2D (it will be in
% units of Hz)
SheppLogan2d_vol.generate_deltaB0('2d_linearIP', [5 0]);
% plot it to see how it looks
figure; imagesc(SheppLogan2d_vol.deltaB0)

% simulate MRI magnitude and phase data for your Shepp-Logan phantom in the
% presence of the deltaB0 field you generated
TE = [0.001 0.0015 0.0020 0.0025 0.0030 ]; % echo times
SheppLogan2d_vol.simulate_measurement(15, TE, 100);

% get magnitude and phase data
magn = SheppLogan2d_vol.getMagnitude;
phase = SheppLogan2d_vol.getPhase;

% save the magnitude and phase data as a nifti file
SheppLogan2d_vol.save('Phase', 'SheppLogan2d_simulated_phase.nii');
SheppLogan2d_vol.save('Magnitude', 'SheppLogan2d_simulated_magnitude.nii');

% generate a complex data set from your magnitude and phase data
compl_vol = magn.*exp(1i*phase);

% calculate the deltaB0 map from the magnitude and phase data using the
% SWFM method (it will be in units of Hz)
[dual_echo_delf] = +imutils.b0.dual_echo(compl_vol(:,:,:,:), TE);
% plot it, it should look the same as the deltaB0 you simulated!
figure; imagesc(dual_echo_delf)

% save your deltaB0 map as a nifti file
nii_vol = make_nii(dual_echo_delf);
save_nii(nii_vol, ['dualechoB0_SheppLogan2d' '.nii'])
58 changes: 58 additions & 0 deletions +imutils/+b0/Csigma.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
function sigma = Csigma(dataVol)
% Calculate background noise for each slice of an image by computing the
% average standard deviation of the signal in four (5 x 5 voxel) corners of
% an input image (which would presumably only contain noise) across all
% echoes
%
% _SYNTAX_
%
% sigma = bkgrnd_noise(dataVol)
%
% _DESCRIPTION_
%
% _INPUT ARGUMENTS_
%
% dataVol
% 4D data set dataVol(x,y,z,nEcho)
%
% _OUTPUTS_
%
% sigma
% standard deviation of the noise for each echo in 'dataVol'


% "default" noise level...


% Prepare noise mask
noiseMask = zeros(size(dataVol,1),size(dataVol,2),size(dataVol,3));

% Pick four corners of 5x5
noiseMask(1:5,1:5,:) = 1;
noiseMask(1:5,end-5:end,:) = 1;
noiseMask(end-5:end,1:5,:) = 1;
noiseMask(end-5:end,end-5:end,:) = 1;

% apply mask to the data
for echo = 1:size(dataVol,4)
noiseData(:,:,:,echo) = dataVol(:,:,:,echo).*noiseMask(:,:,:);
end

%noiseData = reshape(noiseData, size(dataVol,1)*size(dataVol,2), size(dataVol,3), size(dataVol,4));
stdNoise = zeros(1, size(dataVol,4) );

for echo = 1:size(dataVol,4)
if size(dataVol,3)==1
stdNoise(echo) = std(squeeze(noiseData(:,:,:,echo)),0,'all');
else
stdNoise(echo) = std(noiseData(:,:,:,echo),0,'all');
end
end
sigma=stdNoise


% SNR = zeros(1, size(dataVol,4) );
% for echo = 1:size(dataVol,4)
% SNR(echo) = snr(dataVol(:,:,:,echo),noiseData(:,:,:,echo));
% end
% SNR
155 changes: 155 additions & 0 deletions +imutils/+b0/SWFM_Formulecorrige.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
function [b0] = SWFM_Formulecorrige(complVol, TE)
% SWFM computes c0 fieldmaps based on "Selective-Weighted-Field-Mapping"
% method.
%
% _SYNTAX_
%
% [b0] = SWFM(complVol, TE)
%
% _DESCRIPTION_
%
% _INPUT ARGUMENTS_
%
% complVol
% complex 4D data set complVol(x,y,z,t)
% TE
% array of TEs in [s]
%
% _OUTPUTS_
%
% b0
% field map in units of Hz
%

gamma = 267.52218744 * 10^6; % rad*Hz/Tesla
factor= 1./(gamma.*pi.*TE(end));

% create magnitude and phase data volumes
mag_data = abs(complVol);
ph_data = angle(complVol);


X=zeros(size(complVol));
for t=1:length(TE)
X(:,:,:,t)= mag_data(:,:,:,t).*exp(1i*ph_data(:,:,:,t));
end


% SDQWF correspond au "Signal Decay Quality Weighting Factor"
% FM repr�sente la fraction de la magnitude de l'�cho j par rapport � la
% somme des magnitudes.
% FB repr�sente la fraction inverse du bruit de l'�cho j par rapport � la
% somme des bruits du signal.

% Calcul du facteur FM
magtotal=zeros( size(mag_data,1), size(mag_data,2), size(mag_data,3) );
FM=zeros(size(mag_data));
for t=1:length(TE)
magtotal(:,:,:)=magtotal(:,:,:)+mag_data(:,:,:,t);
end
for t=1:length(TE)
FM(:,:,:,t)=mag_data(:,:,:,t)./magtotal;
end

% Calcul du facteur FB pr�sent dans la fonction arctangente � 4 cadrants

sigmatot=0;
FB=zeros(1,size(mag_data,4));
sigma = imutils.b0.Csigma(mag_data);
for t=1:length(TE)
sigmatot = sigmatot + sigma(t);
end
for t =1:length(TE)
FB(t)=sigmatot./sigma(t);
end


% Calcul du facteur SDQWF
SDQWF=zeros(size(mag_data));

for t=1:length(TE)
SDQWF(:,:,:,t)=FM(:,:,:,t).*FB(t);
end

% Calcul du facteur produit X1Xj pr�sent dans la fonction arctangente � 4 cadrants
X1X_T=zeros(size(complVol));
deltaphi=zeros(size(complVol));
for t=1:length(TE)
deltaphi(:,:,:,t)=ph_data(:,:,:,1)-ph_data(:,:,:,t);
X1X_T(:,:,:,t)= mag_data(:,:,:,1).*mag_data(:,:,:,t).*exp(1i.*(-deltaphi(:,:,:,t)));
end

% Calcul du facteur somme pr�sent dans la fonction arctangente � 4 cadrants

% cond1 = zeros( size(complVol) );
% cond2 = zeros( size(complVol) );
% cond3 = zeros( size(complVol) );
% cond4 = zeros( size(complVol) );
% cond5 = zeros( size(complVol) );
%
% for t = 1:length(TE)
% cond1(:,:,:,t) = real(X1X_T(:,:,:,t))>0 ;
% cond2(:,:,:,t) = (real(X1X_T(:,:,:,t))<0 & imag(X1X_T(:,:,:,t))~=0) ;
% cond3(:,:,:,t) = (real(X1X_T(:,:,:,t))>0 & imag(X1X_T(:,:,:,t))==0) ;
% cond4(:,:,:,t) = (real(X1X_T(:,:,:,t))<0 & imag(X1X_T(:,:,:,t))==0) ;
% cond5(:,:,:,t) = (real(X1X_T(:,:,:,t))==0 & imag(X1X_T(:,:,:,t))==0) ;
% end
% for t = 1 : length(TE)
% for i= 1:size(complVol,1)
% for j= 1:size(complVol,2)
% for k= 1: size(complVol,3)
% if cond1(i,j,k,t)==1
% sommable1(i,j,k,t) = (imag(X1X_T(i,j,k,t)).*SDQWF(i,j,k,t))./(sqrt( (imag(X1X_T(i,j,k,t))).^2 + (real(X1X_T(i,j,k,t))).^2 ) +real(X1X_T(i,j,k,t)) );
% elseif cond2(i,j,k,t)==1
% sommable2(i,j,k,t) = ((sqrt((imag( X1X_T(i,j,k,t) ) ).^2 + (real(X1X_T(i,j,k,t))).^2 ) - real(X1X_T(i,j,k,t))) .*SDQWF(i,j,k,t))./(imag(X1X_T(i,j,k,t)));
% end
% end
% end
% end
% end

complexe_imag = zeros( size(complVol) );
complexe_reel = zeros( size(complVol) );

for t = 1:length(TE)
complexe_reel(:,:,:,t) = real(X1X_T(:,:,:,t)) ;
complexe_imag(:,:,:,t) = imag(X1X_T(:,:,:,t)) ;
end

b0=zeros( size(complVol,1) ,size(complVol,2) ,size(complVol,3));
phase=zeros( size(complVol));
indice=1;
for echo = 2 : length(TE)
phase(:,:,:,indice)= +imutils.b0.arctan4quadrant(complexe_imag(:,:,:,echo),complexe_reel(:,:,:,echo),SDQWF(:,:,:,echo));
b0(:,:,:) = b0(:,:,:) + phase(:,:,:,indice);
indice = indice +1;
end
b0=factor.*b0;
b0=b0/(2*pi);
% Sans la fonction arctan4quadrant qui est la fonction arctangente � 4 quadrants
% modifi�e, on aurait obtenu l'algorithme suivant.

% b0=zeros( size(complVol,1) ,size(complVol,2) ,size(complVol,3));
% for t = 2:length(TE)
% for i= 1:size(complVol,1)
% for j= 1:size(complVol,2)
% for k= 1: size(complVol,3)
% if cond1(i,j,k,t)==1
% b0(i,j,k) = b0(i,j,k) + atan(sommable1(i,j,k,t));
% elseif cond2(i,j,k,t)==1
% b0(i,j,k) = b0(i,j,k) + atan(sommable2(i,j,k,t));
% elseif cond3(i,j,k,t)==1
% b0(i,j,k) = b0(i,j,k) + pi/2;
% elseif cond4(i,j,k,t)==1
% b0(i,j,k) = b0(i,j,k) + (-pi/2);
% elseif cond5(i,j,k,t)==1
% b0(i,j,k) = b0(i,j,k);
% end
% end
% end
% end
% end
% b0(:,:,:)= factor.*b0(:,:,:);
% b0 = b0./(2*pi); % [Hz]


26 changes: 26 additions & 0 deletions +imutils/+b0/arctan4quadrant.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
function theta = arctan4quadrant(y,x,SDQWF)
% Fonction arctangente modifi�e et adapt�e au calcul du champ magn�tique
% dans la m�thode SWFM - Selective Weight Field Mapping.
% Theta correspond � la phase.
% y , x et complexdata repr�sentent les donn�es complexes obtenues par
% simulation ou par IRM. SDQWF correspond au Signal Decay Quality Weight Factor.

for i = 1 : size(x,1)
for j = 1 : size(x,2)
for k = 1 : size(x,3)
if x(i,j,k)> 0
theta(i,j,k) = atan( (y(i,j,k).*SDQWF(i,j,k))./(sqrt( y(i,j,k)^2+x(i,j,k)^2 )+ x(i,j,k)));
elseif x(i,j,k)<0 && y(i,j,k)~=0
theta(i,j,k)= atan( ((sqrt( y(i,j,k)^2+x(i,j,k)^2 )- x(i,j,k)).*SDQWF(i,j,k))./y(i,j,k));
elseif x(i,j,k) > 0 && y(i,j,k)==0
theta(i,j,k) = pi/2;
elseif x(i,j,k) < 0 && y(i,j,k)==0
theta(i,j,k) = -pi/2;
elseif x(i,j,k) == 0 && y(i,j,k)==0
theta(i,j,k)=0;
end
end
end
end


16 changes: 0 additions & 16 deletions .gitignore
View file
Open in desktop

This file was deleted.

Binary file added dualechoB0_SheppLogan2d.nii
View file
Open in desktop
Binary file not shown.
11 changes: 0 additions & 11 deletions info.xml
View file
Open in desktop

This file was deleted.