0001 function hh=show_fem( mdl, options)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 switch nargin
0031 case 0;
0032 error('Insufficient parameters for show_fem');
0033 case 1;
0034 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0035 options = [];
0036 end
0037 [img,mdl,opts] = proc_params( mdl, options );
0038
0039 if ~ishold
0040 cla;
0041 end
0042
0043 switch opts.dims
0044 case 2; hh= show_2d(img,mdl,opts);
0045 case 3; hh= show_3d(img,mdl,opts);
0046 otherwise; error('model is not 2D or 3D');
0047 end
0048
0049 if opts.show_elem_numbering
0050 xyzc= interp_mesh(mdl);
0051 placenumbers(xyzc, 7, [0,0,0],'none');
0052 end
0053 if opts.show_node_numbering
0054 xyzc= mdl.nodes;
0055 placenumbers(xyzc, 7, [0.0,0,0.5],[1,1,1]);
0056 end
0057
0058 if nargout == 0; clear hh; end
0059
0060
0061
0062
0063
0064
0065 function fix_axis
0066 ax = gca;
0067 axis tight
0068 set(ax,'DataAspectRatio',[1 1 1]);
0069
0070 [az el] = view;
0071 if el==90
0072
0073 xl=get(ax,'XLim'); yl=get(ax,'Ylim');
0074 dx = diff(xl); dy = diff(yl);
0075
0076 xb = strcmp(get(ax,'XAxisLocation'),'bottom');
0077 yn = strcmp(get(ax,'YDir'),'normal');
0078 if xb == yn
0079 set(ax,'YLim',yl + [0 1]*1e-3*dy);
0080 else
0081 set(ax,'YLim',yl + [-1 0]*1e-3*dy);
0082 end
0083 yl = strcmp(get(ax,'XAxisLocation'),'bottom');
0084 xn = strcmp(get(ax,'XDir'),'normal');
0085 if yl == xn
0086 set(ax,'XLim',xl + [0 1]*1e-3*dx);
0087 else
0088 set(ax,'XLim',xl + [-1 0]*1e-3*dx);
0089 end
0090 else
0091
0092 zl = get(ax,'Zlim'); dz = diff(zl);
0093 if strcmp(get(ax,'Zdir'),'normal')
0094 set(ax,'ZLim',zl + [0 1]*1e-3*dz);
0095 else
0096 set(ax,'ZLim',zl + [-1 0]*1e-3*dz);
0097 end
0098 end
0099
0100
0101 function placenumbers(xyzc, fontsize, colour, bgcolour)
0102 xyzc= xyzc * eye(size(xyzc,2),3);
0103 for i= 1:size(xyzc,1)
0104 text(xyzc(i,1),xyzc(i,2), xyzc(i,3), num2str(i), ...
0105 'HorizontalAlignment','center', ...
0106 'FontSize',fontsize,'Color',colour,'BackgroundColor',bgcolour);
0107 end
0108
0109 function [img,mdl,opts] = proc_params( mdl, options );
0110
0111 opts.do_colourbar=0;
0112 opts.number_electrodes=0;
0113 opts.show_numbering =0;
0114 opts.show_elem_numbering = 0;
0115 opts.show_node_numbering = 0;
0116 if nargin >=2
0117
0118 optionstr= zeros(1,100);
0119 optionstr(1:length(options)) = options;
0120
0121 opts.do_colourbar= optionstr(1);
0122 opts.number_electrodes= optionstr(2);
0123 switch optionstr(3)
0124 case 1; opts.show_elem_numbering = 1;
0125 case 2; opts.show_node_numbering = 1;
0126 case 3; opts.show_elem_numbering = 1;
0127 opts.show_node_numbering = 1;
0128 end
0129 end
0130
0131
0132
0133 if strcmp( mdl(1).type , 'image' )
0134 img= mdl;
0135 mdl= img(1).fwd_model;
0136 else
0137 img = [];
0138 end
0139
0140 opts.transparency_thresh = calc_colours('transparency_thresh');
0141 try
0142 opts.transparency_thresh = img.calc_colours.transparency_thresh;
0143 end
0144
0145 opts.dims = size(mdl.nodes,2);
0146
0147
0148 function hh= show_2d(img,mdl,opts)
0149 hax= gca;
0150 pax= get(hax,'position');
0151 if ~isempty(img)
0152 colours= calc_colours(img, [] );
0153 else
0154 colours= [1,1,1];
0155 end
0156 hh= show_2d_fem( mdl, colours );
0157 show_electrodes_2d(mdl, opts.number_electrodes);
0158
0159 set(hax,'position', pax);
0160 view(0, 90); axis('xy'); grid('off');
0161
0162
0163
0164 if exist('img','var') && opts.do_colourbar;
0165 colours= calc_colours(img, [], opts.do_colourbar);
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178 end
0179 if ~ishold
0180 fix_axis
0181 end
0182
0183
0184
0185
0186 function hh= show_3d(img,mdl,opts)
0187 hh= show_3d_fem( mdl );
0188
0189 if ~isempty(img)
0190 elem_data = get_img_data(img);
0191 show_inhomogeneities( elem_data , mdl, img, opts);
0192 if opts.do_colourbar
0193 calc_colours(img, [], opts.do_colourbar);
0194 end
0195 end
0196
0197
0198
0199 if ~ishold
0200 fix_axis
0201 end
0202 if size(mdl.elems,2) == 3
0203 show_electrodes_surf(mdl, opts.number_electrodes);
0204 else
0205 show_electrodes_3d(mdl, opts.number_electrodes);
0206 end
0207
0208
0209 function show_electrodes_2d(mdl, number_electrodes)
0210 if ~isfield(mdl,'electrode'); return; end
0211
0212 ee= get_boundary( mdl );
0213 ctr_x= mean(mdl.nodes(:,1));
0214 ctr_y= mean(mdl.nodes(:,2));
0215
0216
0217
0218 for e=1:length(mdl.electrode)
0219 if isfield(mdl.electrode(e),'nodes')
0220 elec_nodes= mdl.electrode(e).nodes;
0221
0222 S= 1.00;
0223 vx= (mdl.nodes(elec_nodes,1) - ctr_x)*S;
0224 vy= (mdl.nodes(elec_nodes,2) - ctr_y)*S;
0225
0226 [jnk,idx] = order_loop( [vx,vy] );
0227 elseif isfield(mdl.electrode(e),'pos')
0228 vx = mdl.electrode(e).pos(:,1) - ctr_x;
0229 vy = mdl.electrode(e).pos(:,2) - ctr_y;
0230 idx = 1:length(vx);
0231 else
0232 eidors_msg('show_fem: WARNING: electrode %d has no nodes. Not showing.',e,2);
0233 continue;
0234 end
0235
0236 ecolour = electr_colour( e );
0237 if numel(vx) == 1
0238
0239 line(vx(idx)+ctr_x,vy(idx)+ctr_y, ...
0240 'LineWidth', 2, 'LineStyle','-','Color', ecolour, ...
0241 'Marker','o','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0242 else
0243
0244
0245 line(vx(idx)+ctr_x,vy(idx)+ctr_y, ...
0246 'LineWidth', 3, 'LineStyle','-','Color', ecolour, ...
0247 'Marker','none','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0248 end
0249 if number_electrodes
0250 S= 1.05;
0251 vx= vx*S;
0252 vy= vy*S;
0253 switch number_electrodes
0254 case {1 true}
0255 txt = num2str(e);
0256 case 2
0257 try, txt = mdl.electrode(e).label; end
0258 end
0259 hh= text(mean(vx)+ctr_x, mean(vy)+ctr_y, txt);
0260 set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0261 end
0262 end
0263
0264 function show_electrodes_surf(mdl, number_electrodes)
0265 if ~isfield(mdl,'electrode'); return; end
0266
0267 ee= get_boundary( mdl );
0268 ctr_x= mean(mdl.nodes(:,1));
0269 ctr_y= mean(mdl.nodes(:,2));
0270 ctr_z= mean(mdl.nodes(:,3));
0271
0272
0273 for e=1:length(mdl.electrode)
0274 if isfield(mdl.electrode(e),'pos') && ~isfield(mdl.electrode(e),'nodes')
0275 vx = mdl.electrode(e).pos(:,1) - ctr_x;
0276 vy = mdl.electrode(e).pos(:,2) - ctr_y;
0277 vz = mdl.electrode(e).pos(:,3) - ctr_z;
0278 idx = 1:length(vx);
0279 S = 1.00;
0280 else
0281 elec_nodes= mdl.electrode(e).nodes;
0282
0283 S= 1.00;
0284 vx= (mdl.nodes(elec_nodes,1) - ctr_x);
0285 vy= (mdl.nodes(elec_nodes,2) - ctr_y);
0286 vz= (mdl.nodes(elec_nodes,3) - ctr_z);
0287 end
0288
0289
0290
0291 [jnk,idx] = sort(unwrap(atan2( vy, vx )));
0292
0293 ecolour = electr_colour( e );
0294 if numel(vx) == 1
0295
0296 line(S*vx(idx)+ctr_x,S*vy(idx)+ctr_y, S*vz(idx)+ctr_z, ...
0297 'LineWidth', 2, 'LineStyle','-','Color', ecolour, ...
0298 'Marker','o','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0299 else
0300
0301
0302 line(vx(idx)+ctr_x,vy(idx)+ctr_y, vz(idx)+ctr_z, ...
0303 'LineWidth', 3, 'LineStyle','-','Color', ecolour, ...
0304 'Marker','none','MarkerSize', 6,'MarkerEdgeColor',ecolour);
0305 end
0306 if number_electrodes
0307 S= 1.05;
0308
0309
0310
0311 switch number_electrodes
0312 case {1 true}
0313 txt = num2str(e);
0314 case 2
0315 try, txt = mdl.electrode(e).label; end
0316 end
0317 hh= text(S*mean(vx)+ctr_x, S*mean(vy)+ctr_y,S*mean(vz)+ctr_z,txt);
0318 set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0319 end
0320 end
0321
0322 function show_electrodes_3d(mdl, number_electrodes)
0323
0324 if ~isfield(mdl,'electrode'); return; end
0325
0326 ee= get_boundary( mdl );
0327 for e=1:length(mdl.electrode)
0328 colour= electr_colour( e);
0329
0330 if isfield(mdl.electrode(e),'pos') && ~isfield(mdl.electrode(e),'nodes')
0331 show_electrodes_surf(mdl, number_electrodes);
0332 return
0333 end
0334 elec_nodes= mdl.electrode(e).nodes;
0335
0336
0337 if length(elec_nodes) == 1
0338 vtx= mdl.nodes(elec_nodes,:);
0339 line(vtx(1),vtx(2),vtx(3), ...
0340 'Marker','o','MarkerSize',12, ...
0341 'MarkerFaceColor',colour, 'MarkerEdgeColor', colour);
0342 if number_electrodes
0343 hh= text(vtx(1),vtx(2),vtx(3), num2str(e));
0344 set(hh, 'HorizontalAlignment','center', 'FontWeight','bold');
0345 end
0346 else
0347
0348 map = zeros(length(mdl.nodes),1);
0349 map(mdl.electrode(e).nodes) = 1;
0350 ec = map(ee);
0351 sels= find(all(ec'));
0352
0353 ee= get_boundary( mdl );
0354 paint_electrodes(sels,ee,mdl.nodes,colour,number_electrodes);
0355
0356 if number_electrodes
0357 el_nodes= mdl.nodes(unique(mdl.boundary(sels,:)),:);
0358 switch number_electrodes
0359 case {1 true}
0360 txt = num2str(e);
0361 case 2
0362 try, txt = mdl.electrode(e).label; end
0363 end
0364 hh=text( mean(el_nodes(:,1)), ...
0365 mean(el_nodes(:,2)), ...
0366 mean(el_nodes(:,3)), txt );
0367 set(hh,'FontWeight','bold','Color',[.8,.2,0],'FontSize',8);
0368 end
0369 end
0370 end
0371
0372 function show_inhomogeneities( elem_data, mdl, img, opt)
0373
0374 if size(elem_data,2)>1
0375 q = warning('query','backtrace');
0376 warning('backtrace','off');
0377 warning('EIDORS:FirstImageOnly','show_fem only shows first image');
0378 warning('backtrace',q.state);
0379
0380 end
0381 repaint_inho(elem_data(:,1), 'use_global' , mdl.nodes, mdl.elems, ...
0382 opt.transparency_thresh, img);
0383 if ~exist('OCTAVE_VERSION');
0384 camlight('left');
0385 lighting('none');
0386 end
0387
0388 function paint_electrodes_old(sel,srf,vtx, colour, show_num);
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398 l = srf(sel,1); m = srf(sel,2); n = srf(sel,3);
0399
0400 Xs = [vtx(l,1);vtx(m,1);vtx(n,1)];
0401 Ys = [vtx(l,2);vtx(m,2);vtx(n,2)];
0402 Zs = [vtx(l,3);vtx(m,3);vtx(n,3)];
0403
0404 h=patch(Xs,Ys,Zs, colour);
0405
0406 set(h, 'FaceLighting','none', 'CDataMapping', 'direct' );
0407
0408 function paint_electrodes(sel,srf,vtx, colour, show_num);
0409 if isempty(sel); return; end
0410
0411 [u n m] = unique(srf(sel,:));
0412 fv.vertices = vtx(u,:);
0413 fv.faces = reshape(m,[],3);
0414 h = patch(fv,'FaceColor',colour);
0415
0416 set(h, 'FaceLighting','none', 'CDataMapping', 'direct' );
0417
0418 function hh= show_3d_fem( mdl, options )
0419 if mdl_dim(mdl) == 3
0420 ee= get_boundary( mdl );
0421 elseif mdl_dim(mdl) == 2
0422 ee= mdl.elems;
0423 elseif mdl_dim(mdl) == 1
0424 ee= mdl.elems;
0425 else
0426 error('Simplices are not 2D or 3D in mdl. Unsure how to display');
0427 end
0428 hh= trimesh(ee, mdl.nodes(:,1), ...
0429 mdl.nodes(:,2), ...
0430 mdl.nodes(:,3));
0431 set(hh, 'EdgeColor', [0,0,0]);
0432 axis('image');
0433 hidden('off');
0434
0435 function hh=show_2d_fem( mdl, colours, imgno )
0436 szcolr = size(colours);
0437 if nargin<3;
0438 imgno = 1;
0439 if szcolr(1:2)>1;
0440 eidors_msg('warning: show_fem only shows first image',1);
0441 end
0442 end
0443
0444 if szcolr(1:2) == [1,3]
0445 hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor',colours);
0446 elseif size(colours,1) == num_elems(mdl);
0447 colours = squeeze(colours(:,imgno,:));
0448
0449 hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor','flat', ...
0450 'facevertexcdata',colours,'CDataMapping','direct');
0451 elseif size(colours,1) == num_nodes(mdl);
0452 colours = squeeze(colours(:,imgno,:));
0453
0454 hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor','interp', ...
0455 'facevertexcdata',colours,'CDataMapping','direct');
0456 else
0457 eidors_msg('warning: image elements and mesh do not match. Showing grey',1);
0458 colours= [1,1,1]/2;
0459 hh= patch('Faces',mdl.elems,'Vertices',mdl.nodes, 'facecolor',colours);
0460 end
0461
0462 set(hh, 'FaceLighting','none', 'CDataMapping', 'direct' );
0463
0464 function hh=show_2d_fem_oldmatlab( mdl, colours, imgno )
0465 szcolr = size(colours);
0466 if nargin<3;
0467 imgno = 1;
0468 if szcolr(1:2)>1; eidors_msg('warning: show_fem only shows first image',1); end
0469 end
0470
0471
0472
0473
0474 S= 1;
0475 elem= mdl.elems'; e= size(elem,2);
0476 node= mdl.nodes'; n= size(node,2);
0477 Xs=zeros(3,e);
0478 Xs(:)=mdl.nodes(elem(:),1);
0479 Xs= S*Xs+ (1-S)*ones(3,1)*mean(Xs);
0480 Ys=zeros(3,e); Ys(:)=mdl.nodes(elem(:),2);
0481 Ys= S*Ys+ (1-S)*ones(3,1)*mean(Ys);
0482 Zs = zeros(size(Xs));
0483 if szcolr(1:2) == [1,3]
0484 hh= patch(Xs,Ys,zeros(3,e),colours);
0485 elseif size(colours,1) == e
0486
0487 colours = permute(colours(:,imgno,:),[2,1,3]);
0488 elseif size(colours,1) == n
0489 colours = colours(elem(:), imgno, :);
0490 colours = reshape( colours, size(elem,1), size(elem,2), []);
0491 else
0492 eidors_msg('warning: image elements and mesh do not match. Showing grey',1);
0493 colours= [1,1,1]/2;
0494 end
0495
0496 hh= patch(Xs,Ys,Zs,colours);
0497 set(hh, 'FaceLighting','none', 'CDataMapping', 'direct' );
0498
0499
0500
0501 max_x= max(mdl.nodes(:,1)); min_x= min(mdl.nodes(:,1));
0502 max_y= max(mdl.nodes(:,2)); min_y= min(mdl.nodes(:,2));
0503 axis([ mean([max_x,min_x]) + 0.55*[-1,1]*(max_x-min_x), ...
0504 mean([max_y,min_y]) + 0.55*[-1,1]*(max_y-min_y) ]);
0505
0506
0507
0508 function old_code_3d_plot
0509 domesnum= 0;
0510 donodenum= 0;
0511 dotext=0;
0512
0513 xxx=zeros(4,e); xxx(:)=NODE(1,ELEM(:));
0514 xxx= S*xxx+ (1-S)*ones(4,1)*mean(xxx);
0515 yyy=zeros(4,e); yyy(:)=NODE(2,ELEM(:));
0516 yyy= S*yyy+ (1-S)*ones(4,1)*mean(yyy);
0517 zzz=zeros(4,e); zzz(:)=NODE(3,ELEM(:));
0518 zzz= S*zzz+ (1-S)*ones(4,1)*mean(zzz);
0519 plot3( xxx([1 2 3 1 4 2 3 4],:), ...
0520 yyy([1 2 3 1 4 2 3 4],:), ...
0521 zzz([1 2 3 1 4 2 3 4],:),'k' );
0522 hold on;
0523 xy=NODE(1:2,MES(1,:));
0524 plot3(arrow*xy,arrow*[0 1;-1 0]*xy,ones(8,1)*NODE(3,MES(1,:)),'b')
0525 hold off;
0526 axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ...
0527 [-1.1 1.1]*max(NODE(3,:)) ])
0528
0529
0530 function plot_2d_mesh(NODE,ELEM,el_pos, S, options)
0531
0532
0533
0534
0535
0536 e=size(ELEM,2);
0537 d=size(ELEM,1);
0538 R= zeros(1,e);
0539
0540
0541 arrow= [1.02 0;1.06 .05;1.06 .02;1.1 .02; ...
0542 1.1 -.02; 1.06 -.02;1.06 -.05;1.02 0];
0543
0544
0545 xxx=zeros(3,e); xxx(:)=NODE(1,ELEM(:));
0546 xxx= S*xxx+ (1-S)*ones(3,1)*mean(xxx);
0547 xxx= [xxx;xxx(1,:)];
0548
0549 yyy=zeros(3,e); yyy(:)=NODE(2,ELEM(:));
0550 yyy= S*yyy+ (1-S)*ones(3,1)*mean(yyy);
0551 yyy= [yyy;yyy(1,:)];
0552
0553 if isempty(el_pos)
0554 plot(xxx,yyy,'b');
0555 else
0556 xy= el_pos;
0557 hh=plot([xxx;xxx(1,:)],[yyy;yyy(1,:)],'b', ...
0558 arrow*xy,arrow*[0 1;-1 0]*xy,'r');
0559 set(hh(find(R>0)),'Color',[1 0 0],'LineWidth',2);
0560 set(hh(find(R<0)),'Color',[0 0 1],'LineWidth',2);
0561 end
0562
0563 if options(1)
0564 mesnum= reshape(sprintf(' %-2d',[0:15]),3,16)';
0565 text(NODE(1,MES(1,:))*1.08,NODE(2,MES(1,:))*1.08,mesnum, ...
0566 'Color','green','HorizontalAlignment','center');
0567 end
0568
0569 if options(2)
0570 nodenum= reshape(sprintf('%3d',1:length(NODE)),3,length(NODE))';
0571 text(NODE(1,:),NODE(2,:),nodenum, ...
0572 'Color','yellow','HorizontalAlignment','center','FontSize',14);
0573 end
0574
0575 if options(3)
0576 numeros= reshape(sprintf('%4d',[1:e]),4,e)';
0577
0578 decal= .02;
0579 xcoor=mean(NODE(2*ELEM-1))';
0580 ycoor=mean(NODE(2*ELEM))';
0581 text(xcoor-decal,ycoor,numeros,'FontSize',8, ...
0582 'HorizontalAlignment','center');
0583
0584 end
0585 axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ])
0586
0587 function mes= avg_electrode_posn( mdl )
0588 if ~isfield(mdl,'electrode'); mes=[]; return; end
0589 n_elec= length( mdl.electrode );
0590 nodes = mdl.nodes;
0591 mes= zeros( n_elec, mdl_dim(mdl) )
0592 for i= 1:n_elec
0593 e_nodes= mdl.electrode(i).nodes;
0594 mes(i,:)= mean( nodes(e_nodes,:) , 1 );
0595 end
0596
0597 function colour= electr_colour( e);
0598 if e==1;
0599 colour = [0,.7,0];
0600 elseif e==2
0601 colour = [0,.5,0];
0602 else
0603 colour = [0,.3,0];
0604 end
0605
0606 function ee= get_boundary( mdl )
0607 if isfield(mdl,'boundary')
0608 ee= mdl.boundary;
0609 else
0610
0611 ee = find_boundary( mdl.elems );
0612 end
0613
0614
0615
0616 function do_unit_test
0617 ver= eidors_obj('interpreter_version');
0618 clf
0619
0620 img=calc_jacobian_bkgnd(mk_common_model('a2c0',8));
0621 img.elem_data=rand(size(img.fwd_model.elems,1),1);
0622 subplot(3,4,1); show_fem(img.fwd_model,[0,0,1])
0623 title('regular mesh numbered');
0624
0625 if ~ver.isoctave
0626 imgn = rmfield(img,'elem_data');
0627 imgn.node_data=rand(size(img.fwd_model.nodes,1),1);
0628 subplot(3,4,9); show_fem(imgn)
0629 title('interpolated node colours');
0630 end
0631
0632 img2(1) = img; img2(2) = img;
0633 subplot(3,4,2); show_fem(img,[1]);
0634 title('colours with legend');
0635 subplot(3,4,3); show_fem(img2,[0,1]);
0636 title('colours with legend');
0637 img.calc_colours.mapped_colour = 0;
0638 subplot(3,4,4); show_fem(img,[0,1,1]);
0639 title('RGB colours');
0640 subplot(3,4,4); show_fem(img);
0641 title('RGB colours');
0642
0643 img.elem_data = [1:10];
0644 subplot(3,4,12);show_fem(img);
0645 title('error -> show grey');
0646
0647 if ~ver.isoctave
0648 imgn.calc_colours.mapped_colour = 0;
0649 subplot(3,4,10);show_fem(imgn,[0,1])
0650 title('interpolated node colours');
0651
0652
0653 subplot(3,4,11);hh=show_fem(imgn); set(hh,'EdgeColor',[0,0,1]);
0654 title('with edge colours');
0655 end
0656
0657 img3=calc_jacobian_bkgnd(mk_common_model('n3r2',[16,2]));
0658 img3.elem_data= randn(828,1);
0659 subplot(3,4,5); show_fem(img3.fwd_model)
0660 subplot(3,4,6); show_fem(img3,[1])
0661 subplot(3,4,7); show_fem(img3,[1,1])
0662 subplot(3,4,8); show_fem(img3,[1,1,1])
0663