1 Introduction
This paper addresses the two problems of multimodal image registration and image segmentation within a single framework. Multimodal registration refers to registration of images acquired by different sensor/scanner types. It has been applied to many areas, e.g. medical images, microscopy images, and radar images, to combine information from different modalities and provide more comprehensive understanding about the true object. Image segmentation, the partitioning of an image into meaningful regions, is an important step in image analysis and understanding.
In this work, we focus on multimodal registration and segmentation as applied to scanning electron microscope (SEM) images of materials; the methods to be discussed are equally applicable to other multimodal images that share spatial structure. SEM techniques are widely used in materials science for material characterization, for example detection of defects that may cause fatigue when operating. Segmentation is of interest to map locations of grains, uniform regions occupied by continuous crystal lattices, since grain structure is a principal factor in determining the properties of a polycrystalline material such metallic or ceramic materials [1]. Multimodal registration is desired because different scanning electron modalities carry complementary information [2, 3]. For example, Backscattered Electrons (BSE) provide information about topography and local finescale surface texture [4] while Electron Backscatter Diffraction (EBSD) measures crystal orientation which is useful in locating grains and grain boundaries [5].
Multimodal registration is made challenging by the fact that images from different modalities may have different resolutions, values that lie in different spaces (e.g. scalars vs. vectors), and different levels of distortion. In SEM for instance, these differences are due to different electron beam geometries, sensors, and recording electronics. Furthermore, there is often no complete forward model that jointly characterizes the multimodal signals, nor a transformation model that adequately describes the distortion. In these circumstances, pixellevel registration methods
[6, 7, 8], i.e., those that establish correspondences between pixels, usually resort to interpolation, a somewhat adhoc solution that may bias the resulting images toward excessive smoothness. On the other hand, segmentation of multimodal images, if done independently for each modality using existing methods [9, 10, 5], may fail due to low contrast in some modalities and face difficulties in identifying correspondences between segmented regions from different modalities.In this work, we propose a coercive regionlevel approach to simultaneously register and segment images of different modalities that share similar spatial structure. The algorithm is initialized by segmenting one image by a standard method and coarsely mapping the result onto the other image. Then the two images are registered at the region level and further segmented through alternating minimization of a statisticallybased objective function. There are several advantages of our approach. First, the regionlevel approach is free of pixel value interpolation and its inherent assumptions. Second, it takes advantage of modalities with better discriminative power, improving the overall segmentation result. The approach also preserves region correspondences to facilitate data fusion in [11, 12]. Lastly, both registration and refinement of segmentation are driven by statistical models. In particular, we propose hypothesis tests to detect boundaries that are missed by the initial segmentation due to low contrast.
The paper is organized as follows. In Section 2, we describe the objective function, statistical models for data from different modalities, and optimization methods for the regionlevel registration algorithm. In Section 3, we focus on hypothesis testing for detecting missing boundaries. Section 4 shows experimental results for synthetic and real materials images and compares several different approaches. Section 5 concludes this paper.
2 Algorithm Framework and Models
2.1 Objective Function
We assume that there are two images from different modalities. A pixel location is represented as vector , where is the spatial domain for the th modality. The pixel value at is given by a function . Note that the values and may lie in different spaces. The regionlevel registration problem is to find partitions of , , where each segment is a collection of connected pixel locations and is the number of segments, to minimize the following objective function:
(1) 
where is the intramodal energy function that measures how well the segmentation fits the image data and is the intermodal energy function that coerces the segmentation results to be topologically similar, motivated by the fact that they share the same underlying physical structure. The parameter controls the relative importance of the two terms.
In this paper, we define the intermodal energy to be the number of segment boundaries that are present in one modality but not the other. This number is easily tracked because our algorithm maintains tight correspondences between segments in the two images. More generally, segment structure can be represented by a connected adjacency graph and the intermodal energy can be any function which measures the topological distance between two graphs. The intramodal energy is defined by the statistical models described in the following subsection.
2.2 Statistical Models for Pixel Values
In the materials context, each segment
corresponds to a grain, a continuous crystal lattice. Motivated by this, we assume that the observed values within a segment are similar and can be modeled by i.i.d. random variables following a distribution with the same parameters. In the sequel, the image modality subscript
is suppressed for simplicity. Let the probability density function (PDF) of the distribution for one modality be denoted by
, where represents the parameters specifying the model. The intramodal energy function in (1) given a set of segments is defined as:(2) 
where is the boundary of region with counterclockwise definition and
is the maximumlikelihood (ML) estimate for the parameters of region
. The first term is the negative loglikelihood of observations which penalizes grain inhomogeneity and the second term penalizes the boundary length, where controls the level of smoothness.In this paper, we consider the EBSD and BSE images of one section of a material as our input. Note that other image types can be used directly given properly defined statistical models. For BSE images, since the pixel values are scalars, the intensities in the same grain region are modeled by a univariate Gaussian , where
are the mean and variance of
. Notice that and are unknown parameters to be estimated from image data.For EBSD images, the pixel values characterize the local crystal orientation, which can be represented by Euler angles [13], Rodrigues vectors [14] or quaternions [15]. We choose the unitquaternion representation, i.e. a , the dimensional unit sphere in . This allows use of the von MisesFisher (VMF) distribution in directional statistics [16]
, a natural generalization of the multivariate Gaussian distribution to the sphere
(here ). However, the VMF distribution cannot be used directly due to crystal symmetry, which causes there to be more than one quaternion representation corresponding to a single crystal orientation. The ambiguity in representation may lead to a large artificial diversity of pixels within the same grain, resulting in an oversegmented result. To cope with this problem, we have proposed a VMF mixture distribution model which accounts for symmetry in our previous work [17]. To briefly describe the model, let be a group of quaternion matrices which define the symmetry actions that map a quaternion to its symmetric equivalents. The PDF of the pure VMF distribution is , where , is the mean direction, is the concentration parameter, and is the modified Bessel function of the first kind with order . The density function of the VMF mixture distribution is then given by(3) 
The parameters and
can be estimated from image data through the ExpectationMaximization algorithm derived in
[17].2.3 Optimization
(4)  
(5) 
where is the iteration index. Typically 2–3 iterations suffice.
To initialize the algorithm, the initial segmentation of the first modality, , is obtained by using a suitable image segmentation method. For example, the Voronoibased method in [18] can be applied to EBSD images and the Stabilized Inverse Diffusion Equation method to BSE images [9]. Since EBSD data provides crystal orientation which defines grain regions more accurately, we choose to start with EBSD segmentation in this paper. Next, to account for global misalignment and any resolution difference between the modalities, we determine an affine transformation by treating the material sample as a binary image and registering its outer boundary from one modality to the other using the Elastix toolbox [19]. The transformation is then used to map onto the other modality, yielding the initial segmentation . Note that due to localized distortions between the modalities, the initial segmentation may be misaligned with the image values as shown in Fig.1(a) and therefore needs to be registered.
Optimizing (4) and (5) is done in two steps. The first step is to consider splitting regions in the current segmentation by adding boundaries. In Section 3, we propose a hypothesis testing approach for this purpose based on the statistical model (2). The second step is to register the misaligned boundaries. Due to the fact that adjusting boundary positions does not change the topology of the segment structure, the intermodal energy function is not changed in this step, reducing (4) and (5) to the intramodal energy function alone, which is given by our statistical model. We use the Region Competition algorithm [20] to minimize . This algorithm applies gradient descent to move pixels comprising the boundaries along their respective normal directions. There are two forces driving the movement corresponding to the two terms in (2): the statistics force which comes from the distribution model for the pixel values, and the smoothing force which drives the boundary to have smaller curvature. More details are given in [20].
3 Hypothesis Tests for Missing Boundaries
This section elaborates upon the first step in solving (4) and (5), namely hypothesis testing to determine whether a region should be split into two based on the observed image values. We refer to this as the missing boundary problem. Recall that may come from the initial segmentation result from another modality and may not fit the current image data. Figure 1 shows examples of misalignment and a missing boundary. One can see that both of the situations have multimodal distributions of pixel values within the initially defined regions but only Fig.1(b) shows a missing edge that should be identified. Therefore, a region is declared as having a missing boundary if and only if it satisfies the following two conditions: (1) The pixel values are multimodally distributed. (2) The multimodal distribution is unlikely to be caused by misalignment. We develop two hypothesis tests for the two criteria. The first hypothesis test uses the Generalized Likelihood Ratio Test (GLRT) [21] to test whether the pixel values are multimodally distributed. The second hypothesis test differentiates misalignment from a missing boundary.
3.1 Hypothesis Test for MultiModality
Recall from Section 2.2 that the set of pixel values within a region are modeled by a distribution with unknown parameters , where are the observed pixel values in . Assume there exists a boundary which partitions into two subregions with parameters . The two hypotheses are : region is indeed a single region, i.e. ; and : consists of two regions. The GLRT has the following form:
(6)  
where are the ML estimates of the parameters under the null and alternative hypotheses and is the coefficient in (1). The GLRT can be viewed as a tradeoff between the improvement in the intramodal energy and the penalty of paid in the intermodal energy when inserting a boundary. The boundary length penalty is neglected here for simplicity but can be included easily.
In the following subsections, we derive the GLRT for univariate Gaussian and VMF distributions given the boundary . We only discuss the equal variance (concentration parameter) case due to the paper length constraint. These expressions supply the objective function, denoted as , to be maximized with respect to in (6). We use the Region Growing algorithm [20] to locate the optimal boundary . The algorithm partitions a region starting from two seed pixels and greedily adds neighboring pixels until all pixels in the region are chosen.
3.1.1 Multimodality Test for Univariate Gaussian Distribution
The GLR given boundary for testing mean equality between two Gaussian distributions has the following form [22]:
(7) 
where are the ML estimators of the variances under the null and alternative hypothesis and are the numbers of pixels in regions . The optimization of the boundary then takes the form
(8) 
where is the ML estimate of the mean in .
3.1.2 Multimodality Test for von MisesFisher Distribution
The VMF mixture distribution is reduced to single VMF through transoforming the samples by the symmetry operator towards the mean direction estimated by the EM algorithm. According to the derivation of the ML estimators in [23], has the following form:
(9) 
where and , . The optimization over is
(10) 
The last equality comes from the fact that and are monotonically increasing functions of .
3.2 Hypothesis Test for Misalignment
For regions that pass the previous multimodality test ( declared in (6)), we perform a second hypothesis test to determine whether the multimodal distribution is due to : boundary misalignment, or : a missing boundary. Since in most cases, misalignment causes only a small portion of pixels to differ from the majority, one naive test is to set a threshold on the ratio of the size of the smaller region to the whole region:
(11) 
However, since region size can vary over several orders of magnitude, the same absolute amount of misalignment (in pixels) can result in very different size ratios, making it hard to set a universal threshold. Therefore, we propose an adaptive threshold which incorporates region size. Boundary misalignment is modeled by a displacement in position (see Fig.2), where the displacements are bivariate Gaussian with zero mean and covariance and is the identity matrix. For simplicity, the region is modeled as a circle with radius (Fig. 2), where is the equivalent radius of region
. Based on these assumptions, the test statistic in (
11) can be formulated as the following function of given :The second line follows because is an increasing function. Since the displacement follows a Rayleigh distribution, given the user specified false positive rate , we set
where is the Rayleigh tail distribution. As a result, the threshold is adaptively determined by and the equivalent radius .
4 Experiments
4.1 Boundary Detection Accuracy on Simulated Data
In this section, we compare grain boundary detection performance on simulated EBSD and BSE images using three different approaches: A. Segment the BSE and EBSD images separately by suitable existing segmentation algorithms [10, 18]; B. Segment EBSD and register the boundaries onto BSE using a BSpline deformation model and the mutual information criterion [6, 7]; C. The proposed coercive registration/segmentation algorithm with .
The grain shapes in the testing data are taken from real microscopy images downloaded from BlueQuartz [24] and segmented by their Dream3D toolbox. For each slice, some of the grains are randomly selected and displaced to produce boundary misalignment according to the Gaussian displacement model with (pixels). This creates the ground truth boundaries for evaluation. The pixel values for BSE and EBSD are generated from Gaussian and VMF distributions with random mean and variance/concentration for each grain region.
To evaluate the boundary detection accuracy, we use the “overlapping rate”. Let be the set of boundary pixel locations with boundary width , which is obtained by image dilation with filter disk radius . The overlapping rate is defined as , where are the ground truth and estimated boundary.
Figure 3 shows the overlapping rate of the three approaches for different boundary widths. Independent segmentation has the worst performance since it does not make use of shared substructure between modalities. With Bspline registration, there is some improvement but it is still not satisfactory, especially for small . The proposed coercive registration approach with hypothesis testing is able to accurately register misaligned boundaries and detect missing edges. Therefore, it has much better boundary detection performance.
4.2 Results on Real Microscopy Data
We apply the proposed method to the IN100 data set which contains 170 slices of EBSD and BSE images of a Nibase alloy. Figure 4
shows one registration/segmentation result overlaid on the BSE image. The red lines are the initial boundaries obtained by the EBSD segmentation and affinetransformed to match BSE. The blue lines are the realigned boundaries and the green lines are the missing boundaries detected by the hypothesis tests. The initial red lines are misaligned with the BSE image values but are corrected by our registration algorithm. Using statistical hypothesis tests, we are also able to detect and locate missing boundaries in some grain regions. These results in real data demonstrate that the proposed approach can accurately register boundaries and segment grain regions at the same time.
5 Conclusion
In this work, we proposed a coercive registration/segmentation algorithm for multimodal images. The algorithm alternately utilizes information from one modality to help segment the image in the other modality, resulting in significant performance improvement in both modalities. The proposed hypothesis test based on statistical models of pixel values can accurately detect and locate missing boundaries between regions. Furthermore, our approach identifies and preserves all of the correspondences between regions in different modalities, which is important for fusing information after registration. The experiment results on simulated and real microscopy images show that our approach is able to effectively correct misaligned grain boundaries and detect missing boundaries within grain regions.
References
 [1] DM Shah and DN Duhl, “Effect of Minor Elements on the Deformation Behavior of NickelBase Superalloys,” Superalloys 1988, pp. 693–702, 1988.
 [2] Yun Wang, Hidehiko Kimura, Yoshiaki Akiniwa, and Keisuke Tanaka, “EBSDAFM hybrid analysis of fatigue slip system and crack initiation in polycrystalline metal under cyclic torsional loading,” in MicroNanoMechatronics and Human Science, 2005 IEEE International Symposium on. 2005, pp. 223–228, IEEE.
 [3] Heidi Nordmark, M Di Sabatino, M Acciarri, J Libal, S Binetti, EJ Ovrelid, JC Walmsley, and R Holmestad, “EBIC, EBSD and TEM study of grain boundaries in multicrystalline silicon cast from metallurgical feedstock,” in Photovoltaic Specialists Conference, 2008. PVSC’08. 33rd IEEE. 2008, pp. 1–6, IEEE.
 [4] Joseph Goldstein, Scanning electron microscopy and xray microanalysis, Kluwer Academic/Plenum Publishers, New York, 2003.
 [5] Adam J Schwartz, Mukul Kumar, Brent L Adams, and David P Field, Electron backscatter diffraction in materials science, vol. 2, Springer, 2009.
 [6] Frederik Maes, Andre Collignon, Dirk Vandermeulen, Guy Marchal, and Paul Suetens, “Multimodality image registration by maximization of mutual information,” Medical Imaging, IEEE Transactions on, vol. 16, no. 2, pp. 187–198, 1997.
 [7] D. Rueckert, L. I Sonoda, C. Hayes, D. L G Hill, M. O. Leach, and D.J. Hawkes, “Nonrigid registration using freeform deformations: application to breast MR images,” IEEE Transactions on Medical Imaging, vol. 18, no. 8, pp. 712–721, Aug. 1999.

[8]
Richard Szeliski and James Coughlan,
“Splinebased image registration,”
International Journal of Computer Vision
, vol. 22, no. 3, pp. 199–218, 1997.  [9] HsiaoChiang Chuang, Landis M Huffman, Mary L Comer, Jeff P Simmons, and Ilya Pollak, “An automated segmentation for nickelbased superalloy,” in Image Processing, 2008. ICIP 2008. 15th IEEE International Conference on. 2008, pp. 2280–2283, IEEE.
 [10] R. Nock and F. Nielsen, “Statistical region merging,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 26, no. 11, pp. 1452–1458, Nov. 2004.
 [11] Gemma Piella, “A general framework for multiresolution image fusion: from pixels to regions,” Information fusion, vol. 4, no. 4, pp. 259–280, 2003.
 [12] Tao Wan, Nishan Canagarajah, and Alin Achim, “Segmentationdriven image fusion based on alphastable modeling of wavelet coefficients,” Multimedia, IEEE Transactions on, vol. 11, no. 4, pp. 624–633, 2009.
 [13] David Eberly, “Euler angle formulas,” Geometric Tools, LLC, Technical Report, 2008.
 [14] A Morawiec and DP Field, “Rodrigues parameterization for orientation and misorientation distributions,” Philosophical Magazine A, vol. 73, no. 4, pp. 1113–1130, 1996.
 [15] Simon L. Altmann, “Rotations, quaternions, and double groups,” 2005.
 [16] Kantilal Varichand Mardia and Peter E. Jupp, “Directional statistics,” 1999.
 [17] Y. Chen, D. Wei, G. Newstadt, M. De Graef, J. Simmons, and A. Hero, “Parameter Estimation in Spherical Symmetry Groups,” Signal Processing Letters, IEEE, vol. PP, no. 99, pp. 1–1, 2015.
 [18] Florian Bachmann, Ralf Hielscher, and Helmut Schaeben, “Grain detection from 2d and 3d EBSD data—Specification of the MTEX algorithm,” Ultramicroscopy, vol. 111, no. 12, pp. 1720–1733, 2011.
 [19] Stefan Klein, Marius Staring, Keelin Murphy, Max A Viergever, and Josien PW Pluim, “Elastix: a toolbox for intensitybased medical image registration,” Medical Imaging, IEEE Transactions on, vol. 29, no. 1, pp. 196–205, 2010.
 [20] S.C. Zhu, T. S. Lee, and AL. Yuille, “Region competition: unifying snakes, region growing, energy/Bayes/MDL for multiband image segmentation,” in , Fifth International Conference on Computer Vision, 1995. Proceedings, June 1995, pp. 416–423.
 [21] Jerzy Neyman and Egon S Pearson, On the problem of the most efficient tests of statistical hypotheses, Springer, 1992.

[22]
Theresa K Seize,
“Student’s ttest,”
Southern Medical Journal, vol. 70, no. 11, pp. 1299, 1977.  [23] Inderjit S Dhillon and Suvrit Sra, “Modeling data using directional distributions,” Tech. Rep., Technical Report TR0306, Department of Computer Sciences, The University of Texas at Austin. URL ftp://ftp. cs. utexas. edu/pub/techreports/tr0306. ps. gz, 2003.
 [24] “BlueQuartz Software,” July 2003.
Comments
There are no comments yet.