fwd_solve

 FWD_SOLVE: calculate data from a fwd_model object and an image
 
 fwd_solve can be called as
    data= fwd_solve( img)
 or (deprecated)
    data= fwd_solve( fwd_model, img)

 in each case it will call the fwd_model.solve
                        or img.fwd_model.solve method

 For reconstructions on dual meshes, the interpolation matrix
    is defined as fwd_model.coarse2fine. If required, this takes
    coarse2fine * x_coarse = x_fine

 data      is a measurement data structure
 fwd_model is a fwd_model structure
 img       is an img structure

 Options: (not available on all solvers)
    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
    img.fwd_solve.get_elec_curr= 1 (data.elec_curr = current on electrodes)

 Parameters:
    fwd_model.background = constant conductivity offset added to elem_data
    fwd_model.coarse2fine = linear mapping between img.elem_data and model parameters
    img.params_mapping = function mapping img.elem_data to model parameters

0001 function data = fwd_solve(fwd_model, img)
0002 % FWD_SOLVE: calculate data from a fwd_model object and an image
0003 %
0004 % fwd_solve can be called as
0005 %    data= fwd_solve( img)
0006 % or (deprecated)
0007 %    data= fwd_solve( fwd_model, img)
0008 %
0009 % in each case it will call the fwd_model.solve
0010 %                        or img.fwd_model.solve method
0011 %
0012 % For reconstructions on dual meshes, the interpolation matrix
0013 %    is defined as fwd_model.coarse2fine. If required, this takes
0014 %    coarse2fine * x_coarse = x_fine
0015 %
0016 % data      is a measurement data structure
0017 % fwd_model is a fwd_model structure
0018 % img       is an img structure
0019 %
0020 % Options: (not available on all solvers)
0021 %    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
0022 %    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
0023 %    img.fwd_solve.get_elec_curr= 1 (data.elec_curr = current on electrodes)
0024 %
0025 % Parameters:
0026 %    fwd_model.background = constant conductivity offset added to elem_data
0027 %    fwd_model.coarse2fine = linear mapping between img.elem_data and model parameters
0028 %    img.params_mapping = function mapping img.elem_data to model parameters
0029 %
0030 
0031 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0032 % $Id: fwd_solve.m 6307 2022-04-20 17:45:43Z aadler $
0033 
0034 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end 
0035 
0036 if nargin == 1
0037    img= fwd_model;
0038 else
0039    warning('EIDORS:DeprecatedInterface', ...
0040       ['Calling FWD_SOLVE with two arguments is deprecated and will cause' ...
0041        ' an error in a future version. First argument ignored.']);
0042 end
0043 ws = warning('query','EIDORS:DeprecatedInterface');
0044 warning off EIDORS:DeprecatedInterface
0045 
0046 
0047 fwd_model= img.fwd_model;
0048 
0049 
0050 fwd_model= prepare_model( fwd_model );
0051 
0052 % TODO: This should be handled by the data_mapper
0053 if isfield(img,'params_mapping')
0054 %     fwd_model data is provided using a mapping function
0055     mapping_function= img.params_mapping.function;
0056     img= feval(mapping_function,img);
0057 end
0058 if isfield(fwd_model,'coarse2fine') && isfield(img,'elem_data')
0059    c2f= fwd_model.coarse2fine;
0060    if size(img.elem_data,1)==size(c2f,2)
0061 %     fwd_model data is provided on coarse mesh
0062       img.elem_data = c2f * img.elem_data; 
0063 
0064       if isfield(fwd_model,'background')
0065           img.elem_data = img.elem_data + fwd_model.background; 
0066       end
0067    end
0068 end
0069 
0070 if ~isfield(fwd_model, 'electrode')
0071    error('EIDORS: attempting to solve on model without electrodes');
0072 end
0073 if ~isfield(fwd_model, 'stimulation')
0074    error('EIDORS: attempting to solve on model without stimulation patterns');
0075 end
0076 
0077 solver = fwd_model.solve;
0078 if ischar(solver)
0079     solver = str2func(solver);
0080 end
0081 
0082 copt.fstr = 'fwd_solve';
0083 n_frames = size(img.elem_data,2);
0084 for frame = 1:n_frames;
0085    imgn = img; imgn.elem_data = imgn.elem_data(:,frame,:,:);
0086    copt.cache_obj = imgn;
0087    copt.boost_priority = -2; % fmdl evaluations are low priority
0088    tmp = eidors_cache(solver, {imgn}, copt);
0089    data(frame) = eidors_obj('data',tmp);  % create data object
0090 end
0091 
0092 if isa(fwd_model.solve,'function_handle')
0093     solver = func2str(fwd_model.solve);
0094 else
0095     solver = fwd_model.solve;
0096 end
0097 if strcmp(solver,'eidors_default');
0098     solver = eidors_default('get','fwd_solve');
0099 end
0100 if isfield(fwd_model,'measured_quantity') && ~isfield(data,'measured_quantity')
0101    warning('EIDORS:MeasurementQuantityObliviousSolver',...
0102       ['The solver %s did not handle the requested measurement quantity properly.\n'...
0103        'The results may be incorrect. Please check the code to verify.'], ...
0104        solver);
0105 elseif isfield(fwd_model,'measured_quantity') ... 
0106         && isfield(data,'measured_quantity') ...
0107         && ~strcmp(fwd_model.measured_quantity, data.measured_quantity)
0108    error('EIDORS:MeasurementQuantityDisagreement',...
0109        'The solver %s return measurements as %s, while %s was expected.',...
0110        solver, data.measured_quantity, fwd_model.measured_quantity);
0111 end
0112     
0113 warning on EIDORS:DeprecatedInterface
0114 
0115 
0116 function mdl = prepare_model( mdl )
0117 mdl = mdl_normalize(mdl,mdl_normalize(mdl));
0118 if ~isfield(mdl,'elems');
0119     return;
0120 end
0121 
0122 mdl.elems  = double(mdl.elems);
0123 mdl.n_elem = size(mdl.elems,1);
0124 mdl.n_node = size(mdl.nodes,1);
0125 if isfield(mdl,'electrode');
0126     mdl.n_elec = length(mdl.electrode);
0127 else
0128     mdl.n_elec = 0;
0129 end
0130 
0131 function do_unit_test
0132 img = mk_image(mk_common_model('a2c2',4));
0133 % TODO Try other solvers
0134 img.elem_data([1,5,9,10]) = 2;
0135 A = [1,-1,0,0]; a=A.';
0136 B = [0,0,1,-1]; b=B.';
0137 G = 2; 
0138 P = exp(0.1j);
0139 for solv = 1:4; switch solv
0140     case 1; solver = 'fwd_solve_1st_order';
0141     case 2; solver = 'fwd_solve_2p5d_1st_order';
0142     case 3; solver = 'fwd_solve_higher_order';
0143             img.fwd_model.system_mat = @system_mat_higher_order;
0144             img.fwd_model.approx_type = 'tri3';
0145     case 4; break;
0146     end
0147     for test = 1:4; switch test
0148        case 1; str='reciprocity';
0149                mPat = {A,B}; sPat = {b,a};
0150        case 2; str='Gain';
0151                mPat = {A,A*G}; sPat = {b,b/G};
0152        case 3; str='Phase';
0153                mPat = {A,A*P}; sPat = {b,b*P};
0154        case 4; break
0155        end
0156        img.fwd_model.solve = solver;
0157        img.fwd_model.stimulation  = struct( ...
0158           'meas_pattern', mPat, 'stim_pattern', sPat );
0159         vv= getfield(fwd_solve(img),'meas');
0160         unit_test_cmp([solver,':',str],diff(vv),0,10*eps) 
0161     end
0162 end