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 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

 to change colours, try 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 specifies a set of options
0008 %   options(1) => show colourbar
0009 %   options(2) => show numbering on electrodes
0010 %   options(3) => number elements (==1) or nodes (==2);
0011 %
0012 % for detailed control of colours, use
0013 %    img.calc_colours."parameter" = value
0014 % see help for calc_colours.
0015 %
0016 % for control of colourbar, use img.calc_colours.cb_shrink_move
0017 %
0018 % to change colours, try hh=show_fem(...); set(hh,'EdgeColor',[0,0,1];
0019 
0020 % (C) 2005-2011 Andy Adler. License: GPL version 2 or version 3
0021 % $Id: show_fem.html 2819 2011-09-07 16:43:11Z aadler $
0022 
0023 % TESTS
0024 switch nargin
0025    case 0;
0026      error('Insufficient parameters for show_fem');
0027    case 1;
0028      if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0029      options = [];
0030 end
0031 [img,mdl,opts] = proc_params( mdl, options );
0032 
0033 if ~ishold
0034     cla;
0035 end
0036 
0037 
0038 switch opts.dims
0039    case 2;    hh= show_2d(img,mdl,opts);
0040    case 3;    hh= show_3d(img,mdl,opts);
0041    otherwise; error('model is not 2D or 3D');
0042 end
0043 
0044 if opts.show_numbering
0045    if     bitand( opts.show_numbering, 1 )
0046       xyzc= interp_mesh(mdl);
0047    elseif bitand( opts.show_numbering, 2 )
0048       xyzc= mdl.nodes;
0049    else
0050       error('don''t understand show_numbering value of %d', opts.show_numbering);
0051    end
0052    xyzc= xyzc * eye(size(xyzc,2),3); %convert to 3D
0053    for i= 1:size(xyzc,1)
0054       text(xyzc(i,1),xyzc(i,2), xyzc(i,3), num2str(i), ...
0055             'HorizontalAlignment','center','FontSize',7);
0056    end
0057 end
0058 
0059 if nargout == 0; clear hh; end
0060 
0061 
0062 function [img,mdl,opts] = proc_params( mdl, options );
0063 
0064    opts.do_colourbar=0;
0065    opts.number_electrodes=0;
0066    opts.show_numbering  =0;
0067    if nargin >=2
0068        % fill in default options
0069        optionstr= zeros(1,100);
0070        optionstr(1:length(options)) = options;
0071 
0072        opts.do_colourbar=      optionstr(1);
0073        opts.number_electrodes= optionstr(2);
0074        opts.show_numbering   =  optionstr(3);
0075    end
0076 
0077    % if we have an only img input, then define mdl
0078    if strcmp( mdl(1).type , 'image' )
0079       img= mdl;
0080       mdl= img.fwd_model;
0081    else 
0082       img = [];
0083    end
0084 
0085    opts.dims = size(mdl.nodes,2);
0086 
0087 % 2D Case
0088 function hh= show_2d(img,mdl,opts)
0089    hax= gca;
0090    pax= get(hax,'position');
0091    if ~isempty(img)
0092       colours= calc_colours(img, [] );
0093    else
0094       colours= [1,1,1]; % white elements if no image
0095    end
0096    hh= show_2d_fem( mdl, colours );
0097    show_electrodes_2d(mdl, opts.number_electrodes);
0098 
0099    set(hax,'position', pax);
0100    view(0, 90); axis('xy'); grid('off');
0101 
0102 % IN MATLAB 7 WE NEED TO RERUN THIS BECAUSE STUPID STUPID
0103 % MATLAB WILL RESET THE COLOURBAR EVERY TIME WE RUN PATCH!!!
0104    if exist('img','var') && opts.do_colourbar;
0105       colours= calc_colours(img, [], opts.do_colourbar);
0106       % Matlab is so weird. It puts the first colorbar in the wrong place
0107       %   sometimes ...  (tested with 6.5 and with 7.8)
0108       %   The trick is to never try to move it on the first go
0109       %   OR we reset it and then replace it. STUPID STUPID
0110 
0111      % Here's the magic trick I found. Force a drawnow,
0112      %     then delete and recreate
0113       if ~exist('OCTAVE_VERSION')
0114          drawnow; colorbar('delete');
0115          colours= calc_colours(img, [], opts.do_colourbar);
0116       end
0117    end
0118 
0119    
0120 
0121 % 3D Case
0122 function hh= show_3d(img,mdl,opts)
0123    hh= show_3d_fem( mdl );
0124 
0125    if ~isempty(img)
0126        elem_data = get_img_data(img);
0127        show_inhomogeneities( elem_data , mdl, img);
0128        if opts.do_colourbar
0129            calc_colours(img, [], opts.do_colourbar);
0130        end
0131    end
0132 
0133    show_electrodes_3d(mdl, opts.number_electrodes);
0134 
0135 function show_electrodes_2d(mdl, number_electrodes)
0136     if ~isfield(mdl,'electrode'); return; end
0137 
0138     ee= get_boundary( mdl );
0139     ctr_x= mean(mdl.nodes(:,1));
0140     ctr_y= mean(mdl.nodes(:,2));
0141 
0142 % scale away from model
0143 
0144 for e=1:length(mdl.electrode)
0145     elec_nodes= mdl.electrode(e).nodes;
0146 
0147     S= 1.00;
0148     vx= (mdl.nodes(elec_nodes,1) - ctr_x)*S;
0149     vy= (mdl.nodes(elec_nodes,2) - ctr_y)*S;
0150     % sort nodes around the model (to avoid crossed lines)
0151     [jnk,idx] = sort(unwrap(atan2( vy, vx )));
0152     ecolour = electr_colour( e );
0153     if(length(elec_nodes) == 1)
0154        % Point Electrode Models: put a circle around the node
0155        line(vx(idx)+ctr_x,vy(idx)+ctr_y,  ...
0156             'LineWidth', 2, 'LineStyle','-','Color', ecolour, ...
0157             'Marker','o','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0158     else
0159        % Complete/Shunt Electrode Models (multiple nodes per electrode)
0160        %  put a line along the edges that form the electrode
0161        line(vx(idx)+ctr_x,vy(idx)+ctr_y,  ...
0162             'LineWidth', 3, 'LineStyle','-','Color', ecolour, ...
0163             'Marker','none','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0164     end
0165     if number_electrodes
0166        S= 1.05;
0167        vx= (mdl.nodes(elec_nodes,1) - ctr_x)*S;
0168        vy= (mdl.nodes(elec_nodes,2) - ctr_y)*S;
0169        hh= text(mean(vx), mean(vy), num2str(e));
0170        set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0171     end
0172 end
0173 
0174 function show_electrodes_3d(mdl, number_electrodes);
0175 % show electrode positions on model
0176 if ~isfield(mdl,'electrode'); return; end
0177 
0178 ee= get_boundary( mdl );
0179 for e=1:length(mdl.electrode)
0180     elec_nodes= mdl.electrode(e).nodes;
0181 
0182     colour= electr_colour( e);
0183 
0184     if length(elec_nodes) == 1  % point electrode model
0185         vtx= mdl.nodes(elec_nodes,:);
0186         line(vtx(1),vtx(2),vtx(3), ...
0187             'Marker','o','MarkerSize',12, ...
0188             'MarkerFaceColor',colour, 'MarkerEdgeColor', colour);
0189         if number_electrodes
0190             hh= text(vtx(1),vtx(2),vtx(3), num2str(e));
0191             set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0192         end
0193     else
0194         % find elems on boundary attached to this electrode
0195         nn=ones(size(ee,1),1)*mdl.electrode(e).nodes(:)';
0196         oo=ones(1,size(nn,2));
0197         ec=zeros(size(ee));
0198         for i=1:size(ec,2);
0199            ec(:,i) = any( ee(:,i*oo)==nn, 2);
0200         end
0201         sels= find(all(ec'));
0202 
0203         for u=1:length(sels)
0204             ee= get_boundary( mdl );
0205             paint_electrodes(sels(u), ee, ...
0206                              mdl.nodes, colour, ...
0207                              number_electrodes);
0208         end
0209         if number_electrodes
0210             el_nodes= mdl.nodes(unique(mdl.boundary(sels,:)),:);
0211             hh=text( mean(el_nodes(:,1)), ...
0212                      mean(el_nodes(:,2)), ...
0213                      mean(el_nodes(:,3)), num2str(e) );
0214             set(hh,'FontWeight','bold','Color',[.6,0,0]);
0215         end
0216     end
0217 end
0218 
0219 function show_inhomogeneities( elem_data, mdl, img)
0220 % show
0221 if size(elem_data,2)>1
0222    eidors_msg('warning: show_fem only shows first image',1);
0223 end
0224 repaint_inho(elem_data(:,1), 'use_global' , mdl.nodes, mdl.elems, [], img); 
0225 if ~exist('OCTAVE_VERSION');
0226 camlight('left');
0227 lighting('none'); % lighting doesn't help much
0228 end
0229 
0230 function paint_electrodes(sel,srf,vtx, colour, show_num);
0231 %function paint_electrodes(sel,srf,vtx);
0232 %
0233 % plots the electrodes red at the boundaries.
0234 %
0235 % sel = The index of the electrode faces in the srf matrix
0236 %       sel can be created by set_electrodes.m
0237 % srf = the boundary faces (triangles)
0238 % vtx = The vertices matrix.
0239 
0240 l = srf(sel,1); m = srf(sel,2); n = srf(sel,3);
0241 
0242 Xs = [vtx(l,1);vtx(m,1);vtx(n,1)];
0243 Ys = [vtx(l,2);vtx(m,2);vtx(n,2)];
0244 Zs = [vtx(l,3);vtx(m,3);vtx(n,3)];
0245 
0246 h=patch(Xs,Ys,Zs, colour);
0247 % need 'direct' otherwise colourmap is screwed up
0248 set(h, 'FaceLighting','none', 'CDataMapping', 'direct' );
0249 
0250 function hh= show_3d_fem( mdl, options )
0251    ee= get_boundary( mdl );
0252    hh= trimesh(ee, mdl.nodes(:,1), ...
0253                    mdl.nodes(:,2), ...
0254                    mdl.nodes(:,3));
0255    set(hh, 'EdgeColor', [0,0,0]);
0256    axis('image');
0257    set(gcf,'Colormap',[0 0 0]);
0258    hidden('off');
0259 
0260 function hh=show_2d_fem( mdl, colours, imgno )
0261   szcolr = size(colours);
0262   if nargin<3;
0263       imgno = 1;
0264       if szcolr(1:2)>1;  eidors_msg('warning: show_fem only shows first image',1); end
0265   end
0266   
0267 % el_pos= avg_electrode_posn( mdl );
0268 % plot_2d_mesh(mdl.nodes', mdl.elems', el_pos', .95, [0,0,0]);
0269 
0270   S= 1; %.95; % shrink factor
0271   elem= mdl.elems'; e= size(elem,2);
0272   node= mdl.nodes'; n= size(node,2);
0273   Xs=zeros(3,e);
0274   Xs(:)=mdl.nodes(elem(:),1);
0275   Xs= S*Xs+ (1-S)*ones(3,1)*mean(Xs);
0276   Ys=zeros(3,e); Ys(:)=mdl.nodes(elem(:),2);
0277   Ys= S*Ys+ (1-S)*ones(3,1)*mean(Ys);
0278   Zs = zeros(size(Xs));
0279   if szcolr(1:2) == [1,3]  % no reconstruction  - just use white
0280      hh= patch(Xs,Ys,zeros(3,e),colours);
0281   elseif size(colours,1) == e % defined on elems
0282 % THE STUPID MATLAB 7 WILL RESET THE COLOURBAR WHENEVER YOU DO A PATCH. DAMN.
0283      colours = permute(colours(:,imgno,:),[2,1,3]);
0284   elseif size(colours,1) == n  % multiple images
0285      colours = colours(elem(:), imgno, :);
0286      colours = reshape( colours, size(elem,1), size(elem,2), []);
0287   else
0288     eidors_msg('warning: image elements and mesh do not match. Showing grey',1);
0289     colours= [1,1,1]/2; % Grey to show we're not sure
0290   end
0291 
0292   hh= patch(Xs,Ys,Zs,colours);
0293   set(hh, 'FaceLighting','none', 'CDataMapping', 'direct' );
0294   % FOR NODE RGB MATLAB SCREWS UP THE COLOURS FOR US (ONCE AGAIN). THERE MUST
0295   % BE SOME KIND OF SECRET FLAG@@@
0296 
0297   max_x= max(mdl.nodes(:,1)); min_x= min(mdl.nodes(:,1));
0298   max_y= max(mdl.nodes(:,2)); min_y= min(mdl.nodes(:,2));
0299   axis([ mean([max_x,min_x]) + 0.55*[-1,1]*(max_x-min_x), ...
0300          mean([max_y,min_y]) + 0.55*[-1,1]*(max_y-min_y) ]);
0301 
0302 
0303 
0304 function old_code_3d_plot
0305   domesnum= 0;
0306   donodenum= 0;
0307   dotext=0;
0308 
0309   xxx=zeros(4,e); xxx(:)=NODE(1,ELEM(:));
0310   xxx= S*xxx+ (1-S)*ones(4,1)*mean(xxx);
0311   yyy=zeros(4,e); yyy(:)=NODE(2,ELEM(:));
0312   yyy= S*yyy+ (1-S)*ones(4,1)*mean(yyy);
0313   zzz=zeros(4,e); zzz(:)=NODE(3,ELEM(:));
0314   zzz= S*zzz+ (1-S)*ones(4,1)*mean(zzz);
0315   plot3( xxx([1 2 3 1 4 2 3 4],:), ...
0316          yyy([1 2 3 1 4 2 3 4],:), ...
0317          zzz([1 2 3 1 4 2 3 4],:),'k' );
0318   hold on;
0319   xy=NODE(1:2,MES(1,:));
0320   plot3(arrow*xy,arrow*[0 1;-1 0]*xy,ones(8,1)*NODE(3,MES(1,:)),'b') 
0321   hold off;
0322   axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ...
0323          [-1.1 1.1]*max(NODE(3,:))  ])
0324 
0325 
0326 function plot_2d_mesh(NODE,ELEM,el_pos, S, options)
0327 %  if options(1) -> do mesh num
0328 %  if options(2) -> do node num
0329 %  if options(3) -> do text
0330 %  S=.95; % shrink simplices by this factor
0331 
0332   e=size(ELEM,2);
0333   d=size(ELEM,1);
0334   R= zeros(1,e);
0335 
0336 
0337   arrow= [1.02 0;1.06 .05;1.06 .02;1.1 .02; ...
0338           1.1 -.02; 1.06 -.02;1.06 -.05;1.02 0];
0339   % arrow= [1.06 0;1.18 .15;1.18 .06;1.3 .06; ...
0340   %          1.3 -.06; 1.18 -.06;1.18 -.15;1.06 0];
0341   xxx=zeros(3,e); xxx(:)=NODE(1,ELEM(:));
0342   xxx= S*xxx+ (1-S)*ones(3,1)*mean(xxx);
0343   xxx= [xxx;xxx(1,:)];
0344 
0345   yyy=zeros(3,e); yyy(:)=NODE(2,ELEM(:));
0346   yyy= S*yyy+ (1-S)*ones(3,1)*mean(yyy);
0347   yyy= [yyy;yyy(1,:)];
0348 
0349   if isempty(el_pos)
0350       plot(xxx,yyy,'b');
0351   else
0352       xy= el_pos;
0353       hh=plot([xxx;xxx(1,:)],[yyy;yyy(1,:)],'b', ...
0354               arrow*xy,arrow*[0 1;-1 0]*xy,'r');
0355       set(hh(find(R>0)),'Color',[1 0 0],'LineWidth',2);
0356       set(hh(find(R<0)),'Color',[0 0 1],'LineWidth',2);
0357   end
0358 
0359   if options(1) %domesnum
0360     mesnum= reshape(sprintf(' %-2d',[0:15]),3,16)';
0361     text(NODE(1,MES(1,:))*1.08,NODE(2,MES(1,:))*1.08,mesnum, ...
0362          'Color','green','HorizontalAlignment','center');
0363   end %if domesnum
0364 
0365   if options(2) %donodenum
0366     nodenum= reshape(sprintf('%3d',1:length(NODE)),3,length(NODE))';
0367     text(NODE(1,:),NODE(2,:),nodenum, ...
0368          'Color','yellow','HorizontalAlignment','center','FontSize',14);
0369   end %if domesnum
0370 
0371   if options(3) %dotext
0372     numeros= reshape(sprintf('%4d',[1:e]),4,e)';
0373 %   decal= ( 2-0.2*floor(log10([1:e]')))*sqrt(axis*axis'/4)/100;
0374     decal= .02;
0375     xcoor=mean(NODE(2*ELEM-1))';
0376     ycoor=mean(NODE(2*ELEM))';
0377     text(xcoor-decal,ycoor,numeros,'FontSize',8, ...
0378          'HorizontalAlignment','center');
0379 
0380   end  % if nargin~=0
0381   axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ])
0382 
0383 function  mes= avg_electrode_posn( mdl )
0384    if ~isfield(mdl,'electrode'); mes=[]; return; end
0385    n_elec= length( mdl.electrode );
0386    nodes = mdl.nodes;
0387    mes= zeros( n_elec, mdl_dim(mdl) )
0388    for i= 1:n_elec
0389       e_nodes=  mdl.electrode(i).nodes;
0390       mes(i,:)= mean( nodes(e_nodes,:) , 1 );
0391    end
0392 
0393 function colour= electr_colour( e);
0394     if e==1;
0395        colour = [0,.7,0]; % light green electrode #1
0396     elseif e==2
0397        colour = [0,.5,0]; % mid-green electrode #2
0398     else
0399        colour = [0,.3,0]; % dark green
0400     end
0401 
0402 function ee= get_boundary( mdl )
0403    if isfield(mdl,'boundary')
0404        ee= mdl.boundary;
0405    else
0406        % calc and cache boundary
0407        ee = find_boundary( mdl.elems );
0408    end
0409 
0410 
0411 % TESTS:
0412 function do_unit_test
0413    clf
0414 
0415    img=calc_jacobian_bkgnd(mk_common_model('a2c0',8)); 
0416    img.elem_data=rand(size(img.fwd_model.elems,1),1);
0417    subplot(3,4,1); show_fem(img.fwd_model,[0,0,1]) 
0418    title('regular mesh numbered');
0419 
0420    imgn = rmfield(img,'elem_data');
0421    imgn.node_data=rand(size(img.fwd_model.nodes,1),1);
0422    subplot(3,4,9); show_fem(imgn) 
0423    title('interpolated node colours');
0424 
0425    img2(1) = img; img2(2) = img;
0426    subplot(3,4,2); show_fem(img,[1]);
0427    title('colours with legend');
0428    subplot(3,4,3); show_fem(img2,[0,1]);
0429    title('colours with legend');
0430    img.calc_colours.mapped_colour = 0; % USE RGB colours
0431    subplot(3,4,4); show_fem(img,[0,1,1]);
0432    title('RGB colours');
0433    subplot(3,4,4); show_fem(img);
0434    title('RGB colours');
0435 
0436    imgn.calc_colours.mapped_colour = 0; % USE RGB colours
0437    subplot(3,4,10);show_fem(imgn,[0,1]) 
0438    title('interpolated node colours');
0439 
0440 
0441    subplot(3,4,11);hh=show_fem(imgn); set(hh,'EdgeColor',[0,0,1]);
0442    title('with edge colours');
0443 
0444    imgn.node_data = [1:10];
0445    subplot(3,4,12);show_fem(imgn); %Should show grey
0446    title('error -> show grey');
0447 
0448    img3=calc_jacobian_bkgnd(mk_common_model('n3r2',8));
0449    img3.elem_data= randn(828,1);                       
0450    subplot(3,4,5); show_fem(img3.fwd_model) 
0451    subplot(3,4,6); show_fem(img3,[1])
0452    subplot(3,4,7); show_fem(img3,[1,1])
0453    subplot(3,4,8); show_fem(img3,[1,1,1])
0454