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 may also be specified as
  fwd_model.show_fem_enhanced.field = value

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