[Pybrainsim-activity] SF.net SVN: pybrainsim:[138] trunk/src/HeadModelDipoleSphere.py

[<rg...@us...>]

Revision: 138
          http://pybrainsim.svn.sourceforge.net/pybrainsim/?rev=138&view=rev
Author:   rgoj
Date:     2010-06-27 15:31:45 +0000 (Sun, 27 Jun 2010)

Log Message:
-----------
* Adding a work in progress of a single conducting sphere head model with generators as dipoles at last!

Added Paths:
-----------
    trunk/src/HeadModelDipoleSphere.py

Added: trunk/src/HeadModelDipoleSphere.py
===================================================================
--- trunk/src/HeadModelDipoleSphere.py	                        (rev 0)
+++ trunk/src/HeadModelDipoleSphere.py	2010-06-27 15:31:45 UTC (rev 138)
@@ -0,0 +1,126 @@
+# PyBrainSim Copyright 2009 Roman Goj
+#
+# This file is part of PyBrainSim.
+#
+# PyBrainSim is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# PyBrainSim is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# 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
+__metaclass__ = type # New style classes. Is this necessary?
+
+"""
+This reimplemenation of the HeadModel class models the head as a single 
+conductant sphere containing dipolar sources.
+"""
+
+from HeadModel import HeadModel
+
+class HeadModelDipoleSphere(HeadModel):
+    def __init__(self, head, radius):
+        HeadModel.__init__(self, head)
+        self.head = head
+        self.radius = radius # Radius of the head in [cm]
+        self.leadField = 0;
+    
+    def checkPosition(self,position):
+        # Checking if position is a list of two or five numbers
+        # The position of each registration site is given in spherical coordinates
+        # associated with the center of the spherical head. The radial coordinate of the electrode
+        # is always equal to the radius of the head, hence it doesn't need to be specified.
+        # The only two parameters that need to be specified are the angles theta and phi in this
+        # order.
+        # For each generator, all three cordinates must be given: the depth of the dipole and the
+        # theta and phi angles. Additionaly, angles specyfying the orientation of the dipole in
+        # the spherical coordinates of the dipole must also be given. The theta angle is then the
+        # orientation with respect to the surface of the head - i.e. theta = 0 means a radial
+        # dipole, while theta = pi/2 means a tangential dipole. Only values in this range are
+        # allowed! The specification of the phi angle is considered less important and should be
+        # explained more thoroughly in other documentation.
+        if not isinstance(position, list):
+            return False
+        if len(position) != 2 and len(position) != 5:
+            return False
+        # Check the radial coordinate (depth) of the dipole
+        if len(position) == 5 and (position[0] <= 0 or position[0] >= self.radius):
+            return False
+        # Check the theta angle of the orientation of the dipole
+        if len(position) == 5 and (position[3] < 0 or position[3] > pi/2):
+            return False
+        for i in range(len(position)):
+            if (not isinstance(position[i], float)) and (not isinstance(position[i], int)):
+                return False
+        
+        # If none of the conditions above fail, the position is correct
+        return True
+
+    def calculateLeadField(self):
+        # How many generators and electrodes do we have?
+        nGenerators = len(self.head.generatorSiteList)
+        nElectrodes = len(self.head.registrationSiteList)
+        
+        # Assuming ideal conductivity
+        sigma = 1.0
+
+        # The number of electrodes and generators defines the size of the lead field matrix
+        self.leadField = zeros((nElectrodes,nGenerators))
+        
+        for iElectrode in range(nElectrodes):
+            # Calculating the coordinates of the electrode in the Cartesian coordinates associated with the head
+            electrodeTheta = self.head.registrationSiteList[iElectrode][0]
+            electrodePhi = self.head.registrationSiteList[iElectrode][1]
+            xyzElectrodeHead = zeros(3);
+            xyzElectrodeHead[0] = self.radius * sin(electrodeTheta) * cos(electrodePhi);
+            xyzElectrodeHead[1] = self.radius * sin(electrodeTheta) * sin(electrodePhi);
+            xyzElectrodeHead[2] = self.radius * cos(electrodeTheta);
+            
+            for iGenerator in range(nGenerators):
+                # 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]
+                dipolePhi = self.head.generatorSiteList[iGenerator][2]
+                xyzDipoleHead = zeros(3);
+                xyzDipoleHead[0] = dipoleRadius * sin(dipoleTheta) * cos(dipolePhi);
+                xyzDipoleHead[1] = dipoleRadius * sin(dipoleTheta) * sin(dipolePhi);
+                xyzDipoleHead[2] = dipoleRadius * cos(dipoleTheta);
+                
+                # Calculating the coordinates of the electrode in the coordinates associated with the dipole
+                xyzElectrodeDipole = xyzElectrodeHead - xyzDipoleHead;
+                
+                # Calculating the distance between the dipole and the electrode
+                distance = 0;
+                distance += pow(xyzElectrodeDipole[0],2.0)
+                distance += pow(xyzElectrodeDipole[1],2.0)
+                distance += pow(xyzElectrodeDipole[2],2.0)
+                distance = sqrt(distance)
+
+                self.leadField[iElectrode, iGenerator] = 1.0
+
+    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())
+        
+        # Converting the new recording to list format for compatibility with rest of the code
+        newRecording = list(newRecording[:,0])
+        
+        print newRecording
+        
+        for channel in range(len(recording)):
+            recording[channel].append(newRecording[channel])
+            
+        return recording


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.