slicer_plot.sci

用来实现三维阻抗及光学断层成像重建的matlab程序

function [fc]=slicer_plot(h,BB,vtx,%simp,fc)
[nargout,nargin] = argn(0)
//function [fc] = slicer_plot(h,BB,vtx,simp,fc);
// 
//This function plots a 2D slice of the 3D solution vector BB at z=h.
// 
// 
// 
//h    = The height of the desired solution, max(vtx(:,3))>= h >= min(vtx(:,3)).
//BB   = The caclulated inverse solution
//vtx  = The vertices matrix
//simp = The simplices matrix
//fc   = The edges of the mesh. This is a 2 column matrix required for the 3D plotting. 
//       fc may take some time to be calculated so it is a good idea to save it and use 
//       it thereafter.
 
 
//Extract the nodal solution CC
KK = solution_ext(BB,vtx,%simp);
 
//Need to rescale that to match BB.
 
 
//!  mtlb_max(BB) may be replaced by 
//!     max(BB) if BBis a vector
//!     max(BB,'r') if BBis a matrix
 
//! mtlb_min(BB) may be replaced by 
//!   min(BB) if BBis a vector
//!   min(BB,'r') if BBis a matrix
ranges = (mtlb_max(BB)-mtlb_min(BB))/(max(KK)-min(KK));
 
%v=size(KK)
CC = zeros(%v(1),%v(2));
 
for yu = 1:max(size(KK))
   
  //! CC = mtlb_i(CC,yu,ranges*KK(yu)) may be replaced by CC(yu) = ranges*KK(yu)
  //!   if CC is not a  vector
  //!   or if CC and ranges*KK(yu) are both row or column vectors
  CC = mtlb_i(CC,yu,ranges*KK(yu))
end
 
 
 
if nargin<5 then
  fc = [];
  //Develop the faces           
   
  for f = 1:size(%simp,1)
     
    %v2 = [%simp(f,1),%simp(f,2)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc1 = %v2;
    %v2 = [%simp(f,1),%simp(f,3)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc2 = %v2;
    %v2 = [%simp(f,1),%simp(f,4)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc3 = %v2;
    %v2 = [%simp(f,2),%simp(f,3)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc4 = %v2;
    %v2 = [%simp(f,2),%simp(f,4)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc5 = %v2;
    %v2 = [%simp(f,3),%simp(f,4)]
     
    //! next test may  probably be simplified
    if min(size(%v2))==1 then 
      %v2=-sort(-%v2)
    else
      %v2=-sort(-%v2,'r')
    end
    fc6 = %v2;
     
    fc = [fc;fc1;fc2;fc3;fc4;fc5;fc6];
     
  end
   
   
  //!! Unknown function unique ,the original calling sequence is used
  fc = unique(fc,'rows');
   
end
 
 
vtxp = [];
//Nodes created for the plane
vap = [];
//Value of the node in vtxp
 
for j = 1:size(fc,1)
   
  this_ph = fc(j,:);
  //[nodeA nodeB]
   
  if (max(vtx(this_ph(1),3),vtx(this_ph(2),3))>h)&(min(vtx(this_ph(1),3),vtx(this_ph(2),3))<=h) then
     
    //Section crosses through face
    //If face is crossed through by h1 plane then is eligible to create a plotable node (on the plane).
     
    Pa = this_ph(1);
    Pb = this_ph(2);
     
    xa = vtx(Pa,1);
    ya = vtx(Pa,2);
    za = vtx(Pa,3);
    xb = vtx(Pb,1);
    yb = vtx(Pb,2);
    zb = vtx(Pb,3);
     
     
    xn = xa+(h-za)*(xb-xa)/(zb-za);
    yn = ya+(h-za)*(yb-ya)/(zb-za);
     
    vtxp = [vtxp;[xn,yn]];
     
    //Specifing the value of the new node
     
    la = db23d(xa,ya,za,xn,yn,h);
     
    lb = db23d(xb,yb,zb,xn,yn,h);
     
    vv = CC(Pa)*lb/(la+lb)+CC(Pb)*la/(la+lb);
     
    vap = [vap;vv];
     
  end
  //if
   
end
//for
 
 
//!! Unknown function delaunay ,the original calling sequence is used
tri = delaunay(vtxp(:,1),vtxp(:,2));
 
[vtxp,tri] = delfix(vtxp,tri);
 
X = vtxp(:,1);
Y = vtxp(:,2);
 
//! unknown arg type, using mtlb_length 
Z = h*ones(mtlb_length(vtxp),1);
 
//!! Unknown function trisurf ,the original calling sequence is used
trisurf(tri,X,Y,Z,vap);
 
//!! Unknown function shading ,the original calling sequence is used
shading('interp');
 
//!! Unknown function axis ,the original calling sequence is used
axis([min(vtx(:,1)),max(vtx(:,1)),min(vtx(:,2)),max(vtx(:,2)),min(vtx(:,3)),max(vtx(:,3))]);
 
//You may want to change this setting !!  
 
//! mtlb_min(BB) may be replaced by 
//!   min(BB) if BBis a vector
//!   min(BB,'r') if BBis a matrix
 
//!  mtlb_max(BB) may be replaced by 
//!     max(BB) if BBis a vector
//!     max(BB,'r') if BBis a matrix
 
//!! Unknown function caxis ,the original calling sequence is used
caxis([mtlb_min(BB),mtlb_max(BB)]);
 
 
//!! Unknown function axis ,the original calling sequence is used
axis('square');
 
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// This is part of the EIDORS suite.
// Copyright (c) N. Polydorides 2001
// Copying permitted under terms of GNU GPL
// See enclosed file gpl.html for details.
// EIDORS 3D version 1.0
// MATLAB version 5.3 R11
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%