0001 function [vh,vi,xyr_pt]= simulate_2d_movement( n_sims, fmdl, rad_pr, movefcn )
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
0032
0033
0034
0035
0036
0037
0038
0039 if nargin>=1 && ischar(n_sims) && strcmp(n_sims,'UNIT_TEST'); do_unit_test; return; end
0040
0041 if nargin <1
0042 n_sims = 200;
0043 end
0044
0045 if nargin<2 || isempty(fmdl)
0046 n_circles = 36;
0047 n_elec= 16;
0048 fmdl= mk_fwd_model(n_circles, n_elec);
0049 end
0050
0051 if nargin<3; rad_pr= []; end
0052 if nargin<4; movefcn= 1; end
0053
0054 copt.cache_obj = {n_sims,fmdl, rad_pr, movefcn};
0055 copt.fstr = 'simulate_2d_movement';
0056 [vh,vi,xyr_pt]= eidors_cache(@do_simulate_2d_movement,{n_sims, fmdl, rad_pr, movefcn},copt );
0057
0058
0059 function [vh,vi,xyr_pt]= do_simulate_2d_movement( n_sims, fmdl, rad_pr, movefcn )
0060 if isempty(rad_pr)
0061 radius= 2/3;
0062 rp= .05;
0063 else
0064 radius= rad_pr(1);
0065 rp= rad_pr(2);
0066 end
0067
0068 mv_start = 0;
0069 mv_end = 1;
0070 if isnumeric(movefcn)
0071 if length(movefcn)>=2; mv_start = movefcn(2); end
0072 if length(movefcn)>=3; mv_end = movefcn(3); end
0073 if movefcn(1)==1
0074 movefcn = @rotation_path;
0075 elseif movefcn(1)==2
0076 movefcn = @straight_out;
0077 else
0078 error('value of movefcn not understood');
0079 end
0080 else
0081
0082 end
0083
0084 n_elems= size(fmdl.elems,1);
0085 img= eidors_obj('image','simulate_movement', ...
0086 'fwd_model', fmdl, ...
0087 'elem_data', ones(n_elems,1) );
0088 vh= fwd_solve(img);
0089
0090 if 0
0091 np= 256;
0092 maxxy= max(fmdl.nodes);
0093 minxy= min(fmdl.nodes);
0094 [x,y]=meshgrid( linspace(minxy(1),maxxy(1),np), ...
0095 linspace(minxy(2),maxxy(2),np) );
0096 [eptr,vol]= img_mapper2(fmdl.nodes', fmdl.elems', np, np);
0097 else
0098
0099 mdl_pts = interp_mesh( fmdl, 8);
0100 x= squeeze( mdl_pts(:,1,:) );
0101 y= squeeze( mdl_pts(:,2,:) );
0102 end
0103
0104
0105
0106
0107
0108
0109
0110 target_conductivity= .1;
0111
0112 for i=1:n_sims
0113 f_frac= mv_start + ( (i-1)/n_sims ) * (mv_end - mv_start);
0114 fprintf('simulating %d / %d (f_frac=%0.2f) \n',i,n_sims, f_frac);
0115
0116 [xp,yp]= feval(movefcn, f_frac, radius);
0117
0118 if 0
0119 xyr_pt(:,i)= [xp;-yp;rp];
0120 ff= find( (x(:)-xp).^2 + (y(:)-yp).^2 <= rp^2 )';
0121 obj_n= sparse( eptr(ff)+1,1,1, n_elems+1, 1);
0122 obj_n= full(obj_n(2:end));
0123
0124 img.elem_data= 1 + target_conductivity * (obj_n./vol);
0125 else
0126 xyr_pt(:,i)= [xp;yp;rp];
0127 ff= (x-xp).^2 + (y-yp).^2 <= rp^2;
0128 img.elem_data= 1 + target_conductivity * mean(ff,2);
0129 end
0130
0131 vi(i)= fwd_solve( img );
0132
0133 end
0134
0135
0136 vi= [vi(:).meas];
0137 vh= [vh(:).meas];
0138
0139
0140
0141
0142 function [xp,yp] = rotation_path(f_frac, radius);
0143 xp= radius * cos(f_frac*2*pi);
0144 yp= radius * sin(f_frac*2*pi);
0145
0146 function [xp,yp] = straight_out(f_frac, radius);
0147 xp= radius*f_frac;
0148 yp= 0;
0149
0150
0151
0152
0153 function [EPTR,VOL]= img_mapper2(NODE, ELEM, npx, npy );
0154 xmin = min(NODE(1,:)); xmax = max(NODE(1,:));
0155 xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0156
0157 ymin = min(NODE(2,:)); ymax = max(NODE(2,:));
0158 ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0159
0160 range= max([xrange, yrange]);
0161 [x y]=meshgrid( ...
0162 linspace( xmean - range*0.50, xmean + range*0.50, npx ), ...
0163 linspace( ymean + range*0.50, ymean - range*0.50, npy ) );
0164 v_yx= [-y(:) x(:)];
0165 turn= [0 -1 1;1 0 -1;-1 1 0];
0166 EPTR=zeros(npy,npx);
0167
0168
0169
0170
0171
0172 e= size(ELEM,2);
0173 VOL= zeros(e,1);
0174 for j= 1:e
0175
0176 xy= NODE(:,ELEM(:,j))';
0177
0178
0179 endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0180 & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0181
0182 a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0183 VOL(j)= abs(sum(a));
0184
0185 aa= sum(abs(ones(length(endr),1)*a'+ ...
0186 v_yx(endr,:)*xy'*turn)');
0187 endr( abs( ( VOL(j)-aa ) ./ VOL(j) ) >1e-8)=[];
0188 EPTR(endr)= j;
0189 end
0190
0191
0192 function mdl_2d= mk_fwd_model(n_circles, n_elec)
0193
0194 params= mk_circ_tank(n_circles, [], n_elec);
0195 n_rings= 1;
0196 options= {'no_meas_current','no_rotate_meas','do_redundant'};
0197 [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0198 params.stimulation= st;
0199 params.meas_select= els;
0200 params.solve= 'eidors_default';
0201 params.system_mat= 'eidors_default';
0202 params.jacobian= 'eidors_default';
0203 params.normalize_measurements= 0;
0204 mdl_2d = eidors_obj('fwd_model', params);
0205
0206 function do_unit_test
0207 N_TEST = 5;
0208 imdl = mk_common_model( 'c2c2', 16 );
0209 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST);
0210 subplot(421);
0211 imgs = inv_solve(imdl, vh, vi);
0212 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0213
0214 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [0.3,0.01],1);
0215 subplot(422);
0216 imgs = inv_solve(imdl, vh, vi);
0217 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0218
0219 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [0.9,0.01],2);
0220 subplot(423);
0221 imgs = inv_solve(imdl, vh, vi);
0222 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0223
0224 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [],[1,0.5,0.4]);
0225 subplot(424);
0226 imgs = inv_solve(imdl, vh, vi);
0227 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0228
0229 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [],@test_movefcn);
0230 subplot(425);
0231 imgs = inv_solve(imdl, vh, vi);
0232 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0233
0234 fmdl = mk_common_model('a2c2',16); fmdl= fmdl.fwd_model;
0235 fmdl.nodes = fmdl.nodes*1.5;
0236 [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, fmdl, [],@test_movefcn);
0237 subplot(426);
0238 imgs = inv_solve(imdl, vh, vi);
0239 imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0240
0241 function [xp,yp] = test_movefcn(f_frac, radius);
0242 ff = radius/sqrt(2);
0243 xp= ff*f_frac; yp= ff*f_frac;