0001 function imdl = solve_RM_2Dslice(imdl, sel_fcn)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'); do_unit_test; return; end
0032
0033
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)
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
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
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{:}];