カテゴリ別アーカイブ： MATLAB

Home » カテゴリー別アーカイブ：MATLAB

CONNで前処理したfMRI画像の時系列データを取り出す方法

ある方から、「CONNで前処理したfMRI画像から、デフォルト・モード・ネットワーク(DMN)の時系列データを取り出すにはどうしたらよいですか？」というご質問をいただきました。

鍵となるプログラムは以下の2つです。

conn_matc2nii
spm_summarise

以下に、方法を記載します。

続きを読む →

2023-05-19kiyotakaリンクURL 18 コメント

【FSL】FLICAを用いたマルチモーダル独立成分分析

1. 目的
2. リポジトリ
3. 必要なデータ
4. 注意点
5. 実行
5.1. run_FLICA.m
6. 結果
6.1. index.html
6.2. コンポーネントごとの結果（All_in_one_page）
6.3. モダリティごとの結果（FA, MD, AD, RD, etc.）
7. GBSSへの適応
7.1. index.html
7.2. コンポーネントごとの結果（All_in_one_page）
7.3. モダリティごとの結果（FA, MD, AD, RD, etc.）

続きを読む →

2021-12-15斎藤勇哉リンクURL コメントを残す

【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

2021-12-15斎藤勇哉リンクURL コメントを残す

macOSでのSPM12のコンパイル方法

ある方から、Apple M1のmacでSPMを起動しようとするとspm_check_installation(‘basic’)でエラーが出て起動しないという相談を受けました。

コンパイルしたら問題は解決しました。コンパイル方法を共有します。

ただし、その後、SPMのMLでこのディスカッションに乗ってみたところ、コンパイルは不要だよということも教えていただきました。なので、コンパイルに挑戦してみたい人向けと思ってください。（普通は不要です）

続きを読む →

2021-09-12kiyotakaリンクURL コメントを残す

AC-PC自動設定スクリプト（SPM12対応版）

かつてauto_reorient.mというスクリプトを配布していましたが、これは現在のSPM12で動かなくなってしまいました。理由は単純でSPM12からspm_affregという機能がなくなってしまったからです。もし、過去のSPMからspm_affreg.mを持ってきたら問題なく使えるのですが、この機会に別の手法を考えてみました。

シンプルな方法は、MNI標準脳にCo-registrationすることです。これだけでかなりあいます。
しかし、画像の原点があまりにも違うところに設定されているとエラーが出ることがあります。
そこで、以前、山下先生に教わった方法を採用し、まず、originを画像の中心に設定し、そのうえで、SPM12に搭載されているicbm152.niiにco-registartionするスクリプトを書いてみました。

多くの画像で試してみましたが、それなりにうまくいきますし、処理速度も速いです。

よかったら試してみてください。

acpc_coreg.mをダウンロード（右クリックで保存してください）

Matlabのパスが通っているフォルダにこのファイルを置いていただき、

Matlabから

acpc_coreg

とタイプするだけです。

2017-05-22kiyotakaリンクURL 4 コメント

BashからMatlabスクリプトを実行する方法

先日、ある方と「BashからMatlabを呼び出せないだろうか」という話をしていました。もし、これができたら、シェルスクリプトから、Matlabを呼び出せるので、シェルとMatlabを完全に連携できるわけです。

結論としては、以下でできました。

Short answer

Matlabのスクリプト名を sample_code.m とすると、以下でできます。

$ matlab -nodesktop -nosplash -r 'sample_code; exit'

コツは2つです。

スクリプト名ではなく、コマンドとして指示するため、.mは外す
Matlabから抜けるために exit を追加する

続きを読む →

2017-03-11kiyotakaリンクURL 2 コメント

Matlab R2016a日本語版では、SPM12を起動するときにエラーが出る

タイトル通りなのですが、

Matlab R2016a日本語版では、SPM12を起動しようとするとエラーが出ます。

詳しい現象は、こちらに説明されていますので、下記をご覧ください。

SPM12がMATLAB R2016aで動かない

こちらに解決法が書いてあり、「デスクトップの言語を英語にする」と書いてあります。

Linuxではこれは簡単です。

ターミナルを起動し、

LANG=C matlab &

で英語環境でMatlabが起動します。

Macで同じことができるかと思ったらダメでした。
Macでは、ターミナルでLANG変数を変更しても、Matlab起動時に無視されてしまうとのこと。Mathworks本家が言っているのでどうしようもないです。

Mac プラットフォームのロケール設定

決して賢くはないですが、なんとかする方法がわかりましたので書いておきます。

Matlabを起動する前に、

システム環境設定→言語と地域

で、優先する言語を”English”にします。

そして、そのウィンドウを閉じようとすると、再起動しますかと聞かれますが、再起動しなくてOKです。

そこで、Matlabを起動すると、英語モードで起動します。

そうすると、SPMは問題なく起動します。

終わったあとに、もう一度　System Preferences…　で同様にJapaneseを優先言語にもってきます。

これで終了すれば、日本語に戻ります。

もっとスマートな方法があるかと思いましたが、ありませんでした…。

今後、時間があるときに、SPMのMailing listにバグ報告を出しておこうと思いますが、
一番の対策は、Matlab R2015b以降にはしばらくアップデートしないということになります。

ご参考まで。

2016-05-24kiyotakaリンクURL 7 コメント

A Matlab script to generate ROI masks using an Atlas in SPM12

SPM12 introduces some useful functions such as spm_atlas or new atlas “labels_Neuromorphometrics.” We find the description about labels_Neuromorphometrics in SPM12 Release note.

Maximum probability tissue labels derived from the “MICCAI 2012 Grand Challenge and Workshop on Multi-Atlas Labeling” are available in files tpm/labels Neuromorphometrics.{nii,xml}. These
data were released under the Creative Commons Attribution-NonCommercial (CC BY-NC) with no end date. Users should credit the MRI scans as originating from the OASIS project and the labeled
data as “provided by Neuromorphometrics, Inc. under academic subscription”. These references should be included in all workshop and final publications. See spm templates.man for more details about the generation of this file.

I wanted to generate masks of some regions using this labels_Neuromorphometrics.

Below is the tiny script which generates masks from your preferred atlas.
Running script brings up a file selector. You can choose any atlas you want.
Then it brings up another dialog which lists the region within the atlas. You can choose as many regions as you want, and the scripts generates masks whose file name is the name of the regions.


%generate_masks_from_atlas.m
%This script generate mask files from any atlases you prefer.
%K. Nemoto 25 April 2015

xA=spm_atlas('load');
S=spm_atlas('select',xA);

for i = 1:size(S,2)
    fname=strcat(S{i},'.nii');
    VM=spm_atlas('mask',xA,S{i});
    VM.fname=fname;
    spm_write_vol(VM,spm_read_vols(VM));
end

Download generate_masks_from_atlas.m (right click and save as)

2015-04-25kiyotakaリンクURL 5 コメント

心理のためのMatlabチュートリアル：解答例と解説

2019.10.11: すべての回答ができましたので、アップデートしました。

何人かの方々から、「心理のためのMatlabチュートリアル」の練習の解答がないか問い合わせを受けていました。
回答の完全版ができましたので、公開させていただきます。

心理のためのMatlabチュートリアルの解答例をダウンロード

2014-08-30kiyotakaリンクURL コメントを残す

2024年11月
月	火	水	木	金	土	日
				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