show_fem_enhanced

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

 SHOW_FEM_ENHANCED('refresh') repaints the figure in current axis orientation 
   (only works if options.edge.significant.viewpoint_dependent.callback 
   is true)

 options may be specified by a list (for compatibility purposes)

 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 the following fields (default values
 are indicated in parenthesis)

     options.viewpoint.az (-37.5)
     options.viewpoint.el (30.0)
 
     options.edge.color ([0 0 0])
     options.edge.width (0.5)
     options.edge.significant.show (true)
     options.edge.significant.color ([0 0 0])
     options.edge.significant.width (2)
     options.edge.significant.angle (45)
     options.edge.significant.viewpoint_dependent.show (true)
     options.edge.significant.viewpoint_dependent.color ([0 0 0])
     options.edge.significant.viewpoint_dependent.width (2)
     options.edge.significant.viewpoint_dependent.callback (true)
  
     options.colorbar.show (false)
 
     options.electrode.number.show (false)
     options.electrode.number.font_size (10)
     options.electrode.number.font_weight ('bold')
     options.electrode.number.color ([.6 0 0])
     options.electrode.number.background_color ('none')
     options.electrode.number.scale_position (1.15)
 
     options.element.number.show (false)
     options.element.number.font_size (7)
     options.element.number.font_weight ('normal')
     options.element.number.color ([0 0 0])
     options.element.number.background_color ('none')

     options.node.number.show (false)
     options.node.number.font_size (7)
     options.node.number.font_weight ('normal')
     options.node.number.color ([0 0 0.5])
     options.node.number.background_color ([1 1 1])

 EXAMPLES
   mdl = mk_common_model('c2C2',8);
   img = mk_image(mdl, linspace(-3,3,num_elems(mdl)));
   
 SET PARAMETERS
   img.show_fem.electrode.number.scale_position= 1.01;
   img.show_fem.electrode.number.show =1;
   show_fem(img)

 SET OPTIONS
   opts.colorbar.show = 1;
   show_fem(img.opts)

0001 function hh = show_fem_enhanced(mdl, options)
0002 % SHOW_FEM_ENHANCED: show the EIDORS3D finite element model
0003 % hh = show_fem(mdl, options)
0004 % mdl is an EIDORS3D 'model' or 'image' structure
0005 % hh = handle to the plotted model
0006 %
0007 % SHOW_FEM_ENHANCED('refresh') repaints the figure in current axis orientation
0008 %   (only works if options.edge.significant.viewpoint_dependent.callback
0009 %   is true)
0010 %
0011 % options may be specified by a list (for compatibility purposes)
0012 %
0013 % options specifies a set of options
0014 %   options(1) => show colourbar
0015 %   options(2) => show numbering on electrodes
0016 %   options(3) => number elements (==1) or nodes (==2);
0017 %
0018 % for detailed control of colours, use the following fields (default values
0019 % are indicated in parenthesis)
0020 %
0021 %     options.viewpoint.az (-37.5)
0022 %     options.viewpoint.el (30.0)
0023 %
0024 %     options.edge.color ([0 0 0])
0025 %     options.edge.width (0.5)
0026 %     options.edge.significant.show (true)
0027 %     options.edge.significant.color ([0 0 0])
0028 %     options.edge.significant.width (2)
0029 %     options.edge.significant.angle (45)
0030 %     options.edge.significant.viewpoint_dependent.show (true)
0031 %     options.edge.significant.viewpoint_dependent.color ([0 0 0])
0032 %     options.edge.significant.viewpoint_dependent.width (2)
0033 %     options.edge.significant.viewpoint_dependent.callback (true)
0034 %
0035 %     options.colorbar.show (false)
0036 %
0037 %     options.electrode.number.show (false)
0038 %     options.electrode.number.font_size (10)
0039 %     options.electrode.number.font_weight ('bold')
0040 %     options.electrode.number.color ([.6 0 0])
0041 %     options.electrode.number.background_color ('none')
0042 %     options.electrode.number.scale_position (1.15)
0043 %
0044 %     options.element.number.show (false)
0045 %     options.element.number.font_size (7)
0046 %     options.element.number.font_weight ('normal')
0047 %     options.element.number.color ([0 0 0])
0048 %     options.element.number.background_color ('none')
0049 %
0050 %     options.node.number.show (false)
0051 %     options.node.number.font_size (7)
0052 %     options.node.number.font_weight ('normal')
0053 %     options.node.number.color ([0 0 0.5])
0054 %     options.node.number.background_color ([1 1 1])
0055 %
0056 % EXAMPLES
0057 %   mdl = mk_common_model('c2C2',8);
0058 %   img = mk_image(mdl, linspace(-3,3,num_elems(mdl)));
0059 %
0060 % SET PARAMETERS
0061 %   img.show_fem.electrode.number.scale_position= 1.01;
0062 %   img.show_fem.electrode.number.show =1;
0063 %   show_fem(img)
0064 %
0065 % SET OPTIONS
0066 %   opts.colorbar.show = 1;
0067 %   show_fem(img.opts)
0068 
0069 % (C) 2015 Herve Gagnon. License: GPL version 2 or version 3
0070 % $Id: show_fem_enhanced.m 6342 2022-04-22 11:22:16Z bgrychtol $
0071 
0072 % TESTS
0073 switch nargin
0074     case 0
0075         error('Insufficient parameters for show_fem');
0076     case 1
0077         if ischar(mdl) && strcmp(mdl,'UNIT_TEST')
0078             do_unit_test; 
0079             return; 
0080         end
0081         if ischar(mdl) && strcmp(mdl,'refresh')
0082             refresh_current_axis; 
0083             return; 
0084         end
0085         if (isstruct(mdl) && isfield(mdl, 'show_fem'))
0086             options = mdl.show_fem;
0087         else
0088             options = [];
0089         end
0090     case 2
0091         % No special processing
0092     otherwise
0093         error('Too many parameters for show_fem');
0094 end
0095 
0096 [img, mdl, opts] = proc_params(mdl, options);
0097 
0098 if ~ishold
0099     cla;
0100 end
0101 
0102 mdl = find_sub_elements(mdl);
0103 
0104 % Convert nodes to 3D if necessary.
0105 mdl.nodes = mdl.nodes*eye(size(mdl.nodes, 2), 3);
0106 
0107 hh = draw_all(img, mdl, opts);
0108 
0109 % Setup callback function.
0110 if (opts.edge.significant.viewpoint_dependent.callback)
0111     if ~exist('OCTAVE_VERSION'); %octave's not compativble (4.4.1)
0112     h3d = rotate3d;
0113     set(gca, 'UserData', struct('name', 'show_fem_data', 'img', img, 'mdl', mdl, 'opts', opts, 'handles', hh));
0114     set(h3d, 'ActionPostCallback' , @refresh_current_axis)
0115     end
0116 end
0117 
0118 if (nargout == 0)
0119     clear hh; 
0120 end
0121 
0122 function [img, mdl, opts] = proc_params(mdl, src_opts)
0123 
0124     % Assign default viewpoint option values.
0125     opts.viewpoint.az = -37.5;
0126     opts.viewpoint.el =  30.0;
0127 
0128     % Assign default edge option values.
0129     opts.edge.color   = [0 0 0];
0130     opts.edge.width   = 0.5;
0131     opts.edge.significant.show  = true;
0132     opts.edge.significant.color = [0 0 0];
0133     opts.edge.significant.width = 2;
0134     opts.edge.significant.angle = 45;
0135     opts.edge.significant.viewpoint_dependent.show     = true;
0136     opts.edge.significant.viewpoint_dependent.color    = [0 0 0];
0137     opts.edge.significant.viewpoint_dependent.width    = 2;
0138     opts.edge.significant.viewpoint_dependent.callback = true;
0139  
0140     % Assign default colorbar option values.
0141     opts.colorbar.show = false;
0142 
0143     % Assign default electrode option values.
0144     opts.electrode.number.show             = false;
0145     opts.electrode.number.font_size        = 10;
0146     opts.electrode.number.font_weight      = 'bold';
0147     opts.electrode.number.color            = [.6 0 0];
0148     opts.electrode.number.background_color = 'none';
0149     opts.electrode.number.scale_position   = 1.15;
0150 
0151     % Assign default element option values.
0152     opts.element.number.show             = false;
0153     opts.element.number.font_size        = 7;
0154     opts.element.number.font_weight      = 'normal';
0155     opts.element.number.color            = [0 0 0];
0156     opts.element.number.background_color = 'none';
0157 
0158     % Assign default node option values.
0159     opts.node.number.show             = false;
0160     opts.node.number.font_size        = 7;
0161     opts.node.number.font_weight      = 'normal';
0162     opts.node.number.color            = [0 0 0.5];
0163     opts.node.number.background_color = [1 1 1];
0164 
0165     % Override some defaults in 2D
0166     if mdl_dim( mdl ) == 2
0167        opts.viewpoint.az = 0;
0168        opts.viewpoint.el = 90;
0169        opts.electrode.number.scale_position= 1.05;
0170     end
0171     
0172     if (nargin == 2)
0173         if (isstruct(src_opts))
0174             opts = copy_field(opts, src_opts, 'viewpoint', 'az');
0175             opts = copy_field(opts, src_opts, 'viewpoint', 'el');
0176             opts = copy_field(opts, src_opts, 'edge', 'color');
0177             opts = copy_field(opts, src_opts, 'edge', 'width');
0178             
0179             if (isfield(src_opts, 'edge'))
0180                 opts.edge = copy_field(opts.edge, src_opts.edge, 'significant', 'show');
0181                 opts.edge = copy_field(opts.edge, src_opts.edge, 'significant', 'color');
0182                 opts.edge = copy_field(opts.edge, src_opts.edge, 'significant', 'width');
0183                 opts.edge = copy_field(opts.edge, src_opts.edge, 'significant', 'angle');
0184                 if (isfield(src_opts.edge, 'significant'))
0185                     opts.edge.significant = copy_field(opts.edge.significant, src_opts.edge.significant, 'viewpoint_dependent', 'show');
0186                     opts.edge.significant = copy_field(opts.edge.significant, src_opts.edge.significant, 'viewpoint_dependent', 'color');
0187                     opts.edge.significant = copy_field(opts.edge.significant, src_opts.edge.significant, 'viewpoint_dependent', 'width');
0188                     opts.edge.significant = copy_field(opts.edge.significant, src_opts.edge.significant, 'viewpoint_dependent', 'callback');
0189                 end
0190             end
0191 
0192             opts = copy_field(opts, src_opts, 'colorbar', 'show');
0193 
0194             if (isfield(src_opts, 'electrode'))
0195                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'show');
0196                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'font_size');
0197                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'font_weight');
0198                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'color');
0199                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'background_color');
0200                 opts.electrode = copy_field(opts.electrode, src_opts.electrode, 'number', 'scale_position');                
0201             end
0202             
0203             if (isfield(src_opts, 'element'))
0204                 opts.element = copy_field(opts.element, src_opts.element, 'number', 'show');
0205                 opts.element = copy_field(opts.element, src_opts.element, 'number', 'font_size');
0206                 opts.element = copy_field(opts.element, src_opts.element, 'number', 'font_weight');
0207                 opts.element = copy_field(opts.element, src_opts.element, 'number', 'color');
0208                 opts.element = copy_field(opts.element, src_opts.element, 'number', 'background_color');               
0209             end
0210 
0211             if (isfield(src_opts, 'node'))
0212                 opts.node = copy_field(opts.node, src_opts.node, 'number', 'show');
0213                 opts.node = copy_field(opts.node, src_opts.node, 'number', 'font_size');
0214                 opts.node = copy_field(opts.node, src_opts.node, 'number', 'font_weight');
0215                 opts.node = copy_field(opts.node, src_opts.node, 'number', 'color');
0216                 opts.node = copy_field(opts.node, src_opts.node, 'number', 'background_color');               
0217             end
0218             
0219             
0220         else % Support options from old version of show_fem for compatibility.
0221             % fill in default options
0222             optionstr = zeros(1, 100);
0223             optionstr(1:length(src_opts)) = src_opts;
0224 
0225             opts.colorbar.show         = optionstr(1);
0226             opts.electrode.number.show = optionstr(2);
0227             switch optionstr(3)
0228                 case 1
0229                     opts.element.number.show = 1;
0230                 case 2
0231                     opts.node.number.show = 1;
0232                 case 3
0233                     opts.element.number.show = 1;
0234                     opts.node.number.show = 1;
0235             end
0236         end
0237     end
0238 
0239     % Display first image if several images are available.
0240     if (numel(mdl) > 1)
0241         eidors_msg('warning: show_fem only shows first image',1);
0242         mdl = mdl(1);
0243     end
0244  
0245     % if we have an img input, define mdl
0246     if strcmp(mdl(1).type , 'image' )
0247         img = mdl;
0248         mdl = img.fwd_model;
0249     else
0250         img = [];
0251     end
0252 
0253 function refresh_current_axis(obj, evd)
0254 UserData = get(gca, 'UserData');
0255 if (all(isfield(UserData, {'name', 'img', 'mdl', 'opts'})) && ...
0256         strcmp(UserData.name, 'show_fem_data'))
0257     [az, el] = view(gca);
0258     if (az ~= UserData.opts.viewpoint.az || el ~= UserData.opts.viewpoint.el)
0259         UserData.opts.viewpoint.az = az;
0260         UserData.opts.viewpoint.el = el;
0261         delete(UserData.handles);
0262         UserData.handles = draw_all(UserData.img, UserData.mdl, UserData.opts);
0263         set(gca, 'UserData', UserData);
0264     end
0265 end
0266 
0267 function hh = draw_all(img, mdl, opts)
0268 
0269     hh = draw_fem(img, mdl, opts);
0270 
0271     % Number elements if necessary.
0272     if (opts.element.number.show)
0273         h= draw_numbers(interp_mesh(mdl), opts.element.number);
0274         hh = [hh; h];
0275     end
0276 
0277     % Number nodes if necessary.
0278     if (opts.node.number.show)
0279         h = draw_numbers(mdl.nodes, opts.node.number);
0280         hh = [hh; h];
0281     end
0282 
0283     % Number electrodes if necessary.
0284     if (opts.electrode.number.show && isfield(mdl, 'electrode'))
0285         n_elec = numel(mdl.electrode);
0286         mesh_center = ones(n_elec, 1)*mean(mdl.nodes, 1);
0287 
0288         elec_pos = zeros(n_elec, size(mesh_center, 2));
0289 
0290         for i = 1:n_elec
0291             if (isfield(mdl.electrode(i), 'nodes'))
0292                 elec_pos(i, :) = mean(mdl.nodes(mdl.electrode(i).nodes, :), 1);
0293             elseif (isfield(mdl.electrode(i), 'pos'))
0294                 elec_pos(i, :) = mean(mdl.electrode(i).pos, 1)*eye(size(...
0295                                                         elec_pos(i, :), 2), 3);
0296             end
0297         end
0298 
0299         draw_numbers(opts.electrode.number.scale_position*(elec_pos - mesh_center) + mesh_center, opts.electrode.number);
0300     end
0301 
0302     if (size(mdl.elems, 2) == 4)  
0303         view(opts.viewpoint.az, opts.viewpoint.el);
0304     end
0305 
0306     axis image;
0307     
0308 function hh = draw_fem(img, mdl, opts)
0309     
0310     % Assign colors and alpha to elements.
0311     if (size(mdl.elems, 2) == 4)
0312         if ~isempty(img)
0313             elems = mdl.sub_elements;
0314             electrode_field = 'sub_elements';
0315         else
0316             elems = mdl.boundary;
0317             electrode_field = 'boundary';
0318         end
0319         triangle_alpha = 0.5*ones(size(elems, 1), 1);
0320         triangle_edge_color = 'none';
0321         triangle_edge_width = 0.5;
0322         
0323         % Color mesh edges.
0324         edge_edges = mdl.boundary_edges;
0325         edge_width = ones(size(edge_edges, 1), 1)*opts.edge.width;
0326         edge_color = ones(size(edge_edges, 1), 1)*opts.edge.color;
0327         
0328         if opts.edge.significant.show
0329             if opts.edge.significant.viewpoint_dependent.show 
0330                 if exist('viewmtx')
0331                 % Highlight profile of boundary according to viewing angle.
0332                 T = viewmtx(opts.viewpoint.az, opts.viewpoint.el);
0333                 else % octave doesn't have viewmtx yet
0334                 T = ones(4);
0335                 end
0336                 transformed_boundary_normal_vector = mdl.boundary_normal_vector*T(1:3,1:3)';
0337                 boundary_edges_idx2 = (transformed_boundary_normal_vector(mdl.edge2boundary(:, 1), 3).*transformed_boundary_normal_vector(mdl.edge2boundary(:, 2), 3) <= 0);
0338 
0339                 edge_width(boundary_edges_idx2)    = opts.edge.significant.viewpoint_dependent.width;
0340                 edge_color(boundary_edges_idx2, :) = ones(sum(boundary_edges_idx2), 1)*opts.edge.significant.viewpoint_dependent.color;
0341             end
0342 
0343             % Highlight significant edges where boundary shape bends more than
0344             % 45 degrees.
0345             boundary_edges_idx = dot(mdl.boundary_normal_vector(...
0346                                      mdl.edge2boundary(:, 1), :), ...
0347                                      mdl.boundary_normal_vector(...
0348                                      mdl.edge2boundary(:, 2), :), 2) ...
0349                                      <= cosd(opts.edge.significant.angle);
0350             edge_width(boundary_edges_idx)    = opts.edge.significant.width;
0351             edge_color(boundary_edges_idx, :) = ones(sum(boundary_edges_idx), 1)*opts.edge.significant.color;
0352         end
0353     else
0354         elems = mdl.elems;
0355         electrode_field = 'boundary';
0356         triangle_alpha = ones(size(elems, 1), 1);
0357         triangle_edge_color = opts.edge.color;
0358         triangle_edge_width = opts.edge.width;
0359         
0360         % Highlight mesh boundary.
0361         if opts.edge.significant.show
0362             edge_edges = mdl.boundary;
0363             edge_width = opts.edge.significant.width*ones(size(edge_edges, 1), 1);
0364             edge_color = ones(size(edge_edges, 1), 1)*opts.edge.significant.color;
0365         else
0366             edge_edges = [];
0367         end
0368     end
0369     
0370     if ~isempty(img)
0371         if (size(mdl.elems, 2) == 4) 
0372             img_data = get_img_data(img);
0373             if size(img_data,1) == size(mdl.elems,1)
0374                 triangle_color = calc_colours(mdl.element2sub_element*...
0375                                     img_data, img);
0376             elseif size(img_data,1) == size(mdl.nodes,1)
0377                 triangle_color = calc_colours(img_data, img);
0378             else
0379                 error('wrong size of data');
0380             end
0381 
0382             factor = 3/8;
0383             
0384             alpha_map = zeros(size(colormap, 1), 1);
0385             alpha = linspace(0, 1, round(size(colormap, 1)*factor))';
0386             alpha_map(1:size(alpha, 1)) = alpha(end:-1:1);
0387             alpha_map(end:-1:(end-size(alpha, 1)+1)) = alpha(end:-1:1);
0388             
0389             factor = 3/8;
0390             alpha_map2 = ones(size(colormap, 1), 1);
0391             alpha2 = linspace(0, 1, round(size(colormap, 1)*factor))';
0392             alpha_map2((1:size(alpha2, 1)) + size(colormap, 1)/2) = alpha2;
0393             alpha_map2((end:-1:(end-size(alpha2, 1)+1)) - size(colormap, 1)/2) = alpha2;
0394             
0395             factor = 3/8;
0396             alpha_map3 = ones(size(colormap, 1), 1);
0397             alpha3 = linspace(0, 1, round(size(colormap, 1)*factor))';
0398             idx1 = (size(colormap, 1)/2 - size(alpha3, 1))/2;
0399             alpha_map3((-idx1:idx1) + size(colormap, 1)/2) = 0;
0400             alpha_map3((1:size(alpha3, 1)) + size(colormap, 1)/2 + (size(colormap, 1)/2 - size(alpha3, 1))/2) = alpha3;
0401             alpha_map3((end:-1:(end-size(alpha3, 1)+1)) - size(colormap, 1)/2 - (size(colormap, 1)/2 - size(alpha3, 1))/2) = alpha3;
0402     
0403             triangle_alpha = alpha_map3(triangle_color);
0404             
0405             % Restore transparency for boundary elements;
0406             if size(img_data,1) == size(mdl.elems,1)
0407                 alpha_idx =  triangle_alpha(mdl.boundary_to_sub_element_idx, :) < 0.5;
0408                 triangle_alpha(mdl.boundary_to_sub_element_idx(alpha_idx), :) = 0.5;
0409             else
0410                 alpha_idx = triangle_alpha(mdl.boundary_node_idx, :) < 0.5;
0411                 triangle_alpha(mdl.boundary_node_idx(alpha_idx), :) = 0.5;
0412             end
0413         else
0414             triangle_color = calc_colours(img);
0415         end
0416         triangle_color = squeeze(triangle_color(:, 1, :));
0417     else
0418         triangle_color = ones(size(elems, 1), 1)*[1 1 1];
0419     end
0420     
0421     % Assign colors for electrodes.
0422     marker_position = [];
0423     marker_color    = [];
0424     
0425     if (isfield(mdl, 'electrode'))
0426         for i = 1:numel(mdl.electrode)
0427             if (isfield(mdl.electrode(i), electrode_field) && ...
0428                     ~isempty(mdl.electrode(i).(electrode_field)))
0429                 
0430                 if (size(mdl.elems, 2) == 4)
0431                     triangle_color(...
0432                         mdl.electrode(i).(electrode_field), :) = ...
0433                         ones(numel(mdl.electrode(i).(electrode_field)), ...
0434                         1)*electr_colour(mdl.electrode(i),i, size(triangle_color, 2));
0435 
0436                     triangle_alpha(mdl.electrode(i).(electrode_field), ...
0437                                                                     :) = 1;
0438                     % Highlight electrode edges.
0439                     boundary_edges = [mdl.electrode(i).boundary_inner_edges;
0440                                       mdl.electrode(i).boundary_outer_edges];
0441                     edge_color(boundary_edges, :) = 0.5*ones(size(...
0442                         boundary_edges, 1), 1)*electr_colour(mdl.electrode(i), i, 3);
0443                     edge_width(mdl.electrode(i).boundary_outer_edges) = 1;
0444                 else
0445 
0446                     edge_width(mdl.electrode(i).(electrode_field), :) = 3;
0447                     edge_color(mdl.electrode(i).(electrode_field), :) = ...
0448                         ones(numel(mdl.electrode(i).(electrode_field)), ...
0449                         1)*electr_colour(mdl.electrode(i),i, 3);
0450                 end
0451             else
0452                 if (isfield(mdl.electrode(i), 'nodes'))
0453                     marker_position(end + 1, :) = mean(mdl.nodes(...
0454                                             mdl.electrode(i).nodes, :), 1);
0455                 elseif (isfield(mdl.electrode(i), 'pos'))
0456                     marker_position(end + 1, :) = mean(...
0457                                                mdl.electrode(i).pos, 1)*...
0458                                   eye(size(marker_position(end, :), 2), 3);
0459                 end
0460                 marker_color(end + 1, :) = electr_colour(mdl.electrode(i), i, 3);
0461             end
0462         end
0463     end
0464 
0465     % Draw triangles
0466     hh = draw_triangles(elems, mdl.nodes, ...
0467                       triangle_color, triangle_alpha, ...
0468                       triangle_edge_color, triangle_edge_width);
0469 
0470     % Draw edges if necessary
0471     if (~isempty(edge_edges))
0472         h = draw_edges(edge_edges, mdl.nodes, edge_width, edge_color);
0473         hh = [hh; h];
0474     end
0475                   
0476     % Draw markers if necessary.
0477     if (~isempty(marker_position))
0478         h = draw_markers(marker_position, marker_color, 9);
0479         hh = [hh; h];
0480     end
0481 
0482     if opts.colorbar.show && ~isempty( img )
0483 %       psave = get(gca,'position');
0484         eidors_colourbar( img ); 
0485 %       set(gca,'position',psave); %%% Reset axes after colourbar and move
0486     end
0487     
0488 function hh = draw_numbers(positions, opts)
0489     hh = text(positions(:,1), positions(:,2), positions(:,3), ...
0490         arrayfun(@(x){int2str(x)}, 1:size(positions, 1)), ...
0491         'HorizontalAlignment', 'center', ...
0492         'VerticalAlignment',   'middle', ...
0493         'FontSize',            opts.font_size, ...
0494         'FontWeight',          opts.font_weight, ...
0495         'Color',               opts.color, ...
0496         'BackgroundColor',     opts.background_color);
0497 
0498 function hh = draw_markers(points, colour, marker_size)
0499     [unique_colour, unused, colour_idx] = unique(colour, 'rows');
0500     
0501     for i = 1:size(unique_colour, 1)
0502         points_idx = colour_idx == i;
0503         hh = line(points(points_idx, 1), points(points_idx, 2), ...
0504                   points(points_idx, 3), ...
0505                   'LineStyle', 'none', ...
0506                   'Marker', 'o', 'MarkerSize', marker_size, ...
0507                   'MarkerFaceColor', unique_colour(i, :), ...
0508                   'MarkerEdgeColor', unique_colour(i, :));
0509     end
0510         
0511 function hh = draw_edges(edges, vertices, width_data, color_data)
0512     [unique_width_color_data, unused, width_color_idx] = ... 
0513                                    unique([width_data color_data], 'rows'); 
0514     hh = [];
0515     for i = 1:size(unique_width_color_data, 1)
0516         if (unique_width_color_data(i, 1) > 0)
0517             edge_idx = (width_color_idx == i);
0518             points_list = [];
0519             points_list(1:3:3*size(edges(edge_idx, :), 1), :) = ...
0520                                            vertices(edges(edge_idx, 1), :);
0521             points_list(2:3:3*size(edges(edge_idx, :), 1), :) = ...
0522                                            vertices(edges(edge_idx, 2), :);
0523             points_list(3:3:3*size(edges(edge_idx, :), 1), :) = ...
0524                                        nan(size(edges(edge_idx, :), 1), 3);
0525             h = line(points_list(:, 1), points_list(:, 2), ...
0526                       points_list(:, 3), 'LineWidth', ...
0527                       unique_width_color_data(i, 1), 'LineStyle', '-', ...
0528                       'Color', unique_width_color_data(i, 2:end));
0529             hh = [hh; h];
0530         end
0531     end
0532 
0533 function hh = draw_triangles(faces, vertices, color_data, alpha_data, ...
0534                              edge_color, edge_width)
0535                          
0536     if (size(color_data, 1) == 1)
0537         color_data = ones(size(faces, 1), 1)*color_data;
0538         face_color = 'flat';
0539     elseif (size(color_data, 1) == size(faces, 1))
0540         face_color = 'flat';
0541     elseif (size(color_data, 1) == size(vertices, 1))
0542         face_color = 'interp';        
0543     else
0544         eidors_msg('warning: color data and mesh do not match. Showing grey', 1);
0545         color_data = 0.5*ones(size(faces, 1), 3);
0546         face_color = 'flat';      
0547     end
0548     
0549     alpha_data_mapping = 'none';
0550     if (size(alpha_data, 1) == 1)
0551         alpha_data = ones(size(faces, 1), 1)*alpha_data;
0552         face_color = 'flat';
0553     elseif (size(alpha_data, 1) == size(faces, 1))
0554         face_alpha = 'flat';
0555     elseif (size(alpha_data, 1) == size(vertices, 1))
0556         face_alpha = 'interp';  
0557         alpha_data_mapping = 'scaled';  % why is this needed??
0558     else
0559         eidors_msg('warning: alpha data and mesh do not match. Showing opaque', 1);
0560         alpha_data = 1;
0561         face_alpha = 'flat';        
0562     end
0563 
0564     hh = patch('Faces', faces, 'Vertices', vertices, ...
0565                'AlphaDataMapping', alpha_data_mapping, ...
0566                'CDataMapping', 'direct', ...
0567                'FaceVertexCData', color_data, ...
0568                'FaceVertexAlphaData', alpha_data, ...
0569                'FaceColor', face_color, 'FaceAlpha', face_alpha, ...
0570                'EdgeColor', edge_color, 'FaceLighting', 'flat', ...
0571                'LineWidth', edge_width);
0572 
0573 function colour = electr_colour(elec, e, colormap_width)
0574     if isfield(elec, 'colour')
0575         desired_colour = elec.colour;
0576     else
0577         switch (e)
0578             case 1
0579                 desired_colour = [0 .7 0]; % light green electrode #1
0580             case 2
0581                 desired_colour = [0 .5 0]; % mid-green electrode #2
0582             otherwise
0583                 desired_colour = [0 .3 0]; % dark green
0584         end
0585     end
0586     switch(colormap_width)
0587         case 1
0588             map = colormap;
0589             colormap_idx = find(all(map == ones(size(map, 1), 1)* ...
0590                                                        desired_colour, 2));
0591             if ~isempty(colormap_idx)
0592                 colour = colormap_idx;
0593             else
0594                 map = [map; desired_colour];
0595                 colormap(map);
0596                 colour = size(map, 1);
0597             end
0598         case 3 
0599             colour = desired_colour;
0600         otherwise
0601             error('Invalid colormap width.');
0602     end
0603     
0604 function mdl = find_sub_elements(mdl)
0605     if (~isfield(mdl, 'sub_elements'))
0606         % Number of nodes per elements.
0607         n_nodes_per_elem = size(mdl.elems, 2);
0608         
0609         % Find sub-elements.
0610         combos= nchoosek(1:n_nodes_per_elem, n_nodes_per_elem - 1);
0611         mdl.sub_elements = sort( ...
0612                            reshape(mdl.elems(:, combos')', ...
0613                                    n_nodes_per_elem - 1, []), ...
0614                                  1)';
0615                        
0616         % Vector that associates each sub-element with
0617         % corresponding element.
0618         sub_element_idx = reshape(ones(n_nodes_per_elem, 1) ...
0619                                            *(1:size(mdl.elems, 1)),[] , 1);
0620                                        
0621         sub_element_other_node_idx = reshape((n_nodes_per_elem:-1:1)'...
0622                                      *ones(1, size(mdl.elems, 1)), [] , 1);
0623                                  
0624         sub_element_other_node = mdl.elems(sub2ind(size(mdl.elems), sub_element_idx, sub_element_other_node_idx));
0625                        
0626         % Remove duplicate sub-elements.
0627         % Previously specified "SPORTED" but that is default, and not
0628         %   available in older versions
0629         [mdl.sub_elements, ia, ic] = unique(...
0630                                        mdl.sub_elements, 'rows');
0631                                    
0632          sub_element_other_node = sub_element_other_node(ia, :);                    
0633                                 
0634         % Compute a matrix that converts a property defined on the elements
0635         % to a property defined on the sub-elements.
0636         mdl.element2sub_element = sparse(ic, sub_element_idx, ...
0637                              ones(n_nodes_per_elem*size(mdl.elems, 1), 1));
0638                          
0639         % Normalize each row to one to account for boundary sub-elements
0640         % and internal sub-elements.
0641         mdl.element2sub_element = spdiags(1./sum(...
0642                                      mdl.element2sub_element, 2), ...
0643                                     0, size(mdl.sub_elements, 1), ...
0644                       size(mdl.sub_elements, 1))*mdl.element2sub_element;
0645                   
0646         % Extract boundary from sub-elements.
0647         mdl.boundary = mdl.sub_elements;
0648         boundary_other_node = sub_element_other_node;
0649         mdl.boundary_to_sub_element_idx = (1:size(mdl.sub_elements, 1))';
0650 
0651         sorted_ic = sort(ic);
0652 
0653         mdl.boundary(sorted_ic(diff(sorted_ic) == 0), :) = [];
0654         boundary_other_node(sorted_ic(diff(sorted_ic) == 0), :) = [];
0655         mdl.boundary_to_sub_element_idx(sorted_ic(diff(sorted_ic) == 0)) = [];
0656 
0657         % in case we get vertex data we need the boundary nodes
0658         mdl.boundary_node_idx  = (unique(mdl.boundary(:)));
0659         
0660         
0661         if (size(mdl.boundary, 2) == 3)
0662             % Compute normal vectors.
0663             Node1 = mdl.nodes(mdl.boundary(:, 1), :);
0664             Node2 = mdl.nodes(mdl.boundary(:, 2), :);
0665             Node3 = mdl.nodes(mdl.boundary(:, 3), :);
0666             Node4 = mdl.nodes(boundary_other_node, :);
0667             mdl.boundary_normal_vector = cross(Node1 - Node2, ...
0668                                                Node3 - Node2, 2);
0669  
0670             % Normalize normal vectors.
0671             norm = 1./sqrt(sum(mdl.boundary_normal_vector ...
0672                                          .*mdl.boundary_normal_vector, 2));
0673             mdl.boundary_normal_vector = mdl.boundary_normal_vector ...
0674                                                         .*[norm norm norm];
0675               
0676             % Check if normal is outward. If not, invert direction.
0677             normal_vector_idx = dot(mdl.boundary_normal_vector, ...
0678                                                      Node4 - Node2, 2) > 0;
0679             mdl.boundary_normal_vector(normal_vector_idx, :) = ...
0680                          -mdl.boundary_normal_vector(normal_vector_idx, :);
0681                                                     
0682             % Find boundary edges.
0683             mdl.boundary_edges = sort(reshape(mdl.boundary(:, ...
0684                                             nchoosek(1:3, 2)')', 2, []), 1)';
0685 
0686             % Vector that associates each edge with its
0687             % corresponding two boundary elements.
0688             sub_boundary_edges_idx = reshape(ones(3, 1) ...
0689                                        *(1:size(mdl.boundary, 1)), [] , 1);
0690 
0691             [mdl.boundary_edges, sorted_idx] = sortrows(mdl.boundary_edges);
0692 
0693             mdl.boundary_edges = mdl.boundary_edges(1:2:end, :);
0694 
0695             mdl.edge2boundary = reshape(sub_boundary_edges_idx(...
0696                                                       sorted_idx), 2, [])';               
0697         end
0698 
0699         % Extract electrode boundary if necessary.
0700         if (isfield(mdl, 'electrode'))
0701             for i = 1:numel(mdl.electrode)
0702                 mdl.electrode(i).boundary = find(all(ismember(...
0703                                 mdl.boundary, mdl.electrode(i).nodes), 2));
0704                 mdl.electrode(i).sub_elements = ...
0705                                     mdl.boundary_to_sub_element_idx( ...
0706                                                 mdl.electrode(i).boundary);
0707                 
0708                 % Find electrode edges.
0709                 if (isfield(mdl, 'boundary_edges'))
0710                     edge_candidates = sum(ismember(mdl.edge2boundary, ...
0711                                             mdl.electrode(i).boundary), 2);
0712                     mdl.electrode(i).boundary_inner_edges = find(...
0713                                                      edge_candidates == 2);
0714                     mdl.electrode(i).boundary_outer_edges = find(...
0715                                                      edge_candidates == 1);                                           
0716                 end
0717             end
0718         end
0719     end
0720    
0721 function dest_struct = copy_field(dest_struct, src_struct, parent_field_name, child_field_name)
0722     if (isfield(src_struct, parent_field_name) && ...
0723         isfield(dest_struct, parent_field_name) && ...
0724         isfield(src_struct.(parent_field_name), child_field_name))
0725             dest_struct.(parent_field_name).(child_field_name) = ...
0726                 src_struct.(parent_field_name).(child_field_name);
0727     end
0728 
0729 % TESTS:
0730 function do_unit_test
0731    ver= eidors_obj('interpreter_version');
0732    clf
0733 
0734    img=calc_jacobian_bkgnd(mk_common_model('a2c0',8)); 
0735    img.elem_data= 3*sin(linspace(-2,2,num_elems(img))');
0736    subplot(3,4,1); show_fem_enhanced(img.fwd_model,[0,0,1]) 
0737    title('regular mesh numbered');
0738 
0739 if ~ver.isoctave 
0740    imgn = rmfield(img,'elem_data');
0741    imgn.node_data= 3*sin(linspace(-5,5,num_nodes(img))');
0742    subplot(3,4,9); show_fem_enhanced(imgn) 
0743    title('interpolated node colours');
0744 end
0745 
0746    img2(1) = img; img2(2) = img;
0747    subplot(3,4,2); show_fem_enhanced(img,[1]);
0748    title('colours with legend');
0749    subplot(3,4,3); show_fem_enhanced(img2,[0,1]);
0750    title('colours with legend');
0751    img.calc_colours.mapped_colour = 0; % USE RGB colours
0752    subplot(3,4,4); show_fem_enhanced(img,[0,1,1]);
0753    title('RGB colours');
0754    subplot(3,4,4); show_fem_enhanced(img);
0755    title('RGB colours');
0756 
0757    img.elem_data = [1:10];
0758    subplot(3,4,12);show_fem_enhanced(img); %Should show grey
0759    title('error -> show grey');
0760 
0761 if ~ver.isoctave
0762    imgn.calc_colours.mapped_colour = 0; % USE RGB colours
0763    subplot(3,4,10);show_fem_enhanced(imgn,[0,1]) 
0764    title('interpolated node colours');
0765 
0766 
0767    subplot(3,4,11);hh=show_fem_enhanced(imgn); set(hh(1),'EdgeColor',[0,0,1]);
0768    title('with edge colours');
0769 
0770 end
0771 
0772    img3=calc_jacobian_bkgnd(mk_common_model('n3r2',[16,2]));
0773    img3.elem_data= randn(828,1);                       
0774    subplot(3,4,5); show_fem_enhanced(img3.fwd_model) 
0775    subplot(3,4,6); show_fem_enhanced(img3,[1])
0776    subplot(3,4,7); show_fem_enhanced(img3,[1,1])
0777    subplot(3,4,8); show_fem_enhanced(img3,[1,1,1])