【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

著者情報：斎藤勇哉

順天堂大学医学部大学院医学研究科放射線診断学講座所属
脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS

月	火	水	木	金	土	日
						1
2	3	4	5	6	7	8
9	10	11	12	13	14	15
16	17	18	19	20	21	22
23	24	25	26	27	28	29
30	31