crop_model

 CROP_MODEL: Crop away parts of a fem model

 USAGE #1: crop display to show model internals
   crop_model( axis_handle, fcn_handle );

   fcn_handle ==1 where model is *kept*
 
   if axis_handle==[], then use the current axis
   examples:
     crop_model([],  inline('z==3','x','y','z'))
     crop_model(gca, inline('x+y>0','x','y','z'))
     crop_model([],  @(x,y,z) z<0);
   if the fcn_handle is a string, a function with x,y,z is created
     crop_model(gca, 'x+y>0') % same as previous

 USAGE #2: crop fwd_model to create new fwd_model
   fmdl_new= crop_model( fwd_model, fcn_handle );
 
   example:
   fmdl2= crop_model(fmdl1, @(x,y,z) x+y>0);

  with two parameters output
 [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
     c2f_idx maps each elemen in fmdl_new to fwd_model

0001 function [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
0002 % CROP_MODEL: Crop away parts of a fem model
0003 %
0004 % USAGE #1: crop display to show model internals
0005 %   crop_model( axis_handle, fcn_handle );
0006 %
0007 %   fcn_handle ==1 where model is *kept*
0008 %
0009 %   if axis_handle==[], then use the current axis
0010 %   examples:
0011 %     crop_model([],  inline('z==3','x','y','z'))
0012 %     crop_model(gca, inline('x+y>0','x','y','z'))
0013 %     crop_model([],  @(x,y,z) z<0);
0014 %   if the fcn_handle is a string, a function with x,y,z is created
0015 %     crop_model(gca, 'x+y>0') % same as previous
0016 %
0017 % USAGE #2: crop fwd_model to create new fwd_model
0018 %   fmdl_new= crop_model( fwd_model, fcn_handle );
0019 %
0020 %   example:
0021 %   fmdl2= crop_model(fmdl1, @(x,y,z) x+y>0);
0022 %
0023 %  with two parameters output
0024 % [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
0025 %     c2f_idx maps each elemen in fmdl_new to fwd_model
0026 
0027 % (C) 2006-2008 Andy Adler. License: GPL version 2 or version 3
0028 % $Id: crop_model.m 6030 2019-10-28 06:50:48Z aadler $
0029 
0030 if ischar(axis_handle) && strcmp(axis_handle,'UNIT_TEST'); do_unit_test; return; end
0031 
0032 % TODO (update 2 apr 2012):
0033 %  - make crop_model work for 2D fems
0034 
0035 if isstr(fcn_handle)
0036   fcn_handle = inline(fcn_handle,'x','y','z');
0037 end
0038 
0039 type= isfield(axis_handle, 'type');
0040 if type; type = axis_handle.type; end
0041 
0042 if isempty(axis_handle)
0043    axis_handle= gca;
0044    crop_graphics_model(axis_handle, fcn_handle);
0045 elseif ishandle( axis_handle )
0046    crop_graphics_model(axis_handle, fcn_handle);
0047 elseif strcmp(type, 'fwd_model');
0048    [fmdl,c2f_idx]= crop_fwd_model(axis_handle, fcn_handle);
0049 elseif strcmp(type, 'image');
0050    [fmdl_,c2f_idx]= crop_fwd_model(axis_handle.fwd_model, fcn_handle);
0051    fmdl = axis_handle; % input parameter
0052    fmdl.fwd_model = fmdl_;
0053    fmdl.elem_data = fmdl.elem_data(c2f_idx,:);
0054 %  keyboard
0055 else
0056    error('Not sure how to process first parameter');
0057 end
0058 
0059 % CROP GRAPHICS
0060 function crop_graphics_model(axis_handle, fcn_handle);
0061    kk= get(axis_handle,'Children');
0062    % only the FEM frame
0063    %k=kk( find( kk== min(kk) ));
0064 
0065    for k= kk(:)'
0066       try
0067          x= get(k,'XData');
0068          y= get(k,'YData');
0069          z= get(k,'ZData');
0070          c= get(k,'CData');
0071          idx= ~all( feval(fcn_handle,x,y,z) );
0072          if any(size(c)>[1,1])
0073             set(k,'Xdata', x(:,idx), ...
0074                   'Ydata', y(:,idx), ...
0075                   'Zdata', z(:,idx), ...
0076                   'Cdata', c(:,idx));
0077          else
0078             set(k,'Xdata', x(:,idx), ...
0079                   'Ydata', y(:,idx), ...
0080                   'Zdata', z(:,idx));
0081          end
0082       end
0083    end
0084 
0085 % CROP fwd_model
0086 function [fmdl1,c2f_idx]= crop_fwd_model(fmdl0, fcn_handle);
0087    fmdl1= fmdl0;
0088 
0089 % Find nodes to remove
0090    nodes= fmdl0.nodes;
0091    [n,d]= size(nodes);
0092    n2xyz= eye(d,3); 
0093    xyz= nodes*n2xyz;
0094 % Keep these nodes
0095    idx0=  all( feval(fcn_handle,xyz(:,1), ...
0096                                 xyz(:,2), ...
0097                                 xyz(:,3)),2);
0098 % remove these nodes
0099    fmdl1.nodes(idx0,:) = [];
0100 
0101 % renumber nodes, set unused ones to 0
0102    idx1= zeros(n,1);
0103    idx1(~idx0)= 1:sum(~idx0);
0104 
0105    fmdl1.elems(:) = idx1(fmdl0.elems);
0106    remove= any( fmdl1.elems == 0, 2);
0107    fmdl1.elems(remove,:)= [];
0108 % c2f_idx maps each new elem to its original position
0109    c2f_idx= find(~remove);
0110 
0111    fmdl1.boundary(:) = idx1(fmdl0.boundary);
0112    remove= any( fmdl1.boundary == 0, 2);
0113    fmdl1.boundary(remove,:)= [];
0114 
0115 % renumber nodes, set unused ones to 0
0116 if isfield(fmdl1,'electrode');
0117    n_elecs = length(fmdl1.electrode);
0118    rm_elec_list = zeros(n_elecs,1);
0119    for i=1:n_elecs;
0120       el_nodes= fmdl0.electrode(i).nodes;
0121       el_nodes(:)= idx1(el_nodes);
0122       if any(el_nodes==0);
0123          rm_elec_list(i) = 1;
0124       end
0125       fmdl1.electrode(i).nodes= el_nodes;
0126    end
0127    if any(rm_elec_list)
0128       str = sprintf('%d,', find(rm_elec_list));
0129       eidors_msg('crop_model: removing electrodes (%s)',str(1:end-1),1);
0130       fmdl1.electrode = fmdl1.electrode( find(~rm_elec_list));
0131    end
0132 end
0133 
0134 
0135 function do_unit_test
0136    imdl = mk_common_model('a2c0',8); fmdl= imdl.fwd_model;
0137    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0138    unit_test_cmp('crop_model-a2c0-01',length(fmdl.electrode),5);
0139    unit_test_cmp('crop_model-a2c0-02',size(fmdl.elems),[32,3]);
0140    unit_test_cmp('crop_model-a2c0-03',size(fmdl.nodes),[25,2]);
0141 
0142    fmdl = crop_model(fmdl,'x<0'); % verify it's same
0143    unit_test_cmp('crop_model-str-a2c0-01',length(fmdl.electrode),5);
0144    unit_test_cmp('crop_model-str-a2c0-02',size(fmdl.elems),[32,3]);
0145    unit_test_cmp('crop_model-str-a2c0-03',size(fmdl.nodes),[25,2]);
0146 
0147    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0148    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0149    unit_test_cmp('crop_model-n3r2-01',length(fmdl.electrode),16);
0150    unit_test_cmp('crop_model-n3r2-02',size(fmdl.elems),[342,4]);
0151    unit_test_cmp('crop_model-n3r2-03',size(fmdl.nodes),[128,3]);
0152    fmdl = crop_model(fmdl,inline('z<2','x','y','z'));
0153    unit_test_cmp('crop_model-n3r2-04',length(fmdl.electrode),8);
0154    unit_test_cmp('crop_model-n3r2-05',size(fmdl.elems),[114,4]);
0155    unit_test_cmp('crop_model-n3r2-06',size(fmdl.nodes),[64,3]);
0156 
0157 
0158 
0159 
0160    subplot(331)
0161    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0162    show_fem(fmdl);
0163    subplot(332)
0164    show_fem(fmdl); hh= gca; 
0165    subplot(333)
0166    show_fem(fmdl);
0167    crop_model([],inline('z<2','x','y','z'));
0168    crop_model(hh,inline('x<0','x','y','z'));
0169 
0170    subplot(334)
0171    fmdlc = fmdl;
0172    fmdlc = crop_model(fmdlc,inline('x<0','x','y','z'));
0173    show_fem(fmdlc);
0174 
0175    subplot(335)
0176    img = mk_image(fmdl,1); 
0177    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0178    imgs =  crop_model(img,'y-z<-2.5');
0179    show_fem( imgs );
0180    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0181 
0182    subplot(336)
0183    img = mk_image(fmdl,1); 
0184    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0185    imgs =  crop_model(img,@(x,y,z) y-z<-2.5);
0186    show_fem( imgs );
0187    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0188 
0189 
0190 
0191    subplot(337)
0192    imdl = mk_common_model('d2c2'); fmdl= imdl.fwd_model;
0193    show_fem(fmdl);
0194    subplot(338)
0195    show_fem(fmdl); hh= gca; 
0196    title('expected fail');
0197    subplot(339)
0198    show_fem(fmdl);
0199    crop_model([],inline('y<0','x','y','z'));
0200    crop_model(hh,inline('x<0','x','y','z'));
0201    title('expected fail');
0202    
0203