Open Access Research

High resolution micrograph synthesis using a parametric texture model and a particle filter

Ramakrishna Tipireddy1*, Roger Ghanem1, Somnath Ghosh2 and Daniel Paquet3

Author Affiliations

1 Department of Civil Engineering at University of Southern California, Los Angeles, CA, USA

2 Department of Civil Engineering at Johns Hopkins University, Baltimore, MD, USA

3 Hydro-Québec Research Institute (IREQ), 1800, boul. Lionel-Boulet, Varennes, J3X 1S1, Québec, Canada

For all author emails, please log on.

Integrating Materials and Manufacturing Innovation 2013, 2:2  doi:10.1186/2193-9772-2-2


The electronic version of this article is the complete one and can be found online at: http://www.immijournal.com/content/2/1/2


Received:21 February 2013
Accepted:29 April 2013
Published:10 June 2013

© 2013 Tipireddy et al.; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We present a methodology for synthesizing high resolution micrographs from low resolution ones using a parametric texture model and a particle filter. Information contained in high resolution micrographs is relevant to the accurate prediction of microstructural behavior and the nucleation of instabilities. As these micrographs may be tedious and uneconomical to obtain over an extended spatial domain, we propose a statistical approach for interpolating fine details over a whole computational domain starting with a low resolution prior and high resolution micrographs available only at a few spatial locations. As a first step, a small set of high resolution micrographs are decomposed into a set of multi-scale and multi-orientation subbands using a complex wavelet transform. Parameters of a texture model are computed as the joint statistics of the decomposed subbands. The synthesis algorithm then generates random micrographs satisfying the parameters of the texture model by recursively updating the gray level values of the pixels in the input micrograph. A density-based Monte Carlo filter is used at each step of the recursion to update the generated micrograph, using a low resolution micrograph at that location as a measurement. The process is continued until the synthesized micrograph has the same statistics as those from the high resolution micrographs. The proposed method combines a texture synthesis procedure with a particle filter and produces good quality high resolution micrographs.

Keywords:
Micrograph synthesis; Parametric texture model; Particle filters; Bayesian methods; Wavelet transform

Background

Recent capabilities in the multi-scale modeling and simulation of material behavior have been matched by technologies for sensing and characterizing materials at distinct spatial scales. Several technical challenges ensue in the process of pairing multi-scale models to multi-scale data, one of which is addressed in this paper namely that of characterizing knowledge about multi-scale material structure in a manner that is conducive to probabilistic conditioning and updating.

We will consider experimental data provided through micrographs, which are digital images of the microstructure taken through a microscope, showing a magnified view of the material domain, and allowing digital representation of fine-scale features. In this paper, the terms “image”, “texture” and “micrograph” are used interchangeably. High resolution microstructural information is required as a prerequisite for the multi-scale modeling and analysis of composite structures, and in particular to track the nucleation of instabilities in material behavior. As it remains expensive to experimentally obtain high resolution micrographs throughout a computational domain where coarse-scale behavior is being assessed, procedures for gleaning equivalent information from lower resolution images have demonstrated their value in computational materials science [1,2]. This paper contributes to that body of knowledge by developing a statistical procedure that updates low resolution priors with spatially localized high resolution data using a parametric texture model and a particle filter. Three approaches are commonly used, in light of their simplicity, for generating high resolution images from lower resolution ones. These are based on nearest neighbor, bilinear and cubic interpolation methods [3]. In the nearest neighbor interpolation method, a piecewise constant interpolation maps the low resolution to the high resolution image, typically producing blocky images. Bilinear and cubic interpolation methods smoothen these discontinuities, while still allowing for a blur in the high resolution image. Methods based on wavelet based interpolation with gradient based enhancement (WIGE) and higher-order polynomial interpolation with gradient based enhancement (PIGE) techniques have been proposed and successfully adapted to problems of computational multi-scale material modeling [1,2]. In WIGE, high resolution images are obtained by interpolating low resolution images using wavelet basis with the interpolated images enhanced using gradient based methods calibrated from a few high resolution images. In PIGE, higher order polynomial bases are used for interpolation. In these methods, gradient based enhancement accounts for the position of calibration images relative to the spatial location of the simulated image. Both WIGE and PIGE methods are shown to be very effective in synthesizing high resolution images from low resolution ones with the aid of few calibrating high resolution images. In [4], digital micrographs are simulated by representing the microstructures with two-point correlation function. Characterization and simulation of three-dimensional microstructures is also discussed in [5-8]. In these approaches, the synthesized micrographs are deterministic and do not account for uncertainty in experimentally obtained micrographs. Other common and effective approaches for microstructure modeling include Markov random field (MRF) models whereby a Gibbs distribution is used to represent a microstructure and Gibbs sampling method is used for micrograph synthesis [9,10]. Parameters of the Markov random field model are typically estimated from a small number of calibrating images and the micrographs are synthesized using Markov chain Monte Carlo (MCMC) methods. A large number of parameters is required to produce high resolution features with MRF models and estimating them with credible confidence from few calibrating images remains a challenge. Numerous variations on and enhancements to standard MRF have adapted it to specific situations which have included deterministic relaxation for image classification [11,12], texture segmentation based on wavelets and hidden Markov tree models [13]. A Markov random field model for image segmentation that combines both color and texture features has also been proposed [14] as well as image imputing based on hierarchical Bayesian models to reconstruct image from partially observed multivariate data [15]. Quantification of microstructural variance in the formalism of stochastic processes is developed in [16]. Algorithms based on sequential testing and feature learning have also been developed for face detection applications [17] while edge preserving image restoration with Huber-Markov random field models has been used to restore degraded images [18], and procedures based on Gibbs reaction-diffusion have yielded a theory for learning generic prior models from a set of natural images [19]. Maximum a posteriori (MAP) estimation methods have also been used for image restoration [20] sometimes in conjunction with adaptive parameter estimation [21]. An approach based on filters, random fields and maximum entropy (FRAME) [22] has also been very effective in capturing high resolution features in the image, but remains computationally expensive. Its equivalence with Julesz ensembles has also been investigated [23]. Finally, a parametric texture models based on an over complete complex wavelet transform was also used [24] to model micrographs based on the steerable pyramid decomposition of texture ([25,26]). In this paper, we rely on a similar texture model [24] with parameters obtained from a few high resolution calibrating images and low resolution image at the synthesis location. A particle filter is used to estimate the high resolution micrograph at a location given the low resolution micrograph at that location.

In our work, we use a Bayesian framework where the experimental low resolution micrograph play the role of measurement used to update a prior model synthesized from a very small set of high resolution micrographs. We use a particle filter, namely a density-based Monte Carlo filter to implement the Bayesian formulation. In [27], Kalman filter is used for system identification based on multiple static and dynamic test data. Kalman filter methodology can not be used in our problem, because the state equation in the filtering problem is modeled as an implicit micrograph synthesis algorithm, where as the Kalman filter requires the state equation to be a linear equation. Monte Carlo based particle filters [28-32] are used for system identification of dynamical system. In our paper a density-based Monte Carlo filter proposed in [32] is used. Synthesis of high resolution micrographs from low resolution images with parametric texture model proposed in [24] along with the density-based Monte Carlo filter [32] is demonstrated in our paper. The main contribution of the present work is to compute the parameters of the texture model from high resolution calibrating images and low resolution images and combine particle filter with a texture synthesis model to characterize the high resolution micrograph. The paper is organized as follows. Motivation of the work is described with the problem statement in subsection “Problem statement”. In subsection “Parametric texture model”, the parametric model along with the synthesis algorithm is described. Subsection “Density-based Monte Carlo filter” illustrates the density-based Monte Carlo filter. Numerical experiments are provided in section “Results and discussion”.

Problem statement

Figure 1 shows a low resolution microscopic image of cast aluminum alloy W319. At locations A and B marked on the low resolution image, high resolution images are obtained from scanning electron microscopy.

thumbnailFigure 1. Low resolution micrograph. Low resolution, low magnification digital image of cast aluminum alloy W319, for which high resolution images are available at points A and B.

Figure 2 shows experimentally taken high resolution micrographs using scanning electron microscope at locations A and B. Figure 3 shows magnified low resolution micrographs available at locations A, B and C. Our objective is to simulate high resolution micrographs at all the locations where the high resolution microscopic images are not available. The proposed high resolution micrograph synthesis methodology is demonstrated by generating a high resolution micrograph at location C. The procedure can be repeated to simulate the high resolution micrographs at other locations.

thumbnailFigure 2. High resolution digital image of cast aluminum alloy W319 at a and b.

thumbnailFigure 3. Low resolution digital image of cast aluminum alloy W319 at a, b and c.

Starting with an image x0 that represents the low-resolution micrograph at location C, our aim is to synthesize an image x at the same location, representing the high-resolution micrograph, that is constrained, in some sense, by our knowledge of high resolution data at locations A and B and low resolution data at C. In order to complete the statement of the problem, we must specify the nature of the constraints (i.e. the specific functionals or features that describe our knowledge of the high resolution data). In what follows, a texture model will be used to provide a context in which features of a micrograph are quantified. This model will be constructed from available high resolution data and presumed to be valid everywhere else. Specifically, the texture model will be identified with the joint statistics of the subbands in a particular wavelet decomposition, namely, a steerable pyramid decomposition. Associated with this model is a mapping f that maps x0 into x such that x has the joint statistics corresponding to the model.

The texture model maps an initial micrograph x0 into a high resolution micrograph x. Let <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M1">View MathML</a> be a set of distinct initial micrographs, then the texture model maps each micrograph <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M2">View MathML</a> into a high resolution micrograph xi such that the parameters (joint statistics of the subbands) of each micrograph, xi are same as those of the texture model. However the visible microstructural features of each xi will be different from each other and also from those observed in the low resolution micrograph at that location. We employ a Bayesian framework using a particle filter, such that a high resolution micrograph x is selected from all possible micrographs xi that is also consistent with the local microstructure. Based on maximum-entropy arguments, prior knowledge for the Bayesian framework takes the form of a Gaussian distribution with an interpolated low resolution micrograph as a mean and an assumed variance. The Likelihood is obtained from the texture model f, that maps low resolution image x0 to high resolution image x. To complete the specification of the Bayesian problem, observations consisting of a low resolution micrograph available at location C are used to update the prior.

A parametric texture model based on complex wavelet transform is used to represent the micrograph. High resolution micrograph is selected randomly from all possible micrographs satisfying a set of constraint functions, which will be discussed in section “Julesz conjecture and constraint functions”. The synthesis is performed by recursively projecting the initial guess of the micrograph onto the constraint surfaces using gradient projection methods [24]. Density-based Monte Carlo filter [32] along with the gradient projection algorithms are used to synthesize the high resolution micrograph at location C that is consistent with experimentally obtained low resolution micrograph at the same location. The problem is formulated as an estimation of a high resolution micrograph at location C using high resolution calibrating images available at locations A and B and low resolution image available at location C. The parametrization of the texture model is discussed in section “Parameter estimation sketch”. Texture parameters at location C are estimated from the parameters of images at locations A and B and those of coarse scale image at location C. The high resolution image at location C is simulated by recursively adjusting the pixel values of the iterates to satisfy a set of constraint functions and is updated using density-based Monte Carlo filter with coarse scale image at C as an observation.

Methods

Parametric texture model

A texture is modeled as a real two dimensional random field, X(i,j) defined on a finite lattice, <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M3">View MathML</a>. Here X(i,j) takes gray level values in <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M4">View MathML</a>. In this section we describe the steerable pyramid decomposition, which is a particular decomposition based on the complex wavelet transform of an image. A steerable pyramid acts as a filter and decomposes the image into multi-level, multi-oriented subbands. Joint statistics of the subbands, described in section “Parameter estimation sketch” are used as the parameters of the texture model. Section “Micrograph synthesis procedure” describes the texture synthesis procedure.

Julesz conjecture and constraint functions

The Julesz conjecture ([22,24,33]) states that there exists a set of functions <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M5">View MathML</a>, called constraint functions [24] such that samples drawn from any two random fields are visually indistinguishable under some conditions, if expectations over this set of functions are equal. That is,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M6">View MathML</a>

(1)

where, {ϕk(·)} is the set of constraint functions acting on a random filed and E[ ·] is mathematical expectation. To compute the expectations, practical ergodicity defined in [24] is assumed. Specifically, a homogeneous random field X has the property of practical ergodicity with respect to function <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M7">View MathML</a>, with tolerance ε, and probability p, if and only if the spatial average of ϕ over a realization drawn from X is a good approximation to the expectation of ϕ with sufficiently high probability. More concretely,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M8">View MathML</a>

(2)

In the present case, where realizations refer to sample images X(i,j), the spatial average of the functional ϕ can be approximated as,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M9">View MathML</a>

(3)

where, ⌊·⌋N means that the result is taken modulo N. With these definitions, our objective then is to draw samples from random field X satisfying statistical constraints of the form,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M10">View MathML</a>

(4)

The values ck in these constraints are computed from calibrating image X(i,j) and constitute the parameters of the texture model. The set of all samples satisfying equation 4 is referred to as the Julesz ensemble [22] given by

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M11">View MathML</a>

(5)

Steerable pyramid decomposition

We use a parametric texture model based on over-complete complex wavelet transform [24]. The complex wavelet transform decomposes the calibrating image into multi-scale, multi-oriented subbands. The joint statistics of the decomposed subbands are used as the parameters of the texture model. The decomposition is known as steerable pyramid when the basis functions are directional derivatives. The filters used in this decomposition are polar-separable in the Fourier domain. The traditional orthogonal separable wavelet decomposition suffers from aliasing whereas the steerable pyramid decomposition overcomes this limitation since the support of lowpass filter obeys the Nyquist sampling criterion [24]. A disadvantage with the steerable pyramid decomposition is that the basis is over-complete. A system is said to be over-complete if the the number of basis functions used to represent a vector is larger than required. Translation and rotation invariance of the steerable pyramid is important to represent the oblique orientations properly. Filters used in the steerable pyramid belong to one of three categories, namely low-pass, high-pass or oriented band-pass filters which are characterized by their respective response functions,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M12">View MathML</a>

(6)

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M13">View MathML</a>

and

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M14">View MathML</a>

(8)

where the angular part of Bk is given by,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M15">View MathML</a>

(9)

Here, r and θ are polar coordinates and <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M16">View MathML</a>. L(r) and H(r) are low-pass and high-pass filters respectively and {Bk(r,θ)} are oriented band-pass filters. The steerable pyramid decomposition is initiated by splitting an image into high-pass and low-pass subbands. The low-pass band is further decomposed into a set of orientation bands and the low-pass residual band. The low-pass residual band is sub-sampled by two and the decomposition in the next level is carried out by further decomposing the sub-sampled low-pass residual band into orientation bands and low-pass residual band. The procedure is continued till the image is decomposed into required number of pyramid levels. A system diagram for steerable pyramid is shown in Figure 4.

thumbnailFigure 4. Steerable pyramid diagram. Block diagram for the steerable pyramid [24,25]. Input image is initially split into high-pass and low-pass bands. Low-pass band is further decomposed into a set of oriented subbands and low-pass residual band. The recursive decomposition is achieved by inserting the diagram in dotted lines into the solid circle.

In our work four orientation bands and four pyramids with a total of 18 subbands are used (16 orientation bands, one high-pass and one low pass residual band). Figure 5 shows the real parts of four level and four orientation steerable pyramid subbands and low-pass residual band. High-pass subband is not shown in Figure 5.

thumbnailFigure 5. Real parts of subbands. Real parts of steerable pyramid orientation bands (4 orientations and 4 scales) and low-pass subband of image at location A. High-pass subband is not shown here.

Parameter estimation sketch

Outline for estimating parameters of the texture model from a sample image is described in the following steps:

1. Decompose the high resolution sample image x into high-pass subband, h0, orientation bands, {bs,k,s=1,⋯,4,k=1,⋯,4} and low-pass residual bands, {l0,l1,l2,l3,l4}.

2. Compute a set of joint statistics on the decomposed subbands. Following are the statistics [24] that are considered as parameters of the texture model:

●  Marginal statistics are the statistics of gray-level values of the pixels in the image.

Mean, μ1(x)=E{x(i,j)}.

Variance, μ2(x)=E{(x(i,j)−μ1(x))2}

Skewness,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M17">View MathML</a>

Kurtosis,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M18">View MathML</a>

Range =[ min(x)), max(x)]

●  where, μn(x)=E{(x(i,j)−μ1(x))n},∀n>1. In the parametric texture model, mean μ1(x), variance μ2(x), skewness η(x), kurtosis κ(x), minimum and maximum of the sample image, variance of the high-pass band μ2(h0), skewness {η(l0),η(l1),η(l2),η(l3),η(l4)} and kurtosis, {κ(l0),κ(l1),κ(l2),κ(l3),κ(l4)} of partially reconstructed low-pass images at each scale, {l0,l1,l2,l3,l4} are considered as a set of parameters from marginal statistics.

●  Raw coefficient correlations, are defined as the central samples of auto-correlation of the partially reconstructed low-pass images, Ak(n,m)=E{lk(i,j)lk(⌊i+nN,⌊j+mM)}. In the texture model, auto correlations of the partially reconstructed low-pass images, {l0,l1,l2,l3,l4} at each scale, are used as a set of parameters. These statistics characterize the salient spatial frequencies and the regularity of the texture [24].

●  Coefficient magnitude statistics are defined as the central samples of auto-correlation of magnitude of each subband, cross correlation of magnitude of each subband with those of the other orientations at the same scale, cross correlation of subband magnitudes with those of all orientations at a coarser scale. Cross-correlation of any two subbands bs(i,j) and bt(i,j) is given as Cs,t(n,m)=E{bs(i,j)bt(⌊i+nN,⌊j+mM)}. These represent the structures such as edges, bars and corners in the image.

●  Cross-scale phase statistics are cross-correlation of the real part of the coefficients with both the real and imaginary part of the phase doubled coefficients at all orientations at the next coarser scale. The rate at which the phase changes for fine scale coefficients is twice the rate change for coarse scale coefficients [24]. To compensate for the difference in rate of phase change, the cross correlations are computed between the phase doubled coarse scale coefficients and fine scale coefficients. These statistics distinguish edges from lines.

All the joint statistics described in the above list are used as the parameters of the texture model. Practical ergodicity described in section “Julesz conjecture and constraint functions” is used to compute expectations E[ ·] in the above expressions.

3. Once a set of parameters is computed from high resolution micrograph, a micrograph with that set of parameters can be constructed using the synthesis procedure described in the next section.

Micrograph synthesis procedure

The synthesis procedure proposed in [24] to generate an image with desired parameters is described in the following steps. The parameters of the texture model are computed from calibrating images as joint statistics of the decomposed subbands.

1. The synthesis procedure is initiated by choosing as initial guess an image with mean and variance equal to those of the target image.

2. The initial image is decomposed into multi-scale, multi-oriented subbands and joint statistics of the desired image are imposed recursively using gradient projection method [24].

3. Starting from coarse scale to fine scale, statistical constrains are imposed on low-pass and orientation bands while simultaneously reconstructing the low-pass image.

4. Auto-correlation, skewness and kurtosis of the reconstructed low-pass image at the fine scale are adjusted and resulting low-pass image is added to variance adjusted high-pass image to obtain the fine scale image.

5. The marginal statistics are adjusted on the pixels of the reconstructed image in step-4.

6. The procedure is repeated from step-2 to step-5 until the image with desired statistics is generated.

7. This algorithm synthesizes an image that is statistically equivalent to the sample images from which the parameters are estimated.

Figure 6 shows one iteration step of synthesis algorithm. This algorithm can be used to simulate many statistically equivalent micrographs for a given set of parameters. But the simulated micrographs will not be visually similar to the calibrating micrograph from which the parameters are computed. When a low resolution microscopic image from an experiment is available, data assimilation methods can be used to synthesize a high resolution micrograph that is consistent with the low resolution image. We use the Matlab tools provided at [34] for steerable pyramid decomposition and micrograph synthesis. In the next section a particle filter known as the density-based Monte Carlo filter [32] along with micrograph synthesis procedure is used to simulate high resolution micrographs that are consistent with the corresponding low resolution one.

thumbnailFigure 6. Block diagram of the recursive micrograph synthesis procedure.

Density-based Monte Carlo filter

The density-based Monte Carlo filter is a nonlinear, non-Gaussian filter proposed in [32] for state estimation of a dynamical system. In the current problem, an image is synthesized by an extension of the procedure described in section “Micrograph synthesis procedure”. In the synthesis procedure, an initial guess of the image is assumed and the desired joint statistics are imposed on the initial guess to obtain a new image. Let f be an implicit algorithm that implements steps 2-5 in the synthesis procedure shown in section “Micrograph synthesis procedure”. The parameters of the texture model are repeatedly imposed on the initial image x0 such that the final image has the desired parameters. That is the implicit function f acts on x0 to generate x1, on x1 to generate x2 and in general acts on xk−1 to generate xk. Action of f on xk−1 to generate xk, is shown in Figure 6. We treat this as one function evaluation in the synthesis algorithm. Let x0,x1,⋯,xk denote the images, modeled as random fields, generated at each function evaluation of the synthesis algorithm. These are construed as successive states of a dynamical system where an implicit algorithm takes a dynamical state (in the present case, the state is a random field) xk−1 to xk where xk depends only on xk−1 and the parameters of the texture model. Here, function evaluation steps, k=0,⋯,K are treated analogous to time steps in the dynamical system. Coarse scale images available at synthesis locations are treated as measurements and are used in the measurement equation of the particle filter algorithm to update the estimated image at each synthesis step. The state and measurement equations of the dynamical system can thus be written as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M19">View MathML</a>

(10)

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M20">View MathML</a>

(11)

where, xk is the synthesized image at kth step in the recursive synthesis process, f is an implicit algorithm relating xk to xk−1 as shown in Figure 6 and g is a coarsening operator relating state xk and observation yk. Here, coarsening operator is an averaging operator that averages 8×8 pixels in fine scale image to produce one pixel in the coarse scale image. Errors in the synthesis process and experimental measurements are modeled as additive white Gaussian noise processes ηk and εk respectively. The covariance of the noises are assumed to be known and are indeed parameters of the filter. Let Xk={x0,⋯,xk} be a vector of states and Yk={y0,⋯,yk} be a vector of measurements. In our work, measurements are assumed to consist of a coarse scale image available at location C and does not change with each synthesis step. Hence, we assume, Yk={y0,y1=y0,⋯,yk=y0}. The joint probability density of all the states and the measurements is given by

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M21">View MathML</a>

(12)

where, P(Yk|Xk) and P(Xk) can be written as,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M22">View MathML</a>

(13)

and

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M23">View MathML</a>

(14)

The conditional densities, P(ys|xs) and P(xs|xs−1), are obtained from the measurement and state equations respectively. From equations 13 and 14, the filter density P(xk|Yk) can be written as,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M24">View MathML</a>

(15)

The estimated micrograph, <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M25">View MathML</a> at step k, can be obtained by computing expectation of xk with respect to the filtering density, P(xk|Yk), i.e.,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M26">View MathML</a>

(16)

Evaluating the integrals in equation 16 using Monte Carlo sampling results in the filter estimate <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M27">View MathML</a> at step k. Generating n random samples of x0 and ηs for s={1,⋯,k}, a set of random draws Xi,s can be obtained from the state equation as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M28">View MathML</a>

(17)

The filtering estimate <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M29">View MathML</a> in equation (16) is approximated by drawing samples of xk using Monte Carlo sampling according to,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M30">View MathML</a>

(18)

Using equation 12,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M31">View MathML</a>

where, (ys|xi,s)∼N((gs(xi,s)),cov(εs)) denotes the likelihood and <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M32">View MathML</a> is the filter estimate. Here, Xi,k is a collection of random draws, defined as Xi,k={xi,0,⋯,xi,k}. Equation (18) can be written as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M33">View MathML</a>

(20)

where the weight function wi,k is defined as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M34">View MathML</a>

(21)

such that <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M35">View MathML</a>. For the initial step, when k=0, equal weights are assigned for all the samples <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M36">View MathML</a>, i.e., <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M37">View MathML</a>. Writing <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M38">View MathML</a>, the variance of ek can be written as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M39">View MathML</a>

(22)

and the sample variance can be written as

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M40">View MathML</a>

(23)

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M41">View MathML</a>

(24)

The following limiting conditions as n goes to infinity can be verified [32]

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M42">View MathML</a>

almost surely,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M43">View MathML</a>

in distribution,

<a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M44">View MathML</a>

almost surely.

Following the above constructions, the computational procedure of the density-based Monte Carlo filter is as follows:

1. The random draws of the initial state x0, represented as xi,0 are taken from the initial density P(x0),

2. Given xi,k−1 and ηi,k−1, estimates of Xi,k are obtained from the state equation 10,

3. Given the initial weight wi,0, weight functions wi,k for i=1,⋯,n and k=1,⋯,K are obtained from equation 21,

4. Finally, the filtered state <a onClick="popup('http://www.immijournal.com/content/2/1/2/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.immijournal.com/content/2/1/2/mathml/M45">View MathML</a> is evaluated from equation 20 while the variance of the error in the estimate is computed from equation 24.

Figure 7 shows the block diagram of the micrograph synthesis using parametric texture model and density-based Monte Carlo filter.

thumbnailFigure 7. Block diagram of synthesis procedure with particle filter. Block diagram of the recursive micrograph synthesis procedure along with the density-based Monte Carlo filter.

Results and discussion

The high resolution micrograph synthesis procedure is demonstrated by synthesizing a high resolution micrograph at location C starting with low resolution micrograph at that location and high resolution calibrating micrographs available at locations A and B as shown in Figures 1 and 2. The procedure can be repeated to construct high resolution micrographs at any other location. Two cases are considered to construct high resolution micrographs. In the first case (Figures 8, 9 and 10), a parametric texture model with the synthesis procedure described in section “Micrograph synthesis procedure” is used to construct a high resolution micrograph, whereas in the second case (Figures 11, 12 and 13), the parametric texture model along with the density-based Monte Carlo filter described in section “Density-based Monte Carlo filter” is used to estimate a high resolution micrograph.

thumbnailFigure 8. Micrograph at A synthesized without particle filter. Micrographs from left to right show the experimental high resolution micrograph and the synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location A. (a) Experimental high resolution micrograph at A; (b) Synthesized micrograph at A with a random initial guess a1; (c) Synthesized micrograph at A with a random initial guess a2; (d) Synthesized micrograph at A with a random initial guess a3.

thumbnailFigure 9. Micrograph at B synthesized without particle filter. Micrographs from left to right show the experimental high resolution micrograph and the synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location B. (a) Experimental high resolution micrograph at B; (b) Synthesized micrograph at B with a random initial guess b1; (c) Synthesized micrograph at B with a random initial guess b2; (d) Synthesized micrograph at B with a random initial guess b3.

thumbnailFigure 10. Micrograph at C synthesized without particle filter. Micrographs from left to right show the experimental high resolution micrograph and the synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location C. High resolution micrograph at location C is used only for comparison and is not used in the synthesis process. (a) Experimental high resolution micrograph at C; (b) Synthesized micrograph at C with a random initial guess c1; (c) Synthesized micrograph at C with a random initial guess c2; (d) Synthesized micrograph at C with a random initial guess c3.

thumbnailFigure 11. Micrograph at A synthesized with particle filter. Micrographs from left to right show respectively, experimental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location A. (a) Experimental high resolution micrograph at A; (b) Experimental low resolution micrograph at A; (c) Synthesized high resolution micrograph at A using parametric texture model and particle filter.

thumbnailFigure 12. Micrograph at B synthesized with particle filter. Micrographs from left to right show respectively, experimental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location B. (a) Experimental high resolution micrograph at B; (b) Experimental low resolution micrograph at B; (c) Synthesized high resolution micrograph at B using parametric texture model and particle filter.

thumbnailFigure 13. Micrograph at C synthesized with particle filter. Micrographs from left to right show respectively, experimental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location C. High resolution micrograph at location C is used only for comparison and is not used in the synthesis process. (a) Experimental high resolution micrograph at C; (b) Experimental low resolution micrograph at C; (c) Synthesized high resolution micrograph at C using parametric texture model and particle filter.

When the particle filter is not used, the synthesis procedure does not make use of the low resolution micrographs available at a given location. It uses only the parameters of the texture model obtained from the calibrating micrographs. The synthesis process with distinct initial guesses results in statistically equivalent high resolution micrographs, that have visually distinct microstructural features. To validate the synthesis procedure, high resolution micrographs are synthesized at locations A and B where experimental high resolution micrographs are available. Figure 8(a) shows a high resolution micrograph obtained experimentally at location A and Figures 8(b)–8(d) show synthesized micrographs whose joint statistics are the same as those of 8(a) but with visibly distinct microstructural features. Similarly Figure 9(a) shows an experimental high resolution micrograph at location B and Figures 9(b)–9(d) show corresponding synthesized micrographs. Parameters of the texture model for the synthesis of these micrographs at locations A and B are obtained using high resolution micrographs available at these locations. In case of location C, on the other hand, we assume that an experimental high resolution micrograph is not available and we aim to synthesize one. The parameters of the texture model for C are obtained from a low resolution micrograph at C and high resolution micrographs available at locations A and B. Pixel statistics and low-pass band statistics which account for total gray scale levels are computed from low resolution micrograph at C and the rest of the parameters which determine high resolution features are computed as averages of the parameters from A and B. Figure 10(a) shows experimental high resolution micrograph and Figures 10(b)–10(d) show corresponding synthesized micrographs at location C. Here high resolution micrograph at location C (Figure 10(a)) is used only for comparison and is not used in the synthesis process.

In the second case, high resolution micrographs of the microstructure are estimated in the Bayesian context using the density-based Monte Carlo filter (section “Density-based Monte Carlo filter”) where the texture synthesis process (section “Micrograph synthesis procedure”) is treated as a dynamical system and the experimental low resolution micrograph at a given location is treated as measurement used to update the dynamical state. It is reiterated that the high resolution micrograph at C is synthesized using the low resolution micrograph at C and high resolution micrographs at A and B together with a texture model and a particle filter. The parameters of the texture model for A, B and C are obtained similarly to before. For estimating high resolution micrographs at locations A, B and C, experimental low resolution micrographs at corresponding locations are used as measurements in the filtering algorithm. Low resolution micrographs are also used as the statistical average of the initial probability density function in the algorithm. The algorithm convergence is slowed for other choices of initial average. Measurement and process noises are assumed to be zero mean Gaussian random fields with 0.01 coefficient of variation. In the particle filter algorithm, 25 Monte Carlo samples are used to estimate the high resolution micrographs. Figure 11(a) and Figure 11(b) show high and low resolution micrographs obtained experimentally at location A, Figure 11(c) is the estimated high resolution micrograph at location A using density-based Monte Carlo filter along with the texture synthesis procedure. Figures 12(a)–12(c) show the experimental high resolution, experimental low resolution, estimated high resolution micrographs at location B. At location C, Figure 13(a) and Figure 13(b) are experimental high and low resolution micrographs, Figure 13(c) is the estimated high resolution micrograph. It should be noted that the high resolution features of micrograph at location C (Figure 13(c)), are recovered using information from experimental low resolution micrograph at C and experimental high resolution micrographs at A and B.

It is again worth comparing Figure 10 showing the synthesis of high resolution micrograph using only the parametric texture model with Figure 13 showing high resolution micrograph synthesized using texture model along with the particle filter. It can be clearly observed that the estimated high resolution micrograph with particle filter (Figure 13(c)) provides information about the local microstructure that cannot be gleaned from previous approaches. Resolving this level of detail is paramount in physical problems where bifurcation is initiated by spatially local events.

Conclusions

An original high resolution micrograph synthesis procedure using parametric texture model and particle filter produced very good quality, high resolution micrographs. When the particle filter is not used, the synthesized micrograph has the same joint statistics as those of sample high resolution micrograph, but visible microstructural features are not consistent with those of local microstructure at that location. Using a particle filter, the estimated micrograph acquires the same high resolution features as those of local microstructure at that location. Experimentally obtained micrographs are often corrupted by measurement noise. The synthesis procedure described in this paper, which includes particle filtering, accounts for measurement noise. To synthesize high resolution micrograph at location C, a set of parameters affecting the coarse scale features is obtained from low resolution micrograph at location C and rest of the parameters, influencing the fine scale features are obtained from high resolution calibrating micrographs available at location A and B. If all the parameters are obtained from, coarse scale image at C, then the synthesized image will not have smooth features observed in high resolution image. If the parameters are obtained only from calibrating images A and B, then, the total gray level content and its spatial distribution will be in consistent with that of the experimental micrograph at that location.

Appendix

Here, we provide the pseudo-code for high resolution micrograph synthesis using the parametric texture model and density based Monte-Carlo filter. A few high resolution micrographs are decomposed using steerable pyramid method and the parameters of the texture model are obtained as joint statistics of the decomposed subbands described in “Parameter estimation sketch”.

Algorithm 1 Recursive micrograph syntheis procedure

Algorithm 2 Density based Monte-Carlo filter

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All four authors, RT, RG, SG and DP collaborated thoroughly in the performance of the reported research. The authors integrated their respective professional experience and knowledge of high resolution micrograph synthesis during this work. All authors read and approved the final manuscript.

Acknowledgements

This research was supported by NSF CMMI Award number 0728304.

References

  1. Ghosh S, Valiveti D, Harris SJ, Boileau J (2006) A domain partitioning based pre-processor for multi-scale modelling of cast aluminium alloys. Model Simul Mater Sci Eng 14:1363-1396 Publisher Full Text OpenURL

  2. Valiveti DM, Ghosh S (2007) Morphology based domain partitioning of multi-phase materials: A preprocessor for multi-scale modelling. Int J Numer Meth Engng 69:1717-1754 Publisher Full Text OpenURL

  3. Lehmann TM, Gönner C, KSpitzer (1999) Survey: Interpolation methods in medical image processing. IEEE Trans Med Imaging 18:1049-1075 PubMed Abstract | Publisher Full Text OpenURL

  4. Singh H, Gokhale AM, Mao Y, Spowart JE (2006) Computer simulations of realistic microstructures of discontinuously reinforced aluminum alloy (dra) composites. Acta Materialia 54:2131-2143 Publisher Full Text OpenURL

  5. Groeber M, Ghosh S, Uchic MD, Dimiduk DM (2008a) A framework for automated analysis and simulation of 3d polycrystalline microstructures.: Part 1: Statistical characterization. Acta Materialia 56:1257-1273 Publisher Full Text OpenURL

  6. Groeber M, Ghosh S, Uchic MD, Dimiduk DM (2008b) A framework for automated analysis and simulation of 3d polycrystalline microstructures. part 2: Synthetic structure generation. Acta Materialia 56:1274-1287 Publisher Full Text OpenURL

  7. Rollett AD, Lee SB, Campman R, Rohrer G (2007) Three-dimensional characterization of microstructure by electron back-scatter diffraction. Annu Rev Mater Res 37:627-658 Publisher Full Text OpenURL

  8. Tewari A, Gokhale AM, Spowart JE, Miracle DB (2004) Quantitative characterization of spatial clustering in three-dimensional microstructures using two-point correlation functions. Acta Materialia 52:307-319 Publisher Full Text OpenURL

  9. Geman S, Geman D (1984) Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Pat Anal Mach Intell 6:721-741 OpenURL

  10. Geman S, Graffigne C (1986) Markov random field image models and their applications to computer vision. In: Proceedings of International Congress of mathematicians. Berkeley, California, USA.

    pp 1496–1517

    OpenURL

  11. Berthod M, Kato Z, Yu S, Zerubia J (1996) Bayesian image classification using markov random fields. Image Vis Comput 14:285-295 Publisher Full Text OpenURL

  12. Fan G, Xia XG (2003) Wavelet-based texture analysis and synthesis using hidden markov models. IEEE Trans Circuits Syst I: Fundam Theory Appl 50:106-120 Publisher Full Text OpenURL

  13. Choi H, Baraniuk RG (2001) Multiscale image segmentation using wavelet-domain hidden markov models. IEEE Trans Image Process 10:1309-1321 PubMed Abstract | Publisher Full Text OpenURL

  14. Kato Z, Pong TC (2000) A markov random field image segmentation model for color textured images. Image Vis Comput 24:1103-1114 OpenURL

  15. Dass SC, Nair VN (2003) Edge detection, spatial smoothing, and image reconstruction with partially observed multivariate data. J Am Stat Assoc 98:77-89 Publisher Full Text OpenURL

  16. Niezgoda SR, Yabansu YC, Kalidindi SR (2011) Understanding and visualizing microstructure and microstructure variance as a stochastic process. Acta Materialia 59:6387-6400 Publisher Full Text OpenURL

  17. Fleuret F, Geman D (2001) Coarse-to-fine face detection. Int J Comput Vision 41:85-107 Publisher Full Text OpenURL

  18. Pan R, Reeves SJ (2006) Efficient huber-markov edge-preserving image restoration. IEEE Trans Image Process 15:3728-3735 PubMed Abstract OpenURL

  19. Zhu SC, Mumford D (1997) Prior learning and gibbs reaction-diffusion. IEEE Trans Pattern Anal Mach Intell 19:1236-1250 Publisher Full Text OpenURL

  20. Schultz R, Stevenson RL (1994) A bayesian approach to image expansion for improved definition. IEEE Trans Image Process 3:233-242 PubMed Abstract | Publisher Full Text OpenURL

  21. Zhang X, Lam KM, Shen L (2005) Image magnification based on adaptive mrf model parameter estimation. In: Proceedings of 2005 International Symposium on Intelligent Signal Proceeding and Communication Systems. Hong Kong.

    pp 653–656

    OpenURL

  22. Zhu S, Wu YN, Mumford D (1998) Filters, random fields and maximum entropy (FRAME):towards the unified theory for texture modeling. Int J Comput Vis 27(2):107-126 Publisher Full Text OpenURL

  23. Wu YN, Zhu SC, Liu X (2000) Equivalence of julesz ensembles and FRAME models. Int J Comput Vis 38:247-265 Publisher Full Text OpenURL

  24. Portilla J, Simoncelli E (2000) A parametric texture model based on joint statistics of complex wavelet coefficients. Int J Comput Vis 40(1):49-70 Publisher Full Text OpenURL

  25. Simoncelli EP Freeman W (1995) The steerable pyramid: A flexible architecture for multi-scale derivative computation. In: In Second Int’l Conf. on Image Proc.,. Washington, DC: IEEE Sig. Proc. Society.

    vol III, pp 444–447

    PubMed Abstract OpenURL

  26. Simoncelli EP, Freeman W, Adelson E, Heeger D (2000) Shiftable multi-scale transforms. IEEE Trans Inf Theory 38(2):587-607 OpenURL

  27. Tipireddy R, Nasrellah HA, Manohar CS (2009) A kalman filter based strategy for linear structural system identification based on multiple static and dynamic test data. Probabilistic Eng Mech 24:60-74 Publisher Full Text OpenURL

  28. Doucet A, Godsill S, Andrieu C (2000) On sequential monte carlo sampling methods for bayesian filtering. Stat Comput 10:197-208 Publisher Full Text OpenURL

  29. Ghosh SJ, Manohar CS, Roy D (2008) A sequential importance sampling filter with a new proposal distribution for state and parameter estimation of nonlinear dynamical systems. Proc R Soc A 464:25-47 Publisher Full Text OpenURL

  30. Gordon NJ, Salmond DJ, Smith AFM (1993) Novel approach to nonlinear/non-gaussian bayesian state estimation. IEEE Trans Radar Signal Proc 140:107-113 Publisher Full Text OpenURL

  31. Manohar CS, Roy D (2006) Monte carlo filters for identification of nonlinear structural dynamical systems. Sadhana 31(4):399-427 Publisher Full Text OpenURL

  32. Tanizaki H (1996) Nonlinear filters: estimation and applications,. Berlin, Germany: Springer. OpenURL

  33. Julesz B (1962) Visual pattern discrimination. IRE Trans Info Theory IT-8. 84-92 OpenURL

  34. Simoncelli EP (2009)

    matlabpyrtools code. http://www.cns.nyu.edu/~eero/steerpyr/ webcite