Topoplots

[Tokazama]

@ElectronicTeaCup brought up topoplots at https://github.com/sam81/BrainWave.jl/issues/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.

[AsimHDar]

Still no luck with this, I tried GR.tricont to plot and it looks a bit better but nothing like the final result.

This might be a good lead (but is in python).

[behinger]

I would be interested to help. I have a package to analyse EEG/fMRI data with multiple (mixed) regression unfold.jl (not yet on the public registry).

But I am at a point where I need good visualization tools :)

I know the eeglab topoplot code a bit, it is a mess due to historical reasons. There is a ggplot based implementation in R here and one in MNE here.

I'd like to get to know makie better - can you give an update on the status here and with the NeuroCore/NeuroSignals in general?

[Tokazama]

NeuroSignals is a place holder right now for where I thought it would be nice to have something for electrophysiology and fMRI based signal analysis. I haven't done anything big with it yet because I wanted it to be a much more collaborative effort than all of the other stuff I've been doing.

NeuroCore is SOOOO close to being ready for a new release. My plans are:

• [ ] Finish a bunch of bug fixes for AxisIndices.jl (array interface that is a lot nicer than AxisArrays.jl) - Should be done today
• [ ] Fix up TimeAxes.jl (which extends AxisIndices.jl with a bunch of time related convenience methods) - Hopefully done today but if not will be done tomorrow
• [ ] Submit a new PR for NeuroCore for review (tomorrow or Wednesday).

There's a lot of other stuff to do but that's the best estimations I make for things right now. When I have the NeuroCore PR ready I'll include more details on the overall intention and usage of the package. For now it's probably sufficient to say that I've been creating datatypes with the intention that they can be generically be passed to methods for statistical modeling or plotting.

[behinger]

Ah now I see, the documentation is a bit hidden in the issues - sorry could have checked.

my unfold.jl seems to fit very well into the "statistical modeling for neuroscience" bit. I am happy to adopt my required datastructures to something general (well, it is an AbstractArray and a DataFrame right now, so I guess not much work needed ;))

---

coming back to the issue at hand here, any progress in plotting topos?

[AsimHDar]

@behinger Nice to have you here! From my end I'm still trying to make a simple interpolation based EEG topography plot. My attempts haven't been very successful and every time I open up topoplot.m to consider porting it, my legs get a little shaky. Would be willing to help out if someone takes the lead, I'm doing this as part of a time-frequency analysis course—for which code is only available in MATLAB.

[Tokazama]

@ElectronicTeaCup, if you want to get get started the best thing would be to familiarize yourself with Makie.jl and GeometryBasics.jl. Probably the biggest thing I'll be introducing to NeuroCore relevant to this issue is pairing closely with GeometryBasics.jl, which provides the underlying types for Makie.jl plotting.

We also need to decide exactly what we want for a topoplot. This discussion started with code that mapped cartesian indices to a polar projection (which you should be able to accomplish with CoordinateTransformations.jl). However, how we want to "interpolate" points between electrodes should be discussed. Do we actually want to use interpolation or estimate dipole coordinates and such?

[behinger]

I don't think dipole coordinates are important here. I think how to transform the 3D channel locations to 2D is an important thing and then with what way to interpolate the points. I will have a look at it. For me GeometryBasics.jl seems too basic right now

[Tokazama]

GeometryBasics.jl definitely isn't going to solve the problem, but it's the raw types that Makie uses for plotting so if we can work within compatible types then there's no need to convert between types. For example, I think this is where interpolation is performed on stuff within Makie. I've yet to figure out the full plotting pipeline there yet though.

I've used CoordinateTransformations for this sort of thing in the past and it's pretty straightforward.

[behinger]

Ok, I looked at it a bit yesterday & today.

• I think interpolating manually first and then using Makie is easier - we could also interpolate in 3D and then move to 2D. Best of course would be to let Makie handle all the interpolation. But at the moment I don't think it has that functionality
• Regarding coordinateTransformation: As you can see, Eeglab somehow maps the X,Y,Z Coordinates to 2D coordinates in a different way than just dropping Z. I wonder if such a perspective transformation can be handled by coordinateTransformation.jl?
• I tried using Dierckx (if you want to use ScatteredInterpolation, you could try Gaussian instead of Multiquadratic?). I tried Interpolations.jl but couldn't get the unregular grid function to work. Trouble with Dierckx is that I can't reduce the amount of smoothing and the smoothing parameter is defined in a (for me) unintuitive fashion.

This is what I came up with

Of course the head+nose+ears are missing,

using Dierckx
using Makie
using ColorSchemes

# Load eeglab dataset
include("test/debug_readEEGlab.jl")
filename = "C:/Users/behinger/Downloads/1_P3_shifted_ds_reref_ucbip_hpfilt_ica_corr_cbip_elist.set"
data,srate,evts_df,chanlocs_df = import_eeglab(filename)
##
X= Float64.(chanlocs_df.X[1:end-7]) # cast to Float, remove some irrelevant channels
Y= -Float64.(chanlocs_df.Y[1:end-7]) # Reverse Y, not sure why
D = data[1:end-7,10000] #choose random timepoint

# calculate 2D spline (Dierckx)
spl = Spline2D(Y,X,D,kx=3,ky=3,s=800) # s = smoothing parameter, kx/ky = spline order

# get extrema and extend by 20%
xlim = extrema(X) .+ abs.(extrema(X)).*[-.2, .2]
ylim = extrema(Y).+ abs.(extrema(Y)).*[-.2, .2]

# generate and evaluate a grid
xg = range(xlim[1],stop=xlim[2],  step=0.5)
yg = range(ylim[1],stop=ylim[2], step=0.5)
v = evalgrid(spl,yg,xg) # evaluate the spline at the grid locs

# remove everything in a circle (should be ellipse at some point)
ix = sqrt.([i.^2+j.^2 for i in yg, j in xg]).>90
v[ix] .= NaN
# nicer colormaps
cmap = ColorSchemes.vik;
# generate scene+subscene to keep aspect ratio
scene = Scene(resolution = (1000, 1000))
sub = Scene(scene, IRect(0, 0, 1000, 1000))

heatmap!(sub,yg,xg,v,colormap=cmap, show_axis = false)
contour!(sub,yg,xg,Float64.(v),linewidth=3,colormap=cmap)
scatter!(sub,Y,X,markersize=3,color="white",strokewidth=3,strokecolor="black") # add electrodes
[text!(sub,String(chanlocs_df.labels[i]),position=(Y[i],X[i]),textsize=7) for i in 1:length(X)]; # add labels

[AsimHDar]

That's awesome! Thanks for the informative commenting, it looks great.

[palday]

Regarding the projection from 3D to 2D, see here

[sappelhoff]

just chiming in to advertise some work I did, hoping that it can be helpful: https://github.com/sappelhoff/eeg_positions

in that repo I provide the (Python) code to compute standard EEG electrode positions (10-5, 10-10, 10-20 systems) based on a spherical head. The main algorithm is here: https://github.com/sappelhoff/eeg_positions/blob/863ca5e0160f00f78a6ab944e1f29e8559277bf1/eeg_positions/utils.py#L10

and could be easily ported to Julia. --> that way, users who have only (standard) channel names and no coordinates could also plot topo maps, because the coordinates can be "computed on the fly".

[behinger]

Just to give a brief update. I put everything together in a bunch of functions here: https://github.com/unfoldtoolbox/UnfoldMakie.jl/blob/topoplot/src/topoplot.jl

With some random test-case here: https://github.com/unfoldtoolbox/UnfoldMakie.jl/blob/topoplot/test/topoplot.jl

I still run into some problems when plotting many topoplots, i.e. first topoplot takes ~100ms, 20th something like 5s.

[Tokazama]

I've waited far too long to get this working in hopes of a foundation that won't lead to breaking changes. I think we're going to have to just accept that this package will change frequently along with the plotting ecosystem. @behinger , I sent you an invite. Feel free to start making PRs and will get people working together on this.

[behinger]

I don't think I can invest time right now to do this unfortunately.

There is also Neuroimaging.jl (https://github.com/rob-luke/Neuroimaging.jl) which got a recent development boost by @rob-luke. Potentially, there are synergies?

[Tokazama]

That's totally understandable. Everyone's welcome to make PR's BTW. Every time I make a good neuro related plot that I want to set up here my PI's get excited and give me more work; then I'm not taking time to formalize it into a public code base. Perhaps wider involvement will help mature some of that code.

[palday]

@Tokazama can you init the repo with a README and the MIT license?

[Tokazama]

I gave it the full PkgTemplates treatment and added a ColPrac badge so we're on a common understanding for how to work together here.

If we can really get this going we may still need a .zenodo.json or something, but I'll wait until we get there.

[palday]

@Tokazama thanks! the license is also a prereq for me contributing things as part of my dayjob :smile:

[Tokazama]

@Tokazama thanks! the license is also a prereq for me contributing things as part of my dayjob 😄

Totally get it. I'm sitting on tons of code that I've been strongly discouraged from making public until it results in a publication. Problem is that I'm terrible at writing.

Just an PSA on why this project is a bit more difficult to make push on:

Ideally we have a minimal recipe system and maybe some geometric types so we don't force a specific plotting backend. However, there's GeometryBasics.jl, GeometryTypes.jl, and Meshes.jl. They all are trying to accomplish the same thing. I have deliberately not gotten involved in that battle but others are welcome to try to sort out that mess. In order to avoid horribly breaking changes we probably shouldn't get too dependent on a particular set of user facing types.

[SimonDanisch]

However, there's GeometryBasics.jl, GeometryTypes.jl, and Meshes.jl. They all are trying to accomplish the same thing

It's admittedly bit messy, but I don't think it's huge problem. You can cross out GeometryTypes from that list at least, and GeometryBasics vs Meshes is also not too hard: If you want to do things with Makie, just use GeometryBasics, if you need to algorithms from Meshes, use that ;) It's also relatively easy to convert between those.

If the Makie topoplot recipe is at the right place here, I'll be able to help maintain it as part of this package and maybe contribute some other recipes.

Is there some information anywhere about the general direction of this package?

[palday]

@SimonDanisch the package is just a skeleton at this point, so we would probably just follow your recommendations for organization. I'm assuming just depend on Makie for the main package and maybe use Percy + CairoMakie for testing?

[Tokazama]

I was planning on it being more of a dedicated set of recipes building on related packages for visualizing neuroscience data. My most immediate action plan is:

1. Figure out this issue https://github.com/JuliaArrays/SpatioTemporalTraits.jl/pull/11
2. Some basic traits related to anatomy that will be useful for labeling plots of neuro data https://github.com/JuliaNeuroscience/AnatomicalTraits.jl.
3. Package my old brain visualization code in PRs here

Plotting is not my strong suit but I constantly need it, so I'm happy to work with others to make this thing really work.

[likanzhan]

For your information, here is what I've achieved with the triangulate function in GMT.jl. The code is here.

• Julia rewritten version

• Original Matlab version

[behinger]

There is a new package for TopoPlots in Makie: https://github.com/MakieOrg/TopoPlots.jl