(注意)
ここで公開しているソースコードは自由に改変・再配布してかまいませんが、
内容について作者は一切の責任を負いません。
(1) ファイル mom_inv.zip
をダウンロードして適当なフォルダに展開してください。
(2) [4]のサイトの"MATLABコード"からファイルをダウンロードして適当なフォルダに展開し、
この中から3個の数値ファントムのファイル
（mediumLabelList_Reso1mm_Case1.txt～mediumLabelList_Reso1mm_Case3.txt）
を(1)のフォルダにコピーしてください。
(3) MATLABを起動してエディターウィンドウに(1)のフォルダの
"mom_inv.m"と"input_parameter.m" を開いてください。
（必要に応じて他のファイルも開いてください）
以上でプログラムを実行する環境が整います。
計算を実行する前に、リスト7.1のようにinput_parameter.mの入力データを編集してください。
リスト7.1は表5-2のNo.4の状態です。
リスト7-1 入力データの編集（input_parameter.mの一部）
1 %% input parameters
2 function [Parm, Na, Nf] = input_parameter()
3 %%
4 Parm.restart = 0; % restart (0/1)
5 Parm.H = 8; % ボクセル合体数(整数)
6
7 % 計算条件
8 Parm.maxit = 100; % 最大反復回数
9 Parm.tol = 1e-4; % 収束判定条件
10 Parm.alpha = 0.01; % 正則化パラメータα
11 Parm.c0 = 3.0; % コントラスト初期値c0
12
13 % アンテナ
14 ra = 60e-3; % アンテナ面と原点の距離[m]
15 nxa = 4; % X方向アンテナ数
16 nya = 4; % Y方向アンテナ数
17 nza = 3; % Z方向アンテナ数
18
19 % (1)周囲 (2)脂肪 (3)乳腺 (4)カップ (5)がん (6)皮膚
20 %Parm.epsr = [6.2 7.0 35.0 9.4 50.0 35.0 ];
21 %Parm.sigm = [0.15 0.21 1.05 0 1.5 1.05];
22 Parm.epsr = [5.0 6.5 40.0 9.4 54.0 40.0];
23 Parm.sigm = [0.1 0.2 1.0 0 1.3 1.0];
24 Parm.epsr(1) = Parm.epsr(2); % background=fat
25 Parm.sigm(1) = Parm.sigm(2);
26 Parm.epsr(4) = Parm.epsr(2); % cup=fat
27 Parm.sigm(4) = Parm.sigm(2);
28
29 % 周波数
30 freq = 2.0e9; % 1周波数のとき
31 %freq = linspace(2.0e9, 2.2e9, 2); % 複数周波数のとき
32 Nf = length(freq);
計算を実行するには、コマンドウィンドウで "mom_inv" と実行してください。
コマンドウィンドウに計算経過が表示され、同時に図形出力が行われます。
図形出力はfigure1～figure6から成ります。
それらはファイル "mom_inv.m" のコメントでON/OFFにすることができます。
それぞれの内容は以下の通りです。
本シミュレーターはリスタート機能を持っています。
リスト7-1の4行目を"1"として再度実行します。
これは以下の用途に使用します。
（リスタートが終了したら4行目は"0"に戻してください）
リスト7-2 主プログラムのソースコード(mom_inv.m)
1 %{
2 mom_inv.m (chapter 5_2_1)
3 %}
4
5 clearvars;
6 %close all;
7 tic;
8
9 % input parameters
10 [Parm, Na, Nf] = input_parameter();
11 %figure1(Parm, Na); shg;
12
13 % read numerical phantom
14 medLabFile = 'mediumLabelList_Reso1mm_Case1.txt';
15 [Parm, Epsr, Sigm] = read_phantom(medLabFile, Parm);
16 if Parm.err
17 return;
18 end
19
20 % setup
21 [Parm, N, T10, T20, Robj, Iobj] = setup(Parm, Epsr, Sigm);
22 strtitle = sprintf('H=%d, N=%d, 3N=%d, Na=%d, Na^2=%d, Nf=%d, α=%g, c0=%g', Parm.H, N, 3 * N, Na, Na^2, Nf, Parm.alpha, Parm.c0);
23 fprintf("%s\n", strtitle);
24 I = eye(3 * N);
25
26 % green's function
27 % G : 3N x 3N x Nf
28 % Ga : Na x 3N x Nf
29 [G, Ga] = green(Parm, N, Na, Nf, Robj);
30
31 % E-inc
32 Ei = einc(Parm, N, Na, Nf, Robj); % 3N x Na x Nf
33
34 % forward calculation
35 if (Parm.restart == 0)
36 % solve E
37 E_fwd = zeros(Na^2, Nf);
38 for ifreq = 1 : Nf
39 tau = T10 + (T20 * Parm.joe(ifreq)); % 1 x N
40 S = diag([tau; tau; tau]); % 3N x 3N
41 [Et, E_fwd(:, ifreq)] = solve(ifreq, N, Na, S, G, Ga, Ei);
42 %figure2(Parm, N, Iobj, Et); shg;
43 end
44
45 % initial value for iteration
46 T1 = zeros(N, 1);
47 T2 = zeros(N, 1);
48 for n = 1 : N
49 T1(n) = Parm.c0 * 0.03;
50 T2(n) = Parm.c0;
51 end
52 elseif (Parm.restart == 1)
53 % restart : load
54 load 'mom_inv.mat' E_fwd T1 T2
55 end
56
57 % backward calculation (iteration)
58 fprintf('iter dE/E dεr/εr dσ/σ (εr-εr0)/εr0 (σ-σ0)/σ0 g\n');
59 fiter = [];
60 for iter = 1 : Parm.maxit
61 sumf = zeros(1, 6);
62 for ifreq = 1 : Nf
63 % solve E
64 % Et : 3N x Na
65 % E_bwd : Na^2 x 1
66 % S : 3N x 3N
67 tau = T1 + (T2 * Parm.joe(ifreq));
68 S = diag([tau; tau; tau]);
69 [Et, E_bwd] = solve(ifreq, N, Na, S, G, Ga, Ei);
70
71 % dE : Na^2 x 1
72 dE = E_bwd - E_fwd(:, ifreq);
73
74 % matrix D
75 D = zeros(Na^2, 3 * N); % Na^2 x 3N
76 invISG = inv(I - S * G(:, :, ifreq)); % 3N x 3N
77 for ir = 1 : Na
78 D((ir - 1) * Na + 1 : ir * Na, :) = Ga(:, :, ifreq) * invISG * diag(Et(:, ir)); % Na x 3N
79 end
80
81 % g : regularization factor
82 tr = trace(D' * D);
83 eerr = norm(dE) / norm(E_bwd);
84 g = Parm.alpha * tr / (3 * N) * eerr^2;
85
86 % update T
87 dm = (D' * D + g * I) \ D'; % 3N x Na^2
88 dS = dm * dE; % 3N x 1
89 dS = reshape(dS, N, 3); % N x 3
90 dT = (dS(:, 1) + dS(:, 2) + dS(:, 3)) / 3;
91 T1 = T1 - real(dT);
92 T2 = T2 - imag(dT) / Parm.oe(ifreq);
93
94 % correct T
95 [T1, T2] = correct(T1, T2, Parm);
96
97 % monitor
98 sumf(1) = sumf(1) + eerr;
99 sumf(2) = sumf(2) + norm(imag(dT)) / norm(T2) / Parm.oe(ifreq);
100 sumf(3) = sumf(3) + norm(real(dT)) / norm(T1);
101 sumf(4) = sumf(4) + norm(T2 - T20) / norm(T20);
102 sumf(5) = sumf(5) + norm(T1 - T10) / norm(T10);
103 sumf(6) = sumf(6) + g;
104 end
105 sumf = sumf / Nf;
106
107 % monitor
108 if (iter == 1)
109 fid = fopen('mom_inv.log', 'wt');
110 end
111 fprintf( '%3d %e %e %e %e %e %e\n', iter, sumf);
112 fprintf(fid, '%3d %e %e %e %e %e %e\n', iter, sumf);
113
114 % plot iteration
115 fiter = figure3(strtitle, fiter, sumf); shg;
116 figure4(Parm, N, T10, T20, T1, T2); shg;
117
118 % check convergence
119 if (sumf(2) < Parm.tol) && (sumf(3) < Parm.tol)
120 break;
121 end
122 end
123
124 % plot epsr/sigma contour and histogram
125 figure5(Parm, N, T1, T2, Iobj, Epsr, Sigm);
126 figure6(Parm, T1, T2);
127
128 % output
129 save 'mom_inv.mat' E_fwd T1 T2
130
131 toc
132 %end
133
134 %% solve E
135 function [Et, Es] = solve(ifreq, N, Na, S, G, Ga, Ei)
136 %%
137 I = eye(3 * N);
138 igs = I - G(:, :, ifreq) * S; % 3N x 3N
139 Et = igs \ Ei(:, :, ifreq); % 3N x Na
140 es = Ga(:, :, ifreq) * S * Et; % Na x Na
141 Es = reshape(es, Na^2, 1); % Na^2 x 1
142 end
143
144 %% correct T
145 function [T1, T2] = correct(T1, T2, Parm)
146 %%
147 sgm0 = T1 + Parm.sigm(1);
148 eps0 = T2 + Parm.epsr(1);
149 eps1 = max(1, min(60, eps0));
150 sgm1 = max(0.02 * eps1, min(0.04 * eps1, sgm0));
151 T1 = sgm1 - Parm.sigm(1);
152 T2 = eps1 - Parm.epsr(1);
153 end