solve_RM_2Dslice

 SELECT_RM_SLICE: cut slices out of an inverse model with a 
  reconstruction matrix solving onto a 3D (voxel) grid

 NOTE that the 3D grid must be of the structure created by MK_GRID_MODEL
 or MK_VOXEL_VOLUME.

 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= select_RM_slice(imdl, 1.0);
 Or using 
    sel_fcn = @(z) abs(z-1.0)<0.2; OR
    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ; % two planes
    imdl2= select_RM_slice(imdl, sel_fcn);

 SELECT_RM_SLICE applies the selection criteria to element centres. If no
 elements are captured, increase vert_tol (see below).
 
 Note that the reconstruction planes are between the
 cut planes specified in vopt.zvec

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

 See also MK_GRID_MODEL, MK_VOXEL_VOLUME, SELECT_RM_SLICE

0001 function imdl = select_RM_slice(imdl, sel_fcn)
0002 % SELECT_RM_SLICE: cut slices out of an inverse model with a
0003 %  reconstruction matrix solving onto a 3D (voxel) grid
0004 %
0005 % NOTE that the 3D grid must be of the structure created by MK_GRID_MODEL
0006 % or MK_VOXEL_VOLUME.
0007 %
0008 % Basic usage
0009 % For example, given a 3D GREIT model
0010 %    vopt.zvec= 0.5:0.2:2.5;
0011 %    vopt.imgsz = [32 32];
0012 %    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0013 %    imdl = mk_GREIT_model(imdl, 0.20, [], opt);
0014 %
0015 % One could cut the z=1 slice as
0016 %    imdl2= select_RM_slice(imdl, 1.0);
0017 % Or using
0018 %    sel_fcn = @(z) abs(z-1.0)<0.2; OR
0019 %    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ; % two planes
0020 %    imdl2= select_RM_slice(imdl, sel_fcn);
0021 %
0022 % SELECT_RM_SLICE applies the selection criteria to element centres. If no
0023 % elements are captured, increase vert_tol (see below).
0024 %
0025 % Note that the reconstruction planes are between the
0026 % cut planes specified in vopt.zvec
0027 %
0028 % options:
0029 %  imdl.select_RM_slice.vert_tol = .001;
0030 %       % Vertical tolerance for elems on plane (Default .001)
0031 %  imdl.select_RM_slice.keep3D = 1; %keep 3D aspect of cuts
0032 %
0033 % See also MK_GRID_MODEL, MK_VOXEL_VOLUME, SELECT_RM_SLICE
0034 
0035 % (C) 2015-2024 Andy Adler and Bartlomiej Grychtol
0036 % License: GPL version 2 or version 3
0037 % $Id: solve_RM_2Dslice.m 7087 2024-12-20 17:44:41Z aadler $
0038 
0039 warning('EIDORS:deprecated','SELECT_RM_SLICE is deprecated as of 01-Dec-2024. Use SELECT_RM_SLICE instead.');
0040 
0041 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 % Options
0044 vert_tol = .001; try
0045   vert_tol = imdl.select_RM_slice.vert_tol;
0046 end
0047 
0048 keep3D=0; try
0049   keep3D = imdl.select_RM_slice.keep3D;
0050 end
0051 
0052 ctr = interp_mesh(imdl.rec_model); ct3= ctr(:,3);
0053 
0054 if isnumeric( sel_fcn )
0055    ff = abs(ct3-sel_fcn) < vert_tol;
0056 else
0057     if ischar(sel_fcn) %  we have a string, create a function
0058       sel_fcn = inline(sel_fcn, 'z');
0059     end
0060     ff = feval(sel_fcn,ct3);
0061 end
0062 
0063 % mk_grid_model creates 6 tets per voxel. We want to capture the entire
0064 % voxel.
0065 c2f = imdl.rec_model.coarse2fine;
0066 vox_frac = c2f' * ff / 6; % how much of each voxel got captured
0067 ff = c2f * vox_frac>0; % now capture the entire voxel
0068 
0069 if ~any(ff)
0070    error('Found no elements. Check your input or increase opt.vert_tol.')
0071 end
0072 
0073 imdl.rec_model.coarse2fine(~ff,:) = [];
0074 imdl.rec_model.elems(~ff,:) = [];
0075 
0076 fb = any(imdl.rec_model.coarse2fine,1);
0077 imdl.rec_model.coarse2fine(:,~fb) = [];
0078 imdl.fwd_model.coarse2fine(:,~fb) = [];
0079 imdl.solve_use_matrix.RM(~fb,:) = [];
0080 if isfield(imdl.solve_use_matrix,'PJt')
0081    imdl.solve_use_matrix.PJt(~fb,:) = [];
0082 end
0083 
0084 if keep3D; 
0085    imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0086    return;
0087 end
0088 
0089 % convert to 2D (only makes sense for a single slice)
0090 imdl.rec_model.nodes(:,3) = [];
0091 keep = false(size(imdl.rec_model.elems,1),1);
0092 keep(1:6:end) = true; keep(6:6:end) = true; % intimate knowledge of mk_grid_model
0093 imdl.rec_model.elems = imdl.rec_model.elems(keep,:);
0094 imdl.rec_model.elems(:,4) = [];
0095 imdl.rec_model.coarse2fine = imdl.rec_model.coarse2fine(keep,:);
0096 
0097 nelems = imdl.rec_model.elems;
0098 nnodes = unique(nelems(:));
0099 nnidx = zeros(max(nnodes),1);
0100 nnidx(nnodes) = 1:length(nnodes);
0101 nnodes = imdl.rec_model.nodes(nnodes,:);
0102 nelems = reshape(nnidx(nelems),size(nelems));
0103 imdl.rec_model.nodes = nnodes;
0104 imdl.rec_model.elems = nelems;
0105 imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0106 
0107 function do_unit_test
0108    warning('off','EIDORS:FirstImageOnly');
0109    [vh,vi] = test_fwd_solutions;
0110    % inverse geometry
0111    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0112    fmdl= mystim_square(fmdl);
0113 
0114    vopt.imgsz = [32 32];
0115    vopt.zvec = linspace( 0.75,3.25,6);
0116    vopt.save_memory = 1;
0117 %  opt.noise_figure = 2;
0118    weight = 8.42412109211;
0119    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0120    imdl = mk_GREIT_model(imdl, 0.2, weight, opt);
0121 
0122    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [4280,928]);
0123 %    unit_test_cmp('RM', imdl.solve_use_matrix.RM(1,1:2), ...
0124 %        [7.896475314622105 -3.130412179150207], 1e-8);
0125 
0126 
0127 
0128    img = inv_solve(imdl, vh, vi);
0129    unit_test_cmp('img size', size(img.elem_data), [4280,5]);
0130    [mm,ll] =max(img.elem_data(:,1));
0131 %    unit_test_cmp('img', [mm,ll], ...
0132 %        [0.426631207850641, 1274], 1e-10);
0133 
0134    img.show_slices.img_cols = 1;
0135    slice1 = calc_slices(img,[inf inf 2]);
0136    subplot(231); show_fem(fmdl); title 'fmdl'
0137    subplot(234); show_fem(img);  title '3D'
0138    subplot(232); show_slices(img,[inf,inf,2]); title '3D slice';
0139    got_error = false;
0140    try 
0141       imd2 = select_RM_slice(imdl, 2.2);
0142    catch 
0143       got_error = true;
0144    end
0145    unit_test_cmp('No elem error', got_error, true);
0146    
0147    imd2 = select_RM_slice(imdl, 2.0);
0148 
0149    unit_test_cmp('RM size', size(imd2.solve_use_matrix.RM), [856,928]);
0150 %    unit_test_cmp('RM', imd2.solve_use_matrix.RM(1,1:2), ...
0151 %        [-13.546930204198315   9.664897892799864], 1e-8);
0152 
0153    img = inv_solve(imd2, vh, vi);
0154    unit_test_cmp('img size', size(img.elem_data), [856,5]);
0155    [mm,ll] =max(img.elem_data(:,1));
0156 %    unit_test_cmp('img', [mm,ll], ...
0157 %        [0.031608449353163, 453], 1e-10);
0158    slice2 = calc_slices(img);
0159 
0160    unit_test_cmp('slice Nan', isnan(slice1), isnan(slice2))
0161    slice1(isnan(slice1))= 0;
0162    slice2(isnan(slice2))= 0;
0163    unit_test_cmp('slice value', slice1, slice2, 1e-10)
0164 
0165    
0166    imdl.select_RM_slice.keep3D = 1;
0167    imd3 = select_RM_slice(imdl, 1.0);
0168    img = inv_solve(imd3, vh, vi);
0169    subplot(235); show_fem(img);
0170    
0171    sel_fcn = inline('abs(z-1.0)<0.2','z');
0172    if exist('OCTAVE_VERSION')
0173       sel_fcn = @(z) abs(z-1.0)<0.2; 
0174    end
0175    imd3 = select_RM_slice(imdl, sel_fcn);
0176 
0177    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ;
0178    if exist('OCTAVE_VERSION')
0179       sel_fcn = @(z) abs(z-1.0)<0.2 | abs(z-2.0)<0.2;
0180    end
0181    imd3 = select_RM_slice(imdl, sel_fcn);
0182    img = inv_solve(imd3, vh, vi);
0183    subplot(236); show_fem(img);
0184 
0185    warning('on','EIDORS:FirstImageOnly');
0186 
0187 function fmdl = mystim_square(fmdl);
0188    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0189    [fmdl.stimulation, fmdl.meas_select] =  ...
0190        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0191 
0192 function [vh,vi] = test_fwd_solutions;
0193    posns= linspace(1.0,3.0,5);
0194    str=''; for i=1:length(posns);
0195       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0196       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0197    end
0198    extra{i+1} = str;
0199    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0200    fmdl = mystim_square(fmdl);
0201    
0202    img= mk_image(fmdl,1);
0203    vh = fwd_solve(img);
0204    %vh = add_noise(2000, vh);
0205    for i=1:length(posns);
0206       img= mk_image(fmdl,1);
0207       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0208       vi{i} = fwd_solve(img);
0209    end;
0210    vi = [vi{:}];