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

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005