Image Anal Stereol 2014;33:131-145 doi:10.5566/ias.v33.p131-145 Original Research Paper INTERRUPTED IN-SITU COMPRESSIVE DEFORMATION EXPERIMENTS ON MMC FOAMS IN AN XCT: EXPERIMENTS AND ESTIMATION OF DISPLACEMENT FIELDS KATHARINA LOSCHC,1, KATJA SCHLADITZ1, UTA BALLASCHK2, HARRY BEREK2 AND CHRISTOS G. ANEZIRIS2 1Fraunhofer-Institut f¨ur Techno-und Wirtschaftsmathematik; 2TU Bergakademie Freiberg, Institute of Ceramics, Glass and Construction Materials e-mail: katharina.losch@itwm.fraunhofer.de; katja.schladitz@itwm.fraunhofer.de; uta.ballaschk@ikgb.tu-freiberg.de; harry.berek@ikgb.tu-freiberg.de; aneziris@ikgb.tu-freiberg.de (Received December 13, 2013; revised April 22, 2014; accepted May 26, 2014) ABSTRACT The mechanical properties of a metal-matrix composite foam are investigated by interrupted in-situ compressive deformation experiments within an X-ray computed tomography device (XCT). Each in-situ experiment generates a sequence of reconstructed 3D images of the foam microstructure. From these data, the deformation .eld is estimated by registring the images corresponding to three consecutive steps. To this end, the generic registration framework of the itk software suite is exploited and combined with several image pre­processing steps. Both segmented (binary) images having just two grey values for foreground (strut structure) and background (pore space) and the result of the Euclidean distance transform (EDT) on pore space and solid phase are used. The estimation quality is evaluated based on a sequence of synthetic data sets, where the foam’s microstructure is modelled by a random Laguerre tessellation. For large deformations, a combination of non-rigid registration for the EDT images and partwise-rigid registration on strongly deformed regions of the binary images, yields surprisingly small estimation errors. Keywords: 3D image data, displacement .eld, image registration, volume images. INTRODUCTION The deformation and collapse behaviour of cellular structures is of growing interest for a lot of applications as light weight crash absorbing systems. In the present study, mechanical properties of metal-matrix composite (MMC) foams combining an austenitic TRIP-steel and magnesia partially stabilised zirconia (Mg-PSZ) are investigated. Both components exhibit a martensitic phase transformation – a diffusionless change of the atomic arrangement – generating the potential for improved mechanical properties like strength, strain, and energy absorption. Interrupted in-situ compressive deformation experiments are performed within an XCT device. 3D imaging allows to observe the spatial microstructural deformation as well as to select regions of interest for subsequent metallographic preparation and phase analysis (Berek et al., 2011). Quantitative description of the shear bands and calculation of displacement .elds contribute essentially to comprehension and prediction of the damage progress. Each in-situ experiment generates a sequence of reconstructed 3D images of the foam microstructure. The problem of estimating the deformation .eld from the image data is equivalent to registration of the images corresponding to consecutive experimental steps. Registration – .nding the transformation mapping one image to another such that the features match best – is a well-known and widely studied problem in medical imaging. Thus we use the generic registration framework of “The Insight Segmentation and Registration Toolkit” (itk) software suite (Yoo, 2004; Ibanez et al., 2005; ITK, 2013) in order to estimate the deformation .eld. Previous to the actual registration, the images are denoised, binarised, and a Euclidean distance transform is applied to both pore space and solid component to enhance structural information in spite of the high porosity. Subsequently, the best .tting transform is found iteratively, measuring image similarity using a mutual information metric. The non­rigid transform is modeled by continuous, piecewise polynomial functions (B-splines). In order to be able to estimate the structural failure, a “re.nement step” is applied on areas of poor .t. In order to validate our method, the estimation error is determined for a synthetic data set being very close to the real observations. A realisation of a Laguerre tessellation generated by a random sphere packing serves as model foam (Redenbach et al., 2011). It is dilated, digitised, and virtually deformed. For the thus derived image sequence, the estimated deformation .elds are compared to the true ones. This paper is organised as follows. First, the experimental setup for the in-situ experiments as well as the resulting image data are described. The used image processing and segmentation methods are shortly summarised in the following section. The subsequent section contains an overview of image registration and its main algorithmic ingredients, a description of the itk registration framework, and details of our algorithm. Then the validation using random Laguerre tessellations is explained followed by the section presenting the results. Finally, we discuss the results and give an outlook to future work. EXPERIMENTAL SETUP AND IMAGE DATA Composites with improved mechanical properties in comparison to the single materials can be manufactured by the combination of ductile metals with ceramics. Guo et al. (2003) created metal matrix composites (MMC) based on a low alloyed TRIP-steel (transformation induced plasticity). The incorporation of up to 20 vol. % of yttrium partially stabilised zirconia (Y-PSZ) yielded strengths of about 1 400 up to 2 100 MPa under dynamic impact loading accompanied by an increase in the global hardness of approximately 25%. Indeed, the composite did not show plasticity and failed in a brittle behaviour. We consider a high alloyed stainless Cr-Mn-Ni TRIP-steel reinforced by magnesia partially stabilised zirconia (Mg-PSZ)-particles. The austenitic face centred cubic crystal structure of the steel can be retained to room temperature depending on the chemical composition. Thus, stress and deformation induced transformations to the tetragonal space centred martensitic phase become possible at room temperature. This formation of martensite during plastic deformation of austenite results in a signi.cant increase of ductility (Transformation Induced Plasticity) and strength as shown by Kovalev et al. (2012). Martensitic phase transformations also appear in zirconia. Between melting point (2 680.C) and 2 370.C the cubic crystal structure is stable. Further cooling causes the transformation into the tetragonal modi.cation. Below 1 170.C highly distorted monoclinic crystal structure of zirconia is thermodynamically stable (Kisi and Howard, 1998; Stevens, 1986). The high temperature modi.cations of zirconia can be stabilised to room temperature by the addition of positive ions like Mg2+ and Y3+ . In connection with oxygen vacancies metastable tetragonal zirconia is formed. Under suf.cient mechanical loading partially stabilised zirconia shows stress induced transformation from the tetragonal to the monoclinic modi.cation. This is connected with a volume increase leading to compressive stresses in the surrounding of the zirconia particles. This results in transformation toughening (Garvie et al., 1975; Green et al., 2000). Thus, in the corresponding MMC the effects of martensitic phase transformations within matrix and reinforcement with the particle reinforcement effect lead to improved mechanical performance (Aneziris et al., 2010). The resulting enhanced energy absorption capability is strongly desired for crash absorbers. Open cell foams are manufactured by powder metallurgy using a replica method (Schwartzwalder, 1963). Thereby, polyurethane (PU) foams with 30 pores per inch are applied as templates. TRIP-steel (TLS Bitterfeld, Germany) and 10 or 5 vol.% magnesia stabilised zirconia powders (Saint Gobain, USA) are mixed with organic additives, respectively. The particle size distribution and the chemical composition of the powders are speci.ed in Tables 1-3. The manufacturing technology and the composition of the organic additives are described in detail by Aneziris et al. (2010) and Sieber et al. (2010), respectively. Table 1. Particle size distribution of the powders. More precisely, d10 ,d50, and d90 denote the 10, 50, and 90% quantiles of the particle diameter distribution, respectively. d10 [µm] d50 [µm] d90 [µm] X5CrNi18-10 11.1 34.8 91.25 Mg-PSZ 0.1 1.3 10.8 Table 2. Chemical composition of the X5CrNi18-10 TRIP-steel powder. C Si Mn Nb S Cr wt.% 0.05 0.12 0.96 0.01 0.01 16.4 Ni N Mo Ti Fe wt.% 9.3 0.07 1.8 0.004 bal. Table 3. Chemical composition of the Mg-PSZ powder. MgO SiO2 Al2O3 CaO wt.% 3.4 2.4 0.6 0.2 TiO2 Fe2O3 NaO2 ZrO2 wt.% 0.1 0.1 0.1 bal. Image Anal Stereol 2014;33:131-145 The polyurethane foam is impregnated with the prepared slurry, squeezed out and dried at 90.C for 1 hour. Subsequently, cold spray coatings based on the same mixture composition are applied in order to eliminate defects and to reach the critical wall thickness for the desired mechanical properties. The structure of the MMC foam is mainly determined by the polyurethane foam. After pyrolysis process, the burned out PU foam leaves cavities within the walls of the metal-ceramic foam. It is presumed that these voids have an essential in.uence on the mechanical and thermo-mechanical behaviour of the foam. After removing the organics (debindering) the samples are sintered for 2 hours at 1 350.C in argon atmosphere. For the deformation experiments samples with a size of about 10 x 10 x 16 mm3 are cut. Interrupted in-situ compressive deformation experiments are performed in a computer tomograph (ProCon X-Ray Garbsen, Germany) with a 160 kV transmission X-ray tube. The nominal spatial resolution (pixel size) was about 30 µm after reconstruction, the volume rendering of the MMC foam sample is shown in Fig. 1. Fig. 1. Reconstructed CT image of the MMC foam sample. Sample size 10.2 × 10.2 × 16.45mm3. Volume rendering using MAVI (2005). Two loading plates are positioned in a carbon .bre reinforced polymer tube between the top load cell and the bottom mechanical actuator. The sample is placed between the loading plates. The deformation speed was about 10-3/s at quasi-static and isothermal conditions. Between the deformation steps, the position was held constant in order to obtain CT images of the structures. After relaxation periods (e.g., recovery effects by movement of dislocations), CT scans were taken under constant stress. Fig. 2 shows the time schedule of an XCT in-situ deformation experiment for an MMC foam with 10 vol.% up to 20% average compressive strain corresponding to the initial loaded structure area. The particular stress maxima are points of the respective continuous stress strain curve. Three areas can be separated: Elastic deformation, cracking leading to a plateau and densi.cation involving a signi.cant increase of the stress. Fig. 2. In situ quasi-static compression test in XCT on a MMC foam with 10 vol.% Mg-PSZ. Stress and strain versus time. The combination of test device and computed tomography allows the observation of the damage process in selected areas. Due to the manufacturing process, cell wall thicknesses vary within the foam sample. Vertices and struts are thicker than the cell walls in between. This effect is due to the surface tension of the applied slurries. Furthermore, caused by the spray process, cell walls in the external areas are thicker than in the central area. Berek et al. (2011) use the cell wall distribution for a .rst quantitative description of the structural change under loading. The deformation process during our experiment is characterised by a bottom-up shifting of the structure due to the movement of the lower plate of the compression device. Thanks to the thicker cell walls, the external regions are rather stable, while the central region is deformed. First cracks appear at the weakest cells, which are rather large, feature thin cell walls, and are stretched perpendicular to the deformation axis. Starting from the .rst cracks, deformation spreads to the neighbouring cells – a so-called deformation band forms. These deformation bands can be explained by the dramatic change of the stress distribution in the neighbouring cells after the .rst break of a cell wall. On the other hand, compact structure elements of the foam remain stable up to high degrees of deformation. (a) 2% (b) 6% (c) 10% (d) 13% IMAGE PRE-PROCESSING AND SEGMENTATION Before actually registering the images corresponding to consecutive experimental steps, the original grey value images have to be pre-processed as imaging artifacts like noise, ring artifacts, or global grey value .uctuations can distort the registration results considerably. Thus these sources of error are minimised by .rst .ltering the images using a 5×5×5 median .lter. Subsequently, the solid component is segmented. That is, a new, binary image is created with just two grey values for foreground – the MMC strut system –, and background – the pore space. This segmentation is achieved by combining two standard methods – Otsu’s (1979) global and Sauvola and Pietik¨ainen’s (2000) local grey value thresholdings. For dividing the two classes (foreground and background) according to the pixels’ grey values, Otsu’s method .nds the threshold value so that the intra-class variances are minimised, while the variance between the two classes is maximised. Sauvola’s method chooses the threshold for each pixel x as t(x)= µW (w,x)1 + ks(W (w,x))/R - 1, (e) 16% (f) 20% Fig. 3. Reconstructed CT images of an MMC foam sample with 10 vol.% Mg-PSZ at 2, 6, 10, 13, 16 and 20% deformation. The 2D slices were cut at the same position within the 3D images. A quantitative description of the shear bands and displacement .elds is supposed to improve comprehension and prediction of the damage initiation and progress signi.cantly. Estimation of deformation .elds based on 2D image correlation is not able to capture the close relation between failure mechanisms and features of the microstructure. Hence the following sections are dedicated to an algorithm allowing precise estimation of the 3D deformation .elds based on the obtained sequence of CT image data. with a cubic window W (w,x) of size w and centre x, the mean µ(W (w,x)) of the grey values inside the window, their standard deviation s(W (w,x)), and the dynamic range of the standard deviation R, see Sauvola and Pietik¨ainen (2000). Finally, masking the result of the local thresholding with the one of the global thresholding yields both the outline of the foam structure and the .ne pores within the struts (Ohser and Schladitz, 2009, Fig.4.11). Some grey value .uctuations stem from the locally varying thickness of the material and thus contain spatial structural information. This implicit spatial information is lost by binarising the images. In order to reconstruct it and even gain information about the pore space, we applied the Euclidean distance transform (EDT) on both foreground and background of the binarised images and smoothen the result by a mean .lter. This way, we assigned each pixel its signed shortest Euclidean distance to the other class while keeping a smooth grey value histogram. The resulting image therefore contains signi.cantly more spatial information than the simple black-and-white image: midpoints of struts and cells are bright, the surfaces of the struts are dark. The whole process is visualised in Image Anal Stereol 2014;33:131-145 Fig. 4. All pre-processing and segmentation steps were performed using MAVI (2005). (a) original (b) median .ltered (c) binarised, Otsu (d) binarised, Sauvola (e) (c) masked by (d) (f) smoothed EDT Fig. 4. Pre-processing steps illustrated using 2D slices through the 3D image of the 10% MMC sample. IMAGE REGISTRATION OVERVIEW Image registration can be loosely de.ned as the task to identify corresponding regions in different images that somehow “show the same”. For example, image registration is widely used in a medical context for aligning images obtained by different methods (e.g., computed tomography and ultra-sound) or at different times (e.g., showing the lung in a relaxed and an unrelaxed state). For two images A and B that have to be aligned, this is generally achieved by maximising a function m(A,B(T )) that indicates the well-alignedness of the transformed image B(T ). Here B(T ) denotes the image holding the grey values of B under the spatial transformation de.ned by T . More precisely, the images A and B(T ) are aligned if 1. T is the correct transformation between A and B, and 2. each point of the image domain of A is .rst transformed by T and subsequently assigned the grey value corresponding to this position in B. This results in B(T ), see Fig. 5. Fig. 5. Schematic overview of registration. In order to track sub-pixel displacements, one has to consider an interpolation I(B) of the image. In case of an af.ne transform or a rigid transform (displacement of the whole image by the same vector or rotation), T would be a matrix; however for more complex, non-rigid transformations, theoretically it could be any function. Allowing arbitrary functions for T would obviously result in misalignment since the general problem is ill-posed. As a remedy usually certain assumptions for T such as af.nity or continuity are made and generally only certain types of functions are allowed as possible transforms, depending on the previous knowledge about the image content. Still, image registration is generally prone to misalignment, especially in case of large deformations or poor image quality. Moreover, the limitation to certain functions also limits the transformations that can be observed, also resulting in bad .ts or slow convergence in cases an inapt class of functions was picked. We decided to use grid nodes on the images and interpolated the transformation between them by B-splines. In its objective, image registration is closely related to the method of digital image and volume correlation, which also seeks to align images in order to .nd the deformations between the images. Digital Image Correlation (DIC) and its three-dimensional version Digital Volume Correlation (DVC) seek to obtain the transformation between images by a different method than image registration: Instead of minimising an overall cost function as in image registration, small subvolumes of the original image are .tted into the deformed image. In two dimensions, this method can be applied to a neighbourhood of every single pixel (or every second pixel) and yields usually precise results. Similar to standard image registration, there also arise problems with discontinuities and large deformations, but there have been recent proposals how to elegantly solve these issues (Poissant and Barthelat, 2010; Tang et al., 2012). Still, in three dimensions, DVC is very costly due to the large number of pixels even in relatively small images: Usually, one would only select few subregions and interpolate the transformation between the subregions or rely on parallel computing methods (Gates et al., 2011). FINDING INITIAL PARAMETERS It is generally well-known that the (fast) convergence of image registration and image correlation methods highly depends on a good initial guess of the deformation parameters. Especially in the case of large deformations as we observed them in our data, this is crucial for acceptable results. Albeit precise, DVC would have been computationally too costly to apply to our data due to the high number of pixels we needed in order to resolve discontinuities. However, we used the general idea of .tting only subvolumes to increase precision for guessing initial parameters for a B-spline transform. Therefore, we used rigid registration of “slices” in order to .nd initial parameters for the more sophisticated B-spline transform. The slices were set perpendicular to the direction of loading, this allowed us to initially register parts of the image, which experience little deformation, and then use the parameters of these regions as initial parameters for neighbouring, more strongly deformed slices, see Fig. 6: In this example one would register region A at .rst, since this region experiences only little change compared to the deformed image. Therefore, even a simple algorithm yields good results if we assume the initial parameter to be 0. The .nal parameters obtained by this registration would then be used as initial parameter for region B, which experiences more deformation than region A. After the registration of region B, its .nal parameters would be used as initial parameters for region C, and so forth. By this we make sure that the initial parameter for region E, which experiences the biggest deformation, are close to the real deformation parameters, avoiding local minima in the metric. The resulting parameters are set as initial parameters for the B-spline transform. Despite 30% strain corresponding to deformations as large as 70 pixels in our test set, we still obtain good results using this method. (a) (b) Fig. 6. Schematic overview of the “slices” registration. The next part, the B-spline transformation, makes the assumption of underlying continuity in the deformation which can be considered as legitimate for large parts of the images and generally for small deformations. At last, with the “re.nement step”, we sought out possible discontinuities in the already non-rigidly registered images and applied another registration in their neighbourhoods. Details of the algorithm are explained in the following. THE ITK REGISTRATION FRAMEWORK Fig. 7. Overview of the registration framework implemented in itk, (source: (Yoo, 2004; Ibanez et al., 2005)). Fig. 7 gives an overview over the itk registration framework: Using an optimisation algorithm (“Optimizer”), we minimise a measure (“Metric”) of the images’ difference with respect to the transformation (“Transform”) parameters. Interpolation (“Interpolator”) guarantees registration on sub-pixel level. Image Anal Stereol 2014;33:131-145 DIGITAL VOLUME CORRELATION For .nding initial parameters of the non-rigid registration, we registered subvolumes of the binarised images. For this registration, we used the mean squares metric as a measure of difference. For each subvolume V of A, a transformation tV =[ux,uy,uz]V between image A and B is determined by 2 t^V = min . A(i) - I(B(i +~t)) . i.V In this case, I is the above-mentioned interpolation function needed for detection of sub-pixel displacements. In our case of binarised images, we used a nearest neighbour interpolation, meaning that in order to interpolate grey values of B at non-grid points, these positions get assigned the grey value of the nearest pixel. For minimisation we used a regular step gradient descent algorithm with 300 iterations, a maximal step length of 4.0 and a minimal step length of 0.001. NON-RIGID REGISTRATION For the non-rigid registration setup, we used the algorithm proposed by Mattes et al. (2003), which provided a framework for non-rigid multimodality image registration, using a mutual-information-based measure and modeling the non-rigid deformation on a B-spline grid on the EDT images. They used a multiresolution approach for avoiding local minima and increased resolution. However, as we do not expect our materials to be as soft as their chest tissue and as the foam experiences not only large deformations, but also material failure, we abstained from this multiresolution approach. Instead, we estimated initial parameters for the B-spline grids by rigid registration of slices as described above. The main parts of this algorithm will be shortly described below. We used the implementation from ITK (2013). Difference measure Due to the application of the EDT on the binarised images, as described above, we preserved some spatial information in the images that would have been lost otherwise. However, now we cannot assume anymore that corresponding physical points have the same grey value within the images, leading to the need for a multimodal measure of the images’ well-adjustedness. A correlation-based measure only takes into account linear relationships in the image modalities, which is not suf.cient for our EDT registration. Therefore, we decided to use the above-mentioned mutual information approach from Mattes et al. (2001): The mutual information I(A,B)= -. p(a)log p(a)+ . p(a|b)log p(a|b) aa,b of two images A and B is a measure of the information A contains about B and vice versa. Here, p(a) is the probability of the occurrence of the grey value a in A and p(a|b) is the probability of the occurrence of the grey value a in a pixel in A given that the grey value b appears on the corresponding pixel in B (Pluim et al., 2003). Mattes’s mutual information uses B-spline kernels to smoothen the image grey value histograms. It is de.ned as follows: For the grid X and two images A : X › IA, B : X › IB and a sample S . X, we set p^(iA,iB) I^(A,B)= - .. p^(iA,iB)ln , p^A(iA)p^B(iB) iA.IA iB.IB with A(x,y, z) - A' ß(0) p^(iA,iB)= . . iA - bA (x,y,z).S B(x,y,z) - B' ß (3) · iB - ,bB with A' resp. B' being A’s and B’s grey value minima and bA and bB the grey value ranges of each histogram bin. This is necessary to .t a speci.c number of histogram bins to the grey value histogram. . is a normalizing factor to ensure . p^(iA,iB)= 1 and ß(0) and ß(3) are zero-order and cubic spline kernels, respectively. Following from that we get p^B(iB)= . p^(iA,iB) , iA.IA and A(x,y,z) - A' ß(0) p^A(iA)= . . iA - .bA (x,y,z).S We used a sample size of 30% for the calculation of the mutual information. Interpolation In order to allow sub-pixel registration, we sampled pixel values at non-grid positions of the original image B by a linear interpolation. This interpolation method was chosen heuristically as it performed better than cubic or tricubic B-splines. B-spline deformation In order to obtain a smooth, non-rigid deformation .eld, we used a set of grid nodes which were transformed independently. For the estimation of the transformation between the grid nodes we used tricubic B-splines (Unser, 1999). Each grid node’s initial parameters tD were obtained by the “slice” registration described above. The resulting deformation is x - xD(r) T (tD,x)= .tD(h)ß (3) - h . h .xD The displacement .eld T^ that was found in this step is then used to calculate I(B(T^ )). The gradient magnitude of the displacement .eld was subsequently used to identify the deformation bands. Minimisation Following the proposal from Mattes et al. (2003), we used the limited memory Broyden, Fletcher, Goldfarb and Shannon minimisation for bounded problems (L-BFGS-B optimiser) without bounding constraints (Byrd et al., 1994), which is a quasi-Newton minimisation algorithm. This algorithm is very well-suited for problems with many parameters since it only requires relatively small memory for storing the update steps. We found a cost function convergence factor of 103, a projected gradient tolerance of 10-7 and 60 evaluations suf.cient for a good result. This algorithm is provided by (ITK, 2013), too. REFINEMENT STEP With the setup of non-rigid registration that was described above we were only able to detect continuous deformations due to the assumption of B-spline transformation. However, as the foam structure is expected to break eventually, one has to assume that there are also discontinuities present. These discontinuities would lead to misregistrations on the image data. This thought gave rise to a third step of the algorithm, which we call “re.nement step”. Roughly said, the re.nement follows the following idea: 1. Identify badly registered regions. 2. Divide them into much smaller subregions. 3. Register these subregions separately using the previous continuous approach. This method can capture discontinuities by movement of the subregions w.r.t. each other while keeping computational effort practicable as only a small fraction of the original volume is used. In order to .nd discontinuities of the deformation, we compared the resampled binary image B(T^ ) with A and calculated their squared difference. The resulting image was .ltered with a median .lter in order to consider only larger misregistrations. The resulting binary image was labeled and the labels’ midpoints were subsequently used as central points to split their neighbourhood into eight cubes of size 50 × 50 × 50 (Fig. 8). On each of the resulting cubes we then performed again the non­rigid registration on the corresponding EDT images as described above. By this method, we allowed discontinuities on labels’ midpoints, which correspond to misregistrations, while keeping the continuous deformations on parts of the image that were already well-registered. Initial values for the Bsplines were calculated from the deformation .eld of the non­rigid registration step. The parameters were set to 30 evaluations during the minimisation, using 7 B-spline grids. The cost function convergence factor was set to 105, the projected gradient tolerance to 10-9. Additionally, we used a sample size of 80% to calculate the mutual informations of the cubes. In order to keep the computational times low, we set the size of the median .lter’s mask usually to 5. Fig. 8. Re.nement for each label. VALIDATION USING A SYNTHETIC FOAM STRUCTURE For validating our method, we used a synthetic 230×230×230 pixel data set of a structure resembling the real sample. Lautensack et al. (2008) and Redenbach et al. (2011) showed that dilated random Laguerre tessellations generated by a sphere packing are well suited to model ceramic foams. Here, we used the Laguerre tessellation deduced from a force-biased packing of spheres with gamma distributed volumes. The strut and wall thicknesses were not .t to the real sample in the strict sense, however, the model still showed similar properties regarding the strut thickness and closed facets. See Fig. 9 for a visualisation of the synthetic foam structure. Image Anal Stereol 2014;33:131-145 Fig. 9. Volume rendering of the synthetic foam used for method validation. Now this synthetic structure is virtually deformed. Due to the stiffness of the observed foam, we chose the general sigmoid function as most appropriate displacement in direction of the z-axis for testing our setup, while the synthetic sample is not deformed in x-or y-direction: Tx = Ty = 0 , K - 0.2 Tz = -0.2 - . 1 + 0.5exp(-0.04(z - 130)) The deformation .elds and resulting deformations of the synthetic foam are displayed in Fig. 10. For simulating discontinuities similar to the deformation bands observed in the real data, we chose T x = z * 0.05 * K , T y = 0 , K - 0.2 T z = -0.2 - , 1 + 0.5exp(-0.04(z - 130)) for z < 130 and T x =(z - 230) * 0.05 * K , T y = 0 , K - 0.2 T z = -0.2 - , 1 + 0.5exp(-0.04(z - 130)) else. This results in a displacement .eld similar to the one we would expect for material failure in the foam with both parts moving slightly sidewards to different directions. Fig. 11 shows the resulting deformation .elds and deformed synthetic foam. (a) K = 12 (b) K = 12 (c) K = 70 (d) K = 70 Fig. 10. Continuous deformation .elds and resulting deformation of the synthetic foam structure visualised using ParaView (2002) and MAVI, respectively. (a) K = 12 (b) K = 12 (c) K = 70 (d) K = 70 Fig. 11. Discontinuous deformation .elds and resulting deformation of the synthetic foam structure visualised using ParaView and MAVI, respectively. We used the parameters K = 12,24,56,70, corresponding to the approximate displacements in pixels observed in the experimental data. RESULTS In case of continuous deformations simulated by sigmoid functions, we obtained very good results for each value of K by only using the B-spline deformation and omitting the re.nement step: When comparing the obtained displacement .elds with the original displacement .elds, all large errors occur close to the image borders, they can by eliminated by cropping the image by 5 pixels at each edge. The resulting average errors, mean squared errors (mse), peak signal to noise ratios (psnr) and median errors can be found in Table 4. Table 4. Quality indices of B-spline registration for continuous deformation .elds for different K after cropping by 5 pixels at each side. Displayed are the mean differences between original and estimated deformation .eld for each vector component in pixels. K 12 mse 0.0467 psnr 34.84 medx -0.0123 medz 0.0740 med|.|2 0.1593 24 0.0945 37.80 -0.0068 0.0717 0.2013 56 0.222 41.45 0.001 0.0873 0.2580 70 0.3906 40.94 0.0058 0.1256 0.3157 For discontinuous test .elds we only obtained slightly better results when additionally using the “re.nement step”. Figs. 12,13 show histograms of the deviations between original discontinuous displacement .elds and the displacement .elds obtained by the B-spline algorithm, Figs. 14,15 show the respective histograms when using the “re.nement step”. The mean differences and and mean squared errors (mse) for the B-spline registration can be found in Table 5, the respective values for the re.nement step can be found in Table 6. The quality indices indicate that the re.nement step actually does not improve the quality of the registration. However, with our experimental data we found differences which we consider relevant. This might indicate that our discontinuous test .eld was a too simple model of the deformation of the experimental foam. Both the histograms and the fact that the error’s median is absolutely smaller than its mean, indicate that it is advised to use methods of robust regression for .tting: It is clear from the histograms that the error’s distribution is heavy-tailed and not Gaussian. LOSCH K ET AL: Displacement .elds of foams Table 5. Quality indices of the B-spline registration for discontinuous deformation .elds for different K after cropping by 5 pixels at each side. Displayed are the mean differences between original and estimated deformation .eld for each vector component in pixels. K mse psnr medx medz med|.|2 12 0.2465 27.62 0.0056 0.0760 0.3139 24 0.9066 27.98 -0.0216 0.0898 0.4590 56 4.0243 28.87 -0.0043 0.0925 0.7312 70 6.28072 28.8719 -0.0285 0.123 0.8943 Table 6. Quality indices of registration with re.nement step for discontinuous deformation .elds for different K after cropping by 5 pixels at each side. Displayed are the mean differences between original and estimated deformation .eld for each vector component in pixels. K 12 mse 0.2538 psnr 27.49 medx 0.0071 medz 0.0770 med|.|2 0.3140 24 0.8994 28.02 -0.0198 0.0892 0.4569 56 4.1744 28.71 -0.0014 0.0889 0.7241 70 6.632 28.64 -0.0306 0.1076 0.8937 Image Anal Stereol 2014;33:131-145 Finally, we applied our algorithm to subimages of the image data obtained by the experiments described earlier, performing the same preprocessing steps as described above. We registered two consecutive images. Fig. 16 shows the resulting deformation .elds for the experiments. Due to differences in coordinate systems of the programs used for displaying the data, the foams seem to be compressed from above and not from below as described in the experimental setup. We used the gradient magnitude computation described by Sapiro and Ringach (1996) to mark areas with large local changes in the deformation. This algorithm, which was used in its itk implementation, sets as gradient magnitude the difference of the two largest eigenvalues of the matrix [gi j] with .T . T gi j = · .ui . uj Fig. 14. Histograms of the errors in vector compontents (.x, .y, .z) and the length of the errors (.tot ) for discontinuous test .elds at K = 12 when using the “re.nement step”. for the deformation T (u1,u2,u3). These areas correspond to the above-mentioned deformation bands. The gradient magnitudes shown in Fig. 16 were calculated for the smooth deformation .elds obtained after the B-spline step. DISCUSSION For purely continuous displacement .elds, we were able to obtain a peak signal to noise ration of 40.94 dB and a mean squared error of 0.3906 for strains as large as 30% (for K = 70) when cropping the images by 5 pixels. For the more realistic cases of 5%, 10% and 20% strain, we had mean squared error of 0.0467, 0.0945 and 0.222. Fig. 15. Histograms of the errors in vector compontents (.x, .y, .z) and the length of the errors (.tot ) for discontinuous test .elds at K = 70 when using the “re.nement step”. The failure behaviour of the foam clearly exhibits discontinuities when and where the structure actually fails. When only being based on B-splines, our registration algorithm always yields a smooth deformation .eld T^ and is thus not able to capture the locally discontinuous changes when the .rst struts break. We therefore added a DVC-like procedure to improve deformation .eld estimation in these regions. DVC is computationally far more expensive, thus we applied it only to regions where material failure and discontinuities in the displacement .eld have to be expected. To this end, points of material failure were located by comparing the binarised original image A with the binarised image B(T^ ), both strongly smoothed. Restriction to central regions of the images and labelling of the connected components resulted in a small set of misaligned pixels. Only the surrounding pixels of the labels’ centres were taken into consideration for the third part of our algorithm. Surprisingly, this re.nement of the algorithm did not improve the estimation results signi.cantly on the test .elds. To the contrary, in some cases, the pure B­splines-registration led to more precise results. In order to estimate the quality of the algorithm on the images obtained by the experiment we compared the second binary image (B(T^ ) in the description of the algorithm) to the .rst binary image (A). For this we found a reduction of the amount of errors larger than 2 pixels on the experimental data, see Fig. 17 for a comparison between simple B-spline registration and registration with re.nement step. For example, in case of the resampling of the binary image corresponding to the 6% to 8% compression, the proportion of error pixels (which is the mean grey value distribution in case of binary images) is 0.0304 after the re.nement step and 0.0324 for a B-spline registration. In general, we have to expect errors of 1 pixel size due to the nearest neighbour interpolation. However, after eroding by 2 × 2 × 2 cubes, the proportion of error pixels is 0.0010 after the re.nement step and 0.002 for the B-spline registration, meaning the “re.nement step” essentially bisects the man of larger errors. For other registered images corresponding to a compression step of 2% we observe the same relation between the proportions. This indicates that our discontinuous test .elds were chosen unrealistically since we cannot observe this improvement on the test deformations. Still, the re.nement step is computationally quite expensive and dependent on the image data. The seemingly good results on the experimental data could be rooted in the fact that we do not know the real deformation of the experimental data and that therefore our comparison gives unrealistically good results. More importantly, we also showed that our algorithm is capable of acquiring discontinuous displacements via the re.nement step routine, although it is not sure to which amount since the peak signal to noise ratios and mean squared errors do not improve when compared to the B-spline transformation alone. However, we observe an improvement on our experimental data. Future work would consist in creating benchmarks for realistic discontinuous non­rigid deformations that allow testing of algorithms on realistic data sets, in parallelizing the algorithm for increased speed and precision and the implementation of a three-dimensional “subset splitting” procedure (Poissant and Barthelat, 2010) for better resolution of discontinuities and material failures. Image Anal Stereol 2014;33:131-145 Fig. 16. Resulting estimations of deformation .elds and the gradient magnitudes after B-spline registration, calculated on the experimental data. (a) (B(T^B-spline) - A)2 (b) (B(T^re fined ) - A)2 (c) (B(T^B-spline) - A)2 (d) (B(T^re fined ) - A)2 eroded by a 2 × 2 × 2 cube eroded by a 2 × 2 × 2 cube Fig. 17. Comparison of the misalignments in binary images for different algorithm steps on corresponding slices. White pixels correspond to errors. CONCLUSIONS Image registration yields useful results with respect to both precision and computational ef.ciency not only in medical applications but also for material sciences. By intelligent and cost-ef.cient initialisation of the image registration, we were able to successfully determine the deformation of open MMC foams under load. Analysis of the derived deformation .elds opens new opportunities to relate local failure with potentially distinctive structural features. For instance, according to the estimated deformation .eld, areas of high interest could be identi.ed and prepared for further (destructive) testing, e.g., to investigate the local crystal structure. Moreover, combining realistic models for both microstructure and deformation .elds would enable virtual experiments allowing to study the effects of the cellular structure, the strut and wall thickness distributions, local structural deviations, and the crystal structure separately. ACKNOWLEDGEMENTS The authors from TU Bergakademie Freiberg thank the German Research Foundation (DFG) for funding within the Collaborative Research Center TRIP-Matrix Composites (SFB 799). Mr. M. Hasterok is acknowledged for providing the specimen. REFERENCES Aneziris C, Berek H, Hasterok M, Biermann H, Wolf S, Kr¨ uger L (2010). Novel TRIP-steel/Mg-PSZ composite open cell foam structures for energy absorption. Adv Eng Mater 12:197–204. Berek H, Ballaschk U, Aneziris CG (2011). In situ characterization of internal damage in TRIP-steel/Mg-PSZ composites under compressive stress using X-ray computed tomography. Adv Eng Mater 13:1101–7. Byrd RH, Lu P, Nocedal J, Zhu C (1994). A limited memory algorithm for bound constrained optimization. Tech Rep NAM-08. Northwestern University, Dept Electr Eng Comput Sci. Garvie RC, Hannink RHJ, Pascoe RT (1975). Ceramic steel. Nature 258:703–4. Gates M, Lambros J, Heath M (2011). Towards high performance digital volume correlation. Exp Mech 51:491–507. Green DJ, Hannink RHJ, Swain MV (2000). Transformation toughening of ceramics. Boca Raton: CRC Press. Guo Y, Zhou Y, Li D, Duan X, Lei T (2003). Microstructure and performance of 2Y-PSZ/TRIP steel composites. J Mater Sci Technol 19:137–40. Iba´~nez L, Schroeder W, Ng L, Cates J (2005). The ITK Software Guide, 2nd Ed. Clifton Park, NY, USA: Kitware, Inc. http://www.itk.org/ItkSoftwareGuide.pdf ITK (2013). The insight segmentation and registration toolkit. http://www.itk.org/. Kisi EH, Howard CJ (1998). Crystal structures of zirconia phases and their inter-relation. Key Eng Mater 153-4:1– 36 Kovalev A, Jahn A, Weiß A, Wolf S, Scheller PR (2012). STT and DTT diagrams of austenitic Cr-Mn-Ni as­ ' cast steels and crucial thermodynamic aspects of . - . transformation. Steel Res Int 83:576–83. Lautensack C, Giertzsch M, Godehardt M, Schladitz K (2008). Modelling a ceramic foam using locally adaptable morphology. J Microsc 230:396–404. Mattes D, Haynor DR, Vesselle H, Lewellen TK, Eubank W (2001). Nonrigid multimodality image registration. Proc SPIE 4322:1609–20. Mattes D, Haynor DR, Vesselle H, Lewellen TK, Eubank W (2003). PET-CT image registration in the chest using free-form deformations. IEEE T Med Imaging 22:120– 8. MAVI (2005). Modular algorithms for volume images. http: //www.mavi-3d.de/. Ohser J, Schladitz K (2009). 3D Images of materials structures – processing and analysis. Weinheim: Wiley VCH. Otsu N (1979). A threshold selection method from gray level histograms. IEEE T Syst Man Cyb 9:62–6. ParaView (2002). Parallel visualization application. http: //www.paraview.org/. Image Anal Stereol 2014;33:131-145 Pluim JPW, Maintz JBA, Viergever MA (2003). Mutual information based registration of medical images: a survey. IEEE T Med Imaging . 22:986–1004. Poissant J, Barthelat F (2010). A novel ”subset splitting” procedure for digital image correlation on discontinuous displacement .elds. Exp Mech 50:353–64. Redenbach C, Wirjadi O, Rief S (2011). Modelling a ceramic foam for .ltration simulation. Adv Eng Mater 13:171–7. Sapiro G, Ringach DL (1996). Anisotropic diffusion of multivalued images with application to color .ltering. IEEE T Image Proc 5:1582–6. Sauvola J, Pietik¨ainen M (2000). Adaptive document image binarization. Pattern Recogn 33:225–36. Schwartzwalder K (1963). Method of making porous ceramic articles. US Patent 3090094. Sieber T, M¨uhlich U, Liedke T, Ballaschk U, Berek H, Aneziris CG, Ehinger D, Wolf S, Kr¨uger L (2010). Deformation and failure of open-cell foams made of trip-steel-zro2-composite materials: Experimental observations versus numerical simulations. Steel Res Int 82:1004–16. Stevens R (1986). Zirconia and Zirconia Ceramics. Magnesium Elektron Ltd. Tang Z, Liang J, Ciao Z, Guo C (2012). Large deformation measurement scheme for 3D digital image correlation method. Opt Laser Eng 50:122–30. Unser M (1999). Splines: A perfect .t for signal and image processing. IEEE Signal Proc Mag 16:22–38. Yoo TS (2004). Insight into images: principles and practice for segmentation, registration, and image analysis. Wellesley, MA, USA: AK Peters.