# Example M3: Color mapping to visualize local lattice orientation

The Polyhedredral Template Matching (PTM) function of OVITO allows computing the local lattice orientation for each atom in a (poly)crystal. The computed local orientations are stored by the modifier as quaternions, i.e. as rotations within the fundamental zone, in the particle property named `Orientation`. Each per-particle quaternion can be translated into an RGB color to visualize the local lattice orientation. This can be achieved by inserting a custom Python modifier into the pipeline which translates the output of the PTM modifier into RGB values and stores them in the `Color` particle property.

In the graphical OVITO version, simply insert a new Python modifier and copy/paste the following script into the source code window:

    from ovito.data import *
    import math
    import numpy as np

    def quaternions_to_colors(qs):
        """ Takes a list of quaternions (Nx4 array) and returns a list of
            corresponding RGB colors (Nx3 array) """

        if len(qs.shape) != 2:
            raise RuntimeError("qs must be a 2-dimensional array")

        if qs.shape[1] != 4:
            raise RuntimeError("qs must be a n x 4 dimensional array")

        # Project quaternions into Rodrigues space: rs = (qs.X/qs.W, qs.Y/qs.W, qs.Z/qs.W)
        # Note that the qs.W may be zero for particles for which no lattice orientation
        # could be computed by the PTM modifier.
        rs = np.zeros_like(qs[:,:3])
        np.divide(qs[:,0], qs[:,3], out=rs[:,0], where=qs[:,3] != 0)
        np.divide(qs[:,1], qs[:,3], out=rs[:,1], where=qs[:,3] != 0)
        np.divide(qs[:,2], qs[:,3], out=rs[:,2], where=qs[:,3] != 0)

        # Compute vector lengths rr = norm(rs)
        rr = np.linalg.norm(rs, axis=1)
        rr = np.maximum(rr, 1e-9) # hack

        # Normalize Rodrigues vectors.
        rs[:,0] /= rr
        rs[:,1] /= rr
        rs[:,2] /= rr

        theta = 2 * np.arctan(rr)
        rs[:,0] *= theta
        rs[:,1] *= theta
        rs[:,2] *= theta

        # Normalize values.
        rs += math.radians(62.8)
        rs[:,0] /= 2*math.radians(62.8)
        rs[:,1] /= 2*math.radians(62.8)
        rs[:,2] /= 2*math.radians(62.8)

        return rs

    def modify(frame, data):
        """ The user-defined modifier function """

        # Input:
        orientations = data.particles['Orientation']

        # Output:
        data.particles_.create_property('Color', data=quaternions_to_colors(orientations))