Ubuntu への AnyDesk のインストール

AnyDeskはマルチプラットフォームで動作するリモートデスクトップアプリケーションです。動きがとても軽くて非常に使いやすいです。

本家サイトに ‘AnyDesk DEB repository how-to‘ が紹介されていますが、root で作業するように書かれています。

Ubuntuで普通のユーザーが使う時はちょっとだけ違う方がすんなりインストールできますので紹介します。

  • リポジトリキーを信頼するソフトウェア・プロバイダー・リストに登録する
wget -qO - https://keys.anydesk.com/repos/DEB-GPG-KEY | sudo apt-key add -
  • リポジトリを追加する
sudo sh -c 'echo "deb http://deb.anydesk.com/ all main" > /etc/apt/sources.list.d/anydesk-stable.list'
  • リポジトリを更新し、AnyDeskをインストールする
sudo apt update
sudo apt install anydesk

【FSL】Gray matter-Based Spatial Statistics: GBSS



1. 目的
2. GBSSとは
2.1. データの準備
2.2. 解析パラメータの設定
2.3. gbss_1_T1wpreproc:T1WIの前処理
2.4. gbss_2_DWIwpreproc:拡散定量値の前処理
2.5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影
2.6. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)
3. おまけ


1. 目的

  • Gray matter-Based Spatial Statistics: GBSS

2. GBSSとは

Gray matter-Based Spatial Statistics(GBSS)は、灰白質の統計解析をするための手法。TBSSの灰白質版。

灰白質の中心線(skeleton)に定量値を投影する。通常の脳画像の統計解析では、脳構造の個人差を除外するために空間的「平滑化」を用いる。しかし、平滑化の程度に原則がなく、平滑化をかけては情報があいまいになり、MRIの高空間分解能を生かせないという問題がある。一方、GBSSでは、灰白質の中心線と思われるところにskeletonを生成し、そこに個人ごとの定量値を投影するという手法をとる。これにより、平滑化せずに群間比較をすることができるため、平滑化による問題を回避できるという利点がある。

GBSS解析では、次のような処理をする。

  1. データの準備
  2. 解析パラメータの設定
  3. gbss_1_T1wpreproc:T1WIの前処理
  4. gbss_2_DWIwpreproc:拡散定量値の前処理
  5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影
  6. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)

2.1. データの準備

まず、次のデータを準備する。

  • b0フォルダ:b0 (DWI, b=0 s/mm^2)
  • T1wフォルダ:3D-T1WI
  • 拡散定量値:FA, MD, FW, etc.

フォルダ構造は次のようにする。

ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。各フォルダに入れるファイル名は同じにする必要がある(例:FA/Subj001.nii.gz, T1w/Subj001.nii.gzは同じファイル名)。

.
├── FA
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── FW
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── MD
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── T1w
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
└── b0
    ├── Con0001.nii.gz
    ├── Con0002.nii.gz
    ├── ...
    └── Pat0010.nii.gz

2.2. 解析パラメータの設定

解析に用いるパラメータを設定する。

PEDIRECHOSPACINGepi_regコマンドのパラメータであり、それぞれ位相エンコード方向とecho spacing timeを指定する。SKELETON_THRはスケルトン化する際のしきい値である。

MAP_LIST=$(ls | grep -v b0 | grep -v T1w | grep -v stats)
SUBJ_LIST=$(ls T1w | cut -d . -f1)
PEDIR='-y'                        # Phase encoding direction: "epi_reg" config
ECHOSPACING='0.0380544'           # Echo spacing time (s): "epi_reg" config
SKELETON_THR=0.2                  # Skeleton threshold

2.3. gbss_1_T1wpreproc:T1WIの前処理

T1WIの前処理内容は、次の通り。

  1. 脳頭蓋除去
  2. 灰白質をセグメント
  3. 灰白質を標準脳(MNI152)に位置合わせ
  4. 標準空間上の各被験者の灰白質を平均化

以上の処理内容を実行するために、関数gbss_1_T1wpreprocを定義。

function gbss_1_T1wpreproc() {
    echo "T1w preprocessing..."
    # T1w preprocessing
    for SUBJ in ${SUBJ_LIST}; do
        ## Segment GM
        ### Skull-stripping
        bet T1w/${SUBJ} T1w/${SUBJ}_skull_stripped -f 0.3 -R -S -B
        ### Segment T1w into CSF/GM/WM
        fast -t 1 -g -B -b -p -o T1w/${SUBJ}_fast T1w/${SUBJ}_skull_stripped

        ## Register GM to MNI
        ### Register T1w to MNI
        flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
            -in T1w/${SUBJ}_fast_restore \
            -omat T1w/${SUBJ}_indiv2MNI.mat
        fnirt --in=T1w/${SUBJ}_fast_restore \
            --aff=T1w/${SUBJ}_indiv2MNI.mat \
            --cout=T1w/${SUBJ}_indiv2MNI_warp \
            --config=T1_2_MNI152_2mm
        ### applywarp to GM and move it into MNI
        applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
            --in=T1w/${SUBJ}_fast_pve_1 \
            --warp=T1w/${SUBJ}_indiv2MNI_warp \
            --out=T1w/${SUBJ}_fast_pve_1_inMNI &
    done
    ## Create 4D all_GM image
    ### GM_pve1 in MNI are merged into one volume
    mkdir stats
    fslmerge -t stats/all_GM $(imglob T1w/*_fast_pve_1_inMNI.nii.gz)
    ### Create merged GM_pve1 (all_GM) mask and apply it to all_GM
    fslmaths stats/all_GM -max 0 -Tmin -bin stats/mean_GM_mask -odt char
    fslmaths stats/all_GM -Tmean stats/mean_GM
}

関数が定義できたら、実行する。

gbss_1_T1wpreproc

処理が完了すると、statsフォルダにすべての被験者の灰白質画像がまとめられた「all_GM.nii.gz」とそれを平均化した「mean_GM.nii.gz」が生成される。

2.4. gbss_2_DWIwpreproc:拡散定量値の前処理

拡散定量値の前処理の過程は、次の通り。

  1. b0をT1WIに位置合わせ
  2. b0 to T1WIとT1WI to 標準脳のWarp Fieldを使って、拡散定量値マップを標準空間に移動
  3. 標準空間上にあるすべての被験者の拡散定量値マップをひとつの4D volume dataにまとめる
  4. 平均灰白質マスクmean_GM_maskで、3の拡散定量値マップをマスク処理
function gbss_2_DWIwpreproc() {
    echo "Diffusion MAP preprocessing..."
    # Diffusion MAP preprocessing
    for SUBJ in ${SUBJ_LIST}; do
        ## Register b0 to MNI
        ### Register b0 to T1WI
        epi_reg --epi=b0/${SUBJ} \
            --t1=T1w/${SUBJ} --t1brain=T1w/${SUBJ}_skull_stripped \
            --echospacing=${ECHOSPACING} --pedir=${PEDIR} \
            --wmseg=T1w//${SUBJ}_fast_seg_2.nii.gz \
            --out=b0/${SUBJ}_BBR
        for MAP in ${MAP_LIST}; do
            ### applywarp to diffusion map and move it to T1w based on epi_reg
            flirt -in ${MAP}/${SUBJ} \
                -out ${MAP}/${SUBJ}_inT1 \
                -ref T1w/${SUBJ}_skull_stripped \
                -init b0/${SUBJ}_BBR.mat -applyxfm
            ### applywarp to dMAP in T1w and move it to MNI
            applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
                --in=${MAP}/${SUBJ}_inT1 \
                --warp=T1w/${SUBJ}_indiv2MNI_warp \
                --out=${MAP}/${SUBJ}_inMNI
        done
    done
    ## 4D all_${dMAP}
    for MAP in ${MAP_LIST}; do
        ### Merge dMAP in MNI into one volume, all_${MAP}
        fslmerge -t stats/all_${MAP} $(imglob ${MAP}/*_inMNI.nii.gz)
        ### Masking dMAP in MNI using mean GM mask
        fslmaths stats/all_${MAP} -mas stats/mean_GM_mask stats/all_${MAP}
    done
}

関数が定義できたら、実行する。

gbss_2_DWIwpreproc

処理が完了すると、拡散定量値の種類に応じてstatsフォルダに「all_.nii.gz」が生成される。

以下は、標準空間上における灰白質領域のFA画像である。

2.5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影

拡散定量値をスケルトンに投影するまでの手順は、次の通り。

  1. 平均灰白質画像からスケルトンを生成
  2. 灰白質スケルトンをしきい値処理をしてスケルトンマスクを生成
  3. スケルトンの距離マップ(distance map)を生成
  4. 各拡散定量値マップを、スケルトンに投影
function gbss_3_skelpreproc() {
    echo "Projecting MAP into skeleton..."
    # Skeletonize dMAP
    ## Create GM skeleton
    tbss_skeleton -i stats/mean_GM -o stats/mean_GM_skeleton
    ## threshold GM skeleton
    fslmaths stats/mean_GM_skeleton -thr ${SKELETON_THR} -bin stats/mean_GM_skeleton_mask
    ## Create skeleton distancemap
    fslmaths stats/mean_GM_mask -mul -1 -add 1 -add stats/mean_GM_skeleton_mask stats/mean_GM_skeleton_mask_dst
    distancemap -i stats/mean_GM_skeleton_mask_dst -o stats/mean_GM_skeleton_mask_dst
    for MAP in ${MAP_LIST}; do
        ## Project masked dMAP in MNI into skeleton
        tbss_skeleton -i stats/mean_GM \
            -p ${SKELETON_THR} \
            stats/mean_GM_skeleton_mask_dst \
            ${FSLDIR}/data/standard/LowerCingulum_1mm \
            stats/all_${MAP} \
            stats/all_${MAP}_skeletonised
    done
}

処理が完了すると、statsフォルダに灰白質のスケルトンmean_GM_skeleton.nii.gz、灰白質のスケルトンの距離マップmean_GM_skeleton_mask_dst.nii.gz、拡散定量値が投影されたスケルトンall_<MAP>_skeletonised.nii.gzが保存される。

2.6. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)

まず、GLMのデザインマトリックス(計画行列)とコントラストを設定する。

今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。これらのファイルを昇順に並び替えると、先に健常者10名の画像、次に患者10名の画像の順に並ぶ。

Con0001.nii.gz
Con0002.nii.gz
Con0003.nii.gz
...
Pat0010.nii.gz

次に、GLMのデザインマトリックス(計画行列)とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。

design_ttest2 stats/design 10 10

statsフォルダに、デザインマトリックス(design.mat)とコントラスト(design.con)が生成される。

デザインマトリックス(design.mat)の中身を確認。

/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者ファイルの順番(昇順)に対応する。したがって、これらは対応があるようにしておかなければならない。

cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1

コントラスト(design.con)の中身を確認してみる。

/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。

cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1

デザインマトリックス(計画行列)とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定(permutation test)を実行する。

randomiseコマンドの各オプションは、次の通り。

  • -i:入力画像
  • -m:マスク画像
  • -o :出力画像
  • -o :デザインマトリックス
  • -o :デザインコントラスト
  • -n:並べ替え検定の数
  • –T2:2D最適化を用いたTFCE
  • -x:voxel-wiseのcorrected P値マップ
  • –uncorrp:un-corrected P値マップ
  • -R:統計値マップ
for MAP in ${MAP_LIST}; do
    randomise -i stats/all_${MAP}_skeletonised \
        -o stats/gbss_${MAP} \
        -m stats/mean_GM_mask \
        -d stats/design.mat \
        -t stats/design.con \
        -n 10000 --T2 -V -x --uncorrp -R
done

次のようなファイルが、生成される。

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(gbss_FA_tfce_corrp_tstat1.nii.gz)を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい

fsleyes ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
    stats/mean_GM_skeleton -cm Green \
    stats/gbss_FA_tfce_p_tstat1 -cm Red-Yellow -dr 0.95 1

スケルトンは細いため有意差が見づらい場合がある。そのような時、tbss_fillコマンドが役に立つ。

tbss_fillコマンドの基本的な使い方は、以下の通り。

tbss_fill <P値マップ> <しきい値> <平均FA画像> <出力画像>

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(gbss_FA_tfce_corrp_tstat1.nii.gz)の有意差があった領域のみを0.95でしきい値処理をして抽出し、その領域を膨張させる。

tbss_fill stats/gbss_FA_tfce_p_tstat1 0.95 stats/mean_GM stats/gbss_FA_tfce_p_tstat1_fill

赤く表示されている領域は、健常群が患者群よりも有意(FWE-corrected)にFA値が大きいことを示している。

3. おまけ

大量に定量値があり、それらをすべて検定する場合、有意差があったかどうかをすべて確認するのは大変である。そこで、各定量値画像のP値マップが0.95以上の値を持つかどうかを判定し、有意差があった場合のみ、tbss_fillコマンドを実行する。

for PMAP in $(ls stats/ | grep tfce_corrp); do
    PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
    echo ${PMAP} >>stats/tmp1.txt
    echo ${PMAX} >>stats/tmp2.txt
    if [ $(echo "${PMAX} > 0.95" | bc) == 1 ]; then
        tbss_fill stats/${PMAP} 0.95 stats/mean_GM stats/${PMAP}_fill
    fi
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
    >stats/corrected_P_report.txt
rm stats/tmp*

上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。

検定結果を、ファイル名でソート(上段)したものと、P値でソートしたもの(下段)に分けて保存している。

cat stats/corrected_P_report.txt
gbss_FA_tfce_corrp_tstat1.nii.gz        0.992000
gbss_FA_tfce_corrp_tstat2.nii.gz        0.416000
gbss_FW_tfce_corrp_tstat1.nii.gz        0.361839
gbss_FW_tfce_corrp_tstat2.nii.gz        0.997261
gbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
gbss_MD_tfce_corrp_tstat2.nii.gz    0.985748


gbss_FW_tfce_corrp_tstat2.nii.gz        0.997261
gbss_FA_tfce_corrp_tstat1.nii.gz        0.992000
gbss_MD_tfce_corrp_tstat2.nii.gz    0.985748
gbss_FA_tfce_corrp_tstat2.nii.gz        0.416000
gbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
tbss_FW_tfce_corrp_tstat1.nii.gz        0.361839

【QSM】QSM解析 ~MEDIの使い方~


1. 目的
2. 準備
2.1. MEDI toolboxにパスを通す
2.2. QSM解析実行のコード準備
2.3. Write_DICOM.mを修正
2.4. QSM_3Dを集める
3. 実行
3.1. 作業フォルダ(analyzeフォルダ)へ移動
3.2. runme.mを実行
3.3. Outputファイルの確認
4. 複数の被験者データをまとめて処理したい場合
4.1. run_sequential.m


1. 目的

Cornell大学が開発したQSM software (MEDI)を使って解析

(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)

2. 準備

####MEDI toolboxのダウンロード
こちらからMEDI_toolboxをダウンロード
MEDI_toolboxの中身はこのような感じ。

2.1. MEDI toolboxにパスを通す

MATLABを起動し、パスの設定からサブフォルダも追加をクリックしfunctionsフォルダを選択した後、保存。

2.2. QSM解析実行のコード準備

QSM解析を実行するためのコードは以下。
以下をrunme.m(※任意の名前でOK)で保存し、作業フォルダとなるanalyzeフォルダ(※任意の名前でOK)に保存する。

clear all;clc;close all;

% Set path
% MEDI_set_path

% STEP 1: Import data
[iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('AXL_QSM');

%%%% Generate the Magnitude image %%%%
iMag = sqrt(sum(abs(iField).^2,4));

%%%%% provide a Mask here if possible %%%%%%
if (~exist('Mask','var'))                     
%     Mask = genMask(iField, voxel_size);
    Mask = BET(iMag,matrix_size,voxel_size);
end

%%%%% provide a noise_level here if possible %%%%%%
if (~exist('noise_level','var'))
    noise_level = calfieldnoise(iField, Mask);
end

%%%%% normalize signal intensity by noise to get SNR %%%
iField = iField/noise_level;

% STEP 2a: Field Map Estimation
%%%%%Estimate the frequency offset in each of the voxel using a 
%%%%%complex fitting %%%%
[iFreq_raw N_std] = Fit_ppm_complex(iField);

% STEP 2b: Spatial phase unwrapping %%%%
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);

% STEP 2c: Background Field Removal
%%%% Background field removal %%%%
[RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);

% CSF Mask for zero referencing
R2s=arlo(TE,abs(iField));
Mask_CSF = extract_CSF(R2s, Mask, voxel_size);

% STEP 3: Dipole Inversion
save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
     voxel_size delta_TE CF B0_dir Mask_CSF;

%%%% run MEDI %%%%%
QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');

%%% save to DICOM
% ignore warnings...
Write_DICOM(QSM,files,'QSM')

2.3. Write_DICOM.mを修正

解析後のQSMの保存には、Write_DICOM.mを使用する。(runme.mを実行すれば自動で実行される。)

修正前のファイルでは、dicomwriteでエラーが発生する。
そのため、MEDI_toolbox/functions/Write_DICOM.mを以下のコードにそっくりそのまま書き換える。
やっていることは、dicomwriteのCreateModeをCopyからCreateへチェンジ。

function Write_DICOM(M,files,outdir,opts)
%WRITE_DICOM Summary of this function goes here
%   Detailed explanation goes here


defopts.SeriesDescription = 'QSM';
defopts.SeriesInstanceUID = [];
defopts.SeriesNumber = [];
defopts.Convert = @(x) convert(x);
defopts.Window = 0.500;
defopts.Level = 0;
defopts.flag3D = 0;
% defopts.EchoNumber = [];
% defopts.EchoTime = 0.0;
        
if nargin<4; opts=struct; end
deffields=fieldnames(defopts);
for i=1:length(deffields)
    if ~isfield(opts, deffields{i})
        opts.(deffields{i})=defopts.(deffields{i});
    end
end

if isfield(files,'M')
    filenames=files.M;
elseif isfield(files,'R')
    filenames=files.R;
elseif isfield(files, 'info3D')
    opts.flag3D=1;
else
    error('No filenames (M nor R) found.');
end
    
flag_signed=min(M(:))<0;

if ~opts.flag3D
    if size(M,3) ~= size(filenames,1)
        error([num2str(size(filenames,1)) ' filenames given for ' num2str(size(M,3)) ' slices.']);
    end
end

if isempty(opts.SeriesInstanceUID)
   opts.SeriesInstanceUID=dicomuid;
end
progress='';

mkdir(outdir);

warning('off','images:dicom_add_attr:invalidAttribChar');
load_flag=1;
insert_flag=~opts.flag3D;
for slice=1:size(M,3)
    if load_flag
        if opts.flag3D
            filename=files.info3D.filename;
        else
            filename=filenames{slice,end};
        end
        info = dicominfo(filename);
        info.SeriesDescription = opts.SeriesDescription;
        info.SeriesInstanceUID = opts.SeriesInstanceUID;
        if isempty(opts.SeriesNumber)
            opts.SeriesNumber=info.SeriesNumber*100;
        end
        info.SeriesNumber = opts.SeriesNumber;
        info.SOPInstanceUID = dicomuid;
        info.InstanceNumber = slice;
        if opts.flag3D
            load_flag=0;
        end
    end
    if opts.flag3D
        item=files.info3D.slice2index{slice};
%         info.PerFrameFunctionalGroupsSequence.(item).PlanePositionSequence.Item_1.ImagePositionPatient;
        info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.SOPInstanceUID = dicomuid;
    end
    im = M(:,:,slice);
    if (isfield(info, 'SmallestImagePixelValue'))
        info.SmallestImagePixelValue=opts.Convert(min(im(:)));
    end
    if (isfield(info, 'LargestImagePixelValue'))
        info.LargestImagePixelValue=opts.Convert(max(im(:)));
    end
    if (isfield(info, 'RescaleIntercept'))
        info.RescaleIntercept=0;
    end
    if (isfield(info, 'RescaleSlope'))
        info.RescaleSlope=1;
    end
    info.WindowCenter=opts.Convert(opts.Level);
    info.WindowWidth=opts.Convert(opts.Window);
%     if opts.flag3D
%         info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
%         info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleIntercept=0;
%         info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleSlope=1;
%         info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleSlope=1;
%         info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
%         info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);
%     end
	info.SamplesPerPixel=1;
    info.BitsAllocated=16;
    info.BitsStored=16;
    info.HighBit=15;
    info.PixelRepresentation=flag_signed;
    if size(M,3)==slice
        insert_flag=1;
    end
    if insert_flag
        outfile=fullfile(outdir,sprintf('IM%05d.dcm', slice));
        print_progress(outfile);
        if opts.flag3D
           f=fieldnames(info.PerFrameFunctionalGroupsSequence);
           f=setdiff(f,files.info3D.slice2index,'stable');
           for i=1:length(f)
               info.PerFrameFunctionalGroupsSequence=rmfield(info.PerFrameFunctionalGroupsSequence, f{i});
           end
           for i=1:length(files.info3D.slice2index)
               item=files.info3D.slice2index{i};
               info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.InstanceNumber=1;
               info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
               info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleIntercept=0;
               info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleSlope=1;
               info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleSlope=1;
               info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
               info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);   
           end
           info.NumberOfFrames=length(files.info3D.slice2index);
           sz=size(M);
           M=reshape(opts.Convert(M), sz(1), sz(2), 1, []);
           M=permute(M, [2 1 3 4]);
           if isfield(files, 'slices_added')
               if files.slices_added
                   warning('Removing empty slice at bottom of volume');
                   M=M(:,:,1:end-1);
               end
           end
           %dicomwrite(M,outfile, ...
           %    'CreateMode', 'copy', 'WritePrivate', true, info);
           dicomwrite(M,outfile, ...
               'WritePrivate', true, info);
        else
            if isfield(files, 'slices_added')
                if files.slices_added
                    warning('Removing empty slice at bottom of volume');
                    M=M(:,:,1:end-1);
                end
            end
            %dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
            %    'CreateMode', 'copy', 'WritePrivate', true, info);
            %
            dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
                'WritePrivate', true, info);
        end
    end
end
fprintf('\n');


    function print_progress(arg)
        num=length(progress);
        num=num-numel(regexp(progress, '\\\\'));
        for ii=1:num; fprintf('\b'); end
        progress=['Writing file ' arg];
        progress=regexprep(progress,'\','\\\\');
        fprintf(progress);
    end

    function y=convert(x)
        if flag_signed
            y=int16(x*1000);
        else
            y=uint16(x*1000);
        end
    end
end

2.4. QSM_3Dを集める

各フォルダには、強度画像と位相画像が入っている。
TEを7回変えて撮像している。
国際脳QSMの撮像は、1回撮像に128 slicesなので強度画像と位相画像はそれぞれ、896枚(=128 slices ×7 phase)

(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)

QSM_3Dにある強度画像と位相画像をanalyzeフォルダのrawdata(※必ずrawdata)にまとめて保存。

3. 実行

3.1. 作業フォルダ(analyzeフォルダ)へ移動

MATLABを起動し、analyzeフォルダへ移動。

3.2. runme.mを実行

runme.mをMATLABへDrag & Dropし、QSM解析を実行。

3.3. Outputファイルの確認

runme.mを実行後、

  • QSM
  • results

フォルダが生成される。

QSMフォルダにQSM画像がDICOM形式で保存される。

(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)

4. 複数の被験者データをまとめて処理したい場合

被験者ごとにフォルダを作成し一つのフォルダにまとめます。

さらに、後で紹介するまとめてMEDIを実行するスクリプトも入れておきます。

各被験者フォルダには強度画像と位相画像のDICOMがまとめて入ったrawdataフォルダがあります。

以下のスクリプトを実行すると、すべての被験者のQSM mapが計算できます。

4.1. run_sequential.m

clear all;clc;close all;

% change direct to study folder
selpath = uigetdir('Select the folder including all subject');
cd(selpath);
% list folder
folderlist = dir(selpath);
folderlist = folderlist(~ismember({folderlist.name}, {'.', '..'}));
folderlist = folderlist([folderlist.isdir]);

% run MEDI
for i = 1:length(folderlist)
    disp(['processing '  folderlist(i).name ' ...'])
    main(folderlist(i).name, selpath);
end

% define function of MEDI processings
function main(folder, basepath)
    % initialize
    clearvars -except selpath folderlist folder basepath i
    % move to subject folder
    cd(folder)

    % STEP 1: Import data
    [iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('rawdata');

    %%%% Generate the Magnitude image %%%%
    iMag = sqrt(sum(abs(iField).^2,4));

    %%%%% provide a Mask here if possible %%%%%%
    if (~exist('Mask','var'))                     
    %     Mask = genMask(iField, voxel_size);
        Mask = BET(iMag,matrix_size,voxel_size);
    end

    %%%%% provide a noise_level here if possible %%%%%%
    if (~exist('noise_level','var'))
        noise_level = calfieldnoise(iField, Mask);
    end

    %%%%% normalize signal intensity by noise to get SNR %%%
    iField = iField/noise_level;

    % STEP 2a: Field Map Estimation
    %%%%%Estimate the frequency offset in each of the voxel using a 
    %%%%%complex fitting %%%%
    [iFreq_raw N_std] = Fit_ppm_complex(iField);

    % STEP 2b: Spatial phase unwrapping %%%%
    iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);

    % STEP 2c: Background Field Removal
    %%%% Background field removal %%%%
    [RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);

    % CSF Mask for zero referencing
    R2s=arlo(TE,abs(iField));
    Mask_CSF = extract_CSF(R2s, Mask, voxel_size);

    % STEP 3: Dipole Inversion
    save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
        voxel_size delta_TE CF B0_dir Mask_CSF;

    %%%% run MEDI %%%%%
    QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');

    %%% save to DICOM
    % ignore warnings...
    Write_DICOM(QSM,files,'QSM')

    % back to study folder
    cd(basepath)
end