Neighbor based analysing tool

[Tagebh]

Description of the change

I am currently working on hybrid indexing as part of my master at NTNU, and I also wrote a project thesis on hybrid indexing in the autumn of last year, where I used the hybrid indexing tutorial made på Håkon as a starting point. As part of this work I have made a Kernel Average Misorientation (KAM) tool based on the code from @argerlt in issue #531.

Progress of the PR

• Docstrings for all functions
• Unit tests with pytest for all lines
• Clean code style by running black via pre-commit

I have structured the code so it has the potential to be used for many things involving neighbors (not only KAM).

The Neighbors() function performs the general preperations to do neighbor-based calculations on a xmap or any 2d array really, and can be used on it's own to do whatever you want with neighbors.

The neighbor_misorientation() and KAM_calc() then uses Neighbors() to make a KAM map.

I have also included a slight variation of KAM which I am calling Number of Same Neighbors, which also uses the results from Neighbors() and neighbor_misorientation() but does a slightly different calculation. A user defines what degree of misorientation that counts as a "different" orientation, and then it counts how many neighbors has the same orientation as a central point, where each point in the dataset acts as a central point once.

import numpy as np

from orix.crystal_map import CrystalMap
from orix.quaternion import Orientation

from skimage.morphology._util import _raveled_offsets_and_distances
from skimage.util._map_array import map_array

def Neighbors(data,foot,mask=None):  
    """Performs the general preperations to do neighbor-based calculations on CrystalMaps in a vectorized form
        
    Args:
        data: 2d array. Can be a CrystalMap object
        mask: 2D array with same shape as CrystalMap indicating what points are valid or not eg. Non-indexed points
        foot: The kernel that defines which neighbors to look at. Can be both binary and not
        
    Returns: (maybe should be returned as an object instead?)
        indices: Array where the index for each valid point is repeated as many times as it has valid neighbors
        neighbor_indices: Array with the corosponding valid neighbors
        foot_values: contains how each neighbor in neighbor_indices is weighted defined by a given footprine/kernel
        num_neigbors_valid: the number of neighbors each valid point has
        nodes_valid: The index for valid points with at least one valid neighbor
        nodes_no_valid_neighbors: The index for valid points with no valid neighbors
        nodes_not_valid: The index for non-valid points 
        
    Raises:
        Have to add ValueErrors
    """
    if mask is None:
        if isinstance(data, CrystalMap):
            mask = data.is_indexed.reshape(data.shape)
        else:
            mask = np.full(data.shape, True) #All pointd are valid
    nodes = np.flatnonzero(mask)
    
    pad = np.max(foot.shape)//2
    padded_mask = np.pad(mask, pad, mode='constant', constant_values=False) 
    padded_nodes = np.flatnonzero(padded_mask)

    #The offset are given in 1D for the neighbors defined in 2D
    neighbor_offsets, dist = _raveled_offsets_and_distances(padded_mask.shape, footprint=foot)

    padded_neighbors = padded_nodes[:, np.newaxis] + neighbor_offsets
    neighbors = map_array(padded_neighbors, padded_nodes, nodes) #finds the indeces (in a non-padded format) for the neighbors belonging to valid points
    neighbors_mask = padded_mask.reshape(-1)[padded_neighbors] #sets which of those neigbors are valid
    
    num_neighbors = np.sum(neighbors_mask, axis=1)#the number of neighbors each valid point has
    indices = np.repeat(nodes, num_neighbors) #Array where the index for each valid point is repeated as many times as it has valid neighbors
    neighbor_indices = neighbors[neighbors_mask] #Array with the corosponding valid neighbors

    #Support for non-binary footprints 
    foot_offsets, dist = _raveled_offsets_and_distances(foot.shape, footprint=foot)
    foot_nodes = foot_offsets+foot.size//2
    foot_values = foot.reshape(-1)[np.repeat([foot_nodes],mask.sum(),axis=0)[neighbors_mask]]#contains how each neighbor in neighbor_indices is weighted defined by a given footprine

    valid_neighbor_mask = num_neighbors !=0
    nodes_valid = nodes[valid_neighbor_mask] #The index for valid points with at least one valid neighbor
    nodes_no_valid_neighbors = nodes[~valid_neighbor_mask] #The index for valid points with no valid neighbors
    nodes_invalid = np.flatnonzero(~mask) #The index for non-valid points (same as the mask input, not really neccesary)
    
    num_neighbors_valid = num_neighbors[valid_neighbor_mask]
    
    return(indices, neighbor_indices, foot_values, num_neighbors_valid, nodes_valid, nodes_no_valid_neighbors, nodes_invalid)


def neighbor_misorientation(xmap, indices, neighbor_indices, degree=False, crystall_symmetry=None):
    """Calculates the misorientation angles between a central point and its neighbors
    Args:
        xmap: CrystalMap object
        indices: Array where the index for each valid point is repeated as many times as it has valid neighbors
        neighbor_indices: Array with the corosponding valid neighbors
        degree (bool): True-return results in degree. False-return results in radiens
        crystall_symmetry: used to differentiate between rotation and orientation
    Returns:
        d: misorientation angles given in radiens
        
    Raises:
        Have to add ValueErrors
    """
    if crystall_symmetry is None:
        crystall_symmetry = xmap.phases[0].point_group #This can pehaps lead to problems, maybe make a check for multiple phases
    
    O_central_points = Orientation(xmap.rotations[indices], crystall_symmetry)
    O_neighbors = Orientation(xmap.rotations[neighbor_indices], crystall_symmetry)
    mis_ori = O_central_points.angle_with(O_neighbors)

    if degree:
        mis_ori = mis_ori*180/np.pi
    
    return mis_ori

def KAM_calc(xmap, mis_ori, non_index_value, no_neighbors_value, num_neighbors_valid, nodes_valid, nodes_no_valid_neighbors, nodes_invalid, foot_values=None):
    
    """Makes the KAM map
    Args:
        xmap: CrystalMap object (or any 2d array actually)
        mis_ori: Misorientation angles given in radiens or degrees (from neighbor_misorientation())
        non_index_value (int or float): Kam value given to a non-indexed points in the xmap
        no_neigbors_value (int or float): Kam value given to an indexed point without a single indexed neighbor

        These are from Neighbors()
        num_neigbors_valid: the number of neighbors each valid point has
        nodes_valid: The index for valid points with at least one valid neighbor
        nodes_no_valid_neighbors: The index for valid points with no valid neighbors
        nodes_not_valid: The index for non-valid points 
        foot_values: The weight of different neighbors if non-binary foot is used
        
    Returns:
        kam_map_im: 2D array with same shape as xmap with all the KAM values
        
    Raises:
        Have to add ValueErrors
    """
    if foot_values is None:
        cumulative_miso = np.add.reduceat(mis_ori, np.r_[0, np.cumsum(num_neighbors_valid)[:-1]]) 
    else:
        cumulative_miso = np.add.reduceat(mis_ori*foot_values, np.r_[0, np.cumsum(num_neighbors_valid)[:-1]])
    
    kam_map = np.full(xmap.size, np.nan, dtype=np.float32)
    kam_map[nodes_valid] = cumulative_miso/num_neighbors_valid
    kam_map[nodes_no_valid_neighbors] = no_neighbors_value
    kam_map[nodes_invalid] = non_index_value
    kam_map_im = kam_map.reshape(xmap.shape)
    
    return kam_map_im

def NSN_calc(xmap, mis_ori, non_index_value, no_neighbors_value, lim, num_neighbors_valid, nodes_valid, nodes_no_valid_neighbors, nodes_invalid, foot_values=None):
    """Makes a Number of same neighbors (NSN) map. Each point in the xmap gets assigned the value equal to the number similar 
       oriented neighbors. Where similar is defined by a user set limit 
    Args:
        xmap: CrystalMap object (or any 2d array actually)
        mis_ori: Misorientation angles given in radiens or degrees (from neighbor_misorientation())
        non_index_value (int or float): NSN value given to a non-indexed points in the xmap
        no_neigbors_value (int or float): NSN value given to an indexed point without a single indexed neighbor
        lim (int or float): The limit for when two neighboring point are defined to have different/same orientations

        These are from Neighbors()
        num_neigbors_valid: the number of neighbors each valid point has
        nodes_valid: The index for valid points with at least one valid neighbor
        nodes_no_valid_neighbors: The index for valid points with no valid neighbors
        nodes_not_valid: The index for non-valid points
        foot_values: The weight of different neighbors if non-binary foot is used
        
    Returns:
        NDN_map_im: 2D array with same shape as xmap with all the NDN values
        
    Raises:
        Have to add ValueErrors
    """

    same_neighbors = mis_ori < lim
    if foot_values is None:
        NSN = np.add.reduceat(same_neighbors, np.r_[0, np.cumsum(num_neighbors_valid)[:-1]]) 
    else:
        NSN = np.add.reduceat(same_neighbors*foot_values, np.r_[0, np.cumsum(num_neighbors_valid)[:-1]])
    
    NSN_map = np.full(xmap.size, np.nan, dtype=np.float32)
    NSN_map[nodes_valid] = NSN 
    NSN_map[nodes_no_valid_neighbors] = no_neighbors_value
    NSN_map[nodes_invalid] = non_index_value
        
    NSN_map_im = NSN_map.reshape(xmap.shape)
    
    return NSN_map_im

Here is the functions as .py file
_neighbors.py

The reulting KAM and NDN (Number of different neigbors (NDN) = total neighbors - NSN) map using a 3x3 binary kernel (8 neighbors with equal weighting) could look something like this (using one of the nickel datasets already in kikuchipy):

I have also made a small jupyter notebook that I used to test the functions
test_neighbors.ipynb

Here is an example from that notebook where the amount of equal colored neighbors is counted in a small dummy dataset

import _neighbors as Neighbors #import my library

#make a dummy dataset
test_mask = np.array([[True,True,True,False,True,True,True,],
                      [True,False,True,True,True,True,True],
                      [True,True,True,False,False,False,True,],
                      [True,True,True,False,True,False,False],
                      [True,True,True,False,False,False,True,]])

test_data = np.array([[1,2,1,0,2,1,1],
                      [1,0,2,1,2,2,1],
                      [2,2,2,0,0,0,2],
                      [2,2,2,0,1,0,0],
                      [2,2,2,0,0,0,2]])

#define the footprint and run the preperation function neighbors()
x=3
foot = np.ones([x,x])
indices, neighbor_indices, foot_values, num_neighbors_valid, nodes_valid, nodes_no_valid_neighbors, nodes_not_valid=Neighbors.Neighbors(data = test_data,
                                                                                                                                        foot = foot,
                                                                                                                                        mask = test_mask)
#calculate how many neighbors how the same value/color as the central point
central_points = test_data.reshape(-1)[indices]
neighbors = test_data.reshape(-1)[neighbor_indices]

l = []
for i in range(len(ori)):
    if central_points[i] == neighbors[i]:
        l.append(1)
    else:
        l.append(0)
same_N = np.array(l)

np.r_[0, np.cumsum(num_neighbors_valid)[:-1]]
sum_same_N = np.add.reduceat(same_N*foot_values, np.r_[0, np.cumsum(num_neighbors_valid)[:-1]])

im = np.full(test_data.size, np.nan)
im[nodes_valid] = sum_same_N
im[nodes_no_valid_neighbors] = -1
im[nodes_not_valid] = -2

im = im.reshape(test_data.shape)

#display the footprint used, the dataset and the number of same neighbors

print('Footprint:')
print()
print(foot)

fig, ax = plt.subplots()
rows = test_data.shape[0]
cols = test_data.shape[1]
ax.imshow(test_data, extent=[0,cols, rows, 0], origin='lower')
ax.set_xlim(0, cols)
ax.set_ylim(0, rows)
ax.set_xticks(np.arange(0, cols + 1))
ax.set_yticks(np.arange(0, rows + 1))
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(which='major', color='k', linestyle='-', linewidth=2)

for i in range(rows):
    for j in range(cols):
        ax.text(0.5 + j, 0.5 + i, im[rows - 1 - i, j], ha='center', va='center', color='red', fontsize=12)
        
ax.set_title('Number of same neighbors')

plt.show()

print('purple: data value is 0 (invalid)')
print('teal: data value is 1')
print('yellow: data value is 2')
print()
print('-2: invalid point (value set by user)')
print('-1: valid point without a valid neighbor (value set by user) (can be both green and yellow)')

The output with a 3x3 binary footprint

The output when the footprint is 3x3 and not binary (0.5 in the corners)

I do not have any experience with git hub, so I do not really know how to impliment it, but I hope this seems usefull and that someone wants to pick it up.

Kind regards
Tage

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.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[hakonanes]

Thank you @Tagebh, the functionality in your comment will be very useful to have in orix!

When adding your own commits, can you use the develop branch as the target branch here in the orix repo?

[Tagebh]

Hey @hakonanes

I changed the PR so it commits to the developer branch, and I added the code to orix/crystal_map/_neighbors.py. I also made the PR ready for review

[hakonanes]

Great, thank you!

Just to confirm, is it fine with you if I push to your branch? I plan to push stuff like tests, documentation, and perhaps some restructuring of your file. Alternatively, I can make a PR to the develop branch in your fork, so you can go through my suggested changes before being visible here.

[argerlt]

@hakonanes and @Tagebh , do either of you care if I try tackling this actually? Same comment as @hakonanes, making pushes directly to the original branch.

@hakonanes, if you've already started, ignore this comment, there are plenty of other things for me to do in ORIX.

Also, general question: how do we feel about this pattern of keeping all the neighborhood functions in a _neighbors.py file? I think i like it, but just trying to think of potential pitfalls.

[Tagebh]

@argerlt, as for me you are more than welcome to tackle this problem if you want to. I don't really know if @hakonanes has started or not, so it's probably best to wait for him to answer as well.

I feel like that it's natural to add new neighborhood functions into _neighbors.py

[argerlt]

@Tagebh, sorry this took so long, thank you for your patience.

Some starting notes:

1. To help improve readability and consistency, ORIX tries to stick as close as possible to the numpydoc style guide. In terms of this PR, this means switching docstrings from this:

    """Performs the general preperations to do neighbor-based calculations on CrystalMaps in a vectorized form
        
    Args:
        data: 2d array. Can be a CrystalMap object
        mask: 2D array with same shape as CrystalMap indicating what points are valid or not eg. Non-indexed points
        foot: The kernel that defines which neighbors to look at. Can be both binary and not
        
    Returns: (maybe should be returned as an object instead?)
        indices: Array where the index for each valid point is repeated as many times as it has valid neighbors
        neighbor_indices: Array with the corosponding valid neighbors
        foot_values: contains how each neighbor in neighbor_indices is weighted defined by a given footprine/kernel
        num_neigbors_valid: the number of neighbors each valid point has
        nodes_valid: The index for valid points with at least one valid neighbor
        nodes_no_valid_neighbors: The index for valid points with no valid neighbors
        nodes_not_valid: The index for non-valid points 
        
    Raises:
        Have to add ValueErrors
    """

to something like this:

def Neighbors(data: np.ndarray,
                          foot:np.ndarray,
                          mask:np.ndarray |None=None):  
    """Calculates the arrays necessary for neighbor-based calculations on CrystalMaps in a vectorized form
        
    Parameters
    ---------
        data
            either a 2D array or a crystal map.
        foot
            The kernel defining the local neighborhood.
        mask
            Boolean array of the same shape as `data`, indicating which points are included in the
            neigborhood lookup.
        
    Returns
    -------
        indices
            Array where the index for each valid point is repeated as many times as it has valid neighbors.
        neighbor_indices
            Array with the corosponding valid neighbors
    """

2. Again for sanity purposes, we try to stick with the PEP 8 naming guidelines, specifically that funcitons should be lowercase, and functions intended for private use should have underscores before them. in this case, Neigbors should be _neighbors, KAM_calc should be kam_calc, etc.

3. I know my bad code is the root source of this error, but we try to avoid adding new dependencies when we can and especially dependencies on private functions like skimage.morphology._util._raveled_offsets_and_distances and skimage.util._map_array.map_array (note the underscores indicating private functions and modules)

4. ORIX as a whole also has a strictly enforced style rules to help ensure people a year from now can still read your code. Check the ORIX style guide as well as the guide for setting up a development environment as a starting point, but ill also add you can cheese the system a bit by autofixing your PR by adding the following comment:

pre-commit.ci autofix

[argerlt]

^ If you click on pre-commit.ci's commit above, you can see what rules it enforced.

[argerlt]

@Tagebh, I added some comments and suggested changes as a review, feel free to ask for clarification or suggest alternatives.

Also, once these changes are made, you/we will need to write unit tests and examples.

[hakonanes]

As discussed off GitHub, @argerlt, I think the public API from this PR should be the kernel average misorientation calculation only. How we handle the other functions internally is then of secondary nature, thankfully.

We can look to make more functions public afterwards, if needed.

[hakonanes]

And it's my understanding from @Tagebh that you can make changes to his branch directly.

[hakonanes]

@argerlt, if you don't have time, I can wrap this up. Up to you.

Went from reviewer to writer, so "approving" this feels a bit too much like skipping review all together

[argerlt]

@hakonanes , sorry this took me a bit to get back around to. I swapped the functions around amd resolved my own review comments.

I'll finish up some changes and documentation, then ping you when this is ready for review.

[hakonanes]

Sounds good.