Home > eidors > graphics > matlab > show_fem_enhanced.m

show_fem_enhanced

PURPOSE ^

SHOW_FEM: show the EIDORS3D finite element model

SYNOPSIS ^

function hh = show_fem(mdl, options)

DESCRIPTION ^

 SHOW_FEM: show the EIDORS3D finite element model
 hh = show_fem(mdl, options)
 mdl is 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005