Topoplots

[Tokazama]

@ElectronicTeaCup brought up topoplots at sam81/BrainWave.jl#1.
I'll borrow an example from the Makie docs to give a more concrete place to start. This example provides gradients and contours (see the plot here).

using AbstractPlotting
x = LinRange(-1, 1, 20)
y = LinRange(-1, 1, 20)
z = x .* y'
vbox(
    contour(x, y, z, levels = 50, linewidth =3),
    contour(x, y, z, levels = 0, linewidth = 0, fillrange = true)
)

If we instead do something like the following we can get a more classic topological for neurophysiology (I think there needs to be some adjustments to the second plot's attributes make the contour lines more visible).

scene = Scene()
contour!(scene, x, y, z, levels = 0, linewidth = 0, fillrange = true)
contour!(scene, x, y, z, levels = 50, linewidth =3),

It might be good to look at how GeoMakie accomplishes the different projection.

Edit: One last thought. I'm using CoordinateTransformations.jl in NeuroCore because it's fast, is used in other packages that do spatial registration in Julia, and I think it may get incorporated into the plots ecosystem.

[AsimHDar]

Cool, thanks for starting this off. Btw, I saw that EEG.jl has something that plots the topography too, I'll check that out to see the general strategy. From the readme:

using EEG

t = read_VolumeImage("example.dat")

p = EEG.plot(t, c=:darkrainbow)
p = EEG.plot(t, find_dipoles(t), l = "Peak Activity", c=:black)

[Tokazama]

That's a good place to start. It looks like this is a manipulation on a scatter plot using orthogonal slices. If that's the case then the only thing I'm missing right now to produce this is the dipole estimation, because I already have some Makie code for this layout and these views.

That being said, this might not by the topological plot you want given this is already estimating the dipoles and inferring the signals' location. If you just want a contour plot using the surface electrodes it would require a slightly different approach.

[AsimHDar]

To get a better idea of the problem I was trying to adapt some code from Mike X Cohen's book and 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("../data/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
)

gridX = repeat(interpX,1,100)'
gridY = repeat(collect(interpY),1,100)

# 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);

% let's look at these matrices
figure
subplot(121)
imagesc(gridX)

subplot(122)
imagesc(gridY)


% 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');

[Tokazama]

I won't be at a computer where I can test this until late tonight. Where does this fail?

[AsimHDar]

It produces an erroneous contour plot, which isn't plausible when compared to the MATLAB version.

[AsimHDar]

Julia on the left and Matlab output on the right

[AsimHDar]

Okay I was messing up the Cartesian coordinates before (updated the pasted code) but I still get this

[Tokazama]

I actually think that's a pretty good start. I haven't figured it out yet but I think once that first image is smoothed and we apply a mask the it will be fine. Im still playing with it so if you beat me too it let me know

[Tokazama]

I'm using Makie right now and I'm still not getting anything like I'm suppose to. I think there should be a simpler way of doing all of this. I'm going to work on some related stuff to help with getting this package ready over the next several days. I'll come back to this issue late next week. Also CoordinateTranaformations has spherical and Cartesian conversion

[AsimHDar]

Yeah, me neither. I couldn't get topoplots to work but I think the output looks something like this:

Also CoordinateTranaformations has spherical and Cartesian conversion

Awesome, will check out. Thanks!

[Tokazama]

I'm definitely working on this still. I know it's obnoxious having to wait for such a critical feature but I'm more worried about releasing broken code at the moment.

Hopefully I can get a bunch of satisfactory tests passing in the branch of NeuroCore I'm working on right now and release that this week. It will have a format that is as close as possible to native plotting data types while still preserving all of the nice performance of AbstractArray types for linear algebra.

It seems like the issue with your compared output (julia left matlab right pic) is that you 1) don't have the circle masked out. There needs to be a threshold there that doesn't interpolate beyond the circle's perimeter. The other problem is just how much smoothing is occurring. You have interpolated the grid values but it doesn't look like it smooths this out very much.

The reason I'm not just diving into solving those two problems is because I wonder if it might be better to figure out how to plot things in a polar space instead of cartesian space. This would theoretically solve both problems, but it depends on getting the plotting ecosystem set up for polar plots.

[AsimHDar]

Yeah that makes sense, and even the MATLAB version is hardly a finished product—it looks pretty bare bones. For the nicer ones, such as topoplot, or topomap (which is a 2600 line beast) the code base is quite extensive.

This makes me think: can you seamlessly run python packages like MNE via julia? (I have not been able to do it) (it wasn't fun—just tried again)

[Tokazama]

This makes me think: can you seamlessly run python packages like MNE via julia? (I have not been able to do it)

Yes, but I don't think it will take 2600 lines of code to do this. If my talk gets accepted to JuliaCon I'll give more details, but for now it probably suffices to say that what you are seeing is one of the biggest reasons neuroscience should be done in Julia.

That being said. If you need something right now I'd just use PyCall to access that function.

[AsimHDar]

Yes, but I don't think it will take 2600 lines of code to do this. If my talk gets accepted to JuliaCon I'll give more details, but for now it probably suffices to say that what you are seeing is one of the biggest reasons neuroscience should be done in Julia.

The argument being that the same stuff can be done with less lines?

I don't need it urgently, this is part of an exercise where I'm trying to port code from a book I'm reading (which has MATLAB code) to Julia. Just trying to get the hang of Julia and what it can do, so that I can scrap MATLAB.

[Tokazama]

The argument being that the same stuff can be done with less lines?

I hate to just lay out stuff without citing examples so this sounds a little hand wavy, but nearly all of computational neuroscience is just using what exists in other fields with semantic layer tying it together. There's no reason they needed to write that code except for the fact that python libraries don't play nice with each other. Therefore, they had to reinvent something that surely exists elsewhere. It's not a maintainable way to write code, having large incompatible libraries maintained by experts in only marginally related fields.