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
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