backproj_solve

PURPOSE ^

BACKPROJ_SOLVE inverse solver using backprojection

SYNOPSIS ^

function img= backproj_solve( inv_model, data1, data2)

DESCRIPTION ^

 BACKPROJ_SOLVE 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= backproj_solve( 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.backproj_solve.type = 'naive' (DEFAULT)
    use naive (unfiltered algorithm)
 inv_model.backproj_solve.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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function img= backproj_solve( inv_model, data1, data2)
0002 % BACKPROJ_SOLVE 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= backproj_solve( 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.backproj_solve.type = 'naive' (DEFAULT)
0017 %    use naive (unfiltered algorithm)
0018 % inv_model.backproj_solve.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: backproj_solve.html 2819 2011-09-07 16:43:11Z aadler $
0029 
0030 try
0031    type= inv_model.backproj_solve.type;
0032 catch
0033    type= 'naive';
0034 end
0035 
0036 fwd_model= inv_model.fwd_model;
0037 
0038 % The one_step reconstruction matrix is cached
0039 one_step_inv = eidors_obj('get-cache', inv_model, 'backproj_solve');
0040 if ~isempty(one_step_inv)
0041     eidors_msg('backproj_solve: using cached value', 2);
0042 else
0043     Jbp = calc_backprojection_mask( fwd_model, type );
0044 
0045     one_step_inv= Jbp;
0046 
0047     eidors_obj('set-cache', inv_model, 'backproj_solve', one_step_inv);
0048     eidors_msg('backproj_solve: setting cached value', 2);
0049 end
0050 
0051 dv = calc_difference_data( data1, data2, inv_model.fwd_model);
0052 
0053 sol = one_step_inv * dv;
0054 
0055 % create a data structure to return
0056 img.name= 'solved by backproj_solve';
0057 img.elem_data = sol;
0058 img.fwd_model= fwd_model;
0059 
0060 function Jbp = calc_backprojection_mask( fmdl , type);
0061    % create homog image
0062    himg= eidors_obj('image','', 'fwd_model',fmdl, ...
0063                     'elem_data',ones(size(fmdl.elems,1),1) );
0064    node_v= calc_all_node_voltages( himg );
0065 
0066    pp= aa_fwd_parameters( fmdl );
0067 
0068    elem_v= reshape(node_v( fmdl.elems',:),pp.n_dims+1,[],pp.n_elec);
0069    elem_v= squeeze(mean(elem_v,1));
0070 
0071    meas_v= pp.N2E * node_v;
0072 
0073    Jbp = zeros(pp.n_elem, sum(fmdl.meas_select) );
0074 
0075    idx= 1;
0076    % This will only work for stim_patterns with bipolar injection
0077    for i = 1:pp.n_stim;
0078      meas_pat_i= fmdl.stimulation(i).meas_pattern;
0079      elem_vi= elem_v(:,i);
0080      for j= 1:size(meas_pat_i,1); 
0081         idx_pl= find(meas_pat_i(j,:)>0);
0082         idx_mi= find(meas_pat_i(j,:)<0);
0083         Jbp_idx  = (elem_vi <= meas_v(idx_pl,i)) & ...
0084                    (elem_vi >  meas_v(idx_mi,i));
0085         Jbp(:,idx) = - Jbp_idx;
0086         idx= idx+1;
0087      end
0088    end
0089 
0090    if strcmp(type,'naive')
0091 % do nothing
0092    elseif strcmp(type,'simple_filter')
0093       xe=mean(reshape(fmdl.nodes(fmdl.elems',1),pp.n_dims+1,pp.n_elem));
0094       ye=mean(reshape(fmdl.nodes(fmdl.elems',2),pp.n_dims+1,pp.n_elem));
0095       filt=sqrt(1-xe.^2-ye.^2);
0096       Jbp= Jbp.*( filt'*ones(1,size(Jbp,2)) );
0097    else
0098       error(['dont know what to do with filter type=',type]);
0099    end
0100

Generated on Tue 09-Aug-2011 11:38:31 by m2html © 2005