show_fem

 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 (fontsize can be controlled by
      the fractional part, ie [0,1.012] => 12pt font)
   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
  img.fwd_model.show_fem.alpha_inhomogeneities
      1 is opaque. A value like 0.02 is transparent     

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

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