-
Notifications
You must be signed in to change notification settings - Fork 21
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
f3e34ad
commit f9cca06
Showing
2 changed files
with
125 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,99 @@ | ||
def define_limits(phs_list): | ||
'''Returns the x-y coordinates limits (min, max) | ||
for a list of persistence diagrams | ||
''' | ||
ph = view.plot.collapse(phs_list) | ||
xlims = [min(np.transpose(ph)[0]), max(np.transpose(ph)[0])] | ||
ylims = [min(np.transpose(ph)[1]), max(np.transpose(ph)[1])] | ||
|
||
return xlims, ylims | ||
|
||
|
||
def persistence_image(ph, norm_factor=None, xlims=None, ylims=None): | ||
'''Create the data for the generation of the persistence image. | ||
If norm_factor is provided the data will be normalized based on this, | ||
otherwise they will be normalized to 1. | ||
If xlims, ylims are provided the data will be scaled accordingly. | ||
''' | ||
from scipy import stats | ||
import numpy as np | ||
|
||
if xlims is None: | ||
xlims = [min(np.transpose(ph)[0]), max(np.transpose(ph)[0])] | ||
if ylims is None: | ||
ylims = [min(np.transpose(ph)[1]), max(np.transpose(ph)[1])] | ||
|
||
X, Y = np.mgrid[xlims[0]:xlims[1]:100j, ylims[0]:ylims[1]:100j] | ||
|
||
values = np.transpose(ph) | ||
kernel = stats.gaussian_kde(values) | ||
positions = np.vstack([X.ravel(), Y.ravel()]) | ||
Z = np.reshape(kernel(positions).T, X.shape) | ||
|
||
if norm_factor is None: | ||
norm_factor = np.max(Z) | ||
|
||
Zn = Z / norm_factor | ||
|
||
return Zn | ||
|
||
|
||
def average_ph_image(images_list): | ||
'''Generates a normalized average image | ||
from a list of images. Careful: images should be | ||
at the same scale (x-y) for appropriate comparison. | ||
''' | ||
average_imgs = images_list[0] | ||
|
||
for im in images_list[1:]: | ||
average_imgs = np.add(average_imgs, im) | ||
|
||
average_imgs = average_imgs / len(images_list) | ||
|
||
return average_imgs | ||
|
||
|
||
def img_diff(Z1, Z2, norm=True): | ||
"""Takes as input two images | ||
as exported from the gaussian kernel | ||
plotting function, and returns | ||
their absolute difference: | ||
diff(abs(Z1 - Z2)) | ||
""" | ||
if norm: | ||
Z1 = Z1 / Z1.max() | ||
Z2 = Z2 / Z2.max() | ||
|
||
return Z1 - Z2 | ||
|
||
|
||
def plot_imgs(img, new_fig=True, subplot=111, title='', xlims=None, ylims=None, cmap=cm.jet, vmin=0, vmax=1., masked=False, threshold=0.01): | ||
'''Plots the gaussian kernel of the input image. | ||
''' | ||
from scipy import stats | ||
from view import common | ||
import numpy as np | ||
|
||
if xlims is None: | ||
xlims = (0,100) | ||
if ylims is None: | ||
ylims = (0,100) | ||
|
||
fig, ax = common.get_figure(new_fig=new_fig, subplot=subplot) | ||
|
||
if masked: | ||
img = np.ma.masked_where((threshold > np.abs(img)), img) | ||
|
||
cax= ax.imshow(np.rot90(img), vmin=vmin, vmax=vmax, cmap=cmap, interpolation='bilinear', extent=xlims+ylims) | ||
|
||
kwargs = {} | ||
|
||
kwargs['xlim'] = xlims | ||
kwargs['ylim'] = ylims | ||
kwargs['title'] = title | ||
|
||
plt.colorbar(cax) | ||
|
||
ax.set_aspect('equal') | ||
|
||
return common.plot_style(fig=fig, ax=ax, aspect='equal') |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,26 @@ | ||
import tmd | ||
import view | ||
|
||
pop1 = tmd.io.load_population(directory1) | ||
pop2 = tmd.io.load_population(directory2) | ||
phs1 = [tmd.methods.get_ph_neuron(n, neurite_type='axon') for n in pop1.neurons] | ||
phs2 = [tmd.methods.get_ph_neuron(n, neurite_type='axon') for n in pop2.neurons] | ||
|
||
# Normalize the limits | ||
xlims, ylims = define_limits(phs1 + phs2) | ||
# Create average images for populations | ||
imgs1 = [persistence_image(p, xlims=xlims, ylims=ylims) for p in phs1] | ||
IMG1 = average_ph_image(imgs1) | ||
imgs2 = [persistence_image(p, xlims=xlims, ylims=ylims) for p in phs2] | ||
IMG2 = average_ph_image(imgs2) | ||
|
||
# You can plot the images if you want to create pretty figures | ||
average_figure1 = plot_imgs(IMG1, title='', xlims=xlims, ylims=ylims, cmap=cm.jet) | ||
average_figure2 = plot_imgs(IMG2, title='', xlims=xlims, ylims=ylims, cmap=cm.jet) | ||
|
||
# Create the diffence between the two images | ||
DIMG = img_diff(IMG1, IMG2) # subtracts IMG2 from IMG1 so anything red IMG1 has more of it and anything blue IMG2 has more of it - or that's how it is supposed to be :) | ||
# Plot the difference between them | ||
diff_image = plot_imgs(DIMG, vmin=-1.0, vmax=1.0) # vmin, vmax important to see changes | ||
# Quantify the absolute distance between the two averages | ||
dist = np.sum(np.abs(DIMG)) |