source_dipole_sphere - calculate surface potential on a sphere containing a dipole
Calculates surface voltage, for a single sphere containing a radial
dipole. This function is limited to radial dipoles (a dipole at the
origin can have any orientation).
The sphere is assigned a homogenous, isotropic conductivity,
sigma = 0.33 Siemens/meter, which is commonly used to estimate brain
conductivity.
Also, the distance of dipole source to sink is 1e-12 (1 pico meter).
FV = source_dipole_sphere([dipole_moment],[dipole_origin],...
[dipole_orient],[sphere_radius],[Lmax]);
-------------------------------------------------------------------------
Inputs
dipole_moment - Ampere meters, 10e-9 by default, which approximates a
small patch of cortex (see Hamalainen et al [1993], Reviews of Modern
Physics, 65(2):413-497).
dipole_origin - The XYZ Cartesian coordinates of the midpoint between the
source and sink of the dipole. The orientation of the dipole will be
radial, with source and sink arranged either side of midpoint. See also
the sphere radius specification below.
dipole_orient - Only if the dipole_origin = [0,0,0] (the default) is it
possible to specify the dipole orientation; eg, [1 0 0] is along the
x-axis, [0 1 0] is along the y-axis, etc. If the dipole is located away
from the origin, the orientation is along the radial line from the origin
to the dipole.
sphere_radius - in meters; the default 0.1 approximates a human head.
Note that any specification of dipole origin should be within the sphere
radius.
Lmax - the order of magnitude for the sum of Legendres, which are used
for the analytic solution of potential at a spherical surface. The
default value is 20, determined after examination of values from 1:250
(20 is conservative, as convergence was observed at Lmax < 5).
-------------------------------------------------------------------------
Returns the FV struct, a spherical tesselation
FV =
vertices: [2562x3 double]
faces: [5120x3 double]
Cdata: [2562x1 double] - surface potential (Volts)
To visualize the result, use:
patch('vertices',FV.vertices,'faces',FV.faces,...
'FaceVertexCData',FV.Cdata,...
'Facecolor','interp');0001 function FV = source_dipole_sphere(dipole_moment,dipole_origin,dipole_orient,sphere_radius,LMAX) 0002 0003 % source_dipole_sphere - calculate surface potential on a sphere containing a dipole 0004 % 0005 % Calculates surface voltage, for a single sphere containing a radial 0006 % dipole. This function is limited to radial dipoles (a dipole at the 0007 % origin can have any orientation). 0008 % 0009 % The sphere is assigned a homogenous, isotropic conductivity, 0010 % sigma = 0.33 Siemens/meter, which is commonly used to estimate brain 0011 % conductivity. 0012 % 0013 % Also, the distance of dipole source to sink is 1e-12 (1 pico meter). 0014 % 0015 % FV = source_dipole_sphere([dipole_moment],[dipole_origin],... 0016 % [dipole_orient],[sphere_radius],[Lmax]); 0017 % ------------------------------------------------------------------------- 0018 % Inputs 0019 % 0020 % dipole_moment - Ampere meters, 10e-9 by default, which approximates a 0021 % small patch of cortex (see Hamalainen et al [1993], Reviews of Modern 0022 % Physics, 65(2):413-497). 0023 % 0024 % dipole_origin - The XYZ Cartesian coordinates of the midpoint between the 0025 % source and sink of the dipole. The orientation of the dipole will be 0026 % radial, with source and sink arranged either side of midpoint. See also 0027 % the sphere radius specification below. 0028 % 0029 % dipole_orient - Only if the dipole_origin = [0,0,0] (the default) is it 0030 % possible to specify the dipole orientation; eg, [1 0 0] is along the 0031 % x-axis, [0 1 0] is along the y-axis, etc. If the dipole is located away 0032 % from the origin, the orientation is along the radial line from the origin 0033 % to the dipole. 0034 % 0035 % sphere_radius - in meters; the default 0.1 approximates a human head. 0036 % Note that any specification of dipole origin should be within the sphere 0037 % radius. 0038 % 0039 % Lmax - the order of magnitude for the sum of Legendres, which are used 0040 % for the analytic solution of potential at a spherical surface. The 0041 % default value is 20, determined after examination of values from 1:250 0042 % (20 is conservative, as convergence was observed at Lmax < 5). 0043 % 0044 % ------------------------------------------------------------------------- 0045 % Returns the FV struct, a spherical tesselation 0046 % 0047 % FV = 0048 % vertices: [2562x3 double] 0049 % faces: [5120x3 double] 0050 % Cdata: [2562x1 double] - surface potential (Volts) 0051 % 0052 % To visualize the result, use: 0053 % 0054 % patch('vertices',FV.vertices,'faces',FV.faces,... 0055 % 'FaceVertexCData',FV.Cdata,... 0056 % 'Facecolor','interp'); 0057 % 0058 0059 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:36 $ 0060 0061 % Licence: GNU GPL, no implied or express warranties 0062 % History: 12/2003, Darren.Weber_at_radiology.ucsf.edu 0063 % developed from notes and advice from 0064 % Dr. Tom.Ferree_at_radiology.ucsf.edu 0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0066 0067 0068 % LMAX is a parameter used in the Legendre approximation 0069 if ~exist('LMAX','var'), LMAX = 20; end 0070 if isempty(LMAX), LMAX = 20; end 0071 0072 % create a spherical surface with a radius that approximates the human 0073 % head, 0.1 meters (10 cm) 0074 if ~exist('sphere_radius','var'), sphere_radius = 0.1; end 0075 if isempty(sphere_radius), sphere_radius = 0.1; end 0076 0077 % sphere_tri is an external function 0078 FV = sphere_tri('ico',4,sphere_radius); 0079 0080 % assign a homogenous, isotropic conductivity to this spherical region 0081 sigma = 0.33; % Siemens/meter, common brain value 0082 fprintf('...assigned brain conductivity (0.33 S/m) to volume\n'); 0083 0084 % distance of dipole source to sink, could approximate the cortical 0085 % thickness 0086 dipole_distance = 1e-12; % 1 pico meter 0087 fprintf('...dipole has source and sink %g meters apart\n',dipole_distance); 0088 0089 % dipole moment, Ampere meters (10 nA m) 0090 if ~exist('dipole_moment','var'), dipole_moment = 10e-9; end 0091 if isempty(dipole_moment), dipole_moment = 10e-9; end 0092 fprintf('...dipole moment = %g Ampere meters\n',dipole_moment); 0093 0094 % dipole current, Amperes 0095 dipole_current = dipole_moment / dipole_distance; 0096 fprintf('...dipole current = %g Amperes\n',dipole_current); 0097 0098 % Define the vector that lies at the midpoint between the source and sink. 0099 % The origin is the default. 0100 if ~exist('dipole_origin','var'), dipole_origin = [0,0,0]; end 0101 if isempty(dipole_origin), dipole_origin = [0,0,0]; end 0102 fprintf('...dipole origin = [ %g, %g, %g ] meters\n',dipole_origin); 0103 0104 % find the magnitude (ie, radius) of the dipole (do not use the matlab norm 0105 % function, it doesn't calculate the same value) 0106 dipole_radius = sqrt( sum( dipole_origin .^2 ) ); 0107 0108 % Calculate the radial dipole unit vector, which is used to find the cosine 0109 % of the angle between the dipole and surface locations below 0110 if dipole_radius, 0111 0112 % the dipole is not located at the origin 0113 dipole_unit_vector = dipole_origin / dipole_radius; 0114 else 0115 0116 % the dipole is located at the origin, so the radial orientation is not 0117 % clearly defined and we need another input parameter to define the 0118 % dipole unit vector 0119 0120 % dipole geometry - oriented along z axis 0121 if ~exist('dipole_orient','var'), dipole_orient = [0 0 1]; end 0122 if isempty(dipole_orient), dipole_orient = [0 0 1]; end 0123 0124 % now calculate the dipole unit vector 0125 dipole_unit_vector = dipole_orient / sqrt( sum( dipole_orient .^ 2 )); 0126 0127 end 0128 fprintf('...dipole unit vector = [ %g, %g, %g ]\n',dipole_unit_vector); 0129 0130 % calculate location of source and sink, just for fun! 0131 source = dipole_origin + (dipole_unit_vector * dipole_distance); 0132 sink = dipole_origin - (dipole_unit_vector * dipole_distance); 0133 fprintf('...source vector = [ %+g, %+g, %+g ] meters\n', source); 0134 fprintf('...sink vector = [ %+g, %+g, %+g ] meters\n', sink); 0135 0136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0137 % solve voltage on surface 0138 0139 fprintf('...calculating surface voltage...'); tic 0140 0141 % initialise the surface potentials 0142 FV.Cdata = zeros(length(FV.vertices),1); 0143 0144 if dipole_radius, 0145 0146 DIPOLE_CONSTANT = ( -dipole_moment / (4 * pi * sigma) ) * (1 / (dipole_radius^2) ); 0147 0148 for vertex = 1:length(FV.vertices), 0149 0150 % extract spherical surface coordinates 0151 surface_vector = FV.vertices(vertex,:); 0152 0153 surface_radius = sqrt( sum ( surface_vector .^2 ) ); 0154 0155 if surface_radius < dipole_radius, 0156 error('surface radius is smaller than dipole radius.'); 0157 end 0158 0159 if surface_radius, 0160 surface_unit_vector = surface_vector / surface_radius; 0161 else 0162 surface_unit_vector = [0 0 0]; 0163 warning('surface radius is zero.'); 0164 end 0165 0166 cos_theta = dot(dipole_unit_vector, surface_unit_vector); 0167 0168 summation = 0; 0169 0170 for L = 1:LMAX, 0171 0172 TERM1 = 2 * L + 1; 0173 0174 TERM2 = ( dipole_radius / surface_radius )^(L + 1); 0175 0176 PnX = LegendreP(L,cos_theta); 0177 if length(PnX) > 1, 0178 PnX_sum = sum(PnX(2:end)); 0179 else 0180 PnX_sum = 0; 0181 end 0182 summation = summation + ( TERM1 * TERM2 * PnX_sum ); 0183 0184 end 0185 0186 FV.Cdata(vertex) = DIPOLE_CONSTANT * summation; 0187 0188 end % for vertex 0189 0190 else 0191 0192 % special case where dipole is at the origin. 0193 0194 DIPOLE_CONSTANT = ( -dipole_moment / (4 * pi * sigma) ); 0195 0196 for vertex = 1:length(FV.vertices), 0197 0198 surface_vector = FV.vertices(vertex,:); 0199 0200 surface_radius = sqrt( sum ( surface_vector .^2 ) ); 0201 0202 if surface_radius < dipole_radius, 0203 error('surface radius is smaller than dipole radius.'); 0204 end 0205 0206 if surface_radius, 0207 surface_unit_vector = surface_vector / surface_radius; 0208 else 0209 surface_unit_vector = [0 0 0]; 0210 warning('surface radius is zero.'); 0211 end 0212 0213 cos_theta = dot(dipole_unit_vector, surface_unit_vector); 0214 0215 FV.Cdata(vertex) = DIPOLE_CONSTANT * (3 / (surface_radius^2) ) * cos_theta; 0216 0217 end % for vertex 0218 0219 end 0220 0221 t=toc; fprintf('done (%5.2f sec)\n',t); 0222 0223 return 0224 0225 0226 0227 0228 0229 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0230 % SUBFUNCTIONS... 0231 0232 function p=LegendreP(nmax,x) 0233 0234 % function p=LegendreP(nmax,x) 0235 % Computes all Legendre polynomials from n=0 to n=nmax, 0236 % them in an array from n=1 to nmax+1 0237 % Uses the recursion relation in Numerical Recipes. 0238 % T Ferree 0239 0240 p=zeros(nmax+1,1); 0241 0242 if nmax>=0 0243 p0=1.0; 0244 p(1)=p0; 0245 end 0246 0247 if nmax>=1 0248 p1=x; 0249 p(2)=p1; 0250 end 0251 0252 if nmax>=2 0253 for j=2:nmax 0254 c1=(2.0*j-1.0)/(1.0*j); 0255 c2=(1.0*j-1.0)/(1.0*j); 0256 p(j+1)=c1*x*p(j-1+1)-c2*p(j-2+1); 0257 end 0258 end 0259 0260 return;