Focussed on “ Crack Paths ” A novel algorithm for crack path identification based on infrared images

The present work aims at developing a new approach towards the analysis of the crack path through the access to the information from the infrared (IR) video imaging of the steady fatigue crack growth. We first test some commonly available infrared image filtering techniques, and then we introduce a new robust algorithm to identify the position of the dynamic crack. The proposed procedure performs the local analysis of the temperature distribution around the fatigue crack tip in the 316L stainless steel. The developed MATLAB framework simplifies the routine processing and permits analyzing multiple experiments reasonably easily. Two filtering procedures Gaussian and wavelet shrinkage were implemented and applied successively to denoise the IR images. Both filters produce comparable and good quality results. The new algorithm capable of automatically locating the crack tip position from consecutive filtered IR images is proposed. Remarkably successful results are demonstrated with only small errors caused by materials distortion.


INTRODUCTION
n recent years, numerous energy-based approaches have gained considerable interest in attempts to characterize the fatigue crack growth.The consensus exists that failure is a dissipative process involving competing mechanisms of accumulation and dissipation of elastic energy at the crack tip.A large portion of the dissipated energy is irradiated as heat which is, nowadays, measurable by high-speed IR-cameras.Already in 1920, Griffith hypothesized that: "The "theorem of minimum potential energy" may be extended so as to be capable of predicting the breaking loads of elastic solids, if account is taken of the increase of surface energy which occurs during the formation of cracks".Later, in 1950, Freudenthal concluded that: "The study of mechanical behavior and properties at different levels of aggregation is possible only in terms of a concept which on all levels has the same meaning.This concept is energy."Multiple energy criteria for fatigue life estimation based on the hysteresis energy loss [1][2][3][4] or stored energy [4,5] have been proposed in the literature.A general assumption that plastic-strain hysteresis is a measure of fatigue damage is neither able to provide a satisfactory fit for experimental data [1] nor able to account for different stress ratios R [4].As an alternative, the energy-based approaches, e.g. the volume-based Strain Energy Density (SED) method [6], have been proved capable of accurate predicting fatigue life, thus illustrating the power of local energy-based concepts in application to the fatigue life characterization.Several experimental attempts have been made to measure the dissipated energy in situ [7,8].In this context, the work of Bannikov and Plekhov [5] is particularly worth noticing.Using thermal camera images of a growing crack and calculating the plastic work at the fatigue crack tip from the Hutchinson, Rice and Rosengren solution, a local parameter of stored energy was defined.It was clearly demonstrated that shortly before fracture, stored elastic energy reduced while the heat dissipation energy increased and reached the value of the plastic work as the material at fracture was no longer able to store the mechanical energy in nice agreement with the Orowan energy criterion of ductile fracture [9].At that time, the rate of the plastic work accumulation was lesser than the rate of the heat dissipation energy.This indicated that the stored energy nullified.It was assumed by the Authors that the damage accumulation mechanism changed as a harbinger of the imminent fracture where the role of macroscopic displacements was essential, and where the energy dissipation increased significantly.They proposed the value of the stored energy in the material to be a failure criterion based on the thermodynamic principles.Meneghetti and Ricotta [8] quantified the energy generated in the small area surrounding the crack tip and found an empiric correlation between the crack growth rate and the specific heat energy.The exponential relation between these two parameters remained even beyond the applicability of LEFM in the final crack growth stage.Exploiting the capacity of modern IR cameras, we performed a series of systematic fatigue crack growth tests with different loading ratios.The temperature distribution around the crack tip was continuously monitored by means of an infrared (IR) camera.The necessity of analysing a huge array of stored data and of improving image quality motivated us to develop a versatile data processing framework in the MATLAB environment.This work is focused on automatic pre-processingfiltering/smoothing -procedures required to reach the necessary image quality before the position of the crack tip can be I located along the crack path.Results of the tests revealing the energy associated with the crack growth under varying loading conditions (loading ratio and amplitude) are beyond the scope of the present paper centred around methodological aspects of crack path determination.The obtained IR images were processed aiming at obtaining the heat source function defined in the next session.For this purpose, the noise from the images must be removed.In fact, experimental data are contaminated by electronic noise, material and surfaces inhomogeneity, surface vibrations, etc.Two filtering procedures were tested.The first one applied an anisotropic three-dimensional Gaussian filter to reduce the time-dependent noise.The second approach was based on the wavelet shrinkage denoising strategy.A novel algorithm has been finally developed to follow the evolution of the crack path and to locate the crack tip on the basis of the temperature distribution in function of time T(x,y,t), as is illustrated in Fig. 1.The crack length is assumed to be constant during every recording.Finally, the crack tip position and the path are reconstructed after denoising.

HEAT SOURCE DEFINITION
btaining a 3D temperature map is challenging.The thermal vision cameras provide a 2D distribution map of the temperature at the surface of the tested specimen.Knowledge of the 2D temperature field assumes an obvious simplification, reducing an actual 3D situation to a two-dimensional diffusion problem and neglecting the temperature through the thickness gradients [10,11].The temperature distribution appears in response to the energy dissipation in the bulk due to a variety of interrelated irreversible processes involving, first of all, plastic deformation (mechanical work done in the plastic zone) and crack jumps (elastic energy release due to a new surface opening).Thus, although the temperature provides access to the thermodynamics of the dissipative processes of the crack growth, it is only indirectly related to the material's behaviour.To get a better understanding of the thermodynamics of underlying processes, it is necessary to recover the properties of the heat sources.This is commonly achieved by solving the heat equation for the heat source function   , , Q x y t [8,12].The general structure of this equation in Cartesian (x,y) coordinates reads as , , , , , , , , where is the material density, C is the specific heat capacity, and k is the heat transfer coefficient.
Thereafter, thermo-elastic coupling and dissipative phenomena can be evaluated [13] and correlated with the stress fields.
Since the solution of this equation requires both spatial and temporal derivation, the noise-induced fluctuations of the temperature function   , , T x y t are of crucial importance and have to be alleviated.Since Eqn.( 1) is a differential equation for   , , T x y t , obtaining the smooth ("denoised") form of   , , T x y t is of key importance for the overall success in the further evaluation of   , , Q x y t and for any next steps of the analysis including monitoring of the crack tip position.Two effective fileting methods are presented in the following sections.Both the heat source function over the whole specimen and its temporal evolution in the crack tip affected zone are presented in what follows.Since the aim is to obtain the heat source function characterizing the irradiated heat, comparisons are made between these differently evaluated quantities instead of the directly measured temperature.

EXPERIMENT SET-UP
atigue crack growth tests were performed on SEN(T) specimens shaped by waterjet from the cold-rolled 316L lowcarbon stainless steel.The specimen surface was coated with the high emissivity black coating.Cyclic load with a maximum of 6500N and fatigue ratio R=0.1 was applied using electrodynamic machine Instron Electropuls- 10000.Sine wave control signal was set at 10 Hz frequency.A larger set of load amplitude and fatigue ratio was tested, but the results are outside the scope of this work.The Telops infrared camera TS-IR MW (Noise Equivalent Temperature Difference -NETD is of 20 mK) with 25 mm lens was used after 1h warming up and stabilization.A spatial resolution was 15 m with the 400x512 pixels.The exposure time O F has been fixed to 2500 μs to increase the frame rate to 190 Hz.The Telops' CameraLINK acquisition card was used to perform continuous testing and data streaming from the camera to the personal computer (PC).Recorded data was continually transferred to the PC hard drive.Special care was exercised to minimize a possible thermal noise contribution.Three main sources of errors are known:  The electronic noise which can only be corrected including time in filtering  Inhomogeneities of the high emissivity paint which is easily solved with a spatial filter  Thermal reflections.Various precautions are suggested in the literature to overcome these problems and reduce these detrimental factors.The electronic noise is dependent on the camera parameters and settings.The thermal reflections were minimized by a slight rotation of the camera with respect of the surface normal.The experiments were conducted in a dark room with the ambient temperature kept constant throughout the testing [14].It is believed, admittedly with caution, that most can be considered adiabatic for frequencies above 2 Hz, but other factors such as the coating and system electronics might introduce concern below 10 Hz [15].In fact, the energy emitted per frame might be of the same order of magnitude as the noise level.

IMAGE FILTERING
he importance of noise removal is particularly acute for high-speed imaging when the exposure time is low.One possible method is to subtract the "baseline" reference image from the acquired sequence of frames to erase thermal reflections [16].This method experiences considerable difficulties in the final stage of the crack growth where the material's deformation is not negligible.The motion compensation technique ("lock-in"), of course, helps to improve the results [17].Admittedly the most common method of the noise smoothing is the median filtering and its variations, which are used to remove the spiky noise before applying a Gaussian filter.This method cannot be regarded as noise-specific as it does not consider any particular properties of the noise.To account for the local features of the noise and improve the filtering results, Vandone tried to use the wavelet decomposition [18].This method is quite flexible as it allows for many tuning strategies including optimization of the mother wavelet shape and judicious ways of choosing threshold coefficients.At the first stage, defective pixels are removed.Then the filtering method is adopted.Hereafter is presented a comparison of the standard method with the new filtering proposals for infrared images.

Defective Pixel Removal
Defective pixels can be defined as those being significantly deviating from the average behaviour of their neighbours.Two types of defective pixels are commonly identified in the literature: "dead" pixels having very little sensitivity and "hot" pixels, which retain a very high level of energy.Defective pixels do not cause any problem for Gaussian filtering.However, due to their locality, they may affect the performance of wavelet-based functions.It is therefore convenient to filter them out.The algorithm is based on the work of Giron and Correa [19] who defined a unique value of the threshold for the whole image frame.Instead of using the median of all the pixels of the frame, we found that the mean value produces more accurate results for our application.The image is sectioned into square 50 pixels wide areas, where the filtering applies.A smaller area size increases the computational time without obtaining notably better results.

T
For every area, a threshold is defined as three times the standard deviation of the temperature values.Then, every defective pixel being outside the three-sigma interval is marked, and for each one a small neighbouring area 8 pixels wide is considered.Finally, the mean of this new area is calculated (defective pixels exclusive), and the corresponding value is substituted to the corrupted pixel.

Standard Method: Gaussian filtering
All the above methods as well as the standard procedures such as two-dimensional Gaussian filtering produce good results for the naked eye.Fig. 2 shows a typical Fourier spectrum after filtering of the raw frame.However, to solve the heat equation (Eqn. 1) and to obtain a smooth Q(x,y,t) heat source function, we employed a more elaborated solution.Simply increasing the smoothing parameter of the Gaussian filter is not sufficient since it does not remove the spatial irregularities as can be seen in Fig. 3.The comparisons below are made in terms of the generated heat generation per unit volume Q and the integration in a finite area region to obtain a time varying function, the main objective of our energetic approach.It is therefore convenient to treat all relevant data for Q as dimensionless, considering the scope of the present work is to find the crack tip position and to compare different denoising methods.The high emissivity paint is the source of the error that causes apparent thermal gradients.Fig. 3 shows irregularities in the heat field, and in the crack tip area where the shape should be sinusoidal as the load.Time irregularities are due to electronic noise.All the sequences utilized for comparison comes from a test performed with fatigue ratio R=0.3 and F max =6500N starting from around ΔK = 30 MPa√m.Every 750 cycles a 1s sequence has been recorded, this is a good compromise between hard disk space and discretization of crack growth.As it can be seen in Fig. 6b quite reasonable points have been obtained with this procedure.To enhance the anisotropic three-dimensional Gaussian filtering capabilities for motion compensation, a sequence in the unstable stage III for crack propagation has been chosen.The sequence has been recorded at 108000 cycles, just 2000 cycles before fracture, when a = 20.05mm and ΔK = 123.5 MPa√m.

First Method: anisotropic three-dimensional Gaussian filtering
The limit of using a common Gaussian filter is that it works only in the space coordinates while electronic noise is timedependent.Having decreased the exposure time, electronic noise has an energy comparable with the signal deriving from the environment.A more sophisticated filtering technique is needed.The simple solution is build up on the threedimensional Gaussian filter, where smoothing in the time domain can be performed too.Anisotropic Gaussian smoothing is achieved using the following kernel [20]: The filtered result is achieved by convolving the   x, y,t G Gaussian distribution function with the image sequence.
, , x y t    are the values of sigma (smoothing parameter) for respective spatial variables and time.This method permits smoothing in the time domain with the adjusted t  , which is independent of , and x y   and can be arbitrarily smaller than for the space dimensions.Fig. 4 illustrates the results which are apparently better smoothed and are better following the load function behaviour than those displayed in Fig. 3.

Second method: Wavelet Shrinkage
Wavelet coefficients are obtained from a single prototype wavelet ψ called mother wavelet by scaling and translation.Various kernel pairs can be chosen.We choose Daubechies-10 (db10) as a wavelet family which has a good compromise between filtering and valuable information preservation.A level 5 decomposition has been utilized in the filtering procedures.Under Besov norms, the magnitudes of wavelet coefficients are directly proportional to the irregularity of a given image.When noise is present, such irregularity grows in the wavelet coefficients.Because coefficients at finer scales tend to be the primary carriers of edge information, we set discrete wavelet transformation (DWT) coefficients at zero if their values are below a certain threshold Tr.These coefficients are mostly corresponding to noise.The performance of the denoising algorithm relies heavily on the optimal value of the threshold.We noted that for the particular type of thermographs, as opposed to conventional light photography, the method for estimating the threshold is not too influencing.The computational speed and easiness of implementation are however of major concern.We choose universal thresholding as the preferred method whereby a threshold value is uniquely chosen for all wavelet coefficients.The threshold Tr for the coefficients is defined as follows: where N is the image size, and σ is the noise variance taken at unity.A soft thresholding function has been applied.Fig. 5 shows the results for the heat source function Q after wavelet filtering.
Figure 5: Results for the heat source function using a wavelet decomposition level of 5 and db10 (left).Results are very similar to those obtained with 3D Gaussian filtering in terms of smoothness in the space dimension.Temporal filtering is not implemented in this processing scheme and therefore the time dependence of the Q-function is not perfectly sinusoidal and is not following the load well (right).

CRACK TIP LOCATION
o determine and visualize the crack path, one has to spot the position of the crack tip in every recorded sequence as in Fig. 1.During the stable crack growth, every recorded sequence 1 s long corresponds to a certain position of the crack tip; at this small time scale, the crack length can be considered constant.The thermo-elastic effect causes a positive heat generation during unloading, and a negative one during loading.Temperature is thus increasing in the unloading phase, especially at stress singularity points such as the crack tip.The simplest idea is to find the maximum temperature in every frame and to calculate the mean of the maxima positions.This is of course not sufficient to obtain an accurate estimate of the position of the crack tip.The motion of the specimen during testing, friction on the crack edges and plasticity ahead of the crack tip distort the results and create additional errors, which should be compensated in separate procedures.Besides, outliers must be excluded for the sake of accuracy.For the scope of calculating ΔK and locating a region around the crack tip, the procedure has proven itself valuable and economical to analyse large quantity of data.The proposed workflow is as follows:  The temperature and the size of the frame are normalized to unity.
 Only the values in the 0.2 to 0.8 temperature range are kept for the analysis.In this way, the crack surface friction affecting the mean values because of the increased temperature, and the frames containing the highest deformation of the specimen are excluded.At the same time, as the crack closure is approached in the unloading part of the cycle, the peak temperature point moves from the crack tip to the centre of the plastic zone as has been noted also by Meneghetti et al. [8].The limit of 0.8 has been found empirically to take into account this phenomenon as well as the friction effects. In the spatial domain, the values outside the 0.01-0.99interval for the horizontal dimension and 0.2-0.8 for the vertical dimension are ignored.This eliminates the border effects which can be quite significant. Points where the temperature decreases are excluded.This is the most important filter.It has been observed that during the loading phase, due to the thermo-elastic effect, the temperature can decrease below the mean value.The concept for which the maximum corresponds to maximum temperature would not be applicable anymore.The plot of maximum temperature against the frame number (time) is shown in Fig. 6, and is self-explanatory of the process adopted.Blue circles correspond to the frames along the ascending half-cycle curve, which were kept for the analysis.Comparing sequence frames with the calculated crack tip point, it resulted that the fatigue crack tip position can be determined with an error lower than 10 pixels related to the camera resolution.This error also can be noticed in Fig. 7, highlighting the difference between the actual crack tip position and the calculated crack path.The method proves to be particularly efficient and accurate during the steady Paris stage of the crack propagation whereas the error approaches the upper limit of 10 pixels in the crack initiation stage and final critical crack growth stage.In fact, the blunt notch for crack initiation does not create a sharp temperature profile where the maximum is precisely detectable, while large surface displacements occur in the final rupture stage, and the mean value for the crack tip position has less significance.Fig. 6 shows the calculated crack length against the number of loading cycles.Results are reasonable and can be used for the location of the crack tip position.However, they are not recommended for calculations of the crack growth rate da/dN before additional smoothing.A regression to high order polynomial has been tested and shown efficient.

CONCLUSIONS
urther study should be conducted to compare the calculated crack length and the crack growth rate with independent measurements by optical observations, crack gage or crack tip opening displacement data.At the current state of research, the proposed algorithm is accurate enough to locate the crack tip in order to integrate the heat source function Q over a finite area as a proof of concept.However, due to chosen computationally cheap measurement strategy when the thousand of cycles occurred between every IR-recorded sequence, the crack growth rate da/dN cannot be measured satisfactory from the thermographic data in the present test.Nevertheless, the method itself does have this capacity, provided the records are performed continuously.In addition, we have shown that both three-dimensional anisotropic Gaussian filtering and wavelet shrinkage-based filtering perform appreciably better than the conventional 2D Gaussian filtering in obtaining the heat source field.In fact, due to their flexibility and tenability, they offer a good compromise between information preservation and noise filtering.Fig. 7 shows that nearly coincident crack path is predicted with both procedures.The Gaussian filtering method has the advantage of being more straightforward to implement and providing a more accurate time representation.Wavelets are a bit more complex to handle.They however have a hood potential for removal specific noise features.Overall, the anisotropic Gaussian filtering appears to be a viable tool for this type of application.In conclusion, a versatile and fully functional framework for automatic IR data processing is presented.It has been successively applied to a large set of experiments.Future development of methodological tools (besides calibration tools) should be based on precise motion compensation algorithms.

Figure 1 :
Figure 1: Crack path determination from experiments using the temperature distribution T(x,y,t)

Figure 2 :
Figure 2: Example of a raw frame with the rainbow color-coded temperature (left) and its Fourier spectrum (middle).Red colors indicate the higher temperature points.On the right-hand image, the FFT spectrum derived from the Gaussian filter.

Figure 3 :
Figure 3: The heat source function ahead of the fatigue crack after the standard Gaussian filtering (left) and the total irradiated heat in the process zone (right).Here and in the following figures the values for the heat function are normalized and dimensionless; therefore, only the shape is important.The red circle corresponds to the frame plotted on the left.

Figure 4 :
Figure 4: Results of the heat source function with the three-dimensional Gaussian filtering applied (σx,y=8.8σt=0.88)(left).The generated energy is shown in the neighbourhood of the crack tip area (right).

Figure 6 :
Figure 6: Maximum temperature in the frame versus the frame number in a sequence.Circle marks are the points used to calculate the crack tip position.On the right figure, the crack length is plotted as a function of the number of cycles.Blue lines correspond to the calculated points; black circles are obtained from a high degree polynomial fitting.Inaccuracies occur in derivative calculations.

TFigure 7 :
Figure 7: Crack path determined for the test specimen from IR images (cropped from a 640x512 px, 9.6x7.7 mm).Cross points represents the crack tip position located after 3D Gaussian filtering and Q-function calculation.Diamond marks are calculated with the wavelet shrinkage procedure.