Bug Hunt

C-Day plus 34

For over a week now, the name of the game has been bug hunting. I had a finished first draft since the last blog post, so I’ve been trying to get the output of my custom SSS filter to match the proprietary version with sample data. One issue that took a couple days to track down was a simple but erroneous switch of two angles in a complex spherical coordinate gradient to Cartesian coordinate gradient transformation matrix. I can’t say that this is a new class of obstacles – spherical coordinates have thrown wrench after wrench into my code since different mathematicians regularly define these coordinates in different ways. (Is it just me, or is having seven separately accepted conventions for the spherical coordinate system a bit absurd?) My project crosses a couple domains of mathematics, so wrestling with these different conventions has helped me deeply appreciate the other mathematical concepts that do have a single accepted formulation.

Regardless, weeding out the spherical coordinate issue and a menagerie of other bugs has left me with a filter that produces filtered data that is similar to (but not exactly matching) the proprietary code (see some example output below). Luckily, I do have several checkpoints in the filter’s processing chain and I know the problem is between the last checkpoint and the final output. My mentors have been fantastic so far, and we have a potential bead on the last issue; the weak magnetic signals produced by the brain are measured with two major classes of MEG pickup coils: magnetometers and gradiometers. In a very simple sense, one measures the magnetic field while the other measures the spatial derivative of the magnetic field, and (because of this difference) they provide readings on very different scales that I have yet to normalize. Given some luck, this last patch could fix the issue and yield a working solution to the first half of my GSoC project! (Knock on wood.)

Exemplar data showing raw unfiltered MEG signal and the same data after the benchmark SSS filter and my own custom filtering (top). Difference between benchmark and my custom implementation (bottom). The filter in progress is close, but not quite the same as the benchmark implying there remains some bugs to fix.

Inner workings of the Maxwell filter

C-Day plus 20

This week, I finished a first draft of the Maxwell filtering project. Remember my goal here is to implement an open-source version of this filter that uses physics to separate brain signal from environmental garbage picked up by the MEG sensors. Now comes the fun part of this project: trying to add all the small tweaks required to precisely match the proprietary Maxwell filter, which I cannot access. I’m sure this will devolve into a tedious comparison between what I’ve implemented and the black box version, so here’s to hoping the proprietary code follows the original white papers.

Most of the Maxwell filter work up until this point was focused on enabling the calculation of the multipolar moment space (comprised of the spherical harmonics), which is the foundation of this Maxwell filter. These multipolar moments are the basis set I’ve mentioned earlier that allow the brain signals to be divided into two parts: those coming from within a sphere and those originating from outside a slightly larger sphere (to see this graphically, cf. Fig 6 in Taulu, et al., 2005). In essence, representing the brain signals as a sum of multipolar moments permits the separation of brain signal from external noise sources like power lines, large moving objects, the Earth’s magnetic field, etc. My most recent code actually projects the brain signals onto this multipolar moment space (i.e., representing MEG data as a sum of these moments), and then reconstructs the signal of interest. These calculations are all standard linear algebra. From Taulu and Simola, 2006 (pg 1762): 

Takeaway: The below equations show how Maxwell filtering is accomplished once the appropriate space has been calculated. We take brain signals recorded using MEG, represent them in a different space (the truncated mutlipolar moment space), and then reconstruct the MEG signal to apply the Maxwell filter and greatly reduce the presence of environmental noise.

ϕ represents MEG recordings
S represents the multipolar moment space (each vector is a spherical harmonic)
x represents the ideal weight of each multipolar moment (how much of each basis vector is present)
a hat represents an estimate
in, out, refer to internal spaces (brain signals), and external spaces (noise), respectively
x^T is the inverse of x


In the ideal case, the signal we recorded can also be represented as a weighted combination of our multipolar moments: 
ϕ = S * x
The S matrix contains multipolar moments but only up to a certain complexity (or degree), so it has been truncated. See my first post (end of 3rd paragraph) about why we cut out the very complex signals. 

Since we can break up the multipolar moments and their weights into an internal and external space (comprised of brain signal and noise), this is equivalent to the last equation:
ϕ = [S_in, S_out] * [x_in, x_out]^T

However, we're not in an ideal world, so we need to estimate these mutlipolar moment weights. x is the unknown so isolate it by taking the pseudo inverse of S to solve for an estimate of multipolar moment weights:

S_pinv * ϕ = S_pinv * S * x
S_pinv * ϕ = x_hat 
x_hat = S_pinv * ϕ
or equivalently,
[x_in_hat, x_out_hat]^T = S_pinv * ϕ

With the multipolar weight estimates in hand, we can finally reconstruct our original MEG recordings, which effectively applies the Maxwell filter. Again, since S_in, and S_out have been truncated, they only recreate signals up to a certain spatial complexity to cut out the noise.

ϕ_in_hat = S_in * x_in_hat
ϕ_out_hat = S_out * x_out_hat

The above ϕ matrices are a cleaner version of the brain signal we started with and the world is now a much better place to live in.

Progress and Paris

C-Day plus 13

·        I’ve spent the past couple of weeks have getting my hands dirty with the underlying physics equations behind solving for forward solutions in source imaging (mostly with Mosher et al., 1999 and Hämäläinen &Sarvas, 1989). The forward solution is a matrix that relates each point on the cortex to the MEEG sensors. Once you have the forward solution, you can use some pretty fancy mathematics (including Tikinov regularization) to find a pseudoinverse for this matrix – called the inverse solution. Then you’re able to relate the sensor measurements with estimates of what areas of the brain are active, which is the fundamental motivation for source imaging.

·        One component of my project is focused on modifying this forward solution, so learning something about the way it was originally formulated currently has been a useful endeavor. I also found out from the algorithm’s creator that no material exists to help bridge the gap between the idealized and published equations and the optimized and cryptic code. Therefore, I’ve added quite a few comments and docstring improvements to make this easier for the next programmer wrestling with these equations. Later in the summer, I’m hoping to use this knowledge to find the relationship between the cortical surface and the multipolar moment space I started describing in my last post (see the discussion on SSS and spherical harmonics). This should provide a number of benefits that I’ll discuss when I pick this portion of the work back up.

·        For now, I’m going to try to make some headway on the first aim of my project: Maxwell filtering. Again, the Maxwell filter implemented in SSS is just an elegant way to exclude noise from MEEG recordings using physics (see my earlier post about steam bowls and floating frogs for more description or Taulu 2005 for one of the SSS papers). 

·        Last thing: the heads of the MNE-Python project have generously offered to fly me to Paris for a weeklong coding sprint in July! I’m pretty excited to finally meet all of the MNE-Python crew and learn more about how the Europeans view science and research.