hs_eit_1.m

MATLAB二维电阻抗断层成像算法!用于医学成像,里面包括有限元剖分正问题,及反问题的算法.并且附有网络剖分数据表!

%算例1 仿真试验数据
%1****************************************************************
load meshdata % 两种剖分格式数据.
%load meshdata24
NNode1=max(size(Node1));                      %剖分1的节点数
NElement1=max(size(Element1));                %剖分1的单元数
NNode2=max(size(Node2));                      %剖分2的节点数
NElement2=max(size(Element2));                %剖分2的单元数

g1=reshape([Node1.Coordinate],2,NNode1)';     %剖分1的节点坐标矩阵  
H1=reshape([Element1.Topology],3,NElement1)'; %剖分1的单元拓朴矩阵
g2=reshape([Node2.Coordinate],2,NNode2)';     %剖分2的节点坐标矩阵
H2=reshape([Element2.Topology],3,NElement2)'; %剖分2的单元拓朴矩阵
%****************************************************************

%2****************************************************************
disp('请选择6个单元，用以非均匀仿真计算。。')
Ind=ChooseElements(Node2,Element2,[],10*4); %设定非均匀介质单元6个，剖分2

sigma=1/400*ones(NElement2,1);           %设定均匀介质特性（电导率）
sigma(Ind)=2/400;			             %设定非均匀介质单元（电导率）

clf,Plotinvsol(1./sigma,g2,H2);colorbar,title(['Your resistivity distribution']);drawnow
%图示非均匀介质特性设定情况
disp('按任意键，继续...'),pause
%*********************************************
%计算仿真试验数据
time=cputime;
disp('计算仿真试验数据')
L=16;					         %电极数16
z=0.005*ones(L,1);			     %电极的接触电阻值
[II1,T]=Current(L,NNode2,'tri'); %电流注入模式－矩阵.三角形注入
%[II1,T]=Current(L,NNode2,'adj');
[Agrad,Kb,M,S,C]=FemMatrix(Node2,Element2,z);%有限元系统矩阵构建，剖分2
A=UpdateFemMatrix(Agrad,Kb,M,S,sigma); %封装 %基于电导率的有限元阵更新

[U,p,r]=ForwardSolution(NNode2,NElement2,A,C,T,[],'real'); %FEM计算
Uel=U.Electrode(:);              %有限元结果含（3部分，使用剖分2）      

Agrad1=Agrad*Ind2;               %剖分2转化为剖分1 
%Group some of the element for the inverse computations Integrate over elements of cause mesh

%3****************************************************************
%%                         解逆问题                             %%
% Approximate the best homogenous resistivity.

disp('使用正则化的高斯－牛顿迭代法求解非线性逆问题')
disp('迭代的初始化。。。.')
%******************************************
A=UpdateFemMatrix(Agrad,Kb,M,S,ones(NElement2,1));  
%基于均匀单元的有限元系统矩阵的更新
Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
%基于剖分1的有限元求解
rho0=Uref.Electrode(:)\U.Electrode(:);
%电阻的基值（归一化值）
A=UpdateFemMatrix(Agrad,Kb,M,S,1./rho0*ones(size(sigma)));  
%基于归一化电导率的有限元基值更新
Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
Urefel=Uref.Electrode(:);

rho=rho0*ones(size(Agrad1,2),1);%电阻的初始化值
J=Jacobian(Node2,Element2,Agrad1,Uref.Current,Uref.MeasField, ...
           rho,'real');         %雅可比矩阵的初始化计算，使用剖分1、2
%4****************************************************************

alpha=0.000005;             %正则化因子设定
R=MakeRegmatrix(Element1);  %计算重构矩阵，剖分1？

iter=5;                     %迭代5次即可

disp('迭代中.......')
%5****************************************************************
time=cputime-time
 
for ii=1:iter
 time=cputime;
 rho=rho+(J'*J+alpha*R'*R)\(J'*(Uel-Urefel)-alpha*R'*R*rho);
                                             %更新电阻矩阵
 rhobig=Ind2*rho;                            %转化为剖分1的电阻矩阵
 A=UpdateFemMatrix(Agrad,Kb,M,S,1./rhobig);  %由38行转化的剖分的更新.
 Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
 %基于剖分2的有限元计算
 Urefel=Uref.Electrode(:);
 J=Jacobian(Node2,Element2,Agrad1,Uref.Current,Uref.MeasField,rho,'real');
 %基于剖分2的雅可比矩阵更新
 clf;Plotinvsol(rho,g1,H1);colorbar,title([num2str(ii) '. step']);drawnow;
 %基于剖分1的电阻矩阵的图示
 time=cputime-time
end
%****************************************************************