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

 USAGE #3: crop img to create new img (preserve elem_data)
   img2= crop_model(img1, @(x,y,z) x+y>0);

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 % USAGE #3: crop img to create new img (preserve elem_data)
0028 %   img2= crop_model(img1, @(x,y,z) x+y>0);
0029 
0030 % (C) 2006-2008 Andy Adler. License: GPL version 2 or version 3
0031 % $Id: crop_model.m 6411 2022-11-25 18:50:48Z aadler $
0032 
0033 if ischar(axis_handle) && strcmp(axis_handle,'UNIT_TEST'); do_unit_test; return; end
0034 
0035 % TODO (update 2 apr 2012):
0036 %  - make crop_model work for 2D fems
0037 
0038 if isstr(fcn_handle)
0039   fcn_handle = inline(fcn_handle,'x','y','z');
0040 end
0041 
0042 type= isfield(axis_handle, 'type');
0043 if type; type = axis_handle.type; end
0044 
0045 if isempty(axis_handle)
0046    axis_handle= gca;
0047    crop_graphics_model(axis_handle, fcn_handle);
0048 elseif ishandle( axis_handle )
0049    crop_graphics_model(axis_handle, fcn_handle);
0050 elseif strcmp(type, 'fwd_model');
0051    [fmdl,c2f_idx]= crop_fwd_model(axis_handle, fcn_handle);
0052 elseif strcmp(type, 'image');
0053    [fmdl_,c2f_idx]= crop_fwd_model(axis_handle.fwd_model, fcn_handle);
0054    fmdl = axis_handle; % input parameter
0055    fmdl.fwd_model = fmdl_;
0056    fmdl.elem_data = fmdl.elem_data(c2f_idx,:);
0057 %  keyboard
0058 else
0059    error('Not sure how to process first parameter');
0060 end
0061 
0062 % CROP GRAPHICS
0063 function crop_graphics_model(axis_handle, fcn_handle);
0064    kk= get(axis_handle,'Children');
0065    % only the FEM frame
0066    %k=kk( find( kk== min(kk) ));
0067 
0068    for k= kk(:)'
0069       ktype = get(k,'Type');
0070       switch ktype; case {'patch','line'}
0071           crop_patch(fcn_handle, k)
0072       end % ignore non-patch
0073    end
0074 
0075 function crop_patch(fcn_handle, k)
0076    x= get(k,'XData');
0077    y= get(k,'YData');
0078    try
0079       z= get(k,'ZData');
0080    catch
0081       z= [];
0082    end
0083    idx= ~all( feval(fcn_handle,x,y,z) );
0084 
0085    set_= {'Xdata', x(:,idx), 'Ydata', y(:,idx)};
0086    if ~isempty(z)
0087       set_{end+1}= 'Zdata';
0088       % Can't assign to empty array, so indirect
0089       dd = z(:,idx); 
0090       set_{end+1}= dd;
0091    end
0092 
0093    try
0094       c= get(k,'CData');
0095 % Matlab doesn't document the shape of Cdata!
0096       if size(c,2)==1; c=c'; end
0097       if ~isempty(c);
0098           c = c(:,idx);
0099       end
0100 %  This should work, but weird errors
0101 %     set_= horzcat(set_,{'Cdata',c(:,idx)});
0102       set_{end+1}= 'Cdata';
0103       % Can't assign to empty array, so indirect
0104       set_{end+1}= c;
0105    end
0106    set(k, set_{:});
0107 
0108 % CROP fwd_model
0109 function [fmdl1,c2f_idx]= crop_fwd_model(fmdl0, fcn_handle);
0110    fmdl1= fmdl0;
0111 
0112 % Find nodes to remove
0113    nodes= fmdl0.nodes;
0114    [n,d]= size(nodes);
0115    n2xyz= eye(d,3); 
0116    xyz= nodes*n2xyz;
0117 % Keep these nodes
0118    idx0=  all( feval(fcn_handle,xyz(:,1), ...
0119                                 xyz(:,2), ...
0120                                 xyz(:,3)),2);
0121 % remove these nodes
0122    fmdl1.nodes(idx0,:) = [];
0123 
0124 % renumber nodes, set unused ones to 0
0125    idx1= zeros(n,1);
0126    idx1(~idx0)= 1:sum(~idx0);
0127 
0128    fmdl1.elems(:) = idx1(fmdl0.elems);
0129    remove= any( fmdl1.elems == 0, 2);
0130    fmdl1.elems(remove,:)= [];
0131 % c2f_idx maps each new elem to its original position
0132    c2f_idx= find(~remove);
0133 
0134    fmdl1.boundary(:) = idx1(fmdl0.boundary);
0135    remove= any( fmdl1.boundary == 0, 2);
0136    fmdl1.boundary(remove,:)= [];
0137 
0138 % renumber nodes, set unused ones to 0
0139 if isfield(fmdl1,'electrode');
0140    n_elecs = length(fmdl1.electrode);
0141    rm_elec_list = zeros(n_elecs,1);
0142    for i=1:n_elecs;
0143       el_nodes= fmdl0.electrode(i).nodes;
0144       el_nodes(:)= idx1(el_nodes);
0145       if any(el_nodes==0);
0146          rm_elec_list(i) = 1;
0147       end
0148       fmdl1.electrode(i).nodes= el_nodes;
0149    end
0150    if any(rm_elec_list)
0151       str = sprintf('%d,', find(rm_elec_list));
0152       eidors_msg('crop_model: removing electrodes (%s)',str(1:end-1),1);
0153       fmdl1.electrode = fmdl1.electrode( find(~rm_elec_list));
0154    end
0155 end
0156 
0157 
0158 function do_unit_test
0159    imdl = mk_common_model('a2c0',8); fmdl= imdl.fwd_model;
0160    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0161    unit_test_cmp('crop_model-a2c0-01',length(fmdl.electrode),5);
0162    unit_test_cmp('crop_model-a2c0-02',size(fmdl.elems),[32,3]);
0163    unit_test_cmp('crop_model-a2c0-03',size(fmdl.nodes),[25,2]);
0164 
0165    fmdl = crop_model(fmdl,'x<0'); % verify it's same
0166    unit_test_cmp('crop_model-str-a2c0-01',length(fmdl.electrode),5);
0167    unit_test_cmp('crop_model-str-a2c0-02',size(fmdl.elems),[32,3]);
0168    unit_test_cmp('crop_model-str-a2c0-03',size(fmdl.nodes),[25,2]);
0169 
0170    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0171    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0172    unit_test_cmp('crop_model-n3r2-01',length(fmdl.electrode),16);
0173    unit_test_cmp('crop_model-n3r2-02',size(fmdl.elems),[342,4]);
0174    unit_test_cmp('crop_model-n3r2-03',size(fmdl.nodes),[128,3]);
0175    fmdl = crop_model(fmdl,inline('z<2','x','y','z'));
0176    unit_test_cmp('crop_model-n3r2-04',length(fmdl.electrode),8);
0177    unit_test_cmp('crop_model-n3r2-05',size(fmdl.elems),[114,4]);
0178    unit_test_cmp('crop_model-n3r2-06',size(fmdl.nodes),[64,3]);
0179 
0180 
0181 
0182 
0183    subplot(331)
0184    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0185    show_fem(fmdl);
0186    subplot(332)
0187    show_fem(fmdl); hh= gca; 
0188    subplot(333)
0189    show_fem(fmdl);
0190    crop_model([],inline('z<2','x','y','z'));
0191    crop_model(hh,inline('x<0','x','y','z'));
0192 
0193    subplot(334)
0194    fmdlc = fmdl;
0195    fmdlc = crop_model(fmdlc,inline('x<0','x','y','z'));
0196    show_fem(fmdlc);
0197 
0198    subplot(335)
0199    img = mk_image(fmdl,1); 
0200    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0201    imgs =  crop_model(img,'y-z<-2.5');
0202    show_fem( imgs );
0203    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0204 
0205    subplot(336)
0206    img = mk_image(fmdl,1); 
0207    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0208    imgs =  crop_model(img,@(x,y,z) y-z<-2.5);
0209    show_fem( imgs );
0210    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0211 
0212 
0213 
0214    subplot(337)
0215    imdl = mk_common_model('d2c2');
0216    fmdl= imdl.fwd_model;
0217    img = mk_image(imdl); img.elem_data(1:16) = 1.1;
0218    show_fem(fmdl);
0219    subplot(338)
0220    show_fem(img);
0221    try
0222    crop_model([],@(x,y,z) y<0)
0223    catch
0224    title('expected fail');
0225    end
0226 
0227    subplot(339)
0228    show_fem(fmdl); hh= gca; 
0229    crop_model(hh,inline('x<0','x','y','z'));
0230    
0231