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 (fontsize can be controlled by
      the fractional part, ie [0,1.012] => 12pt font)
   options(3) => number elements (==1) or nodes (==2);
      fontsize can be controlled as above

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

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

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005