Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

orientationPlotter.py 42.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
"""
Library of SPAM functions for plotting orientations in 3D
Copyright (C) 2020 SPAM Contributors

This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
more details.

You should have received a copy of the GNU General Public License along with
this program.  If not, see <http://www.gnu.org/licenses/>.
"""


20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
# 
# Edward Ando 17/11/2011
#
# Attempt at programming rose plots myself.

# Assuming that the contacts are kind of flat, and so the normal vector pointing up from them is correct
# which means using Visilog's Orientation 2 vectors.

# Modified in order to allow vector components to be taken as input directly

# 2012.08.27 - adding ability to output projected components, in order to allow plotting
#   with gnuplot

# Completely new version with objective of:
#   - reading in different formats (3D coordinates (x,y,z,), Spherical Coordinates, Cylindrical)
#   - outputting any other format
#   - Different projections onto the plane (Lambert, Stereo, direct, etc...)
#   - Possibility to colour the negative part of the projected sphere in a different colour
#   - Projection point-by-point or binned with Hugues Talbot and Clara's code, which
#       gives the very convenient cutting of the circle into in equal area parts:
#         Jaquet, C., Andò, E., Viggiani, G., & Talbot, H. (2013).
#         Estimation of separating planes between touching 3D objects using power watershed.
#         In Mathematical Morphology and Its Applications to Signal and Image Processing (pp. 452-463).
#         Springer Berlin Heidelberg.

# Internal data format will be x,y,z sphere, with x-y defining the plane of projection.

# 2015-07-24 -- EA and MW -- checking everything, there were many bugs -- the points were only plotted to an extent of 1 and not radiusMax
#               created plotting function, which allows radius labels to be updated.

# 2016-11-08 -- MW -- binning was still erroneous. an angularBin (lines 348ff) was usually put to the next higher bin, so not rounded correctly
#               rounding now with numpy.rint and if the orientation doesn't belong to the last bin, it has to be put in the first one
#                 as the first bin extends to both sides of 0 Degrees. -> check validation example in lines 227ff

# 2017-06-01 -- MW -- updating for the current numpy version 1.12.1 and upgrading matplotlib to 2.0.2

# 2017-06-26 -- MW -- there still was a small bug in the binning: for the angular bin numpy.rint was used -- replaced it with numpy.floor (line 375)
57
58
59
60
61
62
63
#                     benchmarking points are in lines 217ff.

# 2018-02-19 -- EA -- changes in the new matplot version revealed a problem. Eddy solved it for the spam client, MW modified this one here.

# 2018-04-20 -- MW -- adding relative bin counts -- to normalise the bincounts by the average overall bin count.
#                     enables to plot many states with the same legend!

64
from __future__ import print_function
65
66
67
68
69
70
71
72

import os
import sys
import numpy
import random
import math
import matplotlib
import matplotlib.pyplot
73
import matplotlib.colors as mcolors
74
75
76
77
78
import multiprocessing
import progressbar
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection

nProcessesDefault = multiprocessing.cpu_count()
79
80
81
82
83
84
85
86

VERBOSE = False

# ask numpy to print 0.000 without scientific notation
numpy.set_printoptions(suppress=True)
numpy.set_printoptions(precision=3)


87
def SaffAndKuijlaarsSpiral(N):
88
    """
Edward Andò's avatar
Edward Andò committed
89
90
    There is no analytical solution for putting equally-spaced points on a sphere.
    This spiral gets close.
91

Edward Andò's avatar
Edward Andò committed
92
93
    Parameters
    ----------
94

Edward Andò's avatar
Edward Andò committed
95
96
    N : integer
        Number of points to generate
97

Edward Andò's avatar
Edward Andò committed
98
99
100
101
    Returns
    -------
    orientations : Nx3 numpy array
        Z,Y,X unit vectors of orientations for each point on sphere
102

Edward Andò's avatar
Edward Andò committed
103
104
105
106
    Note
    ----------
    For references, see
    http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
107

Edward Andò's avatar
Edward Andò committed
108
109
    Which in turn was based on
    http://sitemason.vanderbilt.edu/page/hmbADS
110

Edward Andò's avatar
Edward Andò committed
111
112
113
    From
    Rakhmanov, Saff and Zhou: **Minimal Discrete Energy on the Sphere**, Mathematical Research Letters, Vol. 1 (1994), pp. 647-662:
    https://www.math.vanderbilt.edu/~esaff/texts/155.pdf
114

Edward Andò's avatar
Edward Andò committed
115
116
    Also see discussion here
    http://groups.google.com/group/sci.math/browse_thread/thread/983105fb1ced42c/e803d9e3e9ba3d23#e803d9e3e9ba3d23%22%22
117
118
    """
    M = int(N)*2
Edward Andò's avatar
Edward Andò committed
119

120
    s         = 3.6 / math.sqrt(M)
Edward Andò's avatar
Edward Andò committed
121

122
123
    delta_z   = 2 / float(M)
    z         = 1 - delta_z/2
Edward Andò's avatar
Edward Andò committed
124

125
    longitude = 0
Edward Andò's avatar
Edward Andò committed
126

127
    points = numpy.zeros( (N,3) )
Edward Andò's avatar
Edward Andò committed
128

129
130
131
132
133
134
135
136
137
138
139
    for k in range( N ):
        r           = math.sqrt( 1 - z*z )
        points[k,2] = math.cos( longitude ) * r
        points[k,1] = math.sin( longitude ) * r 
        points[k,0] = z
        z           = z - delta_z
        longitude   = longitude + s/r
    return points



140
def _projectOrientations(projection, coordSystem, vector):
141
142
143
144
145
    #           1. type of projection
    #                        2. coordinate system (XYZ, or spherical)
    #                                    3. input 3-unit vector, or array:
    #                                       for x-y-z that order is maintained
    #                                       for spherical: r, theta (inclination), phi (azimuth)
Edward Andò's avatar
Edward Andò committed
146

147
    # Returns: projection_xy, projection_theta_r
Edward Andò's avatar
Edward Andò committed
148

149
150
    projection_xy_local      = numpy.zeros( (2) )
    projection_theta_r_local = numpy.zeros( (2) )
Edward Andò's avatar
Edward Andò committed
151

152
153
154
    if coordSystem == "spherical":
        # unpack vector 
        r, theta, phi = vector
Edward Andò's avatar
Edward Andò committed
155

156
157
158
        x = r * math.sin( theta ) * math.cos( phi )
        y = r * math.sin( theta ) * math.sin( phi )
        z = r * math.cos( theta )
Edward Andò's avatar
Edward Andò committed
159

160
161
162
163
164
165
166
167
168
169
170
171
    else:
        # unpack vector 
        x, y, z = vector
        # we're in cartesian coordinates, (x-y-z mode) Calculate spherical coordinates
        # passing to 3d spherical coordinates too...
        # From: https://en.wikipedia.org/wiki/Spherical_coordinate_system
        #  Several different conventions exist for representing the three coordinates, and for the order in which they should be written.
        #  The use of (r, θ, φ) to denote radial distance, inclination (or elevation), and azimuth, respectively, is common practice in physics, 
        #   and is specified by ISO standard 80000-2 :2009, and earlier in ISO 31-11 (1992).
        r     = numpy.sqrt( x**2 + y**2 + z**2 )
        theta = math.acos( z / r )   # inclination
        phi   = math.atan2( y, x )   # azimuth
Edward Andò's avatar
Edward Andò committed
172
173

    if projection == "lambert":          # dividing by sqrt(2) so that we're projecting onto a unit circle
174
175
        projection_xy_local[0] = x*( math.sqrt(2/(1+z)) )
        projection_xy_local[1] = y*( math.sqrt(2/(1+z)) )
Edward Andò's avatar
Edward Andò committed
176

177
178
179
180
181
        # sperhical coordinates -- CAREFUL as per this wikipedia page: https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection
        #   the symbols for inclination and azimuth ARE INVERTED WITH RESPEST TO THE SPHERICAL COORDS!!!
        projection_theta_r_local[0] = phi
        #                                                  HACK: doing math.pi - angle in order for the +z to be projected to 0,0
        projection_theta_r_local[1] = 2 * math.cos( ( math.pi - theta ) / 2 )
Edward Andò's avatar
Edward Andò committed
182

183
184
185
186
        # cylindrical coordinates
        #projection_theta_r_local[0] = phi
        #projection_theta_r_local[1] = math.sqrt( 2.0 * ( 1 + z ) )

Edward Andò's avatar
Edward Andò committed
187
188


189
190
191
    if projection == "stereo":
        projection_xy_local[0] = x / ( 1 - z )
        projection_xy_local[1] = y / ( 1 - z )
Edward Andò's avatar
Edward Andò committed
192

193
194
195
196
197
        # https://en.wikipedia.org/wiki/Stereographic_projection uses a different standard from the page on spherical coord Spherical_coordinate_system
        projection_theta_r_local[0] = phi
        #                                        HACK: doing math.pi - angle in order for the +z to be projected to 0,0
        #                                                                             HACK: doing math.pi - angle in order for the +z to be projected to 0,0
        projection_theta_r_local[1] = numpy.sin( math.pi - theta ) / ( 1 - numpy.cos( math.pi - theta ) )
Edward Andò's avatar
Edward Andò committed
198

199
200
201
202
203
    if projection == "equidistant":
        # https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection
        # TODO: To be checked, but this looks like it should -- a straight down projection.
        projection_xy_local[0] = math.sin( phi )
        projection_xy_local[1] = math.cos( phi )
Edward Andò's avatar
Edward Andò committed
204

205
206
        projection_theta_r_local[0] = phi
        projection_theta_r_local[1] = numpy.cos( theta - math.pi/2 )
Edward Andò's avatar
Edward Andò committed
207

208
209
210
    return projection_xy_local, projection_theta_r_local


211
212
213
214
215
216
217
218
219
220
221
222
223
224
def plotOrientations(orientations_zyx,
                     projection="lambert",
                     plot="both",
                     binValueMin=None,
                     binValueMax=None,
                     binNormalisation = False,
                     numberOfRings = 9,
                     # the size of the dots in the points plot (5 OK for many points, 25 good for few points/debugging)
                     pointMarkerSize = 8,
                     cmap = matplotlib.pyplot.cm.RdBu_r,
                     title = "",
                     subtitle = {"points":"",
                                 "bins":""},
                     saveFigPath = None ):
225
    """
Edward Andò's avatar
Edward Andò committed
226
227
    Main function for plotting 3D orientations.
    This function plots orientations (described by unit-direction vectors) from a sphere onto a plane.
228

Edward Andò's avatar
Edward Andò committed
229
    One useful trick for evaluating these orientations is to project them with a "Lambert equal area projection", which means that an isotropic distribution of angles is projected as equally filling the projected space.
230

Edward Andò's avatar
Edward Andò committed
231
232
    Parameters
    ----------
233
234
235
236
237
238
        orientations : Nx3 numpy array of floats
            Z, Y and X components of direction vectors.
            Non-unit vectors are normalised.

        projection : string, optional
            Selects different projection modes:
239
                **lambert** : Equal-area projection, default and highly reccommended. See https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

                **equidistant** : equidistant projection

        plot : string, optional
            Selects which plots to show:
                **points** : shows projected points individually
                **bins** : shows binned orientations with counts inside each bin as colour
                **both** : shows both representations side-by-side, default

        title : string, optional
            Plot main title. Default = ""

        subtitle : dictionary, optional
            Sub-plot titles:
                **points** : Title for points plot. Default = ""
                **bins** : Title for bins plot. Default = ""

Edward Andò's avatar
Edward Andò committed
257
258
259
260
261
262
263
        binValueMin : int, optional
            Minimum colour-bar limits for bin view.
            Default = None (`i.e.`, auto-set)

        binValueMax : int, optional
            Maxmum colour-bar limits for bin view.
            Default = None (`i.e.`, auto-set)
264

265
266
267
        binNormalisation : bool, optional
            In binning mode, should bin counts be normalised by mean counts on all bins
            or absoulte counts?
268
269
270
271
272
273
274
275
276
277
278

        cmap : matplotlib colour map, optional
            Colourmap for number of counts in each bin in the bin view.
            Default = ``matplotlib.pyplot.cm.RdBu_r``

        numberOfRings : int, optional
            Number of rings (`i.e.`, radial bins) for the bin view.
            The other bins are set automatically to have uniform sized bins using an algorithm from Jacquet and Tabot.
            Default = 9 (quite small bins)

        pointMarkerSize : int, optional
279
280
281
282
283
284
            Size of points in point view.
            Default = 8 (quite big points)

        saveFigPath : string, optional
            Path to save figure to -- stops the graphical plotting.
            Default = None
285

Edward Andò's avatar
Edward Andò committed
286
287
    Returns
    -------
288
289
        None -- A matplotlib graph is created and show()n

Edward Andò's avatar
Edward Andò committed
290
291
    Note
    ----
292
293
294
295
296
297
298
299
300
301
302
        Authors: Edward Andò, Hugues Talbot, Clara Jacquet and Max Wiebicke
    """
    import matplotlib.pyplot
    # ========================================================================
    # ==== Reading in data, and formatting to x,y,z sphere                 ===
    # ========================================================================
    numberOfPoints = orientations_zyx.shape[0]

    # ========================================================================
    # ==== Check that all the vectors are unit vectors                     ===
    # ========================================================================
303
    if VERBOSE: print( "\t-> Normalising all vectors in x-y-z representation..." ),
304
305
306
307
308

    # from http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors
    norms = numpy.apply_along_axis( numpy.linalg.norm, 1, orientations_zyx )
    orientations_zyx = orientations_zyx / norms.reshape( -1, 1 )

309
    if VERBOSE: print( "done." )
310
311
312
313

    # ========================================================================
    # ==== At this point we should have clean x,y,z data in memory         ===
    # ========================================================================
314
    if VERBOSE: print( "\t-> We have %i orientations in memory."%( numberOfPoints ) )
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335

    # Since this is the final number of vectors, at this point we can set up the 
    #   matrices for the projection.
    projection_xy       = numpy.zeros( (numberOfPoints, 2) )

    # TODO: Check if there are any values less than zero or more that 2*pi
    projection_theta_r  = numpy.zeros( (numberOfPoints, 2) )

    # ========================================================================
    # ==== Projecting from x,y,z sphere to the desired projection          ===
    # ========================================================================
    # TODO: Vectorise this...
    for vectorN in range( numberOfPoints ):
        # unpack 3D x,y,z
        z,y,x = orientations_zyx[ vectorN ]
        #print "\t\txyz = ", x, y, z

        # fold over the negative half of the sphere
        #     flip every component of the vector over
        if z < 0: z = -z; y = -y; x = -x

336
        projection_xy[ vectorN ], projection_theta_r[ vectorN ] = _projectOrientations( projection, "xyz", [x,y,z] )
337
338
339
340
341
342
343

    # get radiusMax based on projection
    #                                    This is only limited to sqrt(2) because we're flipping over the negative side of the sphere
    if projection == "lambert":         radiusMax = numpy.sqrt(2)
    elif projection == "stereo":        radiusMax = 1.0
    elif projection == "equidistant":   radiusMax = 1.0

344
    if VERBOSE: print( "\t-> Biggest projected radius (r,t) = {}".format( numpy.abs( projection_theta_r[:,1] ).max() ) )
345
346
347
348
349
350
351
352

    #print "projection_xy\n", projection_xy
    #print "\n\nprojection_theta_r\n", projection_theta_r


    if plot == "points" or plot == "both":
        fig = matplotlib.pyplot.figure()
        fig.suptitle( title )
353
354
355
356
        if plot == "both":
          ax  = fig.add_subplot( 121, polar=True )
        else:
          ax  = fig.add_subplot( 111, polar=True)
357

358
        ax.set_title( subtitle['points']+"\n" )
359

360
361
362
        # set the line along which the numbers are plotted to 0°
        #ax.set_rlabel_position(0)
        matplotlib.pyplot.axis( ( 0, math.pi*2, 0, radiusMax ) )
363

364
365
366
367
368
369
        # set radius grids to 15, 30, etc, which means 6 numbers (r=0 not included)
        radiusGridAngles = numpy.arange( 15, 91, 15 )
        radiusGridValues = []
        for angle in radiusGridAngles:
            #                        - project the 15, 30, 45 as spherical coords, and select the r part of theta r-
            #               - append to list of radii -
370
            radiusGridValues.append( _projectOrientations( projection, "spherical", [ 1, angle*math.pi/180.0, 0 ] )[1][1] )
371
372
373
374
        #                                       --- list comprehension to print 15°, 30°, 45° ----------
        ax.set_rgrids( radiusGridValues, labels=[ "%02i$^\circ$"%(x) for x in numpy.arange(  15,91,15) ], angle=None, fmt=None )
        ax.plot( projection_theta_r[:,0], projection_theta_r[:,1] , '.', markersize=pointMarkerSize )

375
376
        if plot == "points":
          matplotlib.pyplot.show()
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398


    if plot == "bins" or plot == "both":
        # ========================================================================
        # ==== Binning the data -- this could be optional...                   ===
        # ========================================================================
        # This code inspired from Hugues Talbot and Clara Jaquet's developments.
        # As published in:
        #   Identifying and following particle-to-particle contacts in real granular media: an experimental challenge
        #   Gioacchino Viggiani, Edward Andò, Clara Jaquet and Hugues Talbot
        #   Keynote Lecture
        #   Particles and Grains 2013 Sydney
        #
        # ...The number of radial bins (numberOfRings)
        # defines the radial binning, and for each radial bin starting from the centre, 
        # the number of angular bins is  4(2n + 1)
        # 
        import matplotlib.patches
        #from matplotlib.colors import Normalize
        import matplotlib.colorbar
        import matplotlib.collections

399
400
401
402
403
404
        if plot == "both":
            ax  = fig.add_subplot( 122, polar=True )
        if plot == "bins":
            fig = matplotlib.pyplot.figure()
            ax  = fig.add_subplot( 111, polar=True)

405
        if VERBOSE: print( "\t-> Starting Data binning..." )
406

407
        # This must be an integer -- could well be a parameter if this becomes a function.
408
        if VERBOSE: print( "\t-> Number of Rings (radial bins) = ", numberOfRings )
409
410


411
412
413
        # As per the publication, the maximum number of bins for each ring, coming from the inside out is 4(2n + 1):
        numberOfAngularBinsPerRing = numpy.arange( 1, numberOfRings+1, 1 )
        numberOfAngularBinsPerRing = 4 * ( 2 * numberOfAngularBinsPerRing - 1 )
414

415
        if VERBOSE: print( "\t-> Number of angular bins per ring = ", numberOfAngularBinsPerRing )
416

417
418
        # defining an array with dimensions numberOfRings x numberOfAngularBinsPerRing
        binCounts = numpy.zeros( ( numberOfRings, numberOfAngularBinsPerRing[-1] ) )
419

420
421
422
423
424
425
        # ========================================================================
        # ==== Start counting the vectors into bins                            ===
        # ========================================================================
        for vectorN in range( numberOfPoints ):
            # unpack projected angle and radius for this point
            angle, radius = projection_theta_r[ vectorN, : ]
426

427
428
429
            # Flip over negative angles
            if angle < 0:             angle += 2*math.pi
            if angle > 2 * math.pi:   angle -= 2*math.pi
430

431
432
            # Calculate right ring number
            ringNumber = int(numpy.floor( radius / ( radiusMax / float(numberOfRings) ) ) )
433

434
435
            # Check for overflow
            if ringNumber > numberOfRings - 1:
436
                if VERBOSE: print( "\t-> Point with projected radius = %f is a problem (radiusMax = %f), putting in furthest  bin"%( radius, radiusMax ) )
437
                ringNumber = numberOfRings - 1
438

439
440
            # Calculate the angular bin
            angularBin = int( numpy.floor( ( angle ) / ( 2 * math.pi / float( numberOfAngularBinsPerRing[ ringNumber ] ) ) ) ) + 1
441

442
443
444
445
446
447
            #print "numberOfAngularBinsPerRing", numberOfAngularBinsPerRing[ringNumber] - 1
            # Check for overflow
            #  in case it doesn't belong in the last angularBin, it has to be put in the first one!
            if angularBin > numberOfAngularBinsPerRing[ringNumber] - 1:
                if VERBOSE: print( "\t-> Point with projected angle = %f does not belong to the last bin, putting in first bin"%( angle ) )
                angularBin = 0
448

449
450
            # now that we know what ring, and angular bin you're in add one count!
            binCounts[ ringNumber, angularBin ] += 1
451

452
453
454
        # ========================================================================
        # === Plotting binned data                                             ===
        # ========================================================================
455

456
457
        plottingRadii = numpy.linspace( radiusMax/float(numberOfRings), radiusMax, numberOfRings )
        #print "Plotting radii:", plottingRadii
458
459
460

        #ax  = fig.add_subplot(122, polar=True)
        #matplotlib.pyplot.axis(  )
461
462
        #ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], polar=True)
        bars = []
463

464
465
        # add two fake, small circles at the beginning so that they are overwritten
        #   they will be coloured with the min and max colour
466
467
468
469
470
        #              theta   radius    width
        bars.append( [   0,   radiusMax,    2*math.pi ] )
        bars.append( [   0,   radiusMax,    2*math.pi ] )
        #bars.append( ax.bar( 0,   radiusMax,    2*math.pi, bottom=0.0 ) )
        #bars.append( ax.bar( 0,   radiusMax,    2*math.pi, bottom=0.0 ) )
471

472
473
        # --- flatifiying binned data for colouring wedges                    ===
        flatBinCounts = numpy.zeros( numpy.sum( numberOfAngularBinsPerRing ) + 2 )
474

475
476
477
        # Bin number as we go through the bins to add the counts in order to the flatBinCounts
        # This is two in order to skip the first to fake bins which set the colour bar.
        binNumber = 2
478

479
        # --- Plotting binned data, from the outside, inwards.                 ===
480
481
482
483
484
        if binNormalisation:
            avg_binCount = float(numberOfPoints)/numpy.sum( numberOfAngularBinsPerRing )
            #print "\t-> Number of points = ", numberOfPoints
            #print "\t-> Number of bins   = ", numpy.sum( numberOfAngularBinsPerRing )
            if VERBOSE: print( "\t-> Average binCount = ", avg_binCount )
485

486
487
        for ringNumber in range( numberOfRings )[::-1]:
            deltaTheta    = 360 / float( numberOfAngularBinsPerRing[ringNumber] )
488
            deltaThetaRad = 2 * math.pi / float( numberOfAngularBinsPerRing[ringNumber] )
489

490
491
            # --- Angular bins                                                 ---
            for angularBin in range( numberOfAngularBinsPerRing[ringNumber] ):
492
493
494
495
                # ...or add bars
                #                           theta                             radius                  width
                bars.append( [ angularBin*deltaThetaRad - deltaThetaRad/2.0, plottingRadii[ ringNumber ], deltaThetaRad ] )
                #bars.append( ax.bar( angularBin*deltaThetaRad - deltaThetaRad/2.0, plottingRadii[ ringNumber ], deltaThetaRad, bottom=0.0 ) )
496

497
                # Add the number of vectors counted for this bin
498
499
500
501
                if binNormalisation:
                    flatBinCounts[ binNumber ] = binCounts[ ringNumber, angularBin ]/avg_binCount
                else:
                    flatBinCounts[ binNumber ] = binCounts[ ringNumber, angularBin ]
502

503
504
                # Add one to bin number
                binNumber += 1
505

506
        del binNumber
507

508
        # figure out auto values if they're requested.
Edward Andò's avatar
Edward Andò committed
509
510
        if binValueMin is None: binValueMin = flatBinCounts[2::].min()
        if binValueMax is None: binValueMax = flatBinCounts[2::].max()
511

512
        # Add two flat values for the initial wedges.
Edward Andò's avatar
Edward Andò committed
513
514
        flatBinCounts[0] = binValueMin
        flatBinCounts[1] = binValueMax
515

516
517
        ##                           theta                   radius                          width
        barsPlot = ax.bar( numpy.array( bars )[:,0], numpy.array( bars )[:,1], width=numpy.array( bars )[:,2], bottom=0.0)
518

519
        for binCount,bar in zip(  flatBinCounts, barsPlot ):
Edward Andò's avatar
Edward Andò committed
520
            bar.set_facecolor( cmap(  ( binCount - binValueMin) / float( binValueMax - binValueMin ) ) )
521

522
523
        #matplotlib.pyplot.axis( [ 0, radiusMax, 0, radiusMax ] )
        matplotlib.pyplot.axis( [ 0, numpy.deg2rad(360), 0, radiusMax ] )
524

525
526
527
        #colorbar = matplotlib.pyplot.colorbar( barsPlot, norm=matplotlib.colors.Normalize( vmin=minBinValue, vmax=maxBinValue ) )
        # Set the colormap and norm to correspond to the data for which
        # the colorbar will be used.
528

Edward Andò's avatar
Edward Andò committed
529
        norm = matplotlib.colors.Normalize( vmin=binValueMin, vmax=binValueMax )
530
531
532
533
534
535
536
537

        # ColorbarBase derives from ScalarMappable and puts a colorbar
        # in a specified axes, so it has everything needed for a
        # standalone colorbar.  There are many more kwargs, but the
        # following gives a basic continuous colorbar with ticks
        # and labels.
        ax3 = fig.add_axes([0.9, 0.1, 0.03, 0.8])
        cb1 = matplotlib.colorbar.ColorbarBase( ax3, cmap=cmap, norm=norm )
538

539
540
        # set the line along which the numbers are plotted to 0°
        #ax.set_rlabel_position(0)
541

542
543
544
545
546
547
548
549
550
        # set radius grids to 15, 30, etc, which means 6 numbers (r=0 not included)
        radiusGridAngles = numpy.arange( 15, 91, 15 )
        radiusGridValues = []
        for angle in radiusGridAngles:
            #                        - project the 15, 30, 45 as spherical coords, and select the r part of theta r-
            #               - append to list of radii -
            radiusGridValues.append( _projectOrientations( projection, "spherical", [ 1, angle*math.pi/180.0, 0 ] )[1][1] )
        #                                       --- list comprehension to print 15°, 30°, 45° ----------
        ax.set_rgrids( radiusGridValues, labels=[ "%02i$^\circ$"%(x) for x in numpy.arange(  15,91,15) ], angle=None, fmt=None )
551

552
553
        fig.subplots_adjust(left=0.05,right=0.85)
        #cb1.set_label('Some Units')
554

555
556
557
558
        if saveFigPath is not None:
          matplotlib.pyplot.savefig( saveFigPath )
        else:
          matplotlib.pyplot.show()
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
          
def distributionDensity( F, step = 50, lim = None, color = None, title = None, saveFigPath = None):
    """
    Creates the surface plot of the distribution density of the deviatoric fabric tensor F

    Parameters
    ----------
        F : 3x3 array of floats
            deviatoric fabric tensor. Usually obtained from spam.label.fabricTensor
            
        step : int, optional
            Number of points for the surface plot
            Default = 50
                
        lim : float, optional
            Limit for the axes of the plot
            Default = None
            
        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.Reds
            
        title : str, optional
            Title for the graph
            Default = None
            
        saveFigPath : string, optional
            Path to save figure to.
            Default = None
    
    Returns
    -------
        None -- A matplotlib graph is created and shown
    
    Note
    ----
        see [Kanatani, 1984] for more information on the distribution density function for the deviatoric fabric tensor
    
    """
    #Create array of angles
    theta, phi = numpy.linspace(0, 2 * numpy.pi, step), numpy.linspace(0, numpy.pi, step)
    #Create meshgrid
    THETA, PHI = numpy.meshgrid(theta, phi)
    #Create radius array
    R = numpy.zeros(THETA.shape)
    #Copmute the radius for each angle
    for r in range(0,step,1):    
        for s in range(0,step,1):
            vect = numpy.array((numpy.cos(phi[r]),
                                numpy.sin(phi[r])*numpy.sin(theta[s]),
                                numpy.cos(theta[s])*numpy.sin(phi[r])))
            R[r,s] = (1/(4*numpy.pi))*(1+numpy.dot(numpy.dot(F,vect),vect))
    #Change to cartesian coordinates
    X = R * numpy.sin(PHI) * numpy.cos(THETA)
    Y = R * numpy.sin(PHI) * numpy.sin(THETA)
    Z = R * numpy.cos(PHI)
    #Create figure
618
619
    import matplotlib
    matplotlib.rcParams.update({'font.size': 10})
620
621
622
623
624
625
626
627
628
629
630
631
632
    fig = matplotlib.pyplot.figure()
    ax  = fig.add_subplot( 111, projection = '3d' )
    #Set limits
    if lim == None:
        lim = round(numpy.max(R), 2)
    ax.set_xlim3d(-lim, lim)
    ax.set_ylim3d(-lim, lim)
    ax.set_zlim3d(-lim, lim)
    #Set ticks
    ax.set_xticks((-lim, 0, lim))
    ax.set_yticks((-lim, 0, lim))
    ax.set_zticks((-lim, 0, lim))
    ax.set_aspect('equal', 'box')
633
634
635
636
    # set axis titles
    ax.set_xlabel('X axis')
    ax.set_ylabel('Y axis')
    ax.set_zlabel('Z axis')
637
    #Title
638
639
    if title is not None:
        ax.set_title( str(title) + "\n" )
640
641
642
643
644
645
646
647
    #Colormap
    if color == None:
        cmap = matplotlib.pyplot.get_cmap(matplotlib.pyplot.cm.Reds)
    else:
        cmap = matplotlib.pyplot.get_cmap(color)
    norm = mcolors.Normalize(vmin=0, vmax=Z.max())
    #Plot
    ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1, 
648
649
650
                           #facecolors = cmap(norm(numpy.abs(Z))),
                           #coloring by max extension
                           facecolors = cmap((R-numpy.amin(R))/numpy.amax(R-numpy.amin(R))),
651
652
653
654
655
656
657
658
                           linewidth = 0, antialiased = True, 
                           alpha = 1)
    
    matplotlib.pyplot.tight_layout()
    if saveFigPath is not None:
        matplotlib.pyplot.savefig( saveFigPath )
    else:
        matplotlib.pyplot.show()
659
660
661
662


if __name__ == "__main__":
    #reference       = numpy.loadtxt( orientationsFile, usecols=(2,3,4))
663
    SPHERE_TEST = False
664
665

    if SPHERE_TEST:
666
        print ("\t-> Working with \"equally\" spaced points of a sphere")
667
668
669
670
671
672
        # This overrides the length of the file, and can simply be set to the number of 
        #   points desired
        numberOfPoints = 4000
        orientations_zyx = SaffAndKuijlaarsSpiral( numberOfPoints )

    else:
673
        # Max's benchmarking points for the binning
674
        orientations_zyx = numpy.array( [ 
675
676
677
678
679
680
681
682
683
684
685
                                  [ 1, 0, 0 ],    \
                                  [ -1, 0, 0 ],   \
                                  [ 0.2, 1, 1 ],  \
                                  [ 0.2,-1,-1],   \
                                  [ 0,1,1],       \
                                  [ 1,1,1],       \
                                  [ 0,1,-0.0],       \
                                  [ 0,1,0.01],      \
                                  [ 1,-1,1 ],     \
                                  [ 1.8,0.1,1 ]  ] 
                                )
686
    plotOrientations( orientations_zyx,
687
688
                    projection="lambert",
                    plot="both",
Edward Andò's avatar
Edward Andò committed
689
690
                    binValueMin=None,
                    binValueMax=None,
691
                    numberOfRings=7,
692
693
694
695
696
                    pointMarkerSize=8,
                    title = "Demonstration of orientation plotting",
                    subtitle = { "points":"Points",
                                 "bins":"Bins"}
                    )
Edward Andò's avatar
Edward Andò committed
697

698
699
700
701
702
703
704
705
706
707
    #plotProjection( orientations_zyx,
                    ##projection="lambert",
                    ##plot="both",
                    ##binValues={ "min": None, "max": None},
                    ##numberOfRings=9,
                    ##pointMarkerSize=8,
                    ##title = "Demonstration of orientation plotting",
                    ##subtitle = { "points":"Points",
                                 ##"bins":"Bins"}
                    #)
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
                    
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,
867
                           verbose=True,
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
                           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.
911
        If it is not given, the plot will be shown but not saved.
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
        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:
969
        print('\nspam.helpers.orientationPlotter.plotSphericalHistogram: The input vectors are not 3D')
970
971
972
973
974
975
976
977
978
        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:
979
        print('\nspam.helpers.orientationPlotter.plotSphericalHistogram: Creating the icosphere')
980
981
982
    icoVerts, icoFaces, icoVectors = createIcosphere(subDiv)
    # Bin the data
    if verbose:
983
        print('\nspam.helpers.orientationPlotter.plotSphericalHistogram: Binning the data')
984
985
986
    counts = binIcosphere(orientations, icoVectors, verbose = verbose)
    # Now we are ready to plot
    if verbose:
987
        print('\nspam.helpers.orientationPlotter.plotSphericalHistogram: Plotting')
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
        
    # 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()