Commit 2a3a5426 authored by Ludovica Griffanti's avatar Ludovica Griffanti
Browse files

Upload QSM

parent 97ca8556
function BiobankQSM(path,TE1,TE2)
% Description: Main QSM processing pipeline
%
% Authors: Chaoyue Wang, Benjamin C. Tendler & Karla L. Miller
%
% Copyright 2021 University of Oxford
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
% path - swMRI data directory
% TE1, TE2 - two echo times
% load data
addpath(genpath('./STISuite_V3.0/'));
fn_pha1 = dir([path '/SWI/PHA_TE1/*.nii.gz']);
num_channels = length(fn_pha1);
for j = 1:num_channels
phs1(:, :, :, j) = double(niftiread([path '/SWI/PHA_TE1/' fn_pha1(j).name]));
end
dim = size(phs1);
M1 = zeros(dim);
fn_mag1 = dir([path '/SWI/MAG_TE1/*.nii.gz']);
for j = 1:num_channels
M1(:, :, :, j) = double(niftiread([path '/SWI/MAG_TE1/' fn_mag1(j).name]));
end
phase1 = ((phs1 - 2048) / 2048) * pi;
S1 = M1 .* exp(1i * phase1);
clear M1 phase1 phs1 fn_pha1 fn_mag1;
phs2 = zeros(dim);
fn_pha2 = dir([path '/SWI/PHA_TE2/*.nii.gz']);
for j = 1:num_channels
phs2(:, :, :, j)=double(niftiread([path '/SWI/PHA_TE2/' fn_pha2(j).name]));
end
M2 = zeros(dim);
fn_mag2 = dir([path '/SWI/MAG_TE2/*.nii.gz']);
for j = 1:num_channels
M2(:, :, :, j) = double(niftiread([path '/SWI/MAG_TE2/' fn_mag2(j).name]));
end
phase2 = ((phs2 - 2048) / 2048) * pi;
clear phs2 j fn_mag2 fn_pha2;
S2 = M2 .* exp(1i * phase2);
clear M2 phase2;
% Coil comb:MCPC-3D-S
hip = zeros(dim(1:3));
for iCha = 1:size(S1, 4)
complexDifference = S2(:, :, :, iCha) .* conj(S1(:, :, :, iCha));
hip = hip + complexDifference;
end
clear complexDifference;
hip_abs = abs(hip);
hip_angle = angle(hip);
nii_info = niftiinfo([path '/SWI/SWI_TOTAL_MAG_brain_mask.nii.gz']);
nii_info.Datatype = 'double';
niftiwrite(hip_abs, [path '/SWI/QSM/hip_abs.nii'], nii_info, ...
'Compressed', true);
niftiwrite(hip_angle, [path '/SWI/QSM/hip_angle.nii'], nii_info,...
'Compressed', true);
%prelude unwrapping
setenv('FSLOUTPUTTYPE', 'NIFTI_GZ');
FSLDIR = getenv('FSLDIR');
system([FSLDIR, '/bin/fslmaths ', path, ...
'/SWI/SWI_TOTAL_MAG_brain_mask.nii.gz -thr 0.6 -bin ', path, ...
'/SWI/QSM/QSM_mask.nii.gz']);
system([FSLDIR, '/bin/prelude -a ', path, '/SWI/QSM/hip_abs.nii.gz -p ', ...
path, '/SWI/QSM/hip_angle.nii.gz -u ', path, ...
'/SWI/QSM/hip_uw.nii.gz -m ', path, '/SWI/QSM/QSM_mask.nii.gz']);
unwrappedHip=double(niftiread([path '/SWI/QSM/hip_uw.nii.gz']));
delete([path '/SWI/QSM/hip_abs.nii.gz']);
delete([path '/SWI/QSM/hip_angle.nii.gz']);
delete([path '/SWI/QSM/hip_uw.nii.gz']);
TEs(1)=str2double(TE1);
TEs(2)=str2double(TE2);
scale = TEs(1) / (TEs(2) - TEs(1));
unwrappedHip = unwrappedHip * scale;
hipComplex = exp(1i * unwrappedHip);
size_compl = size(hipComplex);
po = complex(zeros([size_compl(1:3) num_channels], 'double'));
for iCha = 1:num_channels
po_ang=angle(exp(1i * (angle(S1(:, :, :, iCha)) - unwrappedHip)));
po_double=double(abs(S1(:, :, :, iCha)) .* exp(1i * po_ang));
po(:, :, :, iCha) = po_double;
end
clear hip hip_abs hip_angle hipComplex po_ang po_double unwrappedHip scale;
real_smmothed=double(zeros(size(po)));
imag_smmothed=double(zeros(size(po)));
for j=1:num_channels
real_smmothed(:, :, :, j) = imgaussfilt3(real(po(:, :, :, j)), 4, ...
'FilterDomain', 'spatial');
imag_smmothed(:, :, :, j) = imgaussfilt3(imag(po(:, :, :, j)), 4, ...
'FilterDomain', 'spatial');
end
clear po;
po_c = (real_smmothed + 1i * imag_smmothed) ./ abs(real_smmothed + 1i * imag_smmothed);
clear real_smmothed imag_smmothed;
S1_c = S1 .* squeeze(conj(po_c));
clear S1;
combined1 = weightedCombination(S1_c, abs(S1_c));
combined1(~isfinite(combined1)) = 0;
mask = double(niftiread([path '/SWI/QSM/QSM_mask.nii.gz']));
nii = angle(combined1) .* mask;
nii_info = niftiinfo([path '/SWI/SWI_TOTAL_MAG_brain_mask.nii.gz']);
nii_info.Datatype = 'single';
niftiwrite(single(nii), [path '/SWI/QSM/PHASE_TE1.nii'], nii_info, ...
'Compressed',true);
clear nii S1_c;
S2_c = S2 .* squeeze(conj(po_c));
clear S2 po_c;
combined2 = weightedCombination(S2_c, abs(S2_c));
clear S1_c S2_c;
combined2(~isfinite(combined2)) = 0;
nii = angle(combined2) .* mask;
niftiwrite(single(nii), [path '/SWI/QSM/PHASE_TE2.nii'], nii_info, ...
'Compressed',true);
clear hip hip_abs hip_angle nii po_c S1_c S2_c smoothed_weight ;
clear unwrappedHip ans iCha j num_channels size_compl dim;
% extract DICOM info and update mask
phase1 = mask .* angle(combined1);
phase2 = mask .* angle(combined2);
clear combined1 combined2;
dcm_info = dicominfo([path '/DICOM/SWI.dcm']);
ori = dcm_info.ImageOrientationPatient;
Xz = ori(3);
Yz = ori(6);
Zxyz = cross(ori(1:3), ori(4:6));
Zz = Zxyz(3);
H = [-Xz, -Yz, Zz];
voxsz = [dcm_info.PixelSpacing' dcm_info.SliceThickness];
B0 = dcm_info.MagneticFieldStrength;
clear ori Xz Yz Zxyz Zz;
map = phasevariance_nonlin(mask, phase1, 2);
map2 = phasevariance_nonlin(mask, phase2, 2);
dim = size(phase1);
mask(map < 0.6) = 0;
mask(map2 < 0.5) = 0;
for ii = 1:dim(3)
mask(:, :, ii) = bwareaopen(mask(:, :, ii), 300);
mask(:, :, ii) =~ bwareaopen(~mask(:, :, ii), 50);
end
mask = imfill(mask,26,'holes');
% phase unwrap
[uwphase1, ~] = MRPhaseUnwrap(phase1, 'voxelsize', voxsz, 'padsize', [64 64 64]);
[uwphase2, ~] = MRPhaseUnwrap(phase2, 'voxelsize', voxsz, 'padsize', [64 64 64]);
T2s = 40;
W1 = TEs(1) * exp(-TEs(1) / T2s);
W2 = TEs(2) * exp(-TEs(2) / T2s);
phs_comb = double((W1 * uwphase1 + W2 * uwphase2) / (W1 + W2));
TE = (W1 * TEs(1) + W2 * TEs(2)) / (W1 + W2);
clear phase1 uwphase1 uwphase2 W1 W2 TEs T2s combined phase2;
% v-SHARP and Dipole inversion
[dB_vsf, mask_vsf]=V_SHARP(phs_comb, mask, 'voxelsize', voxsz, 'smvsize', 12);
clear mask;
mask_vsf(map < 0.7) = 0;
mask_vsf(map2 < 0.6) = 0;
for ii = 1:dim(3)
mask_vsf(:, :, ii) = bwareaopen(mask_vsf(:, :, ii), 200);
mask_vsf(:, :, ii) =~ bwareaopen(~mask_vsf(:, :, ii), 30);
end
mask_vsf = imfill(mask_vsf,26,'holes');
clear map phs_comb map2;
qsm_iLSQR_vsf = QSM_iLSQR(dB_vsf, mask_vsf, 'TE', TE, 'B0', B0, ...
'H', H, 'padsize', [64 64 64], ...
'voxelsize', voxsz);
niftiwrite(single(qsm_iLSQR_vsf * 1000), [path '/SWI/QSM/QSM.nii'], ...
nii_info, 'Compressed', true);
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment