solve_RM_2Dslice

 SOLVE_RM_2DSLICE: cut slices out of a 3D model
  which has reconstruction matrix calculated

 Basic usage
 For example, given a 3D GREIT model
    vopt.zvec= 0.5:0.2:2.5;
    vopt.imgsz = [32 32];
    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
    imdl = mk_GREIT_model(imdl, 0.20, [], opt);

 One could cut the z=1 slice as
    imdl2= solve_RM_2Dslice(imdl, 1.0);
 Or using 
    sel_fcn = inline('abs(z-1.0)<0.2'); OR
    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ; % two planes
    imdl2= solve_RM_2Dslice(imdl, sel_fcn);
 
 Note that the reconstruction planes are between the
 cut planes specified in vopt.zvec

 options:
  imdl.solve_RM_2Dslice.vert_tol = .001;
       % Vertical tolerance for elems on plane (Default .001)
  imdl.solve_RM_2Dslice.keep3D = 1; %keep 3D aspect of cuts

0001 function imdl = solve_RM_2Dslice(imdl, sel_fcn)
0002 % SOLVE_RM_2DSLICE: cut slices out of a 3D model
0003 %  which has reconstruction matrix calculated
0004 %
0005 % Basic usage
0006 % For example, given a 3D GREIT model
0007 %    vopt.zvec= 0.5:0.2:2.5;
0008 %    vopt.imgsz = [32 32];
0009 %    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0010 %    imdl = mk_GREIT_model(imdl, 0.20, [], opt);
0011 %
0012 % One could cut the z=1 slice as
0013 %    imdl2= solve_RM_2Dslice(imdl, 1.0);
0014 % Or using
0015 %    sel_fcn = inline('abs(z-1.0)<0.2'); OR
0016 %    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ; % two planes
0017 %    imdl2= solve_RM_2Dslice(imdl, sel_fcn);
0018 %
0019 % Note that the reconstruction planes are between the
0020 % cut planes specified in vopt.zvec
0021 %
0022 % options:
0023 %  imdl.solve_RM_2Dslice.vert_tol = .001;
0024 %       % Vertical tolerance for elems on plane (Default .001)
0025 %  imdl.solve_RM_2Dslice.keep3D = 1; %keep 3D aspect of cuts
0026 
0027 % (C) 2015-2018 Andy Adler and Bartlomiej Grychtol
0028 % License: GPL version 2 or version 3
0029 % $Id: solve_RM_2Dslice.m 6038 2019-12-30 22:33:47Z aadler $
0030 
0031 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'); do_unit_test; return; end
0032 
0033 % Options
0034 vert_tol = .001; try
0035   vert_tol = imdl.solve_RM_2Dslice.vert_tol;
0036 end
0037 
0038 keep3D=0; try
0039   keep3D = imdl.solve_RM_2Dslice.keep3D;
0040 end
0041 
0042 ctr = interp_mesh(imdl.rec_model); ct3= ctr(:,3);
0043 
0044 if isnumeric( sel_fcn )
0045    ff = abs(ct3-sel_fcn) < vert_tol;
0046 else
0047     if ischar(sel_fcn) %  we have a string, create a function
0048       sel_fcn = inline(sel_fcn, 'z');
0049     end
0050     ff = feval(sel_fcn,ct3);
0051 end
0052 imdl.rec_model.coarse2fine(~ff,:) = [];
0053 imdl.rec_model.elems(~ff,:) = [];
0054 
0055 fb = any(imdl.rec_model.coarse2fine,1);
0056 imdl.rec_model.coarse2fine(:,~fb) = [];
0057 imdl.fwd_model.coarse2fine(:,~fb) = [];
0058 imdl.solve_use_matrix.RM(~fb,:) = [];
0059 if isfield(imdl.solve_use_matrix,'PJt')
0060    imdl.solve_use_matrix.PJt(~fb,:) = [];
0061 end
0062 
0063 if keep3D; 
0064    imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0065    return;
0066 end
0067 
0068 imdl.rec_model.nodes(:,3) = [];
0069 imdl.rec_model.elems(:,4) = [];
0070 
0071 nelems = imdl.rec_model.elems;
0072 nnodes = unique(nelems(:));
0073 nnidx = zeros(max(nnodes),1);
0074 nnidx(nnodes) = 1:length(nnodes);
0075 nnodes = imdl.rec_model.nodes(nnodes,:);
0076 nelems = reshape(nnidx(nelems),size(nelems));
0077 imdl.rec_model.nodes = nnodes;
0078 imdl.rec_model.elems = nelems;
0079 imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0080 
0081 function do_unit_test
0082    [vh,vi] = test_fwd_solutions;
0083    % inverse geometry
0084    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0085    fmdl= mystim_square(fmdl);
0086 
0087    vopt.imgsz = [32 32];
0088    vopt.zvec = linspace( 0.75,3.25,6);
0089    vopt.save_memory = 1;
0090    opt.noise_figure = 2;
0091    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0092    imdl = mk_GREIT_model(imdl, 0.2, [], opt);
0093 
0094    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [4280,928]);
0095    unit_test_cmp('RM', imdl.solve_use_matrix.RM(1,1:2), ...
0096        [7.896475314622105 -3.130412179150207], 1e-10);
0097 %      [7.896475314707559 -3.130412179315593], 1e-10); % older values
0098 
0099 
0100 
0101    img = inv_solve(imdl, vh, vi);
0102    unit_test_cmp('img size', size(img.elem_data), [4280,5]);
0103    [mm,ll] =max(img.elem_data(:,1));
0104    unit_test_cmp('img', [mm,ll], ...
0105        [0.426631207850641, 1274], 1e-10);
0106 
0107    img.show_slices.img_cols = 1;
0108    slice1 = calc_slices(img);
0109    subplot(231); show_fem(fmdl); title 'fmdl'
0110    subplot(234); show_fem(img);  title '3D'
0111    subplot(232); show_slices(img,[inf,inf,2]); title '3D slice';
0112    imd2 = solve_RM_2Dslice(imdl, 2.0);
0113 
0114    unit_test_cmp('RM size', size(imd2.solve_use_matrix.RM), [856,928]);
0115    unit_test_cmp('RM', imd2.solve_use_matrix.RM(1,1:2), ...
0116        [-13.546930204198315   9.664897892799864], 1e-10);
0117 %      [-13.546930204456816   9.664897892901678], 1e-10);
0118 
0119    img = inv_solve(imd2, vh, vi);
0120    unit_test_cmp('img size', size(img.elem_data), [856,5]);
0121    [mm,ll] =max(img.elem_data(:,1));
0122    unit_test_cmp('img', [mm,ll], ...
0123        [0.031608449353163, 453], 1e-10);
0124    slice2 = calc_slices(img);
0125 
0126    unit_test_cmp('slice Nan', isnan(slice1), isnan(slice2))
0127    slice1(isnan(slice1))= 0;
0128    slice2(isnan(slice2))= 0;
0129    unit_test_cmp('slice value', slice1, slice2, 1e-10)
0130 
0131    
0132    imdl.solve_RM_2Dslice.keep3D = 1;
0133    imd3 = solve_RM_2Dslice(imdl, 1.0);
0134    img = inv_solve(imd3, vh, vi);
0135    subplot(235); show_fem(img);
0136    
0137    sel_fcn = inline('abs(z-1.0)<0.2','z');
0138    imd3 = solve_RM_2Dslice(imdl, sel_fcn);
0139 
0140    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ;
0141    imd3 = solve_RM_2Dslice(imdl, sel_fcn);
0142    img = inv_solve(imd3, vh, vi);
0143    subplot(236); show_fem(img);
0144 
0145 function fmdl = mystim_square(fmdl);
0146    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0147    [fmdl.stimulation, fmdl.meas_select] =  ...
0148        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0149 
0150 function [vh,vi] = test_fwd_solutions;
0151    posns= linspace(1.0,3.0,5);
0152    str=''; for i=1:length(posns);
0153       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0154       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0155    end
0156    extra{i+1} = str;
0157    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0158    fmdl = mystim_square(fmdl);
0159    
0160    img= mk_image(fmdl,1);
0161    vh = fwd_solve(img);
0162    %vh = add_noise(2000, vh);
0163    for i=1:length(posns);
0164       img= mk_image(fmdl,1);
0165       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0166       vi{i} = fwd_solve(img);
0167    end;
0168    vi = [vi{:}];