Loreta localization with MNE-Python

loreta.py

import mne
 
raw = mne.io.read_raw_bti('jun/MEG/t_SEEG3/7/c,rfDC')
raw_meg = raw.pick_channels([ch for ch in raw.ch_names
if ch.startswith('MEG')], copy=True)
 
tr = mne.read_trans('patient/MEG/run7-trans.fif')
src = mne.setup_volume_source_space (
    subject='patient',
    pos=10.0, # mm
    bem='freesurfer/subjects/patient/bem/patient-320-320-320-bem-sol.fif',
    mri='freesurfer/subjects/patient/mri/T1.mgz'
)
fwd = mne.make_forward_solution(raw_meg.info, mri=tr, src=src, bem=bem)

## construct Loreta inverse, following Pascual Marqui et al '94,
# J = T phi; T = (W B' B W)^-1 K' {K (W B' B W)^-1 K'}+
#                \-----------/    \-------------------/
#                      wb                   kw
#
# J - source current
# T - inverse op
# phi - measurement
# B - Laplacian
# K - lead field
# W - normalization / whitener

## setup Laplacian 
# grid points
rr = src[0]['rr']
# which are used
umask = src[0]['inuse'] > 0
rrm = rr[umask]
# build laplacian op for masked grid
rrd = np.sqrt(((rrm - rrm.reshape((-1, 1, 3)))**2).sum(axis=-1))
rrd_ = rrd.reshape(-1)
d = rrd_[rrd_>0.0].min() # bin size
neighbors = -1 * (rrd > (0.999 * d)) * (rrd < (1.001 * d))
self = 6 * (rrd < (0.9 * d))
L = self + neighbors
# mask op where valid, i.e. all 6 neighbors exist
lmask = L.sum(axis=1) == 0
B = L[lmask][:, lmask] / (d ** 2)
# 3 components per location 
B = np.repeat(np.repeat(B, 3, axis=0), 3, axis=1)
 
# keep fwd soln where laplacian is valid
K = fwd['sol']['data']
K = K.reshape((248, -1, 3))[:, lmask].reshape((248, -1))
# loreta normalization
W = np.diag(map(np.linalg.norm, K.T))

# compute inverse op
wb = np.linalg.inv(W.dot(B.T).dot(B).dot(W))
kw = np.linalg.pinv(K.dot(wb).dot(K.T))
T = wb.dot(K.T).dot(kw)

Volume source space grid. In gray, unused points, in red, used but invalid points for Laplacian. In blue, valid points for Laplacian, in cyan principal components of T the inverse operator.