inv_solve

 INV_SOLVE: calculate imag from an inv_model and data
 
 inv_solve can be called as
     img= inv_solve( inv_model, data1, data2)
   if inv_model.reconst_type = 'difference'
 or
     img= inv_solve( inv_model, data )
   if inv_model.reconst_type = 'static'

   if inv_model.reconst_to = 'nodes' then output
      img.node_data has output data
   else   reconst_to = 'elems' (DEFAULT) then output to
      img.elem_data has output data

 in each case it will call the inv_model.solve

 data      is a measurement data structure
 inv_model is a inv_model structure
 img       is an image structure
           or a vector of images is data1 or data2 are vectors

 For difference EIT:
 data1      => difference data at earlier time (ie homogeneous)
 data2      => difference data at later time   (ie inhomogeneous)

 data can be:
   - an EIDORS data object

   - an M x S matrix, where M is the total number
         of measurements expected by inv_model

   - an M x S matrix, where M is n_elec^2
        when not using data from current injection
        electrodes, it is common to be given a full
        measurement set.  For example, 16 electrodes give
        208 measures, but 256 measure sets are common.
        Data will be selected based on fwd_model.meas_select.

 If S > 1 for both data1 and data2 then the matrix sizes must be equal

 Parameters:
   inv_model.inv_solve.select_parameters: indices of parameters to return
                         DEFAULT: return all paramteres
  Scale solution (to correct for amplitude or other defects)
   inv_model.inv_solve.scale_solution.offset
   inv_model.inv_solve.scale_solution.scale

0001 function img = inv_solve( inv_model, data1, data2)
0002 % INV_SOLVE: calculate imag from an inv_model and data
0003 %
0004 % inv_solve can be called as
0005 %     img= inv_solve( inv_model, data1, data2)
0006 %   if inv_model.reconst_type = 'difference'
0007 % or
0008 %     img= inv_solve( inv_model, data )
0009 %   if inv_model.reconst_type = 'static'
0010 %
0011 %   if inv_model.reconst_to = 'nodes' then output
0012 %      img.node_data has output data
0013 %   else   reconst_to = 'elems' (DEFAULT) then output to
0014 %      img.elem_data has output data
0015 %
0016 % in each case it will call the inv_model.solve
0017 %
0018 % data      is a measurement data structure
0019 % inv_model is a inv_model structure
0020 % img       is an image structure
0021 %           or a vector of images is data1 or data2 are vectors
0022 %
0023 % For difference EIT:
0024 % data1      => difference data at earlier time (ie homogeneous)
0025 % data2      => difference data at later time   (ie inhomogeneous)
0026 %
0027 % data can be:
0028 %   - an EIDORS data object
0029 %
0030 %   - an M x S matrix, where M is the total number
0031 %         of measurements expected by inv_model
0032 %
0033 %   - an M x S matrix, where M is n_elec^2
0034 %        when not using data from current injection
0035 %        electrodes, it is common to be given a full
0036 %        measurement set.  For example, 16 electrodes give
0037 %        208 measures, but 256 measure sets are common.
0038 %        Data will be selected based on fwd_model.meas_select.
0039 %
0040 % If S > 1 for both data1 and data2 then the matrix sizes must be equal
0041 %
0042 % Parameters:
0043 %   inv_model.inv_solve.select_parameters: indices of parameters to return
0044 %                         DEFAULT: return all paramteres
0045 %  Scale solution (to correct for amplitude or other defects)
0046 %   inv_model.inv_solve.scale_solution.offset
0047 %   inv_model.inv_solve.scale_solution.scale
0048 
0049 % (C) 2005-2010 Andy Adler. License: GPL version 2 or version 3
0050 % $Id: inv_solve.html 2819 2011-09-07 16:43:11Z aadler $
0051 
0052 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0053 
0054 inv_model= eidors_model_params( inv_model );
0055 opts = parse_parameters( inv_model );
0056 
0057 eidors_msg('inv_solve: %s', inv_model.solve,2);
0058 
0059 if opts.abs_solve
0060    if nargin~=2;
0061       error('only one data set is allowed for a static reconstruction');
0062    end
0063 
0064    imgc= feval( inv_model.solve, inv_model, ...
0065                filt_data(inv_model,data1) );
0066 else
0067    if nargin~=3;
0068       error('two data sets are required for a difference reconstruction');
0069    end
0070 
0071    % expand data sets if one is provided that is longer
0072    data_width= max(num_frames(data1), num_frames(data2));
0073 
0074    fdata1 = filt_data( inv_model, data1, data_width );
0075    fdata2 = filt_data( inv_model, data2, data_width );
0076 
0077    % TODO: Check if solver can handle being called with multiple data
0078    imgc= feval( inv_model.solve, inv_model, fdata1, fdata2);
0079 end
0080 
0081 img = eidors_obj('image', imgc );
0082 % If we reconstruct with a different 'rec_model' then
0083 %  put this into the img
0084 if isfield(inv_model,'rec_model')
0085    img.fwd_model= inv_model.rec_model;
0086 end
0087 
0088 
0089 if opts.select_parameters;
0090    img.elem_data = img.elem_data( opts.select_parameters, :);
0091 end;
0092 
0093 if ~opts.reconst_to_elems
0094   img.node_data= img.elem_data;
0095   img = rmfield(img,'elem_data');
0096 end
0097 
0098 % Scale if required
0099 try; img.elem_data = opts.offset + opts.scale * img.elem_data; end
0100 try; img.node_data = opts.offset + opts.scale * img.node_data; end
0101 
0102 function opts = parse_parameters( imdl );
0103    if  strcmp(imdl.reconst_type,'static') || ...
0104        strcmp(imdl.reconst_type,'absolute')
0105       opts.abs_solve = 1;
0106    elseif strcmp(imdl.reconst_type,'difference')
0107       opts.abs_solve = 0;
0108    else
0109       error('inv_model.reconst_type (%s) not understood', imdl.reconst_type); 
0110    end
0111 
0112    opts.select_parameters = [];
0113    try
0114       opts.select_parameters = imdl.inv_solve.select_parameters;
0115    end;
0116 
0117    opts.reconst_to_elems = 1;
0118    try; if strcmp( imdl.reconst_to, 'nodes' )
0119       opts.reconst_to_elems = 0;
0120    end; end
0121    
0122    opts.scale  = 1;
0123    try; opts.scale = imdl.inv_solve.scale_solution.scale; end
0124 
0125    opts.offset = 0;
0126    try; opts.offset = imdl.inv_solve.scale_solution.offset; end
0127  
0128 
0129 % TODO: this code really needs to be cleaned, but not before eidors 3.4
0130 function nf= num_frames(d0)
0131    if isnumeric( d0 )
0132       nf= size(d0,2);
0133    elseif d0(1).type == 'data';
0134       nf= size( horzcat( d0(:).meas ), 2);
0135    else
0136       error('Problem calculating number of frames. Expecting numeric or data object');
0137    end
0138    
0139 % test for existance of meas_select and filter data
0140 function d2= filt_data(inv_model, d0, data_width )
0141    if ~isnumeric( d0 )
0142        % we probably have a 'data' object
0143 
0144        d1 = [];
0145        for i=1:length(d0)
0146           if strcmp( d0(i).type, 'data' )
0147               d1 = [d1, d0(i).meas];
0148           else
0149               error('expecting an object of type data');
0150           end
0151        end
0152 
0153    else
0154       % we have a matrix of data. Hope for the best
0155       d1 = d0;
0156    end
0157 
0158    d1= double(d1); % ensure we can do math on our object
0159 
0160    if isfield(inv_model.fwd_model,'meas_select') && ...
0161      ~isempty(inv_model.fwd_model.meas_select)
0162       % we have a meas_select parameter that isn []
0163 
0164       meas_select= inv_model.fwd_model.meas_select;
0165       if     size(d1,1) == length(meas_select)
0166          d2= d1(meas_select,:);
0167       elseif size(d1,1) == sum(meas_select==1)
0168          d2= d1;
0169       else
0170          error('inconsistent difference data: (%d ~= %d). Maybe check fwd_model.meas_select',  ...
0171                d2_width, data_width);
0172       end
0173    else
0174       d2= d1;
0175    end
0176 
0177    if nargin==3 % expand to data width
0178       d2_width= size(d2,2);
0179       if d2_width == data_width
0180          % ok
0181       elseif d2_width == 1
0182          d2= d2(:,ones(1,data_width));
0183       else
0184          error('inconsistent difference data: (%d ~= %d)',  ...
0185                d2_width, data_width);
0186       end
0187    end
0188 
0189 % Test code
0190 function do_unit_test
0191    k=0; N=5; nd = 5;
0192 
0193    imdl = mk_common_model('d2c2',16);
0194    imdl = select_imdl( imdl, {'Choose NF=2.5'});
0195    mvx = linspace(-0.8,0.8,nd);
0196    [vh,vi] = simulate_movement(mk_image(imdl), [mvx;0*mvx;0.05+0*mvx]);
0197 
0198    img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0199    k=k+1; subplot(N,1,k); show_slices(img);
0200    unit_test_cmp('inv_solve: 1a', mean(img.elem_data), 3e-3, 1e-3);
0201    unit_test_cmp('inv_solve: 1b',  std(img.elem_data), 1.5e-2, 1e-2);
0202 
0203    vhm= eidors_obj('data','nn','meas',vh);
0204    img= inv_solve(imdl,vhm,vi);
0205    unit_test_cmp('inv_solve: 2a', mean(img.elem_data), 3e-3, 1e-3);
0206    unit_test_cmp('inv_solve: 2b',  std(img.elem_data), 1.5e-2, 1e-2);
0207 
0208    img= inv_solve(imdl,vh*ones(1,nd),vi);
0209    unit_test_cmp('inv_solve: 3a', mean(img.elem_data), 3e-3, 1e-3);
0210    unit_test_cmp('inv_solve: 3b',  std(img.elem_data), 1.5e-2, 1e-2);
0211 
0212    vim= eidors_obj('data','nn','meas',vi);
0213    img= inv_solve(imdl,vhm,vim);
0214    unit_test_cmp('inv_solve: 4a', mean(img.elem_data), 3e-3, 1e-3);
0215    unit_test_cmp('inv_solve: 4b',  std(img.elem_data), 1.5e-2, 1e-2);
0216 
0217    vhm= eidors_obj('data','nn','meas',vh*ones(1,nd));
0218    img= inv_solve(imdl,vhm,vi);
0219    unit_test_cmp('inv_solve: 5a', mean(img.elem_data), 3e-3, 1e-3);
0220    unit_test_cmp('inv_solve: 5b',  std(img.elem_data), 1.5e-2, 1e-2);
0221 
0222    vhm(1)= eidors_obj('data','nn','meas',vh*ones(1,2));
0223    vhm(2)= eidors_obj('data','nn','meas',vh*ones(1,nd-2));
0224    img= inv_solve(imdl,vhm,vi);
0225    unit_test_cmp('inv_solve: 6a', mean(img.elem_data), 3e-3, 1e-3);
0226    unit_test_cmp('inv_solve: 6b',  std(img.elem_data), 1.5e-2, 1e-2);
0227 
0228    vim(1)= eidors_obj('data','nn','meas',vi(:,1:3));
0229    vim(2)= eidors_obj('data','nn','meas',vi(:,4:end));
0230    img= inv_solve(imdl,vhm,vim);
0231    unit_test_cmp('inv_solve: 7a', mean(img.elem_data), 3e-3, 1e-3);
0232    unit_test_cmp('inv_solve: 7b',  std(img.elem_data), 1.5e-2, 1e-2);
0233 
0234    im2 = imdl; im2.inv_solve.select_parameters = 1:5;
0235    img= inv_solve(im2,vh,vi);
0236    unit_test_cmp('inv_solve: 10', size(img.elem_data), [5,nd]);
0237 
0238 
0239    im2 = imdl;
0240    im2.inv_solve.scale_solution.offset = 1;
0241    im2.inv_solve.scale_solution.scale = 2;
0242    img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0243    unit_test_cmp('inv_solve: 20a', mean(img.elem_data), 1.006, 2e-3);
0244    unit_test_cmp('inv_solve: 20b',  std(img.elem_data), 3e-2, 1e-2);
0245    
0246    im2.inv_solve.scale_solution.offset = 0;
0247    d = interp_mesh( imdl.fwd_model); d= sqrt(sum(d.^2,2));
0248    im2.inv_solve.scale_solution.scale = spdiags(1-d,0,length(d),length(d));
0249    img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0250    k=k+1; subplot(N,1,k); show_slices(img);
0251 
0252    im2 = select_imdl(imdl, {'Nodal GN dif'} );
0253    img= inv_solve(im2,vh,vi);
0254    unit_test_cmp('inv_solve: 30a', mean(img.node_data), 3e-3, 1e-3);
0255    unit_test_cmp('inv_solve: 30b',  std(img.node_data), 1.5e-2, 1e-2);
0256 
0257    im2.inv_solve.scale_solution.offset = 1;
0258    im2.inv_solve.scale_solution.scale = 2;
0259    img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0260    unit_test_cmp('inv_solve: 31a', mean(img.node_data), 1.006, 2e-3);
0261    unit_test_cmp('inv_solve: 31b',  std(img.node_data), 3e-2, 1.5e-2);
0262 
0263    im2 = select_imdl(imdl, {'Basic GN abs'} );
0264    img= inv_solve(im2,vi(:,1));
0265    unit_test_cmp('inv_solve: 40a', mean(img.elem_data), 1.004, 1e-3);
0266    unit_test_cmp('inv_solve: 40b',  std(img.elem_data), 1.5e-2, 1e-2);
0267 
0268