Meg Leadfield Calc

---

Calculate the leadfields from MEG data in python

show code

collapse code

# Copyright 2008 Dan Collins
#
 
# This 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 2 of the License, or
# (at your option) any later version.
 
# And 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 Build; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
# A portion of this code was used and adapted from meg_leadfield1.m function distributed with
# the matlab fieldtrip package at fcdonders.
 
# Thanks to the fieldtrip team for use of their code under the GNU General Public License.
 
"""Compute leadfields
pos=R point
pos= channel location vector
ori=signal channel orientation
 
ex. x=leadfield.calclf(grid) ....
    will return the leadfield sum from both upper and lower coils
 
    grid is N X 3 array. Make sure to reshape it if Ndim=1 to by 1X3
    ex. IF....shape(grid) RETURNS (3,) THEN grid=grid.reshape(1,3)
    or just make grid like... grid=array([[0, 6, 3]])"""
 
from numpy import zeros, array, dot, cross, shape
from numpy.linalg import norm 
import datetime
from meg import cos,pos
#import pdbreturn lambda x: x + n
 
class grid():
    def __init__(self,grid):
        pass
 
class getch():
    def __init__(self):
        "give grid points and return leadfields"
        self.channels=pos.sensors()
        chup=self.channels.chu #upper ch position
        chlp=self.channels.chl #lower ch position
        chud=self.channels.chu #upper ch direction
        chld=self.channels.chu #lower ch direction
 
class getlf():
    """grid=voxel location. dimensions: N X 3
    pos=position of channels. returned by pos.py
    ori=orientation of channels. returned by pos.py"""
    def __init__(self, grid, pos, ori):
        #grid.__init__(self, grid)
        pos=pos*100 #convert channel pos into cm from meters
        self.grid=grid;self.pos=pos;self.ori=ori
 
        c=cos.center()
        spherecenter=c.spherecenter
        self.spherecenter=spherecenter
 
        chp=pos-spherecenter
        loc=grid-spherecenter
        self.loc=loc
        gridshape=shape(loc)
        if len(gridshape)==1:
            loc.shape=loc[1,3]
 
        for j in range(0,len(loc)): #for each grid point
            "for each R point (chp) calculate the leadfield"
            self.nchans = len(chp)
            nchans=self.nchans
            R=loc[j,:]
            self.position=chp;position=chp
            self.orientation=ori;orientation=ori
            lftmp = zeros((nchans,3),float);
 
            tmp2 = norm(R);
 
            for i in range(0,nchans):
                t = position[i]
                #r=r
                o = orientation[i]
                #u=u
 
                tmp1 = norm(t);
                tmp3 = norm(t-R);
                tmp4 = dot(t,R);
                tmp5 = dot(t,t-R);
                tmp6 = dot(R,t-R);
                tmp7 = (tmp1*tmp2)**2 - tmp4**2; #% cross(r,R)**2
                #tmp7=tmp7
 
                self.alpha = 1 / (-tmp3 * (tmp1*tmp3+tmp5));
                alpha=self.alpha
                A = 1/tmp3 - 2*alpha*tmp2**2 - 1/tmp1;
                B = 2*alpha*tmp4;
                C = -tmp6/(tmp3**3);
                #A=A;B=B;C=C;
                self.beta = dot(A*t + B*R + C*(t-R), o)/tmp7;
                beta=self.beta
                #i=i
 
                lftmp[i,:] = cross(alpha*o  + beta*t, R);
                #print cross(alpha*o  + beta*r, R)
                #print lftmp[0,:]
            self.lf = 1e-7*lftmp; #% multiply with u0/4pi
            #lf = lftmp; #% multiply with u0/4pi
 
def calclf(grid):
    """calc leadfield script to use pos.py and getlf class to get the sum of upper and lower coils"""
    #gradiometer config
    coilsperch=2;
    channels=pos.sensors()
 
    chlabels={'set1':'upper','set2':'lower'}
    chpos={'set1':'chu','set2':'chl'}
    chdir={'set1':'chudir','set2':'chldir'}
 
    #for i in chlabels:
    lp=channels.chl;ld=channels.chldir;
    up=channels.chu;ud=channels.chudir;
 
    #sample grid
    #grid=array([6.2746,12.5286,-10.0777])
    #grid.shape=(1,3)
    lowerlf=getlf(grid, lp ,ld)
    upperlf=getlf(grid, up ,ud)
 
    leadfieldfinal=lowerlf.lf+upperlf.lf
    return leadfieldfinal

show howto

collapse howto