Trouble making topographical EEG interpolation

Hi everyone. I’m trying to make an EEG topographical plot, something like this:

image

The steps, in general are (without the aesthetics and circular appearance):

  1. Turn polar coordinates to Cartesian coordinates
  2. Interpolate data onto coordinates
  3. …?

So far I’ve tried to keep it pretty simple, but still can’t manage. I was trying to adapt some MATLAB code but just couldn’t manage to implement the steps. Can you see something that I’m missing?

using MAT, Plots, Random, Measures
using FFTW, DSP
using ScatteredInterpolation
using Statistics
using StatsBase

plotly()

# I/O for the EEG mat file, using MAT package
EEG = matread("sampleEEGdata.mat")
EEG = EEG["EEG"]

function pol2cart(th,r,z=1)
    x = r.*cos.(th)
    y = r.*sin.(th)
    return (x, y, z)
end

# Introduction to Topographical Plots

timepoint2plot = 100 # in ms
trial2plot     = Int(sample(1:EEG["trials"]))
color_limit = 20

# Convert time point from ms to index and get relevant EEG data
_, timepointidx = findmin(dropdims(abs.(EEG["times"].-timepoint2plot), dims=1))
timepoint_trial_data = EEG["data"][:,timepointidx, trial2plot][:,1]

# First step is to get X and Y coordinates of electrodes
# These must be converted to polar coordinates
th = (π/180) .*EEG["chanlocs"]["theta"][1,:]
radi = EEG["chanlocs"]["radius"][1,:]
electrode_locs_X, electrode_locs_Y = pol2cart(th, radi)

# Interpolate to get nice surface
interpolation_level = 100 # you can try changing this number, but 100 is probably good
interpX = range(
    minimum(electrode_locs_X),
    maximum(electrode_locs_X),
    length=interpolation_level
)
interpY = range(
    minimum(electrode_locs_Y),
    maximum(electrode_locs_Y),
    length=interpolation_level
)

# testing ScatteredInterpolation.jl
points = [electrode_locs_X' ; electrode_locs_Y']
samples = timepoint_trial_data
n = 100
x = interpX
y = interpY
X = repeat(x, n)[:]
Y = repeat(y', n)[:]
gridPoints = [X Y]'

itp = interpolate(Multiquadratic(), points, samples)
interpolated = evaluate(itp, gridPoints)
gridded = reshape(interpolated, n, n)
p2 = contour(x, y, gridded)

The matlab code was:

% First step is to get X and Y coordinates of electrodes
% These must be converted to polar coordinates
th=pi/180*[EEG.chanlocs.theta];
[electrode_locs_X,electrode_locs_Y] = pol2cart(th,[EEG.chanlocs.radius]);

% interpolate to get nice surface
interpolation_level = 100; % you can try changing this number, but 100 is probably good
interpX = linspace(min(electrode_locs_X),max(electrode_locs_X),interpolation_level);
interpY = linspace(min(electrode_locs_Y),max(electrode_locs_Y),interpolation_level);

% meshgrid is a function that creates 2D grid locations based on 1D inputs
[gridX,gridY] = meshgrid(interpX,interpY);

% now interpolate the data on a 2D grid
interpolated_EEG_data = griddata(electrode_locs_Y,electrode_locs_X,double(squeeze(EEG.data(:,timepointidx,trial2plot))),gridX,gridY);
contourf(interpY,interpX,interpolated_EEG_data,100,'linecolor','none');

I get the following: (Julia plot on the left, MATLAB on the right)

Does anyone have any better approaches on getting the final product?

You might have some luck with GR.jl’s tricont or trisurf functions.

Out of interest I bought one of these ECG sensors https://www.amazon.co.uk/gp/product/B08216YR9H/ref=ppx_yo_dt_b_search_asin_title?ie=UTF8&psc=1

It should work with a Raspberry Pi.
May I ask what the polar plot represents?

Thanks for the suggestion! Getting this with tricont

image

Any ideas how I can fill this?

The plot represents electrical brain activity (Electroencephalography, EEG) on the surface of a human subject during a experimental task. ECG (the sensor you have) in a similar vein, is for detecting electrical signals read out from the heart and is used clinically by medical doctors.

@ElectronicTeaCup My bad. Entirely. I used to work in a medical physics department and I hang my head in shame.
Please tell us more about your interesting work. What is the task the subjects are asked to do?

Hehe, no issues. This particular data is public, which I’m using for text-book based course on time-frequency analyses so I’m not super sure what the experiment was. I’m hoping to apply the same time-frequency methods to MEG recordings that I will be taking during a simple task which involves passively listening to tones. The idea is to get a general idea of how sounds can form a temporally intact memory trace in the auditory cortex using modelling approaches and empirical data.

Do we have an update on this topic?

@likanzhan Sorry about the late reply, but there was some more discussion in this thread, with some functioning code to make something like this:

FWIW, if the human scalp can be approximated by an ellipsoid, then you can try GMT’s greenspline interpolation which uses a generalized Green’s function for spherical surface splines in tension. See one example here and the technical reference here.

1 Like

usually it’s treated as an ellipsoid, although it’s becoming a bit more common to create a unique 3d mesh representing each persons head structure.

1 Like

Thanks. Based on others’ suggestions and the plot_topography function written in Matlab, I think triangulate methods are better. Here are what I have achieved via the triangulate function in GMT. The code I used to do the plotting is here.

  • Julia rewritten version

  • Original Matlab version

Using GMT and the tip from @joa-quim here, you can produce a contour plot like this (OP’s script has random sampling, so input is different):

# GMT GRIDDING (Smith and Wessel [1990] Splines in tension):
using GMT
data = [electrode_locs_X electrode_locs_Y timepoint_trial_data]
x = y = LinRange(-0.65, 0.65, 100)
G = GMT.surface(data, R=(extrema(x)..., extrema(y)...), inc=(step(x), step(y)), T=0.1)

# GMT PLOTTING
circ0 = reduce(vcat, 0.65*[[cos(θ) sin(θ)] for θ in 0:0.1:2π])    # circle
GMT.contourf(G, clip=circ0, fmt=:png, show=true)

NB:
the spline interpolation method used acts along “flat” x-y support coordinates, but as mentioned it is also possible to use spherical surface splines acting along an ellipsoid

If someone puts together a proposal in a PR I’d be happy to review it for NeuroPlots.

@Zach_Christensen Do you think it is potentially useful to others if I push the function I wrote to NeuroPlots?

could you elaborate why you prefer triangulation methods for topoplots?

Cheers, Bene

@behinger There are some discussion concerning their difference here. As we can see from the four examples, there exist some tiny difference between them. My intuition is that it is more reasonable to suppose the hotpot is in the center of the surface, rather than at the edge of the surface. But I’m no expert on either EEG or GMT, please correct me if I’m wrong. Furthermore, we can easily implement different methods, as the case griddata in Matlab, for the users to select.

Sincerely,
Likan

Let’s start by getting something working. This is long overdue. Then we can work on how to facilitate different interpolation methods.

There is a new package TopoPlots.jl that should resolve this issue. If you have ideas to improve, we are happy to hear them!

grafik

4 Likes