select_RM_slice

 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

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