inv_solve_backproj

 INV_SOLVE_BACKPROJ inverse solver using backprojection
 NOTE: This is the beginnings of an attempt to reproduce
  the backprojection algorithm implemented in the
  Sheffield MKI EIT system. It is far from complete.

 If you wish to use the actual algorithm, use the
  function "mk_common_gridmdl('backproj')"

 img= inv_solve_backproj( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 inv_model.inv_solve_backproj.type = 'naive' (DEFAULT)
    use naive (unfiltered algorithm)
 inv_model.inv_solve_backproj.type = 'filtered' (NOT IMPLEMENTED YET)
    ref: Barber DC Brown BH, "fast reconstruction of resistance
         images", clin Phys Physiol Mes, pp 47-54, vol 8,sup A,1987

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= inv_solve_backproj( inv_model, data1, data2)
0002 % INV_SOLVE_BACKPROJ inverse solver using backprojection
0003 % NOTE: This is the beginnings of an attempt to reproduce
0004 %  the backprojection algorithm implemented in the
0005 %  Sheffield MKI EIT system. It is far from complete.
0006 %
0007 % If you wish to use the actual algorithm, use the
0008 %  function "mk_common_gridmdl('backproj')"
0009 %
0010 % img= inv_solve_backproj( inv_model, data1, data2)
0011 % img        => output image (or vector of images)
0012 % inv_model  => inverse model struct
0013 % data1      => differential data at earlier time
0014 % data2      => differential data at later time
0015 %
0016 % inv_model.inv_solve_backproj.type = 'naive' (DEFAULT)
0017 %    use naive (unfiltered algorithm)
0018 % inv_model.inv_solve_backproj.type = 'filtered' (NOT IMPLEMENTED YET)
0019 %    ref: Barber DC Brown BH, "fast reconstruction of resistance
0020 %         images", clin Phys Physiol Mes, pp 47-54, vol 8,sup A,1987
0021 %
0022 % both data1 and data2 may be matrices (MxT) each of
0023 %  M measurements at T times
0024 % if either data1 or data2 is a vector, then it is expanded
0025 %  to be the same size matrix
0026 
0027 % (C) 2007 Andy Adler. License: GPL version 2 or version 3
0028 % $Id: inv_solve_backproj.m 5025 2015-05-26 06:55:43Z bgrychtol $
0029 
0030 try
0031    type= inv_model.inv_solve_backproj.type;
0032 catch
0033    type= 'naive';
0034 end
0035 
0036 fwd_model= inv_model.fwd_model;
0037 
0038 one_step_inv = eidors_cache(@calc_backprojection_mask, {fwd_model, type});
0039 
0040 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0041 
0042 sol = one_step_inv * dv;
0043 
0044 % create a data structure to return
0045 img.name= 'solved by inv_solve_backproj';
0046 img.elem_data = sol;
0047 img.fwd_model= fwd_model;
0048 
0049 function Jbp = calc_backprojection_mask( fmdl , type);
0050    % create homog image
0051    himg= eidors_obj('image','', 'fwd_model',fmdl, ...
0052                     'elem_data',ones(size(fmdl.elems,1),1) );
0053    node_v= calc_all_node_voltages( himg );
0054 
0055    pp= fwd_model_parameters( fmdl );
0056 
0057    elem_v= reshape(node_v( fmdl.elems',:),pp.n_dims+1,[],pp.n_elec);
0058    elem_v= squeeze(mean(elem_v,1));
0059 
0060    meas_v= pp.N2E * node_v;
0061 
0062    Jbp = zeros(pp.n_elem, sum(fmdl.meas_select) );
0063 
0064    idx= 1;
0065    % This will only work for stim_patterns with bipolar injection
0066    for i = 1:pp.n_stim;
0067      meas_pat_i= fmdl.stimulation(i).meas_pattern;
0068      elem_vi= elem_v(:,i);
0069      for j= 1:size(meas_pat_i,1); 
0070         idx_pl= find(meas_pat_i(j,:)>0);
0071         idx_mi= find(meas_pat_i(j,:)<0);
0072         Jbp_idx  = (elem_vi <= meas_v(idx_pl,i)) & ...
0073                    (elem_vi >  meas_v(idx_mi,i));
0074         Jbp(:,idx) = - Jbp_idx;
0075         idx= idx+1;
0076      end
0077    end
0078 
0079    if strcmp(type,'naive')
0080 % do nothing
0081    elseif strcmp(type,'simple_filter')
0082       xe=mean(reshape(fmdl.nodes(fmdl.elems',1),pp.n_dims+1,pp.n_elem));
0083       ye=mean(reshape(fmdl.nodes(fmdl.elems',2),pp.n_dims+1,pp.n_elem));
0084       filt=sqrt(1-xe.^2-ye.^2);
0085       Jbp= Jbp.*( filt'*ones(1,size(Jbp,2)) );
0086    else
0087       error(['dont know what to do with filter type=',type]);
0088    end
0089