Visualize the grain fitting results

[cjh1]

@joelvbernier Is going to provide some code examples.

[joelvbernier]

Ok -- first let me tabulate the anatomy of the grains.out file. There are 21 columns; they represent:

0    grain ID
1    completeness
2    chi^2
3    exp_map_c[0]      
4    exp_map_c[1] 
5    exp_map_c[2]       
6    t_vec_c[0]
7    t_vec_c[1]
8    t_vec_c[2]
9    inv(V_s)[0,0]
10   inv(V_s)[1,1]
11   inv(V_s)[2,2]
12   inv(V_s)[1,2]*sqrt(2)
13   inv(V_s)[0,2]*sqrt(2)
14   inv(V_s)[0,2]*sqrt(2)
15   ln(V_s)[0,0]
16   ln(V_s)[1,1]
17   ln(V_s)[2,2]
18   ln(V_s)[1,2]
19   ln(V_s)[0,2]
20   ln(V_s)[0,1]

The 12 parameters in the slice [:, 3:15] are the grain parameters that are used by the various models and methods, 3 orientation parameters, 3 center-of-mass coordinates (x, y, z), and 6 strain parameters as the scaled components of the inverse stretch tensor in the sample frame (not particularly useful by itself...). For convenience, I convert those last 6 grain parameters to strain tensor components in the sample frame in the last 6 columns.

For a start, we will want to use the COM coords for a 3-d scatter plot (cols [:, 6:9]), where the colormap is selectable from the following scalars:

[:, 0]      completeness ratio
[:, 1]      goodness of fit
[:, 3:6]    grain orientation *
[:, 15]     xx strain
[:, 16]     yy strain
[:, 17]     zz strain
[:, 18]     yz strain
[:, 19]     xz strain
[:, 20]     xy strain
[:, 15:21]  equivalent strain *

the * quantities require calculation. I'll have to get you a colormap generator for orientation -- it depends on the symmetry group of the material (@saransh13: do you have code to do this?). The equivalent strain is a simple calculation:

from hexrd.matrixutil import vecMVToSymm

ngrains = len(grains_table)  # loadtxt on grains.out

eqv_strain = np.zeros(ngrains)
for i, grain in enumerate(grains_table):
    emat = vecMVToSymm(grain[15:21], scale=False)
    eqv_strain[i]= 2.*np.sqrt(np.sum(emat*emat))/3.

[cjh1]

@joelvbernier Thanks for this @johnkit It going to push a PR to refactor the grain fitting code so it returns the results. It would be great if we can get the colormap generator and example plot code soon.

[cjh1]

@joelvbernier ping

[cjh1]

@joelvbernier, @johnkit has pushed the fit grains work, we will have to put this to one side until we have an example of the viz.

[cjh1]

@joelvbernier WiIl add the mplot 3d code for this.

[cjh1]

@joelvbernier To provide HDF5 file and code for visualizing a single grain/spot data.

[johnkit]

@joelvbernier @cjh1 Here is a look at the beginning of a "fit grains results dialog". There is more to do, but I thought I should show everyone where I am heading and give you a chance to comment/redirect.

[joelvbernier]

@johnkit -- so far so good!

One thing I forgot to mention yesterday would be the ability to visualize the Bragg peak ("spot") data associated with a given grain. These are the data structures that get distilled into the spots_<grain_id>.out files for each detector, and are the output of hexrd.instrument.HEDMInstrument.pull_spots() (here). So it looks like I have a kwarg intended for saving the full data structures that isn't hooked up (yet...). I had written an option to dump the results to HDF5, and below is a plotting script that uses the HDF5 data structure to make montage plots for all spots of a given type (the second column of the spots files). Here is an example of the structure currently getting saved to HDF5: And here are the plotting results from the command calling the attached script:

(hexrdgui) vpna-e-128-15-239-171:dye-869-1 joel$ python spot_montage.py --help
usage: spot_montage.py [-h] [-t THRESHOLD] hdf5_archive gvec_id

Montage of spot data for a specifed G-vector family

positional arguments:
  hdf5_archive          hdf5 archive filename
  gvec_id               unique G-vector ID from PlaneData

optional arguments:
  -h, --help            show this help message and exit
  -t THRESHOLD, --threshold THRESHOLD
                        intensity threshold
(hexrdgui) vpna-e-128-15-239-171:dye-869-1 joel$ python -i spot_montage.py results_fd1-q-1/grain_00004_nearest.hdf5 1

First, here is an example script that calls the pull_spots() method from a config file:

# manual call to pull_spots for dumping hdf5
try:
    import dill as cpl
except(ImportError):
    import pickle as cpl

import glob

import numpy as np

import os

import yaml

from hexrd import config
from hexrd import instrument
from hexrd import imageseries
from hexrd.imageseries.omega import OmegaImageSeries

# %%
# =============================================================================
# START USER INPUT
# =============================================================================

cfg_filename = 'config_fd1-q-1.yaml'

grain_id = 0
tol_loop_idx = -1
block_number = 0    # ONLY FOR TIMESERIES CFG
interp = 'nearest'

# =============================================================================
# END USER INPUT
# =============================================================================

# %%
cfg = config.open(cfg_filename)[block_number]

# load instrument
instr = cfg.instrument.hedm

# load plane data
plane_data = cfg.material.plane_data

imsd = cfg.image_series

# %%
tth_tol = cfg.fit_grains.tolerance.tth[tol_loop_idx]
eta_tol = cfg.fit_grains.tolerance.eta[tol_loop_idx]
ome_tol = cfg.fit_grains.tolerance.omega[tol_loop_idx]
npdiv = cfg.fit_grains.npdiv
threshold = cfg.fit_grains.threshold

eta_ranges = np.radians(cfg.find_orientations.eta.range)
ome_period = np.radians(cfg.find_orientations.omega.period)
# %%
grains_file = os.path.join(cfg.analysis_dir, 'grains.out')

grains_table = np.loadtxt(grains_file, ndmin=2)
grain = grains_table[grain_id]
grain_params = grain[3:15]

# %%
# =============================================================================
# Calls to output
# =============================================================================

# HDF5
complvec, results = instr.pull_spots(
    plane_data, grain_params,
    imsd,
    tth_tol=tth_tol,
    eta_tol=eta_tol,
    ome_tol=ome_tol,
    npdiv=npdiv, threshold=threshold,
    eta_ranges=eta_ranges,
    ome_period=ome_period,
    dirname=cfg.analysis_dir, filename='grain_%05d_%s' % (grain_id, interp),
    save_spot_list=False, output_format='hdf5',
    quiet=True, check_only=False, interp=interp)

# %% for referene, this is the call to txt output
# complvec, results = instr.pull_spots(
#     plane_data, grain_params,
#     imsd,
#     tth_tol=tth_tol,
#     eta_tol=eta_tol,
#     ome_tol=ome_tol,
#     npdiv=npdiv, threshold=threshold,
#     eta_ranges=eta_ranges,
#     ome_period=ome_period,
#     dirname=cfg.analysis_dir, filename='spots_%05d.out' % grain_id,
#     save_spot_list=False,
#     quiet=True, check_only=False, interp=interp)

and here it the plotting script with CLI

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 19 15:29:27 2017

@author: bernier2
"""

import argparse

import numpy as np
import h5py
from matplotlib import pyplot as plt

"""
# UNCOMMENT IF YOU HAVE A SANE LATEX ENV AND WANT NICE FIG LABELS
#
# Options
params = {'text.usetex': True,
          'font.size': 14,
          'font.family': 'mathrm',
          'text.latex.unicode': True,
          'pgf.texsystem': 'pdflatex'
          }
plt.rcParams.update(params)
"""

plt.ion()

def montage(X, colormap=plt.cm.inferno, show_borders=True,
            title=None, xlabel=None, ylabel=None,
            threshold=None, filename=None):
    m, n, count = np.shape(X)
    img_data = np.log(X - np.min(X) + 1)
    if threshold is None:
        threshold = 0.
    else:
        threshold = np.log(threshold - np.min(X) + 1)
    mm = int(np.ceil(np.sqrt(count)))
    nn = mm
    M = np.zeros((mm * m, nn * n))

    # colormap
    colormap.set_under('b')

    fig, ax = plt.subplots()
    image_id = 0
    for j in range(mm):
        sliceM = j * m
        ax.plot()
        for k in range(nn):
            if image_id >= count:
                img = np.nan*np.ones((m, n))
            else:
                img = img_data[:, :, image_id]
            sliceN = k * n
            M[sliceM:sliceM + m, sliceN:sliceN + n] = img
            image_id += 1
    # M = np.sqrt(M + np.min(M))
    im = ax.imshow(M, cmap=colormap, vmin=threshold, interpolation='nearest')
    if show_borders:
        xs = np.vstack(
            [np.vstack([[n*i, n*i] for i in range(nn+1)]),
             np.tile([0, nn*n], (mm+1, 1))]
        )
        ys = np.vstack(
            [np.tile([0, mm*m], (nn+1, 1)),
             np.vstack([[m*i, m*i] for i in range(mm+1)])]
        )
        for xp, yp in zip(xs, ys):
            ax.plot(xp, yp, 'c:')
    if xlabel is None:
        ax.set_xlabel(r'$2\theta$', FontSize=14)
    else:
        ax.set_xlabel(xlabel, FontSize=14)
    if ylabel is None:
        ax.set_ylabel(r'$\eta$', FontSize=14)
    else:
        ax.set_ylabel(ylabel, FontSize=14)
    ax.axis('normal')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['left'].set_visible(False)
    cbar_ax = fig.add_axes([0.875, 0.155, 0.025, 0.725])
    cbar = fig.colorbar(im, cax=cbar_ax)
    cbar.set_label(r"$\ln(intensity)$", labelpad=5)
    ax.set_xticks([])
    ax.set_yticks([])
    if title is not None:
        ax.set_title(title, FontSize=18)
    if filename is not None:
        fig.savefig(filename, bbox_inches='tight', dpi=300)
    return M

def plot_gvec_from_hdf5(fname, gvec_id, threshold=0.):
    """
    """
    f = h5py.File(fname, 'r')

    for det_key, panel_data in f['reflection_data'].items():
        for spot_id, spot_data in panel_data.items():
            attrs = spot_data.attrs
            if attrs['hkl_id'] == gvec_id:
                # grab some data
                tth_crd = np.degrees(spot_data['tth_crd'])
                eta_crd = np.degrees(spot_data['eta_crd'])
                intensities = np.transpose(
                        np.array(spot_data['intensities']),
                        (1, 2, 0)
                )

                # make labels
                figname = r'Spot %d, ' % attrs['peak_id'] \
                    + r"detector '%s', " % det_key \
                    + r'({:^3} {:^3} {:^3})'.format(*attrs['hkl'])
                xlabel = r'$2\theta\in(%.3f, %.3f)$' \
                    % (tth_crd[0], tth_crd[-1])
                ylabel = r'$\eta\in(%.3f, %.3f)$' \
                    % (eta_crd[0], eta_crd[-1])

                # make montage
                montage(intensities, title=figname,
                        xlabel=xlabel, ylabel=ylabel,
                        threshold=threshold)
                pass
            pass
        pass
    f.close()
    return

# =============================================================================
# %% CMD LINE HOOK
# =============================================================================

if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Montage of spot data for a specifed G-vector family")

    parser.add_argument('hdf5_archive',
                        help="hdf5 archive filename",
                        type=str)
    parser.add_argument('gvec_id',
                        help="unique G-vector ID from PlaneData",
                        type=int)

    parser.add_argument('-t', '--threshold',
                        help="intensity threshold",
                        type=float, default=0.)

    args = parser.parse_args()

    h5file = args.hdf5_archive
    hklid = args.gvec_id
    threshold = args.threshold

    plot_gvec_from_hdf5(h5file, hklid, threshold=threshold)

Not sure the best way to do this vis -- maybe 3-d isosurfaces around the spots in the full angle space (tth, eta, ome), like a 3-d version of the polarview? The montages could be packaged into a GridSpec as a start.

As far as wiring this in, we could have a right-click menu on the grains table to execute the HDF5-based pull_spots, and just load the results for plotting... Thoughts?

[psavery]

@joelvbernier @saransh13 Hi Joel.

Does the equation for equivalent strain work for stress as well? Like the following:

eqv_stress = np.zeros(ngrains)
for i in range(ngrains):
    m = vecMVToSymm(stress[i], scale=False)
    eqv_stress[i]= 2.*np.sqrt(np.sum(m*m))/3.

[joelvbernier]

For stress it is 3/2 instead of 2/3. There is one other thing, however; the variable m should be the deviator of the stress.

Sent from my iPhone

On Nov 9, 2020, at 18:00, Patrick Avery notifications@github.com wrote:  @joelvbernier @saransh13 Hi Joel.

Does the equation for equivalent strain work for stress as well? Like the following:

eqv_stress = np.zeros(ngrains) for i in range(ngrains): m = vecMVToSymm(stress[i], scale=False) eqv_stress[i]= 2.np.sqrt(np.sum(mm))/3. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[psavery]

@joelvbernier Can you show us how to make m the deviator of the stress?

[donald-e-boyce]

To make the deviator, you subtract the spherical/dilational part so that the trace becomes zero.

m = \sigma - (1/3) trace(\sigma) I

where I is the identity.

On Tue, Nov 17, 2020 at 9:01 PM Patrick Avery notifications@github.com wrote:

@joelvbernier https://github.com/joelvbernier Can you show us how to make m the deviator of the stress?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/HEXRD/hexrdgui/issues/429#issuecomment-729326859, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAIWT6K55PJFSIPZ62VZKWTSQMTGRANCNFSM4PD45J3Q .

[psavery]

@donald-e-boyce Okay, so this is what I have for the code now.

eqv_stress = np.zeros(ngrains)
for i in range(ngrains):
    sigma = vecMVToSymm(stress[i], scale=False)
    deviator = sigma - (1/3) * np.trace(sigma) * np.identity(3)
    eqv_stress[i]= 3 * np.sqrt(np.sum(deviator**2)) / 2

Let me know if it is incorrect in any way.