[Pybrainsim-activity] SF.net SVN: pybrainsim:[144] trunk/src/HeadModelDipoleSphere.py
Status: Planning
Brought to you by:
rgoj
Menu
▾
▴
[Pybrainsim-activity] SF.net SVN: pybrainsim:[144] trunk/src/HeadModelDipoleSphere.py
|
From: <rg...@us...> - 2010-06-30 08:36:37
|
Revision: 144
http://pybrainsim.svn.sourceforge.net/pybrainsim/?rev=144&view=rev
Author: rgoj
Date: 2010-06-30 08:36:31 +0000 (Wed, 30 Jun 2010)
Log Message:
-----------
* Adding implementations of a dipole in a spherical head model
** First, and implementation of a lead field for a homogeneous infinite conductor
** Implementation of formulas from Brody 1973
** Implementation of formulas from Frank 1952 (with problems, however!)
** Implementation can be selected by commenting out appropriate lines of code
** Brody 1973 is selected as the default implementation
** WARNING: All formulas are now being calculated each time a lead field is calculated - this is extremely resource consuming, something should be done about it.
Modified Paths:
--------------
trunk/src/HeadModelDipoleSphere.py
Modified: trunk/src/HeadModelDipoleSphere.py
===================================================================
--- trunk/src/HeadModelDipoleSphere.py 2010-06-27 19:43:32 UTC (rev 143)
+++ trunk/src/HeadModelDipoleSphere.py 2010-06-30 08:36:31 UTC (rev 144)
@@ -15,7 +15,7 @@
# PyBrainSim. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
-from numpy import arange, array, ones, identity, dot, zeros, sin, cos, pi, sqrt
+from numpy import arange, array, ones, identity, dot, zeros, sin, cos, pi, sqrt, sum, arccos
__metaclass__ = type # New style classes. Is this necessary?
"""
@@ -73,6 +73,9 @@
# The number of electrodes and generators defines the size of the lead field matrix
self.leadField = zeros((nElectrodes,nGenerators))
+ self.leadFieldInfinite = zeros((nElectrodes,nGenerators))
+ self.leadFieldBrody1973 = zeros((nElectrodes,nGenerators))
+ self.leadFieldFrank1952 = zeros((nElectrodes,nGenerators))
for iElectrode in range(nElectrodes):
# Calculating the coordinates of the electrode in the Cartesian coordinates associated with the head
@@ -84,6 +87,10 @@
xyzElectrodeHead[2] = self.radius * cos(electrodeTheta);
for iGenerator in range(nGenerators):
+ #
+ # Infinite homogeneous conductor
+ #
+
# Calculating the coordinates of the dipole in the Cartesian coordinates associated with the head
dipoleRadius = self.radius - self.head.generatorSiteList[iGenerator][0]
dipoleTheta = self.head.generatorSiteList[iGenerator][1]
@@ -102,24 +109,129 @@
distance += pow(xyzElectrodeDipole[1],2.0)
distance += pow(xyzElectrodeDipole[2],2.0)
distance = sqrt(distance)
+
+ # Rotation matrix for translating the coordinates of the electrode in the dipole coordinates parallel
+ # to the reference coordinates at the center of the head to dipole coordinates where the dipole is specified
+ # i.e. where the z axis is radial.
+ rotationMatrix = zeros((3,3));
+ dipoleTheta = self.head.generatorSiteList[iGenerator][1]
+ dipolePhi = self.head.generatorSiteList[iGenerator][2]
+ # Row 1
+ rotationMatrix[0,0] = sin(dipolePhi)
+ rotationMatrix[0,1] = -cos(dipolePhi)
+ rotationMatrix[0,2] = 0
+ # Row 2
+ rotationMatrix[1,0] = cos(dipoleTheta) * cos(dipolePhi)
+ rotationMatrix[1,1] = cos(dipoleTheta) * sin(dipolePhi)
+ rotationMatrix[1,2] = -sin(dipoleTheta)
+ # Row 3
+ rotationMatrix[2,0] = sin(dipoleTheta) * cos(dipolePhi)
+ rotationMatrix[2,1] = sin(dipoleTheta) * sin(dipolePhi)
+ rotationMatrix[2,2] = cos(dipoleTheta)
+
+ rotationMatrixHeadToDipole = rotationMatrix
- self.leadField[iElectrode, iGenerator] = 1.0
+ # The Electrode coordinates in the dipole rotated frame of reference and scaleda
+ xyzElectrodeDipoleRotatedScaled = dot(rotationMatrix, xyzElectrodeDipole) / distance
+ # The Orientation vector
+ orientationTheta = self.head.generatorSiteList[iGenerator][3]
+ orientationPhi = self.head.generatorSiteList[iGenerator][4]
+ xyzOrientationDipoleRotated = zeros(3);
+ xyzOrientationDipoleRotated[0] = sin(orientationTheta) * cos(orientationPhi);
+ xyzOrientationDipoleRotated[1] = sin(orientationTheta) * sin(orientationPhi);
+ xyzOrientationDipoleRotated[2] = cos(orientationTheta);
+
+ # The cosine of the angle between the dipole orientation and the electrode
+ cosAngleOrientationElectrode = sum(xyzElectrodeDipoleRotatedScaled * array(xyzOrientationDipoleRotated))
+
+ # Homogeneous conductor without boundaries
+ self.leadFieldInfinite[iElectrode, iGenerator] = cosAngleOrientationElectrode/distance/distance/4/pi/sigma
+
+ #
+ # Bounded spherical conductor
+ # Brody 1973
+ #
+ rCosPhi = dot(xyzElectrodeHead, xyzDipoleHead) / self.radius
+ #print rCosPhi
+
+ fieldVector = zeros(3);
+ for i in range(3):
+ fieldVector[i] = 2*(xyzElectrodeHead[i] - xyzDipoleHead[i])/pow(distance,2.0)
+ fieldVector[i] += (1/pow(self.radius,2.0)) * (xyzElectrodeHead[i] + (xyzElectrodeHead[i] * rCosPhi - self.radius * xyzDipoleHead[i])/(distance + self.radius - rCosPhi))
+ fieldVector[i] = fieldVector[i] / 4 / pi / sigma / distance
+
+ # Rotation matrix for translating the coordinates in the dipole frame of reference to the
+ # coordinates associated with the dipole parallel to the head coordinates.
+ rotationMatrix = zeros((3,3));
+ dipoleTheta = self.head.generatorSiteList[iGenerator][1]
+ dipolePhi = self.head.generatorSiteList[iGenerator][2]
+ # Row 1
+ rotationMatrix[0,0] = sin(dipolePhi)
+ rotationMatrix[0,1] = cos(dipoleTheta) * cos(dipolePhi)
+ rotationMatrix[0,2] = sin(dipoleTheta) * cos(dipolePhi)
+ # Row 2
+ rotationMatrix[1,0] = -cos(dipolePhi)
+ rotationMatrix[1,1] = cos(dipoleTheta) * sin(dipolePhi)
+ rotationMatrix[1,2] = sin(dipoleTheta) * sin(dipolePhi)
+ # Row 3
+ rotationMatrix[2,0] = 0
+ rotationMatrix[2,1] = -sin(dipoleTheta)
+ rotationMatrix[2,2] = cos(dipoleTheta)
+
+ rotationMatrixDipoleToHead = rotationMatrix
+
+ # Rotating Orientation to translated dipole coordinates
+ xyzOrientationDipole = dot(rotationMatrixDipoleToHead, xyzOrientationDipoleRotated)
+
+ self.leadFieldBrody1973[iElectrode, iGenerator] = dot(fieldVector, xyzOrientationDipole)
+
+ #
+ # Bounded spherical conductor
+ # Frank 1952
+ #
+ xyzElectrodeDipoleRotated = dot(rotationMatrixHeadToDipole, xyzElectrodeHead)
+ # The below line causes problems because of arccos!
+ phi = arccos(xyzElectrodeDipoleRotated[0] / xyzElectrodeDipoleRotated[1]) - orientationPhi
+ cosPhi = cos(phi)
+ psi = orientationTheta
+ cosPsi = cos(psi)
+ sinPsi = sin(psi)
+ cosTheta = dot(xyzDipoleHead,xyzElectrodeHead) / sqrt(dot( xyzDipoleHead, xyzDipoleHead )) / sqrt(dot( xyzElectrodeHead, xyzElectrodeHead))
+ sinTheta = sqrt(1-pow(cosTheta,2.0))
+ b = self.radius - self.head.generatorSiteList[iGenerator][0]
+ f = b/self.radius
+ self.leadFieldFrank1952[iElectrode, iGenerator] = cosPsi * (( 1-pow(f,2.0) )/pow( 1+pow(f,2.0)-2*f*cosTheta, 3/2)-1)
+ self.leadFieldFrank1952[iElectrode, iGenerator] += sinPsi*cosPhi/sinTheta* ( ( 3*f - 3*pow(f,2.0)*cosTheta + pow(f,3.0) - cosTheta )/pow(1+pow(f,2.0)-2*f*cosTheta,3/2) + cosTheta )
+ self.leadFieldFrank1952[iElectrode, iGenerator] = self.leadFieldFrank1952[iElectrode, iGenerator] / 4 /pi / sigma / self.radius / b
+
def runHeadModel(self, generatorOutput, recording):
# Calculating the lead field
self.calculateLeadField()
- print self.leadField
- print (array([generatorOutput])).transpose()
-
# Calculating the recording by multyplying the lead field by the generator output
newRecording = dot(self.leadField, (array([generatorOutput])).transpose())
+ newRecordingInfinite = dot(self.leadFieldInfinite, (array([generatorOutput])).transpose())
+ newRecordingBrody1973 = dot(self.leadFieldBrody1973, (array([generatorOutput])).transpose())
+ newRecordingFrank1952 = dot(self.leadFieldFrank1952, (array([generatorOutput])).transpose())
# Converting the new recording to list format for compatibility with rest of the code
newRecording = list(newRecording[:,0])
+ newRecordingInfinite = list(newRecordingInfinite[:,0])
+ newRecordingBrody1973 = list(newRecordingBrody1973[:,0])
+ newRecordingFrank1952 = list(newRecordingFrank1952[:,0])
- print newRecording
-
+ for i in range(len(newRecording)):
+ # All three of the lead field calculation methods implemented above can be used, however
+ # only Brody1973 and Frank1952 give results in a bounded homogeneous conducting sphere,
+ # and of these only Brody1973 gives results for all combinations of angles, etc.
+ # Additionally, due to and arccos in my implementation above, Frank1952 for now gives even
+ # worse results than it should for many angles...
+ newRecording[i] = newRecordingBrody1973[i]
+ #newRecording[i] = newRecordingFrank1952[i]
+ #newRecording[i] = newRecordingInfinite[i]
+ #newRecording[i] = (newRecordingBrody1973[i])/(newRecordingFrank1952[i])
+
for channel in range(len(recording)):
recording[channel].append(newRecording[channel])
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
