Upload New File

BiobankQSM.m

+ 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