Home > eidors > models > solve_RM_2Dslice.m

solve_RM_2Dslice

PURPOSE ^

SOLVE_RM_2DSLICE: cut slices out of a 3D model

SYNOPSIS ^

function imdl = solve_RM_2Dslice(imdl, sel_fcn)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5791 2018-05-22 15:16:02Z 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 
0060 if keep3D; 
0061    imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0062    return;
0063 end
0064 
0065 imdl.rec_model.nodes(:,3) = [];
0066 imdl.rec_model.elems(:,4) = [];
0067 
0068 nelems = imdl.rec_model.elems;
0069 nnodes = unique(nelems(:));
0070 nnidx = zeros(max(nnodes),1);
0071 nnidx(nnodes) = 1:length(nnodes);
0072 nnodes = imdl.rec_model.nodes(nnodes,:);
0073 nelems = reshape(nnidx(nelems),size(nelems));
0074 imdl.rec_model.nodes = nnodes;
0075 imdl.rec_model.elems = nelems;
0076 imdl.rec_model.boundary = find_boundary(imdl.rec_model);
0077 
0078 function do_unit_test
0079    [vh,vi] = test_fwd_solutions;
0080    % inverse geometry
0081    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0082    fmdl= mystim_square(fmdl);
0083 
0084    vopt.imgsz = [32 32];
0085    vopt.zvec = linspace( 0.75,3.25,6);
0086    vopt.save_memory = 1;
0087    opt.noise_figure = 2;
0088    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0089    imdl = mk_GREIT_model(imdl, 0.2, [], opt);
0090 
0091    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [4280,928]);
0092    unit_test_cmp('RM', imdl.solve_use_matrix.RM(1,1:2), ...
0093        [7.896475314882005  -3.130412179417599], 1e-10);
0094 
0095    img = inv_solve(imdl, vh, vi);
0096    unit_test_cmp('img size', size(img.elem_data), [4280,5]);
0097    [mm,ll] =max(img.elem_data(:,1));
0098    unit_test_cmp('img', [mm,ll], ...
0099        [0.426631207850641, 1274], 1e-10);
0100 
0101    img.show_slices.img_cols = 1;
0102    slice1 = calc_slices(img);
0103    subplot(231); show_fem(fmdl); title 'fmdl'
0104    subplot(234); show_fem(img);  title '3D'
0105    subplot(232); show_slices(img,[inf,inf,2]); title '3D slice';
0106    imd2 = solve_RM_2Dslice(imdl, 2.0);
0107 
0108    unit_test_cmp('RM size', size(imd2.solve_use_matrix.RM), [856,928]);
0109    unit_test_cmp('RM', imd2.solve_use_matrix.RM(1,1:2), ...
0110        [-13.546930204647675   9.664897892828284], 1e-10);
0111    img = inv_solve(imd2, vh, vi);
0112    unit_test_cmp('img size', size(img.elem_data), [856,5]);
0113    [mm,ll] =max(img.elem_data(:,1));
0114    unit_test_cmp('img', [mm,ll], ...
0115        [0.031608449353163, 453], 1e-10);
0116    slice2 = calc_slices(img);
0117 
0118    unit_test_cmp('slice Nan', isnan(slice1), isnan(slice2))
0119    slice1(isnan(slice1))= 0;
0120    slice2(isnan(slice2))= 0;
0121    unit_test_cmp('slice value', slice1, slice2, 1e-10)
0122 
0123    
0124    imdl.solve_RM_2Dslice.keep3D = 1;
0125    imd3 = solve_RM_2Dslice(imdl, 1.0);
0126    img = inv_solve(imd3, vh, vi);
0127    subplot(235); show_fem(img);
0128    
0129    sel_fcn = inline('abs(z-1.0)<0.2','z');
0130    imd3 = solve_RM_2Dslice(imdl, sel_fcn);
0131 
0132    sel_fcn = 'abs(z-1.0)<0.2 | abs(z-2.0)<0.2' ;
0133    imd3 = solve_RM_2Dslice(imdl, sel_fcn);
0134    img = inv_solve(imd3, vh, vi);
0135    subplot(236); show_fem(img);
0136 
0137 function fmdl = mystim_square(fmdl);
0138    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0139    [fmdl.stimulation, fmdl.meas_select] =  ...
0140        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0141 
0142 function [vh,vi] = test_fwd_solutions;
0143    posns= linspace(1.0,3.0,5);
0144    str=''; for i=1:length(posns);
0145       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0146       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0147    end
0148    extra{i+1} = str;
0149    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0150    fmdl = mystim_square(fmdl);
0151    
0152    img= mk_image(fmdl,1);
0153    vh = fwd_solve(img);
0154    %vh = add_noise(2000, vh);
0155    for i=1:length(posns);
0156       img= mk_image(fmdl,1);
0157       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0158       vi{i} = fwd_solve(img);
0159    end;
0160    vi = [vi{:}];

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005