moving_tank_objs

 MOVING_TANK_OBJS: create movies of objects moving in tanks
 Usage:
 imgs= moving_tank_objs(data_sel, inv_sel, options)

 imgs = reconstructed images

 data_sel select data_sel to use
   data_sel: 10 => target across tank data
             20 => target around tank data
             21 => target around tank data (10 positions)
             30 => from netgen simulations (spiral)
             32 => from netgen simulations (spirograph)
             40 => show netgen FEM simulation

   if data_sel is a vector of EIT data, then it will
      be used for the reconstruction. The first measurement
      is the homogeneous, and the next are the data

 inv_sel
  2D reconstructions
   inv_sel = 1   => inv_solve_diff_GN_one_step
   inv_sel = 1.1 => inv_solve_diff_GN_one_step (576 elems)
   inv_sel = 2   => inv_solve_diff_kalman

   if inv_sel is a inv_model structure, then use it to
   do reconstructions.

 options (not changed if value is NaN);
   options(1) - noise figure
   options(2) - time_steps
   options(3) - time_weight
 
 Create moving objects and tanks
 $Id: moving_tank_objs.m 3128 2012-06-08 17:04:21Z bgrychtol $

0001 function imgs= moving_tank_objs(data_sel, inv_sel, options)
0002 % MOVING_TANK_OBJS: create movies of objects moving in tanks
0003 % Usage:
0004 % imgs= moving_tank_objs(data_sel, inv_sel, options)
0005 %
0006 % imgs = reconstructed images
0007 %
0008 % data_sel select data_sel to use
0009 %   data_sel: 10 => target across tank data
0010 %             20 => target around tank data
0011 %             21 => target around tank data (10 positions)
0012 %             30 => from netgen simulations (spiral)
0013 %             32 => from netgen simulations (spirograph)
0014 %             40 => show netgen FEM simulation
0015 %
0016 %   if data_sel is a vector of EIT data, then it will
0017 %      be used for the reconstruction. The first measurement
0018 %      is the homogeneous, and the next are the data
0019 %
0020 % inv_sel
0021 %  2D reconstructions
0022 %   inv_sel = 1   => inv_solve_diff_GN_one_step
0023 %   inv_sel = 1.1 => inv_solve_diff_GN_one_step (576 elems)
0024 %   inv_sel = 2   => inv_solve_diff_kalman
0025 %
0026 %   if inv_sel is a inv_model structure, then use it to
0027 %   do reconstructions.
0028 %
0029 % options (not changed if value is NaN);
0030 %   options(1) - noise figure
0031 %   options(2) - time_steps
0032 %   options(3) - time_weight
0033 %
0034 % Create moving objects and tanks
0035 % $Id: moving_tank_objs.m 3128 2012-06-08 17:04:21Z bgrychtol $
0036 
0037 clim= [];
0038 
0039 if nargin<1; data_sel = 10; end
0040 if nargin<2; inv_sel  = 1; end
0041 
0042 if isnumeric(data_sel) & size(data_sel)==[1,1]
0043   switch data_sel
0044     case 10
0045         load montreal_data_1995;
0046         vh= zc_h_demo4;
0047         vi= zc_demo4;
0048         filename= 'target-across';
0049 
0050     case 20
0051         load montreal_data_1995;
0052         vh= zc_h_demo3;
0053         vi= zc_demo3;
0054         filename= 'target-around';
0055 
0056     case 21
0057         load montreal_data_1995;
0058         vh= zc_h_demo3;
0059         vi= zc_demo3(:,1:2:20);
0060         filename= 'target-around10';
0061 
0062     case 22
0063         load iirc_data_2006
0064         vh= real(v_reference);
0065         vi= real(v_rotate(:,1:30));
0066         filename= 'iirc-around-clean';
0067 
0068     case 22.1
0069         load iirc_data_2006
0070         vh= v_reference;
0071         vi= v_rotate(:,11:20);
0072         filename= 'iirc-around-clean';
0073 
0074     case 23
0075         load iirc_data_2006
0076         vh= v_reference;
0077         vi= v_rotate(:,1:30);
0078         [vi,vh]= add_noise(vi,vh,2)
0079         filename= 'iirc-around-noise';
0080 
0081       case 24
0082         load iirc_data_2006
0083         vh= v_reference;
0084         vi= v_rotate;
0085         filename= 'iirc-around-clean';
0086 
0087 
0088     case 30
0089         load netgen_moving_ball
0090         vh= homg_tank;
0091         vi= target_spiral;
0092         filename= 'netgen-spiral';
0093 
0094     case 31
0095         load netgen_moving_ball
0096         vh= homg_tank;
0097         vi= target_spiral(40:60);
0098         filename= 'netgen-spiral-part';
0099 
0100     case 32
0101         load netgen_moving_ball
0102         vh= homg_tank;
0103         vi= target_spirograph_slow;
0104         filename= 'netgen-spirograph-slow';
0105 
0106     case 33
0107         load netgen_moving_ball
0108         vh= homg_tank;
0109         vi= target_spirograph_slow([31:52,72:90]);
0110         randn('seed',20);
0111         for i=1:length(vi)
0112 %          vi(i).meas = vi(i).meas + 25e-6*randn(size(vi(i).meas));
0113         end
0114         filename= 'netgen-spirograph-slow-part';
0115 
0116     otherwise
0117         error('Don''t recognize data_sel');
0118   end
0119 else
0120   vh= data_sel(:,1);
0121   vi= data_sel(:,2:end);
0122   filename= 'user-data';
0123 end
0124 
0125 % noise_figure?
0126 % Calc BR and posn error
0127 
0128 try; if inv_sel.type=='inv_model'
0129    imdl= inv_sel;
0130    inv_sel= 'set'
0131 end; end
0132 
0133 switch inv_sel
0134     case 'set'
0135         eidors_msg('Using imdl=%s',imdl.name,2);
0136 
0137     case 1
0138         imdl= set_basic( 'b2c' );
0139 
0140     case 1.1
0141         imdl= set_basic( 'c2c' );
0142 
0143     case 1.2
0144         imdl= set_basic( 'd2c' );
0145 
0146     case 2
0147         imdl= set_basic( 'b2c' );
0148         imdl.solve= @inv_solve_diff_kalman;
0149 
0150     case 2.1
0151         imdl= set_basic( 'c2c' );
0152         imdl.solve= @inv_solve_diff_kalman;
0153 
0154     case 3
0155         imdl= set_basic( 'b2c' );
0156         imdl.solve= @aa_inv_conj_grad;
0157 
0158     case 4
0159         imdl= set_basic( 'b2c', 'temporal', [3, .6] );
0160 
0161     case 4.1
0162         imdl= set_basic( 'c2c', 'temporal', [3, .6] );
0163 
0164     case 5
0165         imdl= set_basic( 'b2c', 'temporal', [8, 1] );
0166 
0167     case 11
0168         imdl= set_basic( 'b2c' );
0169 
0170         vi = extract_subframes(vi,4);
0171         vi = repack_subframes(vi,4);
0172 
0173     case 11.1
0174         imdl= set_basic( 'c2c' );
0175 
0176         vi = extract_subframes(vi,4);
0177         vi = repack_subframes(vi,4);
0178 
0179     case 12
0180         imdl= set_basic( 'b2c' );
0181         imdl.solve= @inv_solve_diff_kalman;
0182         [imdl.fwd_model.stimulation(:).delta_time]=deal(1);
0183 
0184         vi = extract_subframes(vi,4);
0185 
0186     case 12.1
0187         imdl= set_basic( 'c2c' );
0188         imdl.solve= @inv_solve_diff_kalman;
0189         [imdl.fwd_model.stimulation(:).delta_time]=deal(1);
0190 
0191         vi = extract_subframes(vi,4);
0192 
0193     case 14
0194         imdl= set_basic( 'b2c' );
0195         time_steps= 3;
0196 
0197         imdl.RtR_prior= @prior_time_smooth;
0198         imdl.prior_time_smooth.space_prior= @prior_noser;
0199         imdl.prior_time_smooth.time_weight= .5;
0200         imdl.prior_time_smooth.time_steps=  time_steps;
0201 
0202         imdl.solve= @inv_solve_time_prior;
0203         imdl.inv_solve_time_prior.time_steps=   time_steps;
0204 
0205         [imdl.fwd_model.stimulation(:).delta_time]=deal(1);
0206 
0207         vi = extract_subframes(vi,4);
0208 
0209     otherwise
0210         error(['inv_sel (' inv_sel ') not recognized']);
0211 end
0212 
0213 % options
0214 if nargin>=3
0215    if length(options) >=1 if ~isnan(options(1))
0216       imdl= rmfield(imdl,'hyperparameter');
0217       imdl.hyperparameter.func= @choose_noise_figure;
0218       imdl.hyperparameter.tgt_elems= [1:4];
0219       imdl.hyperparameter.noise_figure= options(1);
0220          filename = sprintf('%s-hp=%g', filename, options(1));
0221    end; end
0222    if length(options) >=2 if ~isnan(options(2))
0223       imdl.prior_time_smooth.time_steps=  options(2);
0224       imdl.inv_solve_time_prior.time_steps=   options(2);
0225       filename = sprintf('%s-ts=%d', filename, options(2));
0226    end; end
0227    if length(options) >=3 if ~isnan(options(3))
0228       imdl.prior_time_smooth.time_weight=  options(3);
0229       filename = sprintf('%s-tw=%4.2f', filename, options(3));
0230    end; end
0231 end
0232 
0233 t=cputime;
0234 imgs= inv_solve(imdl,vh,vi);
0235 fprintf('solve time=%f\n',cputime-t);
0236 animate_reconstructions(filename, imgs);
0237 
0238 % original = [ 1.1 2.1 3.1 4.1 5.1
0239 %              1.2 2.2 3.2 4.2 5.2
0240 %              1.3 2.3 3.3 4.3 5.3
0241 %              1.4 2.4 3.4 4.4 5.4
0242 %              1.5 2.5 3.5 4.5 5.5  ]
0243 % subseq = 2
0244 % output   = [ 1.1 3.1
0245 %              1.2 4.2
0246 %              2.3 4.3
0247 %              2.4 5.4
0248 %              3.5 5.5 ]
0249 %
0250 % Hardcoded for 16 electrode 2D adjacent stimulation
0251 function ve= extract_subframes( vv, subseq)
0252    vv= remove_curr_elecs(vv);
0253 
0254    ve= zeros( size(vv,1), floor(size(vv,2)*subseq/16) );
0255    dst=0; src=0;
0256    for k=0:16*size(ve,2)-1;
0257       pat = (1:13) + 13*rem(k,16);
0258       if 0==rem(k,16);     dst= dst+1; end
0259       if 0==rem(k,subseq); src= src+1; end
0260       ve(pat,dst) = vv(pat, src);
0261    end
0262 
0263 
0264 % original = [ 1.1 2.1 3.1 4.1 5.1
0265 %              1.2 2.2 3.2 4.2 5.2
0266 %              1.3 2.3 3.3 4.3 5.3
0267 %              1.4 2.4 3.4 4.4 5.4
0268 %              1.5 2.5 3.5 4.5 5.5  ]
0269 % subseq = 2
0270 %
0271 % output   = [ 1.1 2.1 2.1 3.1 3.1
0272 %              1.2 2.2 2.2 2.2 3.2
0273 %              1.3 1.3 2.3 2.3 3.3
0274 %              1.4 1.4 2.4 2.4 2.4
0275 %              1.5 1.5 1.5 3.5 3.5 ]
0276 %
0277 % Hardcoded for 16 electrode 2D adjacent stimulation
0278 
0279 function ve = repack_subframes(vv,subseq);
0280    vv= remove_curr_elecs(vv);
0281 
0282    % duplicate first and last
0283    ve= zeros( size(vv,1), floor(size(vv,2)*16/subseq) );
0284    vv= vv(:,[1:end,end]);
0285    dst=0; src=0;
0286    pat = (1:208)';
0287    for i=1:size(ve,2)
0288       for k=rem(subseq*i + (-subseq:-1),16);
0289          pat(13*k + (1:13)) = pat(13*k + (1:13)) + 208;
0290       end
0291       ve(:,i) = vv(pat);
0292    end
0293 
0294 
0295 function vv= remove_curr_elecs(vv)
0296    if isstruct(vv);
0297       vv= [vv(:).meas];
0298    end
0299    [st, els]= mk_stim_patterns(16, 1, '{ad}','{ad}', {}, 10);
0300 
0301    % thow away measurements at current elecs if needed
0302    if size(vv,1)==size(els,1)
0303       vv= vv(els,:);
0304    end
0305 
0306 function imdl= set_basic( shape_str, type, vals )
0307    imdl = mk_common_model( shape_str, 16 );
0308    imdl.RtR_prior= @prior_noser;
0309    imdl.prior_noser.exponent= .5;
0310    imdl.solve= @np_inv_solve;
0311    imdl.hyperparameter.value= 1e-1;
0312 
0313    if nargin==1; return; end
0314 
0315    switch type
0316      case 'temporal'
0317         time_steps= vals(1); 
0318         time_weight= vals(2);
0319 
0320         imdl.RtR_prior= @prior_time_smooth;
0321         imdl.prior_time_smooth.space_prior= @prior_noser;
0322         imdl.prior_time_smooth.time_weight= time_weight;
0323         imdl.prior_time_smooth.time_steps=  time_steps;
0324 
0325         imdl.solve= @inv_solve_time_prior;
0326         imdl.inv_solve_time_prior.time_steps=   time_steps;
0327         imdl.hyperparameter.value= 1e-1;
0328 
0329      otherwise
0330         error('dont understand parameter');
0331    end
0332 
0333 %  Add noise to a data set, snr_mult is # time SNR
0334 function [vi,vh]= add_noise(vi,vh,snr_mult)
0335    ll= size(vi,2);
0336    snr= norm(vi)/norm(vi - vh(:,ones(1,ll)));
0337    rand('seed',50);
0338    vh= vh+snr*snr_mult*randn(size(vh));
0339    vi= vi+snr*snr_mult*randn(size(vi));