slicer_plot_n

function [fc] = slicer_plot_n(h,sol,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)).
sol  = 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. Initially use [fc] = slicer_plot(h,BB,vtx,simp); to plot the slide 
       and calculate fc.

0001 function [fc] = slicer_plot_n(h,sol,vtx,simp,fc);
0002 %function [fc] = slicer_plot_n(h,sol,vtx,simp,fc);
0003 %
0004 %This function plots a 2D slice of the 3D solution vector BB at z=h.
0005 %
0006 %h    = The height of the desired solution, max(vtx(:,3))>= h >= min(vtx(:,3)).
0007 %sol  = The caclulated inverse solution
0008 %vtx  = The vertices matrix
0009 %simp = The simplices matrix
0010 %fc   = The edges of the mesh. This is a 2 column matrix required for the 3D plotting.
0011 %       fc may take some time to be calculated so it is a good idea to save it and use
0012 %       it thereafter. Initially use [fc] = slicer_plot(h,BB,vtx,simp); to plot the slide
0013 %       and calculate fc.
0014 
0015 % $Id: slicer_plot_n.m 3173 2012-06-27 14:48:20Z aadler $
0016 
0017 if nargin < 5
0018 fc = [];
0019 %Develop the faces
0020 
0021 for f=1:size(simp,1)
0022    
0023    fc1 = sort([simp(f,1),simp(f,2)]);
0024    fc2 = sort([simp(f,1),simp(f,3)]);
0025    fc3 = sort([simp(f,1),simp(f,4)]);
0026    fc4 = sort([simp(f,2),simp(f,3)]);
0027    fc5 = sort([simp(f,2),simp(f,4)]);
0028    fc6 = sort([simp(f,3),simp(f,4)]);
0029    
0030    fc = [fc;fc1;fc2;fc3;fc4;fc5;fc6];
0031    
0032 end
0033 fc = unique(fc,'rows');
0034 end
0035 
0036 %(1) Generate the pseudo-triangulation at plane z=h
0037 vtxp = []; %Nodes created for the plane
0038 vap = []; %Value of the node in vtxp
0039 
0040 for j=1:size(fc,1)
0041       this_ph = fc(j,:); %[nodeA nodeB]
0042    
0043    if max(vtx(this_ph(1),3),vtx(this_ph(2),3))> h && ...
0044       min(vtx(this_ph(1),3),vtx(this_ph(2),3))<= h         
0045      
0046   %If the face is crossed through by the plane then
0047   %create a plotable node on the plane.
0048     Pa = this_ph(1); Pb = this_ph(2);
0049     xa = vtx(Pa,1); ya = vtx(Pa,2); za = vtx(Pa,3);
0050     xb = vtx(Pb,1); yb = vtx(Pb,2); zb = vtx(Pb,3);
0051   
0052     xn = xa + (h-za)*(xb-xa)/(zb-za);
0053     yn = ya + (h-za)*(yb-ya)/(zb-za);
0054     vtxp = [vtxp;[xn,yn]];
0055     
0056   end %if
0057 end %for
0058 tri = delaunay(vtxp(:,1),vtxp(:,2));
0059 
0060 [vtxp,tri] = delfix(vtxp,tri);
0061 %The 2D mesh at h is (vtxp,tri)
0062 
0063 %(2) Evaluate the geometric centers gCts of the new siplices tri
0064 gCts = zeros(size(tri,1),2);
0065 for y=1:size(tri,1)
0066     gCts(y,1) = mean(vtxp(tri(y,:),1));
0067     gCts(y,2) = mean(vtxp(tri(y,:),2));
0068 end
0069 
0070 %(3) Initialise the planar solution
0071 sol2D = zeros(size(gCts,1),1);
0072 
0073 %(4) Now trace which simps contain gCts
0074     
0075     
0076  TT = tsearchn(vtx,simp,[gCts,h*ones(size(gCts,1),1)]);       
0077  sol2D = sol(TT);
0078  
0079   
0080 %(5) Plot the planar solution sol2D with patches
0081 % Autoscale each axes to its own scale
0082 % c_img = calc_colours( sol2D(:), [], 1 );
0083 c_img = calc_colours( sol2D(:), [], 1 );
0084 
0085 for q=1:size(tri)
0086    tri_q= tri(q,:);
0087 % need 'direct' otherwise colourmap is screwed up
0088    patch(vtxp(tri_q,1),vtxp(tri_q,2),squeeze(c_img(q,:,:)), ...
0089          'CDataMapping','direct','EdgeColor','none');
0090 end %for q
0091 % colorbar;
0092 
0093 function [vtx_n,simp_n] = delfix(vtx,simp)
0094 %function [vtx_n,simp_n] = delfix(vtx,simp)
0095 %
0096 % Auxiliary function to remove the zero area faces
0097 % produced by Matlab's delaunay triangulation
0098 %
0099 %
0100 %
0101 %vtx  = The vertices matrix
0102 %simp = The simplices matrix
0103 
0104 
0105 simp_n = [];
0106 tri_a = [];
0107 
0108 for kk=1:length(simp)
0109    
0110    this_tri = simp(kk,:);
0111    
0112    xa = vtx(this_tri(1),1); ya = vtx(this_tri(1),2);
0113    xb = vtx(this_tri(2),1); yb = vtx(this_tri(2),2);
0114    xc = vtx(this_tri(3),1); yc = vtx(this_tri(3),2);
0115    
0116    tria = polyarea([xa;xb;xc],[ya;yb;yc]);
0117    tri_a = [tri_a ; tria];
0118    
0119    if tria > 0.00000000001       
0120       simp_n = [simp_n;this_tri];
0121    end
0122    
0123 end
0124 
0125 vtx_n = vtx;
0126 
0127 
0128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0129 % This is part of the EIDORS suite.
0130 % Copyright (c) N. Polydorides 2003
0131 % Copying permitted under terms of GNU GPL
0132 % See enclosed file gpl.html for details.
0133 % EIDORS 3D version 2.0
0134 % MATLAB version 6.1 R13
0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0136