New Feature: 3D CrystalMap Objects

[smason747]

Description of the change

As discussed in #577, this PR is working to re-introduce the capability for ORIX to handle 3D datasets. As an initial step, I wrote a CrystalMap3D class that extends from the standard CrystalMap class and adds a z variable as a new property. In discussing this approach, it was recommended to absorb those changes into CrystalMap itself, rather than introduce a new class. This has not yet been done at the time of submitting the PR, but opens the floor to do so. Along with updating the data structure itself, changes need to be made to the plotting functionality to plot individual slices, either for a default value if not specified, at a chosen depth, on a chosen plane (xy, yz, xz), or with an interactive slider to allow fly-through viewings. Additionally, this could be a helpful time to address any thoughts on the underlying data structure of the CrystalMap class, as well as ensure loading of 3D files can be handled by orix.io.

Progress of the PR

• Docstrings for all functions
• Unit tests with pytest for all lines
• Clean code style by running black via pre-commit
• Move z functionality to 'CrystalMap' class, rather than create new 'CrystalMap3D' class
• Handle 2D plot behaviour for 3D data

Minimal example of the bug fix or new feature

# import orix
from orix.crystal_map import Phase, CrystalMap, PhaseList
from orix.vector import Miller, Vector3d
from orix.quaternion import Orientation, Rotation
from orix.quaternion.symmetry import Oh
from orix import plot
from diffpy.structure import Lattice, Structure, Atom
from orix import io as oxio

from orix.crystal_map.crystal_map_3d import CrystalMap3D


# %%

structures = [
    Structure(
        title='nickel',
        atoms=[Atom('ni', [0] * 3)],
        lattice=Lattice(3.520,3.520,3.520,90,90,90)
    )
]

# %%

with h5py.File('/Users/simon/data/reference/in100/in100.h5', 'r') as hf:

    xmap3d = CrystalMap3D(
        rotations=Rotation(np.asarray(hf['rotations'])),
        phase_id=np.asarray(hf['pid']),
        x=np.asarray(hf['x']),
        y=np.asarray(hf['y']),
        z=np.asarray(hf['z']),
        phase_list=PhaseList(space_groups=[225], structures=structures),
        prop={'ci': np.asarray(hf['ci']),'iq': np.asarray(hf['iq']),
              'detector_signal': np.asarray(hf['ds']), 
              'fit': np.asarray(hf['fit'])}
              )

# %%

print(xmap3d)
# Out[114]: 
# Phase      Orientations         Name  Space group  Point group  Proper point group     Color
#    -1    668117 (15.2%)  not_indexed         None         None                None         w
#     0   3738607 (84.8%)       nickel        Fm-3m         m-3m                 432  tab:blue
# Properties: ci, iq, detector_signal, fit
# Scan unit: px

print(xmap3d.shape)
# Out[115]: (116, 201, 189)

xmap3d[100,:,:].plot('ci')
xmap3d[:,100,:].plot('ci')
xmap3d[:,:,100].plot('ci')

For reviewers

• The PR title is short, concise, and will make sense 1 year later.
• New functions are imported in corresponding __init__.py.
• New features, API changes, and deprecations are mentioned in the unreleased
  section in CHANGELOG.rst.
• Contributor(s) are listed correctly in __credits__ in orix/__init__.py and in
  .zenodo.json.

[argerlt]

@hakonanes:

I was looking at this PR and trying to fix the errors, and I think crystal_map.create_coordinate_arrays is transposed from the numpy-like behavior. Can you confirm?

What it does (as per the docstring):

>>> create_coordinate_arrays((2, 3))
({'x': array([0, 1, 2, 0, 1, 2]), 'y': array([0, 0, 0, 1, 1, 1])}, 6)

This is doing the tiling column-major (Fortran style) as opposed to row-major (Numpy style). I would expect the opposite, so that it behaved like np.meshgrid:

[x.flatten() for x in np.meshgrid([0,1],[0,1,2])]
[array([0, 1, 0, 1, 0, 1]), array([0, 0, 1, 1, 2, 2])]

The difference is pretty trivial for orix right now, but it matters for how the 3d ebsd is added. Here is what I think it should be corrected to, when extended to three dimensions:

    >>> from orix.crystal_map import create_coordinate_arrays
    >>> create_coordinate_arrays((3, 2))
    ({'x': array([0, 1, 2, 0, 1, 2]), 'y': array([0, 0, 0, 1, 1, 1])}, 6)
    >>> create_coordinate_arrays((2, 3))
    ({'x': array([0, 1, 0, 1, 0, 1]), 'y': array([0, 0, 1, 1, 2, 2])}, 6)
    >>> create_coordinate_arrays((2, 3, 3), (1, 2, 3))
    ({'x': array([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]),
      'y': array([0, 0, 2, 2, 4, 4, 0, 0, 2, 2, 4, 4, 0, 0, 2, 2, 4, 4]),
      'z': array([0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6])},
     18)

[argerlt]

@smason747, check out this branch when you get time and let me know what you think.