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 6811 2024-04-24 20:20:11Z bgrychtol $
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, c2f_idx] = remove_nodes(fmdl0, idx0, 'keep');
0123 
0124 
0125 function do_unit_test
0126 
0127    imdl = mk_common_model('a2c0',8); fmdl= imdl.fwd_model;
0128    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0129    unit_test_cmp('crop_model-a2c0-01',length(fmdl.electrode),8);
0130    unit_test_cmp('crop_model-a2c0-02',size(fmdl.elems),[32,3]);
0131    unit_test_cmp('crop_model-a2c0-03',size(fmdl.nodes),[25,2]);
0132 
0133    fmdl = crop_model(fmdl,'x<0'); % verify it's same
0134    unit_test_cmp('crop_model-str-a2c0-01',length(fmdl.electrode),8);
0135    unit_test_cmp('crop_model-str-a2c0-02',size(fmdl.elems),[32,3]);
0136    unit_test_cmp('crop_model-str-a2c0-03',size(fmdl.nodes),[25,2]);
0137 
0138    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0139    fmdl = crop_model(fmdl,inline('x<0','x','y','z'));
0140    unit_test_cmp('crop_model-n3r2-01',length(fmdl.electrode),32);
0141    unit_test_cmp('crop_model-n3r2-02',size(fmdl.elems),[342,4]);
0142    unit_test_cmp('crop_model-n3r2-03',size(fmdl.nodes),[128,3]);
0143    fmdl = crop_model(fmdl,inline('z<2','x','y','z'));
0144    unit_test_cmp('crop_model-n3r2-04',length(fmdl.electrode),32);
0145    unit_test_cmp('crop_model-n3r2-05',size(fmdl.elems),[114,4]);
0146    unit_test_cmp('crop_model-n3r2-06',size(fmdl.nodes),[64,3]);
0147 
0148 
0149 
0150 
0151    subplot(331)
0152    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0153    show_fem(fmdl);
0154    subplot(332)
0155    show_fem(fmdl); hh= gca; 
0156    subplot(333)
0157    show_fem(fmdl);
0158    crop_model([],inline('z<2','x','y','z'));
0159    crop_model(hh,inline('x<0','x','y','z'));
0160 
0161    subplot(334)
0162    fmdlc = fmdl;
0163    fmdlc = crop_model(fmdlc,inline('x<0','x','y','z'));
0164    show_fem(fmdlc);
0165 
0166    subplot(335)
0167    img = mk_image(fmdl,1); 
0168    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0169    imgs =  crop_model(img,'y-z<-2.5'); % warning about # 11
0170    show_fem( imgs );
0171    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0172 
0173    subplot(336)
0174    img = mk_image(fmdl,1); 
0175    load('datareal.mat','A'); img.elem_data(A)= 1.1;
0176    imgs =  crop_model(img,@(x,y,z) y-z<-2.5);
0177    show_fem( imgs );
0178    unit_test_cmp('crop img',find( imgs.elem_data>1),[476;479; 482; 485])
0179 
0180 
0181 
0182    subplot(337)
0183    imdl = mk_common_model('d2c2');
0184    fmdl= imdl.fwd_model;
0185    img = mk_image(imdl); img.elem_data(1:16) = 1.1;
0186    show_fem(fmdl);
0187    subplot(338)
0188    show_fem(img);
0189    try
0190    crop_model([],@(x,y,z) y<0)
0191    catch
0192    title('expected fail');
0193    end
0194 
0195    subplot(339)
0196    show_fem(fmdl); hh= gca; 
0197    crop_model(hh,inline('x<0','x','y','z'));
0198 
0199