show_fem_enhanced

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

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