Description of aa_fwd

0001 function param = aa_fwd_parameters( fwd_model )
0002 % AA_FWD_PARAMETERS: data= aa_fwd_solve( fwd_model, image)
0003 % Extract parameters from a 'fwd_model' struct which are
0004 % appropriate for Andy Adler's EIT code
0005 %   param.n_elem     => number of elements
0006 %   param.n_elec     => number of electrodes
0007 %   param.n_node     => number of nodes (vertices)
0008 %   param.n_stim     => number of current stimulation patterns
0009 %   param.n_elec     => number of electrodes
0010 %   param.n_dims     => dimentions (2= 2D, 3=3D)
0011 %   param.n_meas     => number of measurements (total)
0012 %   param.boundary   => FEM boundary
0013 %   param.NODE       => vertex matrix
0014 %   param.ELEM       => connection matrix
0015 %   param.QQ         => Current into each NODE
0016 %   param.VOLUME     => Volume (or area) of each element
0017 %   param.normalize  => difference measurements normalized?
0018 %   param.N2E        => Node to electrode converter
0019 %
0020 % If the stimulation patterns has a 'interior_sources' field,
0021 %   the node current QQ, is set to this value for this stimulation.
0022 
0023 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: aa_fwd_parameters.html 2819 2011-09-07 16:43:11Z aadler $
0025 
0026 param = eidors_obj('get-cache', fwd_model, 'aa_fwd_parameters');
0027 
0028 if ~isempty(param)
0029    eidors_msg('aa_fwd_parameters: using cached value', 3);
0030    return
0031 end
0032 
0033 param = calc_param( fwd_model );
0034 
0035 eidors_obj('set-cache', fwd_model, 'aa_fwd_parameters', param);
0036 eidors_msg('aa_fwd_parameters: setting cached value', 3);
0037 
0038 % perform actual parameter calculation
0039 function pp= calc_param( fwd_model );
0040 
0041 pp.NODE= fwd_model.nodes';
0042 pp.ELEM= fwd_model.elems';
0043 
0044 n= size(pp.NODE,2);        %NODEs
0045 d= size(pp.ELEM,1);        %dimentions+1
0046 e= size(pp.ELEM,2);        %ELEMents
0047 try
0048    p = length(fwd_model.stimulation );
0049 catch 
0050    p = 0;
0051 end
0052 try
0053    n_elec= length( fwd_model.electrode );
0054 catch
0055    n_elec= 0;
0056 end
0057 
0058 % calculate element volume and surface area
0059 pp.VOLUME=zeros(e,1);
0060 ones_d = ones(1,d);
0061 d1fac = prod( 1:d-1 );
0062 if d > size(pp.NODE,1)
0063    for i=1:e
0064        this_elem = pp.NODE(:,pp.ELEM(:,i)); 
0065        pp.VOLUME(i)= abs(det([ones_d;this_elem])) / d1fac;
0066    end
0067 elseif d == 3 % 3D nodes in 2D mesh
0068    for i=1:e
0069        this_elem = pp.NODE(:,pp.ELEM(:,i)); 
0070        d12= det([ones_d;this_elem([1,2],:)])^2;
0071        d13= det([ones_d;this_elem([1,3],:)])^2;
0072        d23= det([ones_d;this_elem([2,3],:)])^2;
0073        pp.VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0074    end
0075 elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0076    for i=1:e
0077        this_elem = pp.NODE(:,pp.ELEM(:,i)); 
0078        d12= det([ones_d;this_elem([1],:)])^2;
0079        d13= det([ones_d;this_elem([2],:)])^2;
0080        d23= det([ones_d;this_elem([3],:)])^2;
0081        pp.VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0082    end
0083 else
0084    error('mesh size not understood when calculating volumes')
0085 end
0086 
0087 if isfield(fwd_model,'boundary')
0088     bdy = double( fwd_model.boundary ); % double because of stupid matlab bugs
0089 else
0090     bdy = find_boundary(fwd_model.elems);
0091 end
0092 
0093 % Matrix to convert Nodes to Electrodes
0094 % Complete electrode model for all electrodes
0095 %  N2E = sparse(1:n_elec, n+ (1:n_elec), 1, n_elec, n+n_elec);
0096 %  pp.QQ= sparse(n+n_elec,p);
0097 
0098 cem_electrodes= 0; % num electrodes part of Compl. Elec Model
0099 N2E = sparse(n_elec, n+n_elec);
0100 pp.QQ= sparse(n+n_elec,p);
0101 
0102 for i=1:n_elec
0103     elec_nodes = fwd_model.electrode(i).nodes;
0104     if length(elec_nodes) ==1 % point electrode (maybe inside body)
0105        N2E(i, elec_nodes) = 1;
0106     elseif length(elec_nodes) ==0
0107        error('zero length electrode specified');
0108     else
0109        bdy_idx= find_electrode_bdy( bdy, [], elec_nodes);
0110 
0111        if ~isempty(bdy_idx) % CEM electrode
0112           cem_electrodes = cem_electrodes+1;
0113           N2E(i, n+cem_electrodes) =1;
0114        else % point electrodes
0115             % FIXME: make current defs between point electrodes and CEMs compatible
0116           [bdy_idx,srf_area]= find_electrode_bdy( fwd_model.boundary, ...
0117                          fwd_model.nodes, elec_nodes);
0118           N2E(i, elec_nodes) = srf_area/sum(srf_area);
0119        end
0120     end
0121 end
0122 N2E = N2E(:, 1:(n+cem_electrodes));
0123 pp.QQ= pp.QQ(1:(n+cem_electrodes),:);
0124 
0125 
0126 n_meas= 0; % sum total number of measurements
0127 
0128 if p>0
0129    stim = fwd_model.stimulation;
0130 end
0131 
0132 for i=1:p
0133     src= 0;
0134     try;  src = src +  N2E'* stim(i).stim_pattern; end
0135     try;  src = src +  stim(i).interior_sources;   end
0136     if all(size(src) == [1,1]) && src==0
0137        error('no stim_patterns or interior_sources provided for pattern #%d',i);
0138     end
0139     
0140     pp.QQ(:,i) = src;
0141     n_meas = n_meas + size(stim(i).meas_pattern,1);
0142 end
0143 
0144 
0145 % pack into a parameter return list
0146 pp.n_elem   = e;
0147 pp.n_elec   = n_elec;
0148 pp.n_node   = n;
0149 pp.n_stim   = p;
0150 pp.n_dims   = d-1;
0151 pp.n_meas   = n_meas;
0152 pp.N2E      = N2E;
0153 pp.boundary = bdy;
0154 
0155 if isfield(fwd_model,'normalize_measurements')
0156    pp.normalize = fwd_model.normalize_measurements;
0157 else
0158    pp.normalize = 0;
0159 end
aa_fwd_parameters

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
