-
Notifications
You must be signed in to change notification settings - Fork 6
/
bst_berg.m
149 lines (141 loc) · 6.3 KB
/
bst_berg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
function [mu_berg,lam_berg,f] = bst_berg(R,sigma,nmax,J,mu_berg_init,lam_berg_init)
% BST_BERG: Parameter calaculation for EEG multilayer spherical forward model.
%
% USAGE: [mu_berg,lam_berg,f] = bst_berg(R,sigma,nmax,J,mu_berg_init,lam_berg_init)
% [mu_berg,lam_berg] = bst_berg(R,sigma,nmax,J)
% [mu_berg,lam_berg] = bst_berg(R,sigma,nmax)
% [mu_berg,lam_berg] = bst_berg(R,sigma)
%
% DESCTIPTION:
% This function computes the Berg Eccentricity and Magnitude parameters associated with
% a series Approximiation of a Single Dipole in a Multilayer Sphere -by-
% multiple dipoles in a single shell. Zhang: Eq (1i,2i,1i',5i'' and 7i).
%
% This function has been generalized to compute Berg parameters for an approximation
% based on a total of J (user-specified) dipoles (J=3 is recommended)
%
% (Ref: Z. Zhang "A fast method to compute surface potentials generated by dipoles
% within multilayer anisotropic spheres" (Phys. Med. Biol. 40, pp335-349,1995)
%
% INPUTS (Required):
% R : Radii(in meters) of sphere from
% INNERMOST to OUTERMOST NL x 1
% sigma: conductivity from INNERMOST to OUTERMOST NL x 1
%
% INPUTS (Optional):
% nmax : # of terms used in Legendre Expansion used to
% "fit" Berg Parameters (Default: 100) scalar
% J : Number of Berg Dipoles (Default: 3) scalar
% mu_berg_init: User specified initial value for Berg eccentricity
% factors (default values shown below used otherwise) J x 1
% lam_berg_init: User specified initial value for Berg magnitude
% factors (default values shown below used otherwise) J x 1
%
% where: NL = # of sphere layers; J = # of Berg Dipoles
%
% OUTPUTS:
% mu_berg: Computed Value of Berg eccentricity factors J x 1
% lam_berg: Computed Value of Berg magnitude factors J x 1
% f: Legendre Expansion Weights used to fit Berg Parameters nmax x 1
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: John J.Ermer, 1999
% ----------------------------- Script History ---------------------------------
% JJE 5-May-1999 Creation
% JJE 21-Jun-1999 Addded Legendre Expansion weights as optional output (JE)
% JJE 30-Sep-1999 Fixed Logic to pass Initial Berg Parameters (JE)
% SB Nov-2001 Replaced fmins syntax (obsolete in Matlab R12) by fminsearch
% ------------------------------------------------------------------------------
%%% THIS PART CHECKS INPUT PARAMETERS FOR THEIR DIMENSION AND VALIDITY %%%
%
NL = length(R); % # of concentric sphere layer
%
%%% THIS PART CHECKS FOR THE PRESENCE OF OPTIONAL PARAMETERS AND ASSIGNS %%%
%%% DEFAULT VALUES IF NECESSARY %%%
%
if nargin < 3 % Check # terms to see if number of Legendre expansion terms is specified
nmax = 200; % Default # of Legendre Terms
end
%
if nargin < 4 % Check # of Berg Paramters to be generated
J=3; % Default # of Berg Parameters
end
%
if (nargin < 5) % Check if user specified initial value for Berg Paramters
mu_berg_init = (1/(J+2))*(1:1:J); % Default Value for Initial Eccentricity Parameters
lam_berg_init = 0.2*ones(1,J); % Default Value for Initial Magnitude Parameters
else
if length(mu_berg_init)~=length(lam_berg_init)
error('Initial Berg Eccen and Mag parameters are of Different Lengths!!!')
elseif length(mu_berg_init) ~= J
error('Initial Berg Param Size does not match number of dipoles specified!!!')
end
end
%
%%% THIS PART COMPUTES THE WEIGHTS fn ASSOCIATED WITH A LEGENDRE EXPANSION %%%%%%%%%%%%%%%
%%% THE WEIGHTS fn DEPEND ONLY ON THE MULTISPHERE RADII AND CONDUCTIVITY %%%%%%%%%%%%%%%%%
%%% (This portion for computing fn derived from gainp_sph.m by CCH, Aug/20/1995)
%
Re_mag = R(NL); % Radius of outermost layer (Sensor distance from origin)
if NL==1
f=ones(1,nmax);
else
%
for k = 1:NL-1
s(k) = sigma(k)/sigma(k+1);
end
a = Re_mag./R;
ainv = R/Re_mag;
sm1 = s-1;
twonp1 = 2 * (1:nmax) + 1;
twonp1 = twonp1(:);
f = zeros(nmax,1);
%
for n = 1:nmax
np1 = n+1;
Mc = eye(2);
for k = 2:NL-1
Mc = Mc*[n+np1*s(k), np1*sm1(k)*a(k)^twonp1(n);...
n*sm1(k)*ainv(k)^twonp1(n) , np1+n*s(k)];
end
Mc(2,:) = [n*sm1(1)*ainv(1)^twonp1(n) , np1+n*s(1)]*Mc; % Compute only components of interest
Mc = Mc/(twonp1(n))^(NL-1);
f(n) = n/(n*Mc(2,2)+np1*Mc(2,1));
end
end
%
%%%% End of Calculation of weights fn %%%%%%%
%
%%%% THIS PART COMPUTES THE BERG PARAMETERS ASSOCIATED WITH THE SPHERE RADII AND %%%%%%%%%%%%%
%%%% AND CONDUCTIVITY BY MINIMIZING THE FUNCTION 5I" SPECIFIED IN ZHANG %%%%%%%%%%%%%%%%%%%%%%
%
% OPTIONS=zeros(1,18); % Reset all options flags
% OPTIONS(1:3) = [0 1.e-9 1.e-9]; % Set options termination search criteria
% OPTIONS(14) = J*1000; % Set max number of steps
if(version('-release')<13)
OPTIONS = optimset('MaxFunEvals',J*1000,'MaxIter',J*1000,'TolFun',[0 1.e-9 1.e-9]);
else
OPTIONS = optimset('MaxFunEvals',J*1000,'MaxIter',J*1000,'TolFun', 1e-9);
end
[berg_out,fval,outputs]= fminsearch(@zhang_fit,[mu_berg_init lam_berg_init(2:J)],OPTIONS,R,f);
%num_steps = outputs.iterations;
mu_berg = berg_out(1:J);
lam_berg(2:J) = berg_out(J+1:2*J-1);
lam_berg(1) = f(1) - sum(lam_berg(2:J));
f = f';