Home > eidors > graphics > matlab > show_fem.m

show_fem

PURPOSE ^

SHOW_FEM: show the EIDORS3D finite element model

SYNOPSIS ^

function hh=show_fem( mdl, options)

DESCRIPTION ^

 SHOW_FEM: show the EIDORS3D finite element model
 hh=show_fem( mdl, options )
 mdl is a EIDORS3D 'model' or 'image' structure
 hh= handle to the plotted model

 options may be specified by a list

 options specifies a set of options
   options(1) => show colourbar
   options(2) => show numbering on electrodes
   options(3) => number elements (==1) or nodes (==2);

 for detailed control of colours, use
    img.calc_colours."parameter" = value
 see help for calc_colours.

 for control of colourbar, use img.calc_colours.cb_shrink_move

 parameters
     fwd_model.show_fem.linecolour

 to change line properties:
      hh=show_fem(...); set(hh,'EdgeColor',[0,0,1];

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function hh=show_fem( mdl, options)
0002 % SHOW_FEM: show the EIDORS3D finite element model
0003 % hh=show_fem( mdl, options )
0004 % mdl is a EIDORS3D 'model' or 'image' structure
0005 % hh= handle to the plotted model
0006 %
0007 % options may be specified by a list
0008 %
0009 % options specifies a set of options
0010 %   options(1) => show colourbar
0011 %   options(2) => show numbering on electrodes
0012 %   options(3) => number elements (==1) or nodes (==2);
0013 %
0014 % for detailed control of colours, use
0015 %    img.calc_colours."parameter" = value
0016 % see help for calc_colours.
0017 %
0018 % for control of colourbar, use img.calc_colours.cb_shrink_move
0019 %
0020 % parameters
0021 %     fwd_model.show_fem.linecolour
0022 %
0023 % to change line properties:
0024 %      hh=show_fem(...); set(hh,'EdgeColor',[0,0,1];
0025 
0026 % (C) 2005-2011 Andy Adler. License: GPL version 2 or version 3
0027 % $Id: show_fem.m 5786 2018-05-21 13:37:25Z aadler $
0028 
0029 % TESTS
0030 switch nargin
0031    case 0;
0032      error('Insufficient parameters for show_fem');
0033    case 1;
0034      if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0035      options = [];
0036 end
0037 [img,mdl,opts] = proc_params( mdl, options );
0038 
0039 if ~ishold
0040     cla;
0041 end
0042 
0043 switch opts.dims
0044    case 2;    hh= show_2d(img,mdl,opts); % 2D nodes
0045    case 3;    hh= show_3d(img,mdl,opts); % 3D nodes
0046    otherwise; error('model is not 2D or 3D');
0047 end
0048 
0049 if opts.show_elem_numbering
0050    xyzc= interp_mesh(mdl);
0051    placenumbers(xyzc, 7, [0,0,0],'none');
0052 end
0053 if opts.show_node_numbering
0054    xyzc= mdl.nodes;
0055    placenumbers(xyzc, 7, [0.0,0,0.5],[1,1,1]);
0056 end
0057 
0058 if nargout == 0; clear hh; end
0059 
0060 % this is moved to show_2d and show_3d
0061 % if ~ishold
0062 %    fix_axis
0063 % end
0064 
0065 function fix_axis
0066    ax = gca;
0067    axis tight
0068    set(ax,'DataAspectRatio',[1 1 1]);
0069    % expand the axis limits slightly
0070    [az el] = view;
0071    if el==90
0072       % for 2d plots we need x and y
0073       xl=get(ax,'XLim'); yl=get(ax,'Ylim');
0074       dx = diff(xl); dy = diff(yl);
0075       % make sure to expand away from the origin
0076       xb = strcmp(get(ax,'XAxisLocation'),'bottom');
0077       yn = strcmp(get(ax,'YDir'),'normal');
0078       if xb == yn
0079          set(ax,'YLim',yl + [0 1]*1e-3*dy);
0080       else
0081          set(ax,'YLim',yl + [-1 0]*1e-3*dy);
0082       end
0083       yl = strcmp(get(ax,'XAxisLocation'),'bottom');
0084       xn = strcmp(get(ax,'XDir'),'normal');
0085       if yl == xn
0086          set(ax,'XLim',xl + [0 1]*1e-3*dx);
0087       else
0088          set(ax,'XLim',xl + [-1 0]*1e-3*dx);
0089       end
0090    else
0091       % for 3d plots it's enough to adjust z
0092       zl = get(ax,'Zlim'); dz = diff(zl);
0093       if strcmp(get(ax,'Zdir'),'normal')
0094          set(ax,'ZLim',zl + [0 1]*1e-3*dz);
0095       else
0096          set(ax,'ZLim',zl + [-1 0]*1e-3*dz);
0097       end
0098    end
0099    
0100 
0101 function placenumbers(xyzc, fontsize, colour, bgcolour)
0102    xyzc= xyzc * eye(size(xyzc,2),3); %convert to 3D
0103    for i= 1:size(xyzc,1)
0104       text(xyzc(i,1),xyzc(i,2), xyzc(i,3), num2str(i), ...
0105             'HorizontalAlignment','center', ...
0106             'FontSize',fontsize,'Color',colour,'BackgroundColor',bgcolour);
0107    end
0108 
0109 function [img,mdl,opts] = proc_params( mdl, options );
0110 
0111    opts.do_colourbar=0;
0112    opts.number_electrodes=0;
0113    opts.show_numbering  =0;
0114    opts.show_elem_numbering = 0;
0115    opts.show_node_numbering = 0;
0116    if nargin >=2
0117        % fill in default options
0118        optionstr= zeros(1,100);
0119        optionstr(1:length(options)) = options;
0120 
0121        opts.do_colourbar=      optionstr(1);
0122        opts.number_electrodes= optionstr(2);
0123        switch optionstr(3)
0124           case 1; opts.show_elem_numbering = 1;
0125           case 2; opts.show_node_numbering = 1;
0126           case 3; opts.show_elem_numbering = 1;
0127                   opts.show_node_numbering = 1;
0128        end
0129    end
0130 
0131    
0132    % if we have an only img input, then define mdl
0133    if strcmp( mdl(1).type , 'image' )
0134       img= mdl;
0135       mdl= img(1).fwd_model;
0136    else 
0137       img = [];
0138    end
0139    
0140    opts.transparency_thresh = calc_colours('transparency_thresh');
0141    try
0142        opts.transparency_thresh = img.calc_colours.transparency_thresh;
0143    end
0144    
0145    opts.dims = size(mdl.nodes,2);
0146 
0147 % 2D Case
0148 function hh= show_2d(img,mdl,opts)
0149    hax= gca;
0150    pax= get(hax,'position');
0151    if ~isempty(img)
0152       colours= calc_colours(img, [] );
0153    else
0154       colours= [1,1,1]; % white elements if no image
0155    end
0156    hh= show_2d_fem( mdl, colours );
0157    show_electrodes_2d(mdl, opts.number_electrodes);
0158 
0159    set(hax,'position', pax);
0160    view(0, 90); axis('xy'); grid('off');
0161 
0162 % IN MATLAB 7 WE NEED TO RERUN THIS BECAUSE STUPID STUPID
0163 % MATLAB WILL RESET THE COLOURBAR EVERY TIME WE RUN PATCH!!!
0164    if exist('img','var') && opts.do_colourbar;
0165       colours= calc_colours(img, [], opts.do_colourbar);
0166       % Matlab is so weird. It puts the first colorbar in the wrong place
0167       %   sometimes ...  (tested with 6.5 and with 7.8)
0168       %   The trick is to never try to move it on the first go
0169       %   OR we reset it and then replace it. STUPID STUPID
0170 
0171      % Here's the magic trick I found. Force a drawnow,
0172      %     then delete and recreate
0173 
0174      % Because of a change in matlab colorbar somewhere around
0175      % 2012, none of this compensation code works any more. We
0176      % don't really understand what to do anymore, but this
0177      % seems best, for now ...
0178    end
0179    if ~ishold
0180       fix_axis
0181    end
0182 
0183    
0184 
0185 % 3D Case
0186 function hh= show_3d(img,mdl,opts)
0187    hh= show_3d_fem( mdl );
0188 
0189    if ~isempty(img)
0190        elem_data = get_img_data(img);
0191        show_inhomogeneities( elem_data , mdl, img, opts);
0192        if opts.do_colourbar
0193            calc_colours(img, [], opts.do_colourbar);
0194        end
0195    end
0196    % need to adjust the axis limits before we plot electrodes
0197    % who have thickness in all dimensions and thus cause 'axis tight' to
0198    % mess up the z limits
0199    if ~ishold
0200       fix_axis
0201    end
0202    if size(mdl.elems,2) == 3
0203       show_electrodes_surf(mdl, opts.number_electrodes);
0204    else
0205       show_electrodes_3d(mdl, opts.number_electrodes);
0206    end
0207 
0208 
0209 function show_electrodes_2d(mdl, number_electrodes)
0210     if ~isfield(mdl,'electrode'); return; end
0211 
0212     ee= get_boundary( mdl );
0213     ctr_x= mean(mdl.nodes(:,1));
0214     ctr_y= mean(mdl.nodes(:,2));
0215 
0216 % scale away from model
0217 
0218 for e=1:length(mdl.electrode)
0219     if isfield(mdl.electrode(e),'nodes')
0220         elec_nodes= mdl.electrode(e).nodes;
0221         
0222         S= 1.00;
0223         vx= (mdl.nodes(elec_nodes,1) - ctr_x)*S;
0224         vy= (mdl.nodes(elec_nodes,2) - ctr_y)*S;
0225         % sort nodes around the model (to avoid crossed lines)
0226         [jnk,idx] = order_loop( [vx,vy] );
0227     elseif isfield(mdl.electrode(e),'pos')
0228         vx = mdl.electrode(e).pos(:,1) - ctr_x;
0229         vy = mdl.electrode(e).pos(:,2) - ctr_y;
0230         idx = 1:length(vx);
0231     else
0232        eidors_msg('show_fem: WARNING: electrode %d has no nodes. Not showing.',e,2);
0233        continue;
0234     end
0235         
0236     ecolour = electr_colour( e );
0237     if numel(vx) == 1
0238        % Point Electrode Models: put a circle around the node
0239        line(vx(idx)+ctr_x,vy(idx)+ctr_y,  ...
0240             'LineWidth', 2, 'LineStyle','-','Color', ecolour, ...
0241             'Marker','o','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0242     else
0243        % Complete/Shunt Electrode Models (multiple nodes per electrode)
0244        %  put a line along the edges that form the electrode
0245        line(vx(idx)+ctr_x,vy(idx)+ctr_y,  ...
0246             'LineWidth', 3, 'LineStyle','-','Color', ecolour, ...
0247             'Marker','none','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0248     end
0249     if number_electrodes
0250        S= 1.05;
0251        vx= vx*S;
0252        vy= vy*S;
0253        switch number_electrodes
0254           case {1 true}
0255              txt = num2str(e);
0256           case 2
0257              try, txt = mdl.electrode(e).label; end
0258        end
0259        hh= text(mean(vx)+ctr_x, mean(vy)+ctr_y, txt);
0260        set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0261     end
0262 end
0263 
0264 function show_electrodes_surf(mdl, number_electrodes)
0265     if ~isfield(mdl,'electrode'); return; end
0266 
0267     ee= get_boundary( mdl );
0268     ctr_x= mean(mdl.nodes(:,1));
0269     ctr_y= mean(mdl.nodes(:,2));
0270     ctr_z= mean(mdl.nodes(:,3));
0271 % scale away from model
0272 
0273 for e=1:length(mdl.electrode)
0274     if isfield(mdl.electrode(e),'pos') && ~isfield(mdl.electrode(e),'nodes')
0275         vx = mdl.electrode(e).pos(:,1) - ctr_x;
0276         vy = mdl.electrode(e).pos(:,2) - ctr_y;
0277         vz = mdl.electrode(e).pos(:,3) - ctr_z;
0278         idx = 1:length(vx);
0279         S = 1.00;
0280     else
0281         elec_nodes= mdl.electrode(e).nodes;
0282         
0283         S= 1.00;
0284         vx= (mdl.nodes(elec_nodes,1) - ctr_x);
0285         vy= (mdl.nodes(elec_nodes,2) - ctr_y);
0286         vz= (mdl.nodes(elec_nodes,3) - ctr_z);
0287     end
0288     
0289     % sort nodes around the model (to avoid crossed lines)
0290     % TODO: figure out what to do in different directions
0291     [jnk,idx] = sort(unwrap(atan2( vy, vx )));
0292     
0293     ecolour = electr_colour( e );
0294     if numel(vx) == 1
0295        % Point Electrode Models: put a circle around the node
0296        line(S*vx(idx)+ctr_x,S*vy(idx)+ctr_y, S*vz(idx)+ctr_z,  ...
0297             'LineWidth', 2, 'LineStyle','-','Color', ecolour, ...
0298             'Marker','o','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0299     else
0300        % Complete/Shunt Electrode Models (multiple nodes per electrode)
0301        %  put a line along the edges that form the electrode
0302        line(vx(idx)+ctr_x,vy(idx)+ctr_y, vz(idx)+ctr_z, ...
0303             'LineWidth', 3, 'LineStyle','-','Color', ecolour, ...
0304             'Marker','none','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0305     end
0306     if number_electrodes
0307        S= 1.05;
0308 %        vx= (mdl.nodes(elec_nodes,1) - ctr_x)*S;
0309 %        vy= (mdl.nodes(elec_nodes,2) - ctr_y)*S;
0310 %        vz= (mdl.nodes(elec_nodes,3) - ctr_z)*S;
0311        switch number_electrodes
0312           case {1 true}
0313             txt = num2str(e);
0314           case 2
0315              try, txt = mdl.electrode(e).label; end
0316        end
0317        hh= text(S*mean(vx)+ctr_x, S*mean(vy)+ctr_y,S*mean(vz)+ctr_z,txt);
0318        set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0319     end
0320 end
0321 
0322 function show_electrodes_3d(mdl, number_electrodes)
0323 % show electrode positions on model
0324 if ~isfield(mdl,'electrode'); return; end
0325 
0326 ee= get_boundary( mdl );
0327 for e=1:length(mdl.electrode)
0328     colour= electr_colour( e);
0329     
0330     if isfield(mdl.electrode(e),'pos') && ~isfield(mdl.electrode(e),'nodes')
0331         show_electrodes_surf(mdl, number_electrodes);
0332         return
0333     end
0334     elec_nodes= mdl.electrode(e).nodes;
0335     
0336     
0337     if length(elec_nodes) == 1  % point electrode model
0338         vtx= mdl.nodes(elec_nodes,:);
0339         line(vtx(1),vtx(2),vtx(3), ...
0340             'Marker','o','MarkerSize',12, ...
0341             'MarkerFaceColor',colour, 'MarkerEdgeColor', colour);
0342         if number_electrodes
0343             hh= text(vtx(1),vtx(2),vtx(3), num2str(e));
0344             set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0345         end
0346     else
0347         % find elems on boundary attached to this electrode
0348         map = zeros(length(mdl.nodes),1);
0349         map(mdl.electrode(e).nodes) = 1;
0350         ec = map(ee);
0351         sels= find(all(ec'));
0352         
0353         ee= get_boundary( mdl );
0354         paint_electrodes(sels,ee,mdl.nodes,colour,number_electrodes);
0355         
0356         if number_electrodes
0357             el_nodes= mdl.nodes(unique(mdl.boundary(sels,:)),:);
0358             switch number_electrodes
0359                case {1 true}
0360                   txt = num2str(e);
0361                case 2
0362                   try, txt = mdl.electrode(e).label; end
0363             end
0364             hh=text( mean(el_nodes(:,1)), ...
0365                      mean(el_nodes(:,2)), ...
0366                      mean(el_nodes(:,3)), txt );
0367             set(hh,'FontWeight','bold','Color',[.8,.2,0],'FontSize',8);
0368         end
0369     end
0370 end
0371 
0372 function show_inhomogeneities( elem_data, mdl, img, opt)
0373 % show
0374 if size(elem_data,2)>1
0375    q = warning('query','backtrace');
0376    warning('backtrace','off');
0377    warning('EIDORS:FirstImageOnly','show_fem only shows first image');
0378    warning('backtrace',q.state);
0379 %    eidors_msg('warning: show_fem only shows first image',1);
0380 end
0381 repaint_inho(elem_data(:,1), 'use_global' , mdl.nodes, mdl.elems, ...
0382     opt.transparency_thresh, img); 
0383 if ~exist('OCTAVE_VERSION');
0384 camlight('left');
0385 lighting('none'); % lighting doesn't help much
0386 end
0387 
0388 function paint_electrodes_old(sel,srf,vtx, colour, show_num);
0389 %function paint_electrodes(sel,srf,vtx);
0390 %
0391 % plots the electrodes red at the boundaries.
0392 %
0393 % sel = The index of the electrode faces in the srf matrix
0394 %       sel can be created by set_electrodes.m
0395 % srf = the boundary faces (triangles)
0396 % vtx = The vertices matrix.
0397 
0398 l = srf(sel,1); m = srf(sel,2); n = srf(sel,3);
0399 
0400 Xs = [vtx(l,1);vtx(m,1);vtx(n,1)];
0401 Ys = [vtx(l,2);vtx(m,2);vtx(n,2)];
0402 Zs = [vtx(l,3);vtx(m,3);vtx(n,3)];
0403 
0404 h=patch(Xs,Ys,Zs, colour);
0405 % need 'direct' otherwise colourmap is screwed up
0406 set(h, 'FaceLighting','none', 'CDataMapping', 'direct' );
0407 
0408 function paint_electrodes(sel,srf,vtx, colour, show_num);
0409    if isempty(sel); return; end  % Not required after matlab 2014
0410 
0411    [u n m] = unique(srf(sel,:));
0412    fv.vertices = vtx(u,:);
0413    fv.faces = reshape(m,[],3);
0414    h = patch(fv,'FaceColor',colour);
0415    % need 'direct' otherwise colourmap is screwed up
0416    set(h, 'FaceLighting','none', 'CDataMapping', 'direct' );
0417 
0418 function hh= show_3d_fem( mdl, options )
0419    if mdl_dim(mdl) == 3  % volume simplices
0420       ee= get_boundary( mdl );
0421    elseif mdl_dim(mdl) == 2  % volume simplices
0422       ee= mdl.elems;
0423    elseif mdl_dim(mdl) == 1  % just lines between elements. Resistors?
0424       ee= mdl.elems;
0425    else
0426       error('Simplices are not 2D or 3D in mdl. Unsure how to display');
0427    end
0428    hh= trimesh(ee, mdl.nodes(:,1), ...
0429                    mdl.nodes(:,2), ...
0430                    mdl.nodes(:,3));
0431    set(hh, 'EdgeColor', [0,0,0]);
0432    axis('image');
0433    hidden('off');
0434 
0435 function hh=show_2d_fem( mdl, colours, imgno )
0436   szcolr = size(colours);
0437   if nargin<3;
0438       imgno = 1;
0439       if szcolr(1:2)>1; 
0440           eidors_msg('warning: show_fem only shows first image',1);
0441       end
0442   end
0443 
0444   if szcolr(1:2) == [1,3]  % no reconstruction  - just use white
0445      hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor',colours);
0446   elseif size(colours,1) == num_elems(mdl);
0447      colours = squeeze(colours(:,imgno,:));
0448 
0449      hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor','flat', ...
0450                'facevertexcdata',colours,'CDataMapping','direct'); 
0451   elseif size(colours,1) == num_nodes(mdl);
0452      colours = squeeze(colours(:,imgno,:));
0453      % 'interp' not supported in octave
0454      hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor','interp', ...
0455                'facevertexcdata',colours,'CDataMapping','direct'); 
0456   else
0457     eidors_msg('warning: image elements and mesh do not match. Showing grey',1);
0458     colours= [1,1,1]/2; % Grey to show we're not sure
0459     hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor',colours);
0460   end
0461 
0462   set(hh, 'FaceLighting','none', 'CDataMapping', 'direct' );
0463 
0464 function hh=show_2d_fem_oldmatlab( mdl, colours, imgno )
0465   szcolr = size(colours);
0466   if nargin<3;
0467       imgno = 1;
0468       if szcolr(1:2)>1;  eidors_msg('warning: show_fem only shows first image',1); end
0469   end
0470   
0471 % el_pos= avg_electrode_posn( mdl );
0472 % plot_2d_mesh(mdl.nodes', mdl.elems', el_pos', .95, [0,0,0]);
0473 
0474   S= 1; %.95; % shrink factor
0475   elem= mdl.elems'; e= size(elem,2);
0476   node= mdl.nodes'; n= size(node,2);
0477   Xs=zeros(3,e);
0478   Xs(:)=mdl.nodes(elem(:),1);
0479   Xs= S*Xs+ (1-S)*ones(3,1)*mean(Xs);
0480   Ys=zeros(3,e); Ys(:)=mdl.nodes(elem(:),2);
0481   Ys= S*Ys+ (1-S)*ones(3,1)*mean(Ys);
0482   Zs = zeros(size(Xs));
0483   if szcolr(1:2) == [1,3]  % no reconstruction  - just use white
0484      hh= patch(Xs,Ys,zeros(3,e),colours);
0485   elseif size(colours,1) == e % defined on elems
0486 % THE STUPID MATLAB 7 WILL RESET THE COLOURBAR WHENEVER YOU DO A PATCH. DAMN.
0487      colours = permute(colours(:,imgno,:),[2,1,3]);
0488   elseif size(colours,1) == n  % multiple images
0489      colours = colours(elem(:), imgno, :);
0490      colours = reshape( colours, size(elem,1), size(elem,2), []);
0491   else
0492     eidors_msg('warning: image elements and mesh do not match. Showing grey',1);
0493     colours= [1,1,1]/2; % Grey to show we're not sure
0494   end
0495 
0496   hh= patch(Xs,Ys,Zs,colours);
0497   set(hh, 'FaceLighting','none', 'CDataMapping', 'direct' );
0498   % FOR NODE RGB MATLAB SCREWS UP THE COLOURS FOR US (ONCE AGAIN). THERE MUST
0499   % BE SOME KIND OF SECRET FLAG@@@
0500 
0501   max_x= max(mdl.nodes(:,1)); min_x= min(mdl.nodes(:,1));
0502   max_y= max(mdl.nodes(:,2)); min_y= min(mdl.nodes(:,2));
0503   axis([ mean([max_x,min_x]) + 0.55*[-1,1]*(max_x-min_x), ...
0504          mean([max_y,min_y]) + 0.55*[-1,1]*(max_y-min_y) ]);
0505 
0506 
0507 
0508 function old_code_3d_plot
0509   domesnum= 0;
0510   donodenum= 0;
0511   dotext=0;
0512 
0513   xxx=zeros(4,e); xxx(:)=NODE(1,ELEM(:));
0514   xxx= S*xxx+ (1-S)*ones(4,1)*mean(xxx);
0515   yyy=zeros(4,e); yyy(:)=NODE(2,ELEM(:));
0516   yyy= S*yyy+ (1-S)*ones(4,1)*mean(yyy);
0517   zzz=zeros(4,e); zzz(:)=NODE(3,ELEM(:));
0518   zzz= S*zzz+ (1-S)*ones(4,1)*mean(zzz);
0519   plot3( xxx([1 2 3 1 4 2 3 4],:), ...
0520          yyy([1 2 3 1 4 2 3 4],:), ...
0521          zzz([1 2 3 1 4 2 3 4],:),'k' );
0522   hold on;
0523   xy=NODE(1:2,MES(1,:));
0524   plot3(arrow*xy,arrow*[0 1;-1 0]*xy,ones(8,1)*NODE(3,MES(1,:)),'b') 
0525   hold off;
0526   axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ...
0527          [-1.1 1.1]*max(NODE(3,:))  ])
0528 
0529 
0530 function plot_2d_mesh(NODE,ELEM,el_pos, S, options)
0531 %  if options(1) -> do mesh num
0532 %  if options(2) -> do node num
0533 %  if options(3) -> do text
0534 %  S=.95; % shrink simplices by this factor
0535 
0536   e=size(ELEM,2);
0537   d=size(ELEM,1);
0538   R= zeros(1,e);
0539 
0540 
0541   arrow= [1.02 0;1.06 .05;1.06 .02;1.1 .02; ...
0542           1.1 -.02; 1.06 -.02;1.06 -.05;1.02 0];
0543   % arrow= [1.06 0;1.18 .15;1.18 .06;1.3 .06; ...
0544   %          1.3 -.06; 1.18 -.06;1.18 -.15;1.06 0];
0545   xxx=zeros(3,e); xxx(:)=NODE(1,ELEM(:));
0546   xxx= S*xxx+ (1-S)*ones(3,1)*mean(xxx);
0547   xxx= [xxx;xxx(1,:)];
0548 
0549   yyy=zeros(3,e); yyy(:)=NODE(2,ELEM(:));
0550   yyy= S*yyy+ (1-S)*ones(3,1)*mean(yyy);
0551   yyy= [yyy;yyy(1,:)];
0552 
0553   if isempty(el_pos)
0554       plot(xxx,yyy,'b');
0555   else
0556       xy= el_pos;
0557       hh=plot([xxx;xxx(1,:)],[yyy;yyy(1,:)],'b', ...
0558               arrow*xy,arrow*[0 1;-1 0]*xy,'r');
0559       set(hh(find(R>0)),'Color',[1 0 0],'LineWidth',2);
0560       set(hh(find(R<0)),'Color',[0 0 1],'LineWidth',2);
0561   end
0562 
0563   if options(1) %domesnum
0564     mesnum= reshape(sprintf(' %-2d',[0:15]),3,16)';
0565     text(NODE(1,MES(1,:))*1.08,NODE(2,MES(1,:))*1.08,mesnum, ...
0566          'Color','green','HorizontalAlignment','center');
0567   end %if domesnum
0568 
0569   if options(2) %donodenum
0570     nodenum= reshape(sprintf('%3d',1:length(NODE)),3,length(NODE))';
0571     text(NODE(1,:),NODE(2,:),nodenum, ...
0572          'Color','yellow','HorizontalAlignment','center','FontSize',14);
0573   end %if domesnum
0574 
0575   if options(3) %dotext
0576     numeros= reshape(sprintf('%4d',[1:e]),4,e)';
0577 %   decal= ( 2-0.2*floor(log10([1:e]')))*sqrt(axis*axis'/4)/100;
0578     decal= .02;
0579     xcoor=mean(NODE(2*ELEM-1))';
0580     ycoor=mean(NODE(2*ELEM))';
0581     text(xcoor-decal,ycoor,numeros,'FontSize',8, ...
0582          'HorizontalAlignment','center');
0583 
0584   end  % if nargin~=0
0585   axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ])
0586 
0587 function  mes= avg_electrode_posn( mdl )
0588    if ~isfield(mdl,'electrode'); mes=[]; return; end
0589    n_elec= length( mdl.electrode );
0590    nodes = mdl.nodes;
0591    mes= zeros( n_elec, mdl_dim(mdl) )
0592    for i= 1:n_elec
0593       e_nodes=  mdl.electrode(i).nodes;
0594       mes(i,:)= mean( nodes(e_nodes,:) , 1 );
0595    end
0596 
0597 function colour= electr_colour( e);
0598     if e==1;
0599        colour = [0,.7,0]; % light green electrode #1
0600     elseif e==2
0601        colour = [0,.5,0]; % mid-green electrode #2
0602     else
0603        colour = [0,.3,0]; % dark green
0604     end
0605 
0606 function ee= get_boundary( mdl )
0607    if isfield(mdl,'boundary')
0608        ee= mdl.boundary;
0609    else
0610        % calc and cache boundary
0611        ee = find_boundary( mdl.elems );
0612    end
0613 
0614 
0615 % TESTS:
0616 function do_unit_test
0617      ver= eidors_obj('interpreter_version');
0618    clf
0619 
0620    img=calc_jacobian_bkgnd(mk_common_model('a2c0',8)); 
0621    img.elem_data=rand(size(img.fwd_model.elems,1),1);
0622    subplot(3,4,1); show_fem(img.fwd_model,[0,0,1]) 
0623    title('regular mesh numbered');
0624 
0625 if ~ver.isoctave 
0626    imgn = rmfield(img,'elem_data');
0627    imgn.node_data=rand(size(img.fwd_model.nodes,1),1);
0628    subplot(3,4,9); show_fem(imgn) 
0629    title('interpolated node colours');
0630 end
0631 
0632    img2(1) = img; img2(2) = img;
0633    subplot(3,4,2); show_fem(img,[1]);
0634    title('colours with legend');
0635    subplot(3,4,3); show_fem(img2,[0,1]);
0636    title('colours with legend');
0637    img.calc_colours.mapped_colour = 0; % USE RGB colours
0638    subplot(3,4,4); show_fem(img,[0,1,1]);
0639    title('RGB colours');
0640    subplot(3,4,4); show_fem(img);
0641    title('RGB colours');
0642 
0643    img.elem_data = [1:10];
0644    subplot(3,4,12);show_fem(img); %Should show grey
0645    title('error -> show grey');
0646 
0647 if ~ver.isoctave
0648    imgn.calc_colours.mapped_colour = 0; % USE RGB colours
0649    subplot(3,4,10);show_fem(imgn,[0,1]) 
0650    title('interpolated node colours');
0651 
0652 
0653    subplot(3,4,11);hh=show_fem(imgn); set(hh,'EdgeColor',[0,0,1]);
0654    title('with edge colours');
0655 end
0656 
0657    img3=calc_jacobian_bkgnd(mk_common_model('n3r2',[16,2]));
0658    img3.elem_data= randn(828,1);                       
0659    subplot(3,4,5); show_fem(img3.fwd_model) 
0660    subplot(3,4,6); show_fem(img3,[1])
0661    subplot(3,4,7); show_fem(img3,[1,1])
0662    subplot(3,4,8); show_fem(img3,[1,1,1])
0663

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005