Hierarchical Smoothing

Description

Thie repo contains multiple implementations of the hierarchical smoothing algorithm applicable to voxelated meshes of interface networks. Primarily directed at users of microstructure analysis software, specifically DREAM.3D.

WARNING

It is strongly recommended to compile the C++ version and call it in a Matlab/Octave routine (wrappers are available with the code). The Python and Matlab source codes are basically prototyped versions, and haven't been updated in years. They are definitely slow, not at all optimized and possibly buggy. What Python and Matlab do in a few hours for a large, high-resolution microstructure, C++ can do in a few seconds. I can help with the compilation on a Linux machine; please go through previous issues or open a new one if you run into trouble.

Motivation

1. Faithfulness to mesoscopically smooth interfaces
2. Minimum user interference in deciding smoothing parameters
3. Ability to automate over an entire polycrystal volume
4. Proper treatment of topological features like grain boundary interiors, triple lines and quad points

Reference (citations are greatly appreciated!)

S. Maddali, S. Ta'asan, R. M. Suter, Topology-faithful nonparametric estimation and tracking of bulk interface networks, Computational Materials Science 125, 328-340 (2016).

Software requirements

Import/Export of data from HDF5

1. DREAM.3D v6.2 or higher (to import and export ASCII data) OR
2. HDFView to export from existing DREAM.3D data file

Matlab

1. Matlab R2013a or higher (extensive use of the triangulation object)

Python

1. NumPy (>=1.11.2)
2. SciPy (>=0.18.1)

C++

1. Eigen - a templated C++ library for linear algebra.
2. libIGL - a set of very useful algorithms and routines written using Eigen.

Tutorials

These tutorials do not cover the import and export of the required microstructure data relative to DREAM.3D. Included is a sample data set (in the examples/ex2 directory) containing the mesh of a 794-grain microstructure volume. The following steps illustrate the form of the input and implementation of the central HierarchicalSmooth routine. In each tutorial, the working directory is assumed to be that which contains the corresponding source code.

• Matlab tutorial
• Python tutorial
• C++ tutorial

Matlab usage

• The input array xDat contains a set of N points in 3D in the form of a 3xN array. These represent voxelated sample points of the GB network.

  xdat = load( '../../examples/ex2/SharedVertexList.txt' );
  xdat = xdat';

• The input array tri contains the Delaunay triangulation of the points in xDat in the form of an Mx3 integer array. A +1 offset is applied in order to accommodate the difference between DREAM.3D's 0-indexing and Matlab's 1-indexing.

  tri = 1 + load( '../../examples/ex2/SharedTriList.txt' );

• The Quick Surface mesh filter in DREAM.3D assigns a unique integer ID to every grain in the volume. The following input array serves to identify the grain boundary that a particular triangular element of array tri belongs to, in the form of an Mx2 integer array. The integers in each row denote the grains on either side of that particular triangular patch of grain boundary.

  fl = load( '../../examples/ex2/FaceLabels.txt' );

• By default DREAM3D assigns the 'grain ID' of -1 to the top and bottom surfaces of the imaged volumes, and 0 to the sides respectively. These buffer regions are included in the array fl as well, for example the row [ 30, -1 ] in fl denotes the 'grain boundary' between the grain 30 and the top/bottom surface of the volume. These are unnecessary inputs to the smoothing routine and need to be filtered out to save time.

  f = find( any( fl==-1 | fl==0, 2 ) );
  fl( f, : ) = [];
  tri( f, : ) = [];

• The final required input is the classification of the surface mesh nodes as belonging to surface interiors, triple lines or quad points.

  ntype = load( '../../examples/ex2/NodeType.txt' );
  ntype = ntype';

• The next line runs the smoothing routine with these arrays as input. This could take a long time even for a few hundreds of grains; the Matlab code is not optimized for speed. The log file generated by default is Smooth.Default.log and can be monitored in real time.

  xsmooth = HierarchicalSmooth( xdat, tri, fl, ntype );

  Optional arguments for HierarchicalSmooth are described in the script file. Finally, the script Test.m plots a comparison between the smoothed and unsmoothed versions of a grain of the user's choice.

Python usage

The Python functions were written to mirror the corresponding functions in the Matlab source code as much as possible. Specifically, the primary function HierarchicalSmooth and its auxiliary functions are written in a file HierarchicalSmooth.py which is loaded as a module. A bare-bones version of Matlab's indispensible ismember function is implemented in Base.py. Further, the basic functionality of Matlab's triangulation object is implemented in Triangulation.py. The following steps are the essence of the test script Src/Python/TestSuite/SmoothVolume.py.

• First the required namespaces are loaded:

  import numpy as np
  import HierarchicalSmooth as hs

• The relative path of the input data files is stored in a string:

  filePrefix = '../../examples/ex2'

• The following input arrays are loaded and cast into the appropriate types and dimensions:

  • Sample vertices as a $3 \times N$ array of floats:

    P = np.loadtxt( filePrefix + '/SharedVertexList.txt', dtype=float ).T

  • Mesh triangulation as a $M \times 3$ integer array. Note that unlike for Matlab, there's no need to account for the zero offset.

    tri = np.loadtxt( filePrefix + '/SharedTriList.txt', dtype=int )

  • $M \times 2$ integer array indicating the ID of the grain on either side of each triangular patch. Note that the direction in going from grain 1 to grain 2 corresponds to the right-handed winding sequence of the points of the mesh element. This is the default behavior of Dream.3d and is critical in obtaining the connectivity of the free boundary nodes of each interface.

    nFaces = np.loadtxt( filePrefix + '/FaceLabels.txt', dtype=int )

  • An integer array of length $N$ denoting the type of each node in the array P. According to Dream.3d convention, values of 2, 3 and 4 correspond respectively to interior sample points, triple junctions and quad points in the volume interior, and 12, 13, 14 corresponds to the same topological features on the volume edges.

    nType = np.loadtxt( filePrefix + '/NodeType.txt', dtype=int )

• Finally, the smoothing is achieved by running the command:

  PS, bPS = hs.HierarchicalSmooth(
  xPoints=P, tri=tri,
  nFaceLabels=nFaces,
  nNodeType=nType
  )

  Note that this does not ignore the surfaces on the exterior of the volume.

• Owing to extensive use of dictionary lookups, the Python code is faster by far than the Maltab code. Performance varies with system specifics; the 794 grains in the examples/ex2 data set were smoothed in about 1.5 hours on a laptop running Ubuntu 16.04 with 4GB RAM and a dual core hyper-threaded 2.40 GHz Intel i3 CPU

• There are optional arguments to HierarchicalSmooth that deal with the internal interval bisection threshold and the maximum number of iterations, as well as the text log file to which to redirect sys.stdout if needed.

C++ usage

The C++ implementation of HierarchicalSmooth is by far the fastest and most efficient of the three available on this repo, and it is highly recommended that you use this instead of the much older and slightly buggy Matlab and Python code. It is built on top of Eigen, which is also the linear algebra package of choice for DREAM.3D.

This section describes how to build the compiled libraries and also how to generate wrappers that allow them to be used in Matlab and Octave. This has been tested for 64-bit Linux systems. The source code is found in the Src/Cpp directory and includes a simple hand-written Makefile to generate the shared and static libraries. The procedure for building the project is as follows:

• Download or link the Eigen and libIGL source codes into the working directory with the HierarchicalSmooth source code and Makefile, under the names Eigen and libigl respectively.

• Run the makefile from the command line: make. his should generate a static library libhsmooth.a and a shared (dynamically linked) library libhsmooth.so.

• If you are using Matlab, run the CreateMatlabMex.m script from the Matlab shell. This should generate a compiled binary HierarchicalSmoothMatlab.mexa64 on 64-bit Linux machines. This function can be used in the same manner as the Matlab implementation of HierarchicalSmooth described above, except for the slightly different argument list:

  xsmooth = HierarchicalSmoothMatlab( tri, xdat, fl, ntype );

• Don't forget to remove the buffer mesh elements on the faces of the box before hand:

  f = find( any( fl > 0, 2 ) );
  tri = tri( f, : );
  fl = fl( f, : );

  ... and it will run enormously faster than the Matlab code.

• If you are using Octave, the procedure is very similar, except you should run the CreateOctaveMex.m script from the Octave shell. This should generate a binary HierarchicalSmoothOctave.mex on 64-bit Linux machines. This is implemented in the same manner as the Matlab binary.

  xsmooth = HierarchicalSmoothOctave( tri, xdat, fl, ntype );

• The array dimensions expected by the Matlab/Octave code are:

  • tri: N ⨉ 3
  • xdat: 3 ⨉ N
  • fl: N ⨉ 2
  • ntype: N ⨉ 1

  Failing this Eigen will complain, and you might experience a Matlab/Octave crash, because currently there is no mechanism to handle Eigen errors. A good way to deal with this is to call the compiled MEX binary inside a custom Matlab/Octave function, within which all the necessary error handling may be done.

Acknowledgements

1. Anthony Rollett (Dept of Materials Science and Engineering, Carnegie Mellon University)
2. David Menasche (Dept of Physics, Carnegie Mellon University)
3. Michael Jackson (Bluequartz Software)