1. 静電容量断層撮影（ECT）の基本

静電容量断層撮影（ECT）は，配管や容器内部の誘電率分布を非侵襲的に可視化する技術です。一般に8～16個の電極が円周上に設置され，電極間の静電容量を計測し，逆問題を解くことで断面画像を再構成します。

1.1 ECTシステムの構成要素

• 電極アレイ：配管の周囲に等間隔に配置された電極
• データ収集系：電極間の容量を高速に測定
• 画像再構成アルゴリズム：容量測定値から誘電率分布を推定

1.2 ART（代数再構成法）の原理

ARTは線形方程式系 Ax = b を反復的に解く手法です。
A：感度行列（システム行列），x：再構成画像ベクトル（誘電率分布），b：測定容量ベクトル。

2. ECT-ART MATLABシミュレーションコード（改変版）

以下に，8電極ECTシステムに対するART再構成のMATLAB実装を示す。コードはオリジナルから変数名・構造を変更し，処理を簡略化している。

% ECT-ART デモンストレーション（改訂版）
% 目的：8電極系におけるART再構成を実演
clear; close all; clc;
fprintf('=== ECT ART 再構成シミュレーション開始 ===\n');

%% 1. パラメータ設定
N = 32;                % 画像解像度（N×N）
R = 1.0;               % 配管半径 [m]
nElec = 8;             % 電極数
maxIter = 50;          % 最大反復数
relax = 0.25;          % 緩和係数（学習率）
tol = 1e-4;            % 収束許容誤差
epsAir = 1.0;          % 空気の誘電率
epsMat = 3.0;          % 材料の誘電率

%% 2. 感度行列の生成（簡略モデル）
fprintf('感度行列を計算中...\n');
thElec = linspace(0, 360, nElec+1);
thElec = thElec(1:end-1);  % 電極角度 [deg]

% 独立測定ペア数：nElec*(nElec-1)/2
pairs = [];
for i = 1:nElec
    for j = i+1:nElec
        pairs = [pairs; i j];
    end
end
nMeas = size(pairs,1);
fprintf('電極数: %d, 測定数: %d\n', nElec, nMeas);

% グリッド生成
[xG, yG] = meshgrid(linspace(-R, R, N));
rG = sqrt(xG.^2 + yG.^2);
mask = rG <= R;
nPix = sum(mask(:));

% 感度行列 S [nMeas x nPix]
S = zeros(nMeas, nPix);
for m = 1:nMeas
    i = pairs(m,1); j = pairs(m,2);
    thI = deg2rad(thElec(i)); thJ = deg2rad(thElec(j));
    pI = [R*cos(thI), R*sin(thI)];
    pJ = [R*cos(thJ), R*sin(thJ)];
    idx = 0;
    for ix = 1:N
        for iy = 1:N
            if mask(ix,iy)
                idx = idx+1;
                pp = [xG(ix,iy), yG(ix,iy)];
                dI = norm(pp - pI);
                dJ = norm(pp - pJ);
                sens = 1/(dI+dJ+eps);
                angF = abs(cos(thI-atan2(pp(2),pp(1)))) * ...
                       abs(cos(thJ-atan2(pp(2),pp(1))));
                S(m,idx) = sens * angF;
            end
        end
    end
end
S = S / max(S(:));  % 正規化

%% 3. 真の誘電率分布（テストモデル）の生成
fprintf('テストモデルを生成...\n');
trueImg = zeros(N,N);
% ここでは二つの気泡（双発泡）を例とする
c1 = [-R/3, R/4]; c2 = [R/3, -R/4]; bR = R/4;
for i = 1:N
    for j = 1:N
        if mask(i,j)
            p = [xG(i,j), yG(i,j)];
            if norm(p-c1) < bR || norm(p-c2) < bR
                trueImg(i,j) = epsMat;
            else
                trueImg(i,j) = epsAir;
            end
        end
    end
end

% ベクトル化（配管内のみ）
bTrue = trueImg(mask);
bTrue = bTrue(:);

%% 4. 容量測定値のシミュレーション（順問題）
fprintf('容量測定値を計算...\n');
C = S * bTrue;
noiseLvl = 0.02;
C = C .* (1 + noiseLvl*randn(size(C)));

%% 5. ART反復再構成
fprintf('ART再構成開始（最大 %d 反復, 緩和係数 %.2f）\n', maxIter, relax);

% LBP（線形逆投影）で初期推定
xLBP = S' * C;
xLBP = xLBP / max(xLBP) * (epsMat-epsAir) + epsAir;
x = xLBP;  % ART初期値

errHist = zeros(maxIter,1);
xHist = zeros(length(x), maxIter);

for iter = 1:maxIter
    xOld = x;
    % ランダムな測定順序
    order = randperm(nMeas);
    for k = 1:nMeas
        m = order(k);
        s = S(m,:)';
        pred = s' * x;
        res = C(m) - pred;
        denom = s' * s;
        if denom > eps
            x = x + relax * (res/denom) * s;
        end
    end
    % 物理的制約
    x = max(epsAir, min(epsMat, x));
    
    errHist(iter) = norm(x - xOld)/norm(xOld);
    xHist(:,iter) = x;
    if mod(iter,10)==0
        fprintf('  反復 %d/%d, 相対誤差: %.6f\n', iter, maxIter, errHist(iter));
    end
    if errHist(iter) < tol
        fprintf('  反復 %d で収束\n', iter);
        break;
    end
end
nIter = find(errHist>0,1,'last');
if isempty(nIter); nIter = maxIter; end
errHist = errHist(1:nIter);
xHist = xHist(:,1:nIter);

%% 6. 可視化
fprintf('可視化を実行...\n');
% ベクトルから画像へ変換
reconImg = zeros(N,N);  % ART結果
lbpImg = zeros(N,N);   % LBP結果
idx = 0;
for i=1:N
    for j=1:N
        if mask(i,j)
            idx=idx+1;
            reconImg(i,j)=x(idx);
            lbpImg(i,j)=xLBP(idx);
        else
            reconImg(i,j)=NaN;
            lbpImg(i,j)=NaN;
        end
    end
end

figure('Position',[100,100,1000,600]);
subplot(2,3,1);
imagesc(xG(1,:), yG(:,1), trueImg); axis equal; title('真の分布');
colorbar; caxis([epsAir epsMat]); colormap(jet);
hold on;
for e=1:nElec
    th=deg2rad(thElec(e)); 
    plot(R*cos(th),R*sin(th),'wo','MarkerSize',8,'MarkerFaceColor','r');
end
hold off;

subplot(2,3,2);
imagesc(xG(1,:), yG(:,1), lbpImg); axis equal; title('LBP再構成');
colorbar; caxis([epsAir epsMat]);

subplot(2,3,3);
imagesc(xG(1,:), yG(:,1), reconImg); axis equal; 
title(sprintf('ART再構成 (%d回反復)', nIter));
colorbar; caxis([epsAir epsMat]);

subplot(2,3,4);
semilogy(1:nIter, errHist, 'b-','LineWidth',2); grid on;
xlabel('反復回数'); ylabel('相対誤差'); title('収束曲線');
hold on; yline(tol,'r--'); legend('誤差','許容値'); hold off;

subplot(2,3,5);
errMap = abs(reconImg - trueImg);
errMap(~mask)=NaN;
imagesc(xG(1,:), yG(:,1), errMap); axis equal; title('再構成誤差');
colorbar; colormap(hot);

% 定量的指標
mseLBP = mean((lbpImg(mask)-trueImg(mask)).^2);
mseART = mean((reconImg(mask)-trueImg(mask)).^2);
psnrLBP = 10*log10((epsMat-epsAir)^2/mseLBP);
psnrART = 10*log10((epsMat-epsAir)^2/mseART);
ssimLBP = ssim_val(lbpImg(mask), trueImg(mask));
ssimART = ssim_val(reconImg(mask), trueImg(mask));
fprintf('\n--- 評価指標 ---\n');
fprintf('LBP: MSE=%.4f, PSNR=%.2f dB, SSIM=%.4f\n', mseLBP, psnrLBP, ssimLBP);
fprintf('ART: MSE=%.4f, PSNR=%.2f dB, SSIM=%.4f\n', mseART, psnrART, ssimART);

function s = ssim_val(a,b)
    C1 = (0.01*(max(a)-min(a)))^2;
    C2 = (0.03*(max(a)-min(a)))^2;
    mu1=mean(a); mu2=mean(b);
    s1=var(a); s2=var(b);
    c=cov(a,b); s12=c(1,2);
    s = (2*mu1*mu2+C1)*(2*s12+C2)/((mu1^2+mu2^2+C1)*(s1+s2+C2));
end
fprintf('\n=== シミュレーション完了 ===\n');

3. アルゴリズムの詳細

3.1 感度行列の計算

各画素の各電極対に対する感度は，画素から両電極への距離の逆数に角度依存係数を乗じて近似する。

sens = 1/(dI+dJ+eps);
angF = abs(cos(thI-θ_pix)) * abs(cos(thJ-θ_pix));
S(m, pix) = sens * angF;

3.2 ARTの更新式（測定毎）

pred = s' * x;         % 順投影
res = C(m) - pred;     % 残差
x = x + λ * (res/||s||^2) * s;  % 更新

3.3 正則化付きART（改良版）

function x = art_reg(S, C, lambda, alpha, maxIter)
    [M,N] = size(S);
    x = zeros(N,1);
    L = 2*eye(N) - diag(ones(N-1,1),1) - diag(ones(N-1,1),-1); % 二次差分
    for iter=1:maxIter
        xOld = x;
        order = randperm(M);
        for m=order
            s = S(m,:)';
            pred = s'*x;
            res = C(m)-pred;
            denom = s'*s + alpha * (L(m,:)*L(m,:)');
            if denom>eps
                x = x + lambda * (res/denom) * s;
            end
        end
        x = max(0,x);
        if norm(x-xOld)/norm(xOld) < 1e-4; break; end
    end
end

4. 複数流動パターンでの性能比較

以下のコードでは，5種類の代表的な流動パターン（コア流，層流，二気泡，偏心流，環状流）に対しART再構成を行い，定量的に比較する。

patterns = {'コア流','層流','二気泡','偏心流','環状流'};
results = cell(5,4); % MSE, PSNR, SSIM, 反復数

for pIdx = 1:5
    [trueImg, patName] = genPattern(pIdx, N, R, epsAir, epsMat);
    bTrue = trueImg(mask);
    C = S * bTrue(:) + 0.02*randn(nMeas,1);
    x = artRecon(S, C, relax, maxIter, tol);
    reconImg = vec2img(x, mask, N);
    mseVal = mean((reconImg(mask)-trueImg(mask)).^2);
    psnrVal = 10*log10((epsMat-epsAir)^2/mseVal);
    ssimVal = ssim_val(reconImg(mask), trueImg(mask));
    results{pIdx,1}=mseVal; results{pIdx,2}=psnrVal; results{pIdx,3}=ssimVal;
    results{pIdx,4}=nIter;
    % 図示（省略）
end
% 表形式で出力
fprintf('%-10s %-10s %-10s %-10s\n','パターン','MSE','PSNR(dB)','反復数');
for i=1:5
    fprintf('%-10s %-10.4f %-10.2f %-10d\n', patterns{i}, results{i,1}, results{i,2}, results{i,4});
end

5. 実装上の注意と高速化

• メモリ節約：感度行列が疎である場合，sparse型で保持するとメモリ効率が良い。
• 並列化：各測定の更新は独立ではないが，parforで測定ループを並列化可能（注意：収束特性が変わる場合がある）。
• 適応的緩和係数：反復初期は大きく，後半は小さくすることで収束が安定する。

% 適応的緩和係数の例
lambda0 = 0.5; lambdaMin = 0.01; decay = 0.95;
for iter = 1:maxIter
    lam = max(lambda0 * decay^(iter-1), lambdaMin);
    % ... 更新処理
end

6. 実システムへの応用に関する注意

• システム校正：空管時および満管時の容量測定値で感度行列を校正する必要がある。
• ノイズ対策：ハードウェアLPFとソフトウェア（メディアンフィルタなど）の併用が有効。
• リアルタイム性が要求される場合は，LBPを初期値としてART反復回数を10回程度に制限する。

7. まとめ

本稿では，静電容量断層撮影（ECT）における代数再構成法（ART）のMATLAB実装例を紹介した。コードは教育的な目的のため簡略化されているが，実際の応用には感度行列の精度向上，正則化，およびシステム校正が重要となる。