Adaptive Neutron Radiography Correlation
for Simultaneous Imaging of Moisture Transport
and Deformation in Hygroscopic Materials
S.J. Sanabria & C. Lanvermann & F. Michel & D. Mannes &
P. Niemz
Received: 13 March 2014 /Accepted: 18 September 2014 /Published online: 21 October 2014
#
Society for Experimental Mechanics 2014
Abstract Neutron radiography is a key non-destructive test-
ing technology for the investigation of moisture transport in
materials. However, quantitative moisture measurements in
swelling materials are currently challenging due to the lack
of referencing between moist and dry state radiographs. A
novel adaptive texture correlation algorithm is presented to
simultaneously image inhomogeneous moisture distributions
and moisture-induced strain fields. The proposed method
provides a valuable tool for the study of time- and position-
dependent hygromechanical interactions. Moreover, it re-
quires no modification of existing neutron installations. The
method was validated against gravimetric moisture content
and optic surface deformation measurements. Its applicability
was demonstrated for two actual topics in wood science, the
investigation of moisture gradients within the growth ring
microstructure and the study of moisture transport processes
in wood-fiber composites. The algorithm can be widely used
to characterize hygroscopic materials with heterogeneous
texture, as frequently found in wood constructions, food in-
dustry, engineering and soil science.
Keywords Non-destructive neutron imaging
.
Heterogeneous
texture correlation
.
Moisture-induced deformation strains
.
Wood composite materials
.
Mechano-sorptive rheology
Introduction
Neutronimaging(NI)isanestablishedmethodtonon-
destructively study the dynamics of hydrogenous fluids, such
as bound (hygroscopic) or free water molecules, in geomate-
rials (rocks, soils), engineered media (ceramics, polymers,
metals, concrete, fuel cells, heat pipes, engines…) and biolog-
ical composites (wood, food, natural fiber-reinforced compos-
ites..), among others. Compared to X-ray and Nuclear
Magnetic Resonance (NMR) imaging, NI is highly sensitivity
to hydrogen while being transparent to most solid-state mate-
rials, including heavy metals. Moreover, it allows dynamic
full-field imaging in two (radiography) and three
(tomography) dimensions with tenth-of-millimeter spatial res-
olution and sub-second measurement times [1–5].
Moisture intake in porous media is often accompanied by
swelling, which is especially relevant for biological compos-
ites [6]. Wood is a classical example, volumetric swelling
strains >10 % are not rare within its large hygroscopic range,
which goes up to 30 % moisture content from oven-dry state
until the fiber saturation point (FSP), due to the hydrophilic
behavior of the lignocellulosic cell walls. If not accounted for,
swelling strain gradients can impair dimensional stability and
induce eigenstresses, leading to cracking and delamination,
which reduce the serviceability of wooden constructions [7].
Natural fiber-reinforced composites have drawn widespread
attention in the last decades due to their low cost and biode-
gradability, accompanied by high strength and insulation
Electronic supplementary material The online version of this article
(doi:10.1007/s11340-014-9955-2) contains supplementary material,
which is available to authorized users.
S.J. Sanabria (*)
:
C. Lanvermann
:
F. Michel
:
P. Niemz
Institute for Building Materials ETH Zurich, Stefano-Franscini-Platz
6, CH-8093 Zurich, Switzerland
e-mail: ssanabria@ethz.ch
C. Lanvermann
e-mail: lanvermannchr@ethz.ch
F. Michel
e-mail: frmichel@ethz.ch
P. Niemz
e-mail: niemzp@ethz.ch
D. Mannes
Neutron Imaging and Activation Group, Paul Scherrer Institute,
CH-5232 PSI, Villigen, Switzerland
e-mail: david.mannes@psi.ch
Experimental Mechanics (2015) 55:403–415
DOI 10.1007/s11340-014-9955-2
properties [8]. However, their highly hygromorphic behavior
is challenging. For example, the inner stresses induced in the
hot pressing of wood-based fiber composites lead to larger
thickness swelling than in raw wood material [9]. The cou-
pling between moisture transport and deformation is also a
key element to describe food drying, clay consolidation and
sorption in polymer electrolyte fuel cells. Although sophisti-
cated hygromechanical models are available, they often lack
of accompanying experimental data. In this context, the si-
multaneous full-f ield im aging of moist ure and moist ure-
induced deformation is necessary [6, 10–14].
Quantitative moisture measurements with neutrons also
require hygro-expansion data. The gravimetric moisture con-
tent ω=m
h
/m
o
is defined in terms of the mass m
h
of water with
respect to the mass of the dry solid phase m
o
, which requires
the local referencing of the neutron radiographs in moist T
state to the dry state T
o
. Moreover, the density of the solid
phase is lower in moist than in dry stat e due to sample
swelling. Thus the in-plane swelling strain fields ε
ii
with
respect to the reference dry state (Lagrangian descrip-
tion) need to be included in the calculation of ω
(Equaton 1)[11, 15]:
ω ¼ 1 þ ε
XX
ðÞ1 þ ε
YY
ðÞρ
h
z
h
ρ
o
l
z
ðÞ
−1
ð1Þ
with l
z
the sample thickness, ρ
h
and ρ
o
the densities of water
(10
3
kg m
-3
) and dry solid phase substance, respectively, and
z
h
the effective water column thickness (Equaton 2):
Z
h
¼ − Σ=ρðÞ
h
ρ
h
−1
lnT−lnT
o
1 þ ε
XX
ðÞ
−1
1 þ ε
YY
ðÞ
−1
hi
ρ
o
¼ −lnT
o
Σ=ρðÞ
o
l
z
−1
ð2Þ
(Σ/ρ)
h
and (Σ/ρ)
o
are respectively the mass-attenuation coef-
ficients in water and solid phases. For the energy spectrum of
the thermal neutron beamline NEUTRA at the Paul Scherrer
Institute (Villigen, Switzerland) [41], where the presented
investigations were carried out, (Σ/ρ)
o
=0.18 m
−2
kg
−1
for
wood (Table 1,[42]) and (Σ/ρ)
h
=0.35 m
−2
kg
−1
for water
[43]. Operating Equation 1 shows that neglecting hygroscopic
strains of 10 % in wood leads to unacceptable quantitative
moisture errors over 5 % (absolute error), as shown in
Equation 3:
ω ¼ ω
jε¼0
1 þ ε
XX
þ ε
YY
½
þ Σ=ρðÞ
o
Σ=ρðÞ
h
−1
ε
XX
þ ε
YY
þ ε
XX
ε
YY
½≈ω
jε¼0
þ 0:52 ε
XX
þ ε
YY
½ ð3Þ
Finally, the edge misalignment induced by hygro-
expansion leads to unbounded errors and non-physical mois-
ture values at image discontinuities, which impedes or makes
the study of moisture transport through multi-layered com-
posites and surface coatings challenging [16, 17].
These shortcomings have so far been alleviated with a
global pre-registering of neutron radiographs based on sample
edge detection [18–21], which reduces edge misalignment,
but does not account for inner deformation gradients. Full-
field strain measurements with neutrons have been achieved
in some crystalline materials by using a tunable monochro-
matic beam, for which Bragg cut-off wavelengths provide
lattice spacing [1]. Yet a combination with moisture measure-
ments has not been reported. NI was recently combined with
optical deformation measurements [15]. Although accurate
high-resolution ε
ii
measurements were obtained, the setup
requires sample speckling and an involved synchronized cam-
era installation. Moreover, only surface deformation is mea-
sured, which leads to significant edge misalignment artifacts
in volumetric samples.
In this work, a novel Adaptive Neutron Radiography
Correlation (ANRC) algorithm is presented, which uses the
texture information in neutron images to locally estimate the
strain fields. This m ethod is exclusively based on the
postprocessing of radiographs, therefore requiring no modifi-
cation of existing neutron installations, and providing average
deformation estimates over the sample thickness. The new
approach is validated with gravimetric and optic methods. Its
potential and applicability are demonstrated for two actual
research topics in wood science:
a) Gradients in moisture and hygroscopic swelling in soft-
wood growth rings [15, 16, 22–25]
b) Investigation of moisture transport through swelling
wood-fiber composites [3, 17, 18, 26]
Adaptive Neutron Radiography Correlation (ANRC)
The starting point for the ANRC algorithm is a set of neutron
radiographies of the test samples, which are experimentally
acquired at well-controlled dry and moist states. State of the
art corrections are first applied to compensate for source and
detector inhomogeneities (CCD dark current, median filter,
intensity normalization, flat field correction, spectral effects),
background scattering (black body calibration) and sample
scattering (Monte-Carlo simulations) [43].
Figure 1 summarizes the building blocks of the ANRC
algorithm. The inner core performs Digital Image
Correlation (DIC) processing. Subsets of the reference image
I
ref
(dry state) are automatically searched at specific positions r
of the deformed image I
test
(moist state), resulting in local
estima tes of the deformation vector u(r), from which the
strains fields ε
ij
are differentiated. The DIC method is well-
known and performs best for artificial random speckle
404 Exp Mech (2015) 55:403–415
patterns, where unique subset correspondence is satisfied.
Typical DIC parameters are subset size, search function
and region and deformation order [27]. Many heteroge-
neous materials show enough natural texture to track u
without need of artific ial sample s peckling, the DIC
process is then specifically termed Texture Correlation
(TC) [28]. Neutron radiographs of wood composites are
assumedtoshowenoughtextureforTC,whichis
empirically confirmed by the high correlation statistics
observed in the investigated test samples (Table 1).
The proposed ANRC processing is thus a TC algorithm,
which is optimized and expanded to robustly extract defor-
mation information from NI radiographs while ignoring
“shadow regions”, where only poor correlation statistics are
available. With this purpose, the DIC core implements a zero-
order search (rigid subsets). A zero-normalized cross-correla-
tion function I
corr
(Equation 4) is calculated at each integer
pixel position within a search window W for a subset size Δ.
The peak maximum r provides the displacement vector u,
which is refined by applying bicubic interpolation to I
corr
:
I
corr
b
x;
b
y
¼
X
x;y∈Δ
Tx−
b
x; y−
b
y
−T
hi
T
o
x; yðÞ−T
o
hi
X
x;y∈Δ
Tx−
b
x; y−
b
y
−T
hi
2
X
x;y∈Δ
T
0
x; yðÞ−T
0
hi
2
0
B
B
B
@
1
C
C
C
A
−1
u ¼ arg max I
corr
b
x;
b
y∈W
ð4Þ
In order to minimize false correlation hits, the search region
is adaptively adjusted. The algorithm starts searching in a
large reset window W
R
. Subsequent search positions r
i
are
ordered in a raster scan of the reference image. For each of
them, the search is reduced to a shorter window W
S
cen-
tered at an initial deformation estimate û, which is
calculated based on the linear ex trapolation of N
û
neighboring pixels. An error contro l block ensures the
well-posing of the estimated u by controlling its: a)
correlation: the co rrelation coefficient r is above a
minimum threshold r
min
,b)uniqueness: u is within
image bou nds, is not an e dge pixel of I
corr
, an d does
not exceed a maximum deformation value u
max
,c)
continuity: u does not increase by more than Δu
max
in each iteration. In the c ase of an error of type b), a
second trial is admitted with a larger window W
L
before discarding u. Only non-error u values (e=0)
are used for the adaptive search and the continui ty
error check. Ill-posed u (e≠ 0) increment the error
count c
e
up to a threshold N
e
. Then the reset control
block starts operation, deactivating the adaptive search
and resetting the search window to W
R
. The algorithm
manages errors along individual dimensions of t he de-
formation vector u
j
by only storing well-posed u
j
.The
loss of continuity upon reset is additionally avoided by
managing an additional reset count c
R
and not
accepting u values at first reset positions (c
R
=1) for
updating the deformation estimate û, by using at reset
positions û instead of u for c ontinuity control, and by
filtering out isolated u pixels.
Table 1 Summary of neutron attenuation properties, ANRC settings and correlation statistics for the investigated test materials
ρ
o
(kg m
-3
)(Σ/ρ)
o
(m
-2
kg) d (m) Δ (mm)
b
s
ρ
/ρ (%)
c
CNR (a.u.)
c
r (a.u.) e=1 (%) u
max
(mm)
I: Optical, softwood growth
rings
a
(Fig. 2,4.1.1)
n.a. (not applicable) n.a. n.a. [0.2, 0.2] n.a. 4.6 0.96 [0.7, 0.7] [1.5, 1.0]
II: Neutron, softwood growth
rings
a
(Fig. 3–6,4.1.2)
362 0.18 0.005 [2.9, 2.9] 34 5.1 0.98 [0.5, 16] [0.8, 1.7]
III: Neutron, wood-fiber
composites
a
(Fig. 7–8,4.2)
Isonat chanvre 44 0.27 0.04 [4.5, 9] 15 2.5 0.92 [0, 0] [0.2, 0.04]
Pavaflex 48 0.23 0.04 [4.5, 9] 16 2.5 0.91 [0, 0] [0.2, 0.08]
Pavatherm 135 0.19 0.04 [4.5, 9] 6 1.2 0.76 [0, 0] [0.1, 0.04]
Aerogel Vliesmatte 147 0.10 0.04 [4.5, 9] 9 1.6 0.88 [0, 0] [0.3, 0.2]
Isoroof-natur-KN 225 0.18 0.04 [4.5, 9] 6 1.1 0.64 [1.6, 1.4] [0.6, 0.2]
a
The pixel sizes for I, II and III are 23, 77 and 145 μm, respectively
b
Fixed algorithm settings (px): W
S
=[3, 3], W
L
=[5, 5], Δu
max
=[3, 3], N
e
=5, N
û
=[10, 100]
c
s
ρ
is the standard deviation of the dry material density ρ
o
within the correlated subsets. CNR=s
Δ
/s
n
,withs
Δ
and s
n
respectively for correlation subset
and background image noise. s
Δ
=s
ρ
(Σ/ρ)
o
l
Z
exp[−ρ
o
(Σ/ρ)
o
l
Z
] for the neutron images T, following from Equation 2.s
n
=0.02 for all radiographs
Exp Mech (2015) 55:403–415 405
Continuous deformation fields are calculated by fitting the
estimated deformation vectors u to cubic smoothing splines
f(r
i
) (Equation 5), which trade-off between data fidelity and
smoothness [29]:
X
i¼1
N
r
i
ju
i
j
− f r
i
j
2
þ 1−pðÞ
Z
D
2
jftðÞj
2
dt p ¼ 1 þ h
3
=6
−1
ð5Þ
The weights of the error function are adjusted to the peak
correlation coefficients r. The resolution h, which con-
trols the smoothness, is adjusted with the subset size in
the range h=Δ/3…Δ. The fitted spline functions are
finally differentiated to compute the strain fields ε
ij
.In
order to minimize noise sensitivity, only first order
derivatives are included in ε
ij
=0.5(∂u
i
/∂x
j
+∂u
j
/∂x
i
).
The settings of ANRC for each case study are summarized
in Table 1. The subset size Δ is a trade-off between well-
posing and lateral resolution of u. A threshold Δ is typically
found, above which correlation statistics do not significantly
improve. The ANRC algorithm was implemented in Matlab®
(R2011b, The Mathworks Inc., Natick, MA, USA) and is
provided in the supplementary materials for free use and
development (Supp. Mat.).
Experiments
Moisture and Swelling Gradients in Softwood Growth Rings
at Hygroscopic Equilibrium
Optical surface deformation measurements
The strain field ε
ij
calculation with ANRC was first quantita-
tively validated against a commercial DIC tool (VIC 2D 2009,
Correlation Solutions Inc., USA) for high-resolution optical
images of speckled softwood (Fig. 2). The DIC measurements
were previously used to investigate wood cross-grain hygro-
expansion at the growth ring structural level [15, 24]. The goal
was to identify swelling strain gradients between the alternating
bands of low-density earlywood (EW)– thin-wall wood cells
with large internal lumens grown in spring – and high-density
latewood (LW) – thick-wall cells with small-sized lumens grown
in summer (Fig. 2a). Four cuboid samples (40×40×5 mm
3
)
were cut fr om a Norwa y spruce (Picea abies Karst.) stem, their
edges were well-aligned with respe ct to the main wood material
axes (radial R×tangential T×longitudinal L). The samples were
sprayed with a black/white speckle pattern by means of an
airbrush gun with a nozzle size of 0.2 mm. They were then
subjected in a portable climatic chamber to a moisture sorption–
desorption cycle with relative humidity (RH) adjusted in 15 %
steps between 0 and 95 %. For each nominal RH state, the
samples were conditioned to equilibrium moisture content and
optical images were acquired with a CCD camera (Fig. 2b). The
Fig. 1 Flow diagram of ANRC algorithm
406 Exp Mech (2015) 55:403–415
oven-dry state (RH=0 %) was used as reference. Next, the in-
plane swelling strains were calculated with ANRC and com-
pared to the commercial DIC results (Fig. 2c). The subset size Δ
was 9x9 pixels
2
for both methods, for a pixel size of 23 μm.
Finally, the strains were transformed with respect to the material
axes (R, T) and plotted in function of the relative growth ring
position φ=0…1, with 0 and 1 respectively defining EW and
LW regions (Fig. 2d) . The coordinate transformation method is
further detailed in [15, 24].
Simultaneous neutron imaging of moisture content
and moisture-induced deformation
The moisture and strain gradients within the growth rings
were investigated at equilibrium conditions with NI at the
thermal neutron beam line NEUTRA at the Paul Scherrer
Institute in Villigen, Switzerland. The experimental setup
consists of a neutron beam collimator, a neutron detector and
a portable climate chamber (Fig. 3a) and has been detailed in
[15]. The neutron radiographs (Fig. 3b) show a pixel size of
145 μm and a lateral resolution of ~0.7 mm (minimum de-
tectable line spacing). The deformation vectors u and strain
fields ε
ij
at each RH step were computed with ANRC and a
subset size Δ of 21×21 pixels
2
(Fig. 4). The moisture content
distribution ω was calculated with Equation 1 and plotted in
function of the growth ring position (Fig. 5).
The average ω values were validated gravimetrically (Fig. 6).
At each RH state the weight of the sample was recorded with
0.1 mg accuracy with a precision scale. The local moisture and
strain fields were validated with the combined NI and optical
surface deformation measurement setup of [15]. With this pur-
pose, a stereo-vision installation (Fig. 3a) consisting of two
radiation-shielded CCD cameras and cold-light illumination
captured optical images from the sample laterals without
perturbing the neutron line of sight. The surface deformation
vectors u were then calculated with commercial DIC (VIC3D,
Correlated Solutions). As in 3.1.1, a subset size Δ of 9×9 pixels
2
was used, the resolution was here lower due to the optical
constraints (pixel size 77 μm). The optical images (Fig. 3c)
and the neutron radiographs were referenced with fiducial
markers. The calculated deformation vectors were then used to
aligntheneutronradiographsandtocomputeω.
Transient moisture diffusion and moisture-induced swelling
in wood-fiber composites
The applicability of the ANRC algorithm to simultaneously
investigate dynamic moisture transport and the large coupled
deformations in biological composites was demonstrated for a
selection of multi-layer building insulation materials. Three
commercial wood fiberboards (Pavaflex, Pavatherm and
Isoroof-natur-KN from Pavatex SA, Fribourg, Switzerland),
a hemp fiber composite (Isonat chanvre from Valnaturel SA,
Saxon, Switzerland) and a nanoporous aerogel (Aerogel
Vliesmatte from Aspen Aerogels Inc., Northborough, MA,
USA) were tested (Table 1). The material properties are de-
tailed in [26, 30]. The test objects (40×40 mm
2
section) were
typical three-layer insulation combinations, with two 20 mm
thick Isoroof outer layers and a 35 mm (Isonat, Pavatherm,
Pavaflex) or a 5 mm thick (Aerogel) middle layer. Following
Fig. 2 Quantitative validation of
ANRC processing on high-reso-
lution CCD camera images. a)
Test Norway spruce samples. The
growth ring microstructure is il-
lustrated at cellular scale with LM
and at subcellular scale with
ESEM. b) Optical surface images
at dry (RH=0 %) and moist
(RH=95 %) states, correlation
subsets I
ref
, I
test
andsearchfunc-
tions I
corr
are highlighted. c)
shows hygroscopic strain fields
calculated from b) with both
ANRC and a commercial DIC
tool, which in d) are averaged and
plotted in function of the growth
ring position φ
Exp Mech (2015) 55:403–415 407
