Neighbor based analysing tool by Tagebh · Pull Request #642 · pyxem/orix

Tagebh · 2026-03-23T21:10:49Z

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