Commit b7c3d626 authored by Gustavo Pinzon's avatar Gustavo Pinzon
Browse files

sphericalHistogram for vectorial data

parent 25734c8f
Pipeline #68669 passed with stages
in 27 minutes and 35 seconds
......@@ -71,6 +71,11 @@ import math
import matplotlib
import matplotlib.pyplot
import matplotlib.colors as mcolors
import multiprocessing
import progressbar
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
nProcessesDefault = multiprocessing.cpu_count()
VERBOSE = False
......@@ -700,3 +705,356 @@ if __name__ == "__main__":
##subtitle = { "points":"Points",
##"bins":"Bins"}
#)
def createIcosphere(subDiv):
"""
This function creates a icosphere (convex polyhedron made from triangles) starting from an icosahedron (polyhedron with 20 faces) and then making subdivision on each triangle.
The number of faces is 20*(4**subDiv).
Taken from https://sinestesia.co/blog/tutorials/python-icospheres/
Parameters
----------
subDiv : integer
Number of times that the initial icosahedron is divided.
Suggested value: 3
Returns
-------
icoVerts: numberOfVerticesx3 numpy array
Coordinates of the vertices of the icosphere
icoFaces: numberOfFacesx3 numpy array
Indeces of the vertices that compose each face
icoVectors: numberOfFacesx3
Vectors normal to each face
Note
----------
From
https://sinestesia.co/blog/tutorials/python-icospheres/
"""
# 1. Internal functions
middle_point_cache = {}
def vertex(x, y, z):
""" Return vertex coordinates fixed to the unit sphere """
length = numpy.sqrt(x**2 + y**2 + z**2)
return [i / length for i in (x,y,z)]
def middle_point(point_1, point_2):
""" Find a middle point and project to the unit sphere """
# We check if we have already cut this edge first
# to avoid duplicated verts
smaller_index = min(point_1, point_2)
greater_index = max(point_1, point_2)
key = '{0}-{1}'.format(smaller_index, greater_index)
if key in middle_point_cache:
return middle_point_cache[key]
# If it's not in cache, then we can cut it
vert_1 = icoVerts[point_1]
vert_2 = icoVerts[point_2]
middle = [sum(i)/2 for i in zip(vert_1, vert_2)]
icoVerts.append(vertex(middle[0], middle[1], middle[2]))
index = len(icoVerts) - 1
middle_point_cache[key] = index
return index
# 2. Create the initial icosahedron
# Golden ratio
PHI = (1 + numpy.sqrt(5)) / 2
icoVerts = [
vertex(-1, PHI, 0),
vertex( 1, PHI, 0),
vertex(-1, -PHI, 0),
vertex( 1, -PHI, 0),
vertex(0, -1, PHI),
vertex(0, 1, PHI),
vertex(0, -1, -PHI),
vertex(0, 1, -PHI),
vertex( PHI, 0, -1),
vertex( PHI, 0, 1),
vertex(-PHI, 0, -1),
vertex(-PHI, 0, 1),
]
icoFaces = [
# 5 faces around point 0
[0, 11, 5],
[0, 5, 1],
[0, 1, 7],
[0, 7, 10],
[0, 10, 11],
# Adjacent faces
[1, 5, 9],
[5, 11, 4],
[11, 10, 2],
[10, 7, 6],
[7, 1, 8],
# 5 faces around 3
[3, 9, 4],
[3, 4, 2],
[3, 2, 6],
[3, 6, 8],
[3, 8, 9],
# Adjacent faces
[4, 9, 5],
[2, 4, 11],
[6, 2, 10],
[8, 6, 7],
[9, 8, 1],
]
# 3. Work on the subdivisions
for i in range(subDiv):
faces_subDiv = []
for tri in icoFaces:
v1 = middle_point(tri[0], tri[1])
v2 = middle_point(tri[1], tri[2])
v3 = middle_point(tri[2], tri[0])
faces_subDiv.append([tri[0], v1, v3])
faces_subDiv.append([tri[1], v2, v1])
faces_subDiv.append([tri[2], v3, v2])
faces_subDiv.append([v1, v2, v3])
icoFaces = faces_subDiv
# 4. Compute the normal vector to each face
icoVectors = []
for tri in icoFaces:
# Get the points
P1 = numpy.array(icoVerts[tri[0]])
P2 = numpy.array(icoVerts[tri[1]])
P3 = numpy.array(icoVerts[tri[2]])
# Create two vector
v1 = P2 - P1
v2 = P2 - P3
v3 = numpy.cross(v1,v2)
norm = vertex(*v3)
icoVectors.append(norm)
return icoVerts, icoFaces, icoVectors
def plotSphericalHistogram(orientations,
subDiv,
reflection=True,
maxVal=None,
verbose=False,
lim=None,
color=None,
title=None,
saveFigPath=None):
"""
Generates a spherical histogram for vectorial data, binning the data into regions defined by the faces of an icosphere (convex polyhedron made from triangles).
The icosphere is built from starting from an icosahedron (polyhedron with 20 faces) and then making subdivision on each triangle.
The number of faces is 20*(4**subDiv).
Parameters
----------
orientations : Nx3 numpy array
Vectors to be plotted
subDiv : integer, optional
Number of times that the initial icosahedron is divided.
Default: 3
reflection : bool, optional
If true, the histogram takes into account the reflection of the vectors
Default = True.
maxVal : int, optional
Maximum colour-bar limits for bin view.
Default = None (`i.e.`, auto-set)
verbose : bool, optional
Print the evolution of the plot
Defautl = False
color : colormap class, optional
Colormap class from matplotlib module
See 'https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html' for options
Example : matplotlib.pyplot.cm.viridis
Default = matplotlib.pyplot.cm.viridis_r
title : str, optional
Title for the graph
Default = None
saveFigPath : string, optional
Path to save figure to, including the name and extension of the file.
Default = None
Returns
-------
None -- A matplotlib graph is created and shown
"""
# Internal function for binning data into the icosphere faces
def binIcosphere(data, icoVectors, verbose):
# Create counts array
counts = numpy.zeros((len(icoVectors)))
global computeAngle
def computeAngle(i):
# Get the orientation vector
orientationVect = data[i]
# Exchange Z and X position - for plotting
orientationVect = [orientationVect[2], orientationVect[1], orientationVect[0]]
# Create the result array
angle = []
for i in range(len(icoVectors)):
# Compute the angle between them
angle.append(numpy.arccos(numpy.clip(numpy.dot(orientationVect, icoVectors[i]), -1, 1)))
# Get the index
minIndex = numpy.argmin(angle)
return minIndex
# Create progressbar
if verbose:
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(data))
pbar.start()
finishedOrientations = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcessesDefault) as pool:
for returns in pool.imap_unordered(computeAngle, range(len(data))):
# Update the progressbar
finishedOrientations += 1
if verbose:
widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedOrientations, len(data)))
pbar.update(finishedOrientations)
# Get the results
index = returns
# Add the count
counts[index] += 1
return counts
# Get number of points
numberOfPoints = orientations.shape[0]
# Check that they are 3D vectors
if orientations.shape[1] != 3:
print('\nspam.helpers.orientationPlotter.pllotSphericalHistogram: The input vectors are not 3D')
return
# from http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors
norms = numpy.apply_along_axis(numpy.linalg.norm, 1, orientations )
orientations = orientations / norms.reshape( -1, 1 )
# Check if we can reflect the vectors
if reflection:
orientations = numpy.vstack([orientations, -1*orientations])
# Create the icosphere
if verbose:
print('\nspam.helpers.orientationPlotter.pllotSphericalHistogram: Creating the icosphere')
icoVerts, icoFaces, icoVectors = createIcosphere(subDiv)
# Bin the data
if verbose:
print('\nspam.helpers.orientationPlotter.pllotSphericalHistogram: Binning the data')
counts = binIcosphere(orientations, icoVectors, verbose = verbose)
# Now we are ready to plot
if verbose:
print('\nspam.helpers.orientationPlotter.pllotSphericalHistogram: Plotting')
# Create the figure
fig = matplotlib.pyplot.figure()
ax = fig.gca(projection='3d')
if color is None:
cmap = matplotlib.pyplot.cm.viridis_r
else:
cmap = color
norm = matplotlib.pyplot.Normalize(vmin=0, vmax=1)
if maxVal is None:
maxVal = numpy.max(counts)
# Loop through each of the faces
for i in range(len(icoFaces)):
# Get the corresponding radius
radii = counts[i] / maxVal
if radii != 0:
# Get the face
face = icoFaces[i]
# Get the vertices
P1 = numpy.asarray(icoVerts[face[0]])
P2 = numpy.asarray(icoVerts[face[1]])
P3 = numpy.asarray(icoVerts[face[2]])
# Extend the vertices as needed by the radius
P1 = radii*P1 / numpy.linalg.norm(P1)
P2 = radii*P2 / numpy.linalg.norm(P2)
P3 = radii*P3 / numpy.linalg.norm(P3)
# Combine the vertices
vertices = numpy.asarray([numpy.array([0,0,0]), P1, P2, P3])
# Add the points to the scatter3D
ax.scatter3D(vertices[:, 0], vertices[:, 1], vertices[:, 2], s = 0)
# Create each face
face1 = numpy.array([numpy.array(vertices[0]), numpy.array(vertices[1]), numpy.array(vertices[2])])
face2 = [numpy.array(vertices[0]), numpy.array(vertices[1]), numpy.array(vertices[3])]
face3 = [numpy.array(vertices[0]), numpy.array(vertices[3]), numpy.array(vertices[2])]
face4 = [numpy.array(vertices[3]), numpy.array(vertices[1]), numpy.array(vertices[2])]
# Plot each face!
ax.add_collection3d(Poly3DCollection(face1, facecolors=cmap(norm(radii)), linewidths=0.5, edgecolors='k'))
ax.add_collection3d(Poly3DCollection(face2, facecolors=cmap(norm(radii)), linewidths=0.5, edgecolors='k'))
ax.add_collection3d(Poly3DCollection(face3, facecolors=cmap(norm(radii)), linewidths=0.5, edgecolors='k'))
ax.add_collection3d(Poly3DCollection(face4, facecolors=cmap(norm(radii)), linewidths=0.5, edgecolors='k'))
# Extra parameters for the axis
ax.set_box_aspect([1,1,1])
matplotlib.pyplot.xlim(-1.1,1.1)
matplotlib.pyplot.ylim(-1.1,1.1)
ax.set_zlim(-1.1,1.1)
ax.view_init(25, 45)
# Set the colorbar
norm = matplotlib.colors.Normalize(vmin=0, vmax=maxVal)
sm = matplotlib.pyplot.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
matplotlib.pyplot.colorbar(sm, label="Counts")
# Remove the ticks labels and lines
ax = matplotlib.pyplot.gca()
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.zaxis.set_ticklabels([])
for line in ax.xaxis.get_ticklines():
line.set_visible(False)
for line in ax.yaxis.get_ticklines():
line.set_visible(False)
for line in ax.zaxis.get_ticklines():
line.set_visible(False)
# Title
if title is not None:
ax.set_title( str(title) + "\n" )
matplotlib.pyplot.tight_layout()
if saveFigPath is not None:
matplotlib.pyplot.savefig( saveFigPath )
else:
matplotlib.pyplot.show()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment