eeg_interp_sph_spline - Spherical Spline Laplacian of Potential
Useage: [S,Li] = eeg_lap_sph_spline(voltage,X,Y,Z,Npoints)
where: 'voltage' is an EEG/ERP measurement at time t from
electrode positions (X,Y,Z) on a scalp surface. All
input arrays are the same size, assumed (Nelec x 1).
The origin of X,Y,Z is assumed (0,0,0).
S => the x,y,z points on a hemisphere
Li => spherical spline Laplacian of voltages
S is generated by 'elec_sphere_points' with
Rpoints=Npoints and Epoints=16.
Notes: This function calculates the spherical spline Laplacian
of Perrin et al (1989). Electroenceph. & Clin.
Neurophysiology, 72: 184-187.0001 function [S,Li] = eeg_interp_sph_spline(V,X,Y,Z,Npoints) 0002 0003 % eeg_interp_sph_spline - Spherical Spline Laplacian of Potential 0004 % 0005 % Useage: [S,Li] = eeg_lap_sph_spline(voltage,X,Y,Z,Npoints) 0006 % 0007 % where: 'voltage' is an EEG/ERP measurement at time t from 0008 % electrode positions (X,Y,Z) on a scalp surface. All 0009 % input arrays are the same size, assumed (Nelec x 1). 0010 % The origin of X,Y,Z is assumed (0,0,0). 0011 % 0012 % S => the x,y,z points on a hemisphere 0013 % Li => spherical spline Laplacian of voltages 0014 % 0015 % S is generated by 'elec_sphere_points' with 0016 % Rpoints=Npoints and Epoints=16. 0017 % 0018 % Notes: This function calculates the spherical spline Laplacian 0019 % of Perrin et al (1989). Electroenceph. & Clin. 0020 % Neurophysiology, 72: 184-187. 0021 % 0022 0023 % $Revision: 1.2 $ $Date: 2005/07/12 21:49:43 $ 0024 0025 % Licence: GNU GPL, no implied or express warranties 0026 % History: 08/01 Darren.Weber_at_radiology.ucsf.edu, with 0027 % mathematical and matlab advice from 0028 % Dr. Murk Bottema (Flinders University of SA) 0029 % 0030 % 08/01 Needs testing & verification! 0031 % With large electrode arrays, the result should be 0032 % regularized (not sure why). 0033 % 0034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0035 0036 0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0038 % Check for correct size & orientation of X,Y,Z,V 0039 [x1,x2] = size(X); 0040 [y1,y2] = size(Y); 0041 [z1,z2] = size(Z); 0042 [v1,v2] = size(V); 0043 if ~and(isequal(x1,y1,z1,v1),isequal(x2,y2,z2,v2)) 0044 error('ERROR: all X, Y, Z, V must be size (Nx1)'); 0045 end 0046 if x1 < x2, X = X'; [x1,x2] = size(X); end 0047 if y1 < y2, Y = Y'; [y1,y2] = size(Y); end 0048 if z1 < z2, Z = Z'; [z1,z2] = size(Z); end 0049 if v1 < v2, V = V'; [v1,v2] = size(V); end 0050 if ~(isequal(x2,y2,z2,v2,1)) 0051 error('ERROR: all X, Y, Z, V must be size (Nx1)'); 0052 end 0053 0054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0055 % Computations derived from Perrin et al. (1989) 0056 0057 % G * C = V, where: 0058 % 0059 % G (NxN) spherical spline 0060 % C (Nx1) spline coefficients 0061 % V (Nx1) voltage potentials at (X,Y,Z) electrode locations 0062 % 0063 % solve for C = V * Sp' (C = V\Sp; (see "help slash")) 0064 0065 % First calc G, where G is a function of the cosine of 0066 % the angle (theta) between electrode point vectors 0067 A = [ X Y Z ]; 0068 COS = cosines(A,A); 0069 G = spheric_spline(COS); 0070 C = spline_coefficients(G,V); % spline coefficients (Co,C1,...,Cn) 0071 Co = C(1); 0072 Ci = C(2:end); 0073 0074 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0075 % Obtain interpolated potentials at S (eq.1, Perrin et al., 1989) 0076 % V(s) = Co + sum( Ci * g(x) ) 0077 0078 % get spherical electrode radius 0079 [r,x,y,z] = elec_sphere_fit(X,Y,Z,0,0,0,0); 0080 % Generate spherical points for interpolation 0081 [x,y,z] = elec_sphere_points(16,Npoints,r); 0082 S = [x y z]; 0083 0084 fprintf('...Spherical Interpolation Progress:\n'); 0085 rows = 80; % progress indicator 0086 0087 for[p] = 1:length(S) 0088 0089 B = S(p,:); 0090 Cos = cosines(A,B); %(1xN) 0091 0092 Gx = spheric_spline(Cos); %(1xN) 0093 0094 CiGx = Ci .* Gx'; %(1xN) 0095 0096 Vi(p,1) = Co + sum(CiGx); 0097 0098 fprintf('.'); % progress indicator 0099 if([p] == rows ) fprintf('\n'); rows = rows + 80; end 0100 end 0101 0102 return 0103 0104 0105 0106 0107 0108 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0109 % Solve eq. 3 Perrin et al. (1989) 0110 % g(COS) = 1/4pi * sum[n=1:inf] (( (2*n+1)/(n^m * (n+1)^m) ) * Pn(COS)); 0111 0112 function [Gx] = spheric_spline(Cosine) 0113 0114 m = 4; N = 7; % gives accuracy of 10^-6 0115 0116 P = legendre(N,Cosine); 0117 %P = LEGENDRE(N,X) computes the associated Legendre functions 0118 %of degree N and order m = 0, 1, ..., N, evaluated for each element 0119 %of X. 0120 %In general, P has one more dimension than X. 0121 %Each element P(m+1,i,j,k,...) contains the associated Legendre 0122 %function of degree N and order m evaluated at X(i,j,k,...). 0123 0124 ndim = ndims(P); 0125 switch ndim 0126 case 2, P = P(2:N+1,:); 0127 case 3, P = P(2:N+1,:,:); 0128 case 4, P = P(2:N+1,:,:,:); 0129 otherwise 0130 end 0131 0132 k = (1/4 * pi); 0133 0134 for n = 1:(N), Series(n,1) = (2*n + 1) / (n^(m-1) * (n+1)^(m-1)); end 0135 0136 if min(size(Cosine)) == 1, Gx = k * ( Series' * P ); 0137 else 0138 for i = 1:length(Cosine), Gx(i,:) = k * ( Series' * P(:,:,i) ); 0139 end 0140 end 0141 return 0142 0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0144 % Solve eq. 2 Perrin et al. (1989) 0145 function [C] = spline_coefficients(Gx,V) 0146 0147 % add ones to first row & column of Gx 0148 tmp = ones(size(V)); Gx = [tmp Gx]; 0149 tmp = [ones(size(V))' 1]'; Gx = [tmp Gx']'; 0150 0151 Gx(1,1) = 0; % according to Murk 0152 %Gx(1,1) = 1; % according to EMSE, Greenblatt 0153 0154 CoV = [0 V']'; 0155 0156 C = Gx\CoV; 0157 return 0158 0159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0160 function [Cos] = cosines(A,B) 0161 0162 for a = 1:size(A,1), Aa = A(a,:); A_len = norm(Aa); 0163 for b = 1:size(B,1), Bb = B(b,:); B_len = norm(Bb); 0164 0165 Cos(a,b) = dot(Aa,Bb) / (A_len * B_len); 0166 end 0167 end 0168 return