prior_covar

 PRIOR_COVAR image prior with distance-based interelement covar
 This is a simplification of prior_exponential_covar.m
 Reg= prior_covar( inv_model )
 Reg        => output regularization 
 inv_model  => inverse model struct
 P_type--prior type
 1: elements are globally correlated
 2: elements within/without electrode rings are correlated to elements in same region.
 3: only elements within electrode rings are correlated.

0001 function Reg= prior_covar( inv_model )
0002 % PRIOR_COVAR image prior with distance-based interelement covar
0003 % This is a simplification of prior_exponential_covar.m
0004 % Reg= prior_covar( inv_model )
0005 % Reg        => output regularization
0006 % inv_model  => inverse model struct
0007 % P_type--prior type
0008 % 1: elements are globally correlated
0009 % 2: elements within/without electrode rings are correlated to elements in same region.
0010 % 3: only elements within electrode rings are correlated.
0011 
0012 % (C) 2007, Tao Dai and Andy Adler. Licenced under the GPL Version 2
0013 % $Id: prior_covar.m 5112 2015-06-14 13:00:41Z aadler $
0014 
0015 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0016 
0017 params = {inv_model.fwd_model, inv_model.fourD_prior.P_type};
0018 Reg = eidors_cache(@prior_covar_calc, params, 'prior_covar');
0019 
0020 
0021 function Reg = prior_covar_calc(ff, P_type);
0022    % get average x,y,z of each element
0023    nel = size(ff.elems,1);
0024    eta = .1;%attenuation factor. eta is large when elems're spatially highly correlated
0025 
0026    dist = zeros(nel);
0027 
0028    z1 = ff.nodes(ff.electrode(1).nodes,3);%upper electrode ring
0029    z2 = ff.nodes(ff.electrode(end).nodes,3);%lower electrode ring
0030    zmax = max(z1,z2);
0031    zmin = min(z1,z2);
0032 
0033    for dim = 1: size(ff.nodes,2);
0034        coords = reshape(ff.nodes(ff.elems,dim),[size(ff.elems)]);%coord of four vertex of each elem
0035        m_coords = mean( coords,2);%calc center coord
0036        difm = m_coords*ones(1,nel);
0037        difm = difm - difm';
0038        dist= dist + difm.^2;%
0039    end
0040    dist = sqrt(dist);%elements distance matrix
0041    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0042    dist = dist/max(dist(:));
0043 
0044    Reg = exp(-dist ./ eta);% 3-D elements correlations matrix.
0045        if dim == 3
0046 %          temp = double((m_coords<max(z1,z2))&(m_coords>min(z1,z2)));%regions of interests
0047 %          H = temp*temp';
0048            above = m_coords>zmax;
0049            below = m_coords<zmin;
0050            switch P_type %Prior type
0051                case 2 % elements of ROI are not correlated with elements of non_ROI
0052 %          temp = find( ~above & ~below);
0053 %          Reg(temp,temp) = 1;
0054 %                  temp = double(m_coords>max(z1,z2));
0055 %                  H = H + temp*temp';
0056 %                  above = find(above);
0057 %                  H(above,above) = 1;
0058 %                  temp = double(m_coords<min(z1,z2));
0059 %                  H = H + temp*temp';
0060 %                  temp = m_coords<zmin;
0061 %                  below = find(below);
0062 %                  H(below,below) = 1;
0063 %                  H = eta*(H+1e-6);
0064                    Reg(above,~above) = 0;
0065                    Reg(below,~below) = 0;
0066                case 3 % elements are only correlated within ROI, non-ROI elements are highly attenuated
0067 %                  H = eta*(H+1e-6);
0068                    temp = find( above | below);
0069                    Reg(temp,temp) = 0;
0070                otherwise
0071                    error('Prior type (%d) has not yet been defined',P_type);
0072            end
0073        end
0074 
0075 
0076 function do_unit_test
0077    imdl = mk_common_model('b3cr',[16,3]);
0078    imdl.fourD_prior.P_type = 2;
0079 %tic;
0080    Reg = prior_covar( imdl );
0081    unit_test_cmp('#1:', diag(Reg), ones(size(Reg,1),1));
0082    unit_test_cmp('#2:', Reg(102:103,100:103), [ ...
0083    0.011818547266014   0.367556336897157   1.000000000000000   0.396442626219240; ...
0084    0.008985742075647   0.160077522568152   0.396442626219240   1.000000000000000], 1e-6);
0085 
0086 %toc;
0087 
0088 %tic;
0089    imdl.fourD_prior.P_type = 3;
0090    Reg = prior_covar( imdl );
0091    test = zeros(size(Reg,1),1); test(769:6912) = 1;
0092    unit_test_cmp('#1:', diag(Reg), test);
0093    unit_test_cmp('#2:', Reg(802:803,800:803), [ ...
0094    0.086725242542203   0.132995997778770   1.000000000000000   0.417609404846561; ...
0095    0.157824974662458   0.293750928546438   0.417609404846561   1.000000000000000], 1e-10);
0096 %toc;