Processing of high-speed videos of shock wave boundary layer interactions

This paper presents a framework for processing high-speed videos recorded during gas experiments in a shock tube. The main objective is to study boundary layer interactions of reflected shock waves in an automated way, based on image processing. The shock wave propagation was recorded at a frame rate of 500,000 frames per second with a Kirana high-speed camera. Each high-speed video consists of 180 frames, with image size [768×924\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$768 \times 924$$\end{document}] pixels. An image processing framework was designed to track the wave front in each image and thereby estimate: (a) the shock position; (b) position of triple point; and (c) shock angle. The estimated shock position and shock angle were then used as input for calculating the pressure exerted by the shock. To validate our results, the calculated pressure was compared with recordings from pressure transducers. With the proposed framework, we were able to identify and study shock wave properties that occurred within less than 300μsec\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$300\, \upmu \hbox {sec}$$\end{document} and to track evolveness over a distance of 100 mm. Our findings show that processing of high-speed videos can enrich, and give detailed insight, to the observations in the shock experiments.


Introduction
The introduction of charge-coupled device (CCD) and complementary metal-oxide semiconductor (CMOS) technology revolutionized high-speed photography during 1980s and 1990s. Today, there are numerous types of high-speed cameras that operate at a frame rate of more than a million images per second with resolution of one mega pixel. This development enables researchers to capture fast phenomena like the shock wave propagation and gas explosion in a video as a series of images [1]. Based on the advanced image and computer technology, these images now give an alternative source to estimate the shock wave characteristics like shock speed, pressure etc. [2,3]. However, in early days, images were mainly used for visualization of the phenomena [4]. From an image processing perspective, the ability to extract the desired information automatically was limited.
Image processing is an interdisciplinary research field, where the aim is to extract some desired information from B Samee Maharjan samee.maharjan@usn.no 1 University of South-Eastern Norway, Kjølnes ring 56, 3918 Porsgrunn, Norway images. It is widely used for object detection, object location, classification, segmentation and motion detection, among others. Based on improved computer power, sensor development and algorithmic progress, the applications have been broaden into various other fields, for example security system, road safety, document enhancing and gas dynamics. The study of a shock wave generally called shock on the other hand is one of the important research area in gas dynamics [5]. Over the years, the methods for estimating shock characteristics have evolved from local pressure sensor recordings and numerical simulations [4,6] to advanced simulation techniques [7]. The use of image processing in this particular field was evolved during early 2000s. One of the main advantages of using high-speed video/images over a traditional sensorbased approach is the possibility of extracting large amount of information within a small-scale experiment (< 100 mm). It makes it possible to extract information with intervals < 1 mm, meaning that several images may be recorded in between two pressure transducers. In [2,3], a rather simple image processing method using intensity difference successfully determined the shock front position in the images. These images were captured at low frame rate; thus, they are of a high quality. However, in case of higher frame rates, the quality of images degrades. Thus, improved image process-ing methods are needed. Some previously developed image processing framework to process similar high-speed videos and determining shock front position can be read in [8,9]. The filtering methods suggested in [8,9] have a multistep spatial filtering unit which increases the computing time. Furthermore, in some of the noiser images, the corresponding segmented images contain more background noise.
In this paper, we propose an image processing framework based on a fast preprocessing unit in combination with a robust tracking algorithm. The proposed framework processes images from high-speed videos, recorded during shock wave boundary layer interaction (SWBLI) experiments. The tracked wave fronts were further analyzed to estimate (a) the shock position; (b) position of triple point; and (c) shock angle. The pressure exerted by the reflected shock was calculated by using the incident shock speed, the shock angle and shock polare. For validation, the estimated pressure was compared with the pressure measurement from the transducer.
The rest of the paper is organized as follows. A description of the experimental setup and the high-speed videos is described in Sect. 2. The designed image processing framework along with the method to approximate the shock wave characteristics is described in Sect. 3. The estimated results are presented and discussed in Sect. 4 followed by the conclusions in Sect. 5.

Experimental setup
The experiments were conducted in the Graduate Aeronautical Laboratories at California Institute of Technology (GALCIT) detonation tube [2]. The tube is a 7.6 m long, 280 mm inner-diameter and equipped with a 152.4 mm wide test section and two polycarbonate windows to provide optical access.. The optical window constructed in the test section makes it able to capture shock wave boundary layer interaction when the incident shock is reflected back by the reflecting wall. For details of experimental setup, refer to [2]. The planar incident shock in the CO 2 gas was driven by a detonating slugg of C 2 H 2 − O 2 . The initial pressure was 13 kPa, and Mach number of the incident shock was 2.4.

Shock wave boundary layer interaction
A schematic diagram of a test section is illustrated in Fig. 1. As shown in Fig. 1a, the first event is a planar incident shock propagating toward the reflecting end wall. A boundary layer is simultaneously formed at the tube wall behind the propagating incident shock as seen in Fig. 1b. When this shock hits the reflecting end wall, the reflected shock then propagates Fig. 1 Schematic representation of the shock wave propagating in a shock tube; a shock tube with incident shock; b formation of a boundary layer (B-L) behind the incident shock at the tube wall; c shock tube after the incident shock reflected off the end wall; d reflected shock interacts with the (B-L) and formed a triple shock configuration at the tube wall backward illustrated in Fig. 1c and interacts with the boundary layer. A Mach stem structure with a triple point is formed as a result of the interaction which is shown in Fig. 1d. The phenomena is known as a bifurcation of a reflected shock wave. Mark reported a classic study of boundary layer interaction of reflected shocks in a shock tube in 1958 [4]. The distorted reflected shock close to the wall preceding the normal reflected shock is commonly called an oblique shock or the foot. The third shock behind these two shocks is known as rear shock. The point where all three shocks meet is what called a triple point, and the angle that the oblique shock makes with the boundary is a shock angle. The triple point configuration is also known as a Mach stem structure; further description of bifurcation of reflected shock and Mach stem can be found in [4,5].

Imaging technology
The two imaging techniques used for capturing the wave propagation, namely shadowgraph and schlieren [10]. They translate the phase speed difference of a light passing through the medium, into different intensities in a viewing plane (image). The use of these imaging technologies in the flow visualization can be read in [11]. For the review of recent developments in shadowgraph and schlieren techniques, refer to [1].

High-speed videos
A high-speed camera especially introduced for capturing a fast phenomena like the compressible gas flow named Kirana 1 was used to record the SWBLI experiments. Even though the camera was operated at the frame rate of 500,000 frames per second, the SWBLI is an extremely fast phenomenon, which occurs in less than a millisecond. Thus each highspeed video consists of 180 images of size [768×924] pixels which is about [81 × 97] mm or [60 × 72] mm, depending on the scaling factor of 9.528 or 12.902 pixels/mm, respectively. In this paper, we processed two experiments, one done with schlieren technique (Exp. 2558) and another with shadowgraph (Exp. 2516). A sequence of images from the high-speed video (Exp. 2516) is shown in Fig. 2. The images in the first row show the propagation of the incident shock from left to right toward the reflecting end, and the images in the second row show the propagation of the reflected shock from right to left. The formation of a triple point structure is clearly visible in the second row.

Image processing
In short, an image processing framework was designed to perform these following tasks: 1. Read each image from a high-speed video. 2. Reduce background noise while preserving desired information.
1 https://www.specialised-imaging.com/products/video-cameras/ kirana. 3. Track/locate the front of the shock wave in images with a visual front. 4. In case of reflected shock, segment the tracked front into a normal shock and an oblique shock. 5. Calculate the shock angle and the height of the triple point.
Due to the presence of the camera window in the upper part in all images (see Fig. 2), only the lower 400 rows were considered for the rest of this paper. For the human eyes, the fronts are easily captured in the images, and however, to get a computer program to perform the same task is not trivial. To elaborate this, a 1-D intensity profile of a single row is given in Fig. 3. This figure shows the ambiguity of the solution. The framework to process high- speed videos and to overcome the above-mentioned problem is described in the following subsections.

Image filtering
Background noise in the images comes mainly from the filming technology, i.e., the camera, the experiment equipments and ongoing changes in gas properties. The filtering of an image in the current framework consists of a background subtraction followed by a low-pass frequency filtering. For background subtraction, a background image was generated from each high-speed video by averaging all images without a front. The number of images without fronts varies within each video, and sometimes there might be several images and sometimes just one. The background image was constructed with a prior study of the images and manually choosing the images without shock fronts. The constructed background image was then successively subtracted from all the images with a visual front. The result of background subtraction in one of the images from each experiment is shown in Fig. 4.
Filtering of an image can be defined as an operation in which the value of the any output pixel is determined by the combination of the values of the pixels in the neighborhood of the corresponding input pixel. Convolution is a filtering algorithm which takes a weighted sum of the neighboring input pixels, and the matrix defining the weights is known as convolutional kernel. Filtering using the convolution of a kernel matrix with the original image in the spatial domain is operationally costly, when the image size is too large. In For filtering in frequency domain, at first the original image (spatial) was transferred into frequency domain by using fast Fourier transform (FFT). Then, multiplied the transferred image with the filter function and thereafter, retransformed back to the spatial domain by using inverse Fourier transform. A low-pass filter (LPF) attenuates high frequencies greater than a cutoff frequency, resulting a smoother image in the spatial domain. A low-pass filter of order 3, having a cutoff frequency of 50 Hz, was used for filtering the images. The frequency spectrum of Fig. 4b can be seen in Fig. 5a and Fig. 5b is the result after multiplication of Fig. 5a with LPF, whereas Fig. 5c is the final filtered image. Further details regarding fast Fourier transform and filtering of images in frequency domain can be read in [12]. Figure 5d shows the result of using 'prewitt' edge detection method from MATLAB image processing toolbox corresponding to Fig. 5c. Some of the available edge detection methods in various processing toolboxes were able to detect the edges; however, they did not provide the required clarity of front position. This suggested further processing of the filtered images in order to track the exact front position.

Segmentation
The term image segmentation refers to the partition of an image into a set of regions that cover the entire image. For our application, this implies separating the shock front from the background. To do so, all images were normalized to intensity values in range [0-1]. The images were segmented into a background (black) and the shock wave (white) based on a thresholding technique suggested by Otsu [13]. The Otsu's algorithm assumes that the image contains two classes of pixels following bi-modal histogram (foreground pixels and background pixels). The algorithm then calculates the Fig. 6 The result of segmentation (after low-pass filtering) performed in; a Fig. 4a; b Fig. 4b; c an example of segmentation where filtering did not suppress all noise; d a template consists of two 5 × 5 matrices for the considered pixel optimum threshold separating the two classes with minimum variance. For example, a threshold value given by Otsu's method corresponding to Fig. 5c was 0.25. All the pixels with an intensity level below 0.25 are labeled as background, and all pixel values equal to or greater than 0.25 are labeled as foreground. The result of applying Otsu's algorithm (after low-pass filtering) in Fig. 4a and b is shown in Fig. 6a and b, respectively.

Front tracking using template matching
To get an algorithm to track the front in an image like Fig. 6a,b is rather simple, just slide through all pixels and mark the first white pixel (from left) in each row. However, the output from the segmentation step, i.e., Otsu's algorithm, does not always work that perfect. The amount of noise varies from image to image, video to video, and can affect the segmentation result shown in Fig. 6c. Thus, to overcome these situations, a template matching technique was designed for automatic track of the fronts. A template matching technique is a well-known method for image classification and segmentation. In short, with this approach, the aim is to find the position/object in the scene which match best with a predefined template [14]. For the present work, a template which consists of two 5 × 5 matrix is created as shown in Fig. 6d. As it can be seen in Fig. 6d, the matrices are shifted such that a single template can be used for tracking both the normal (straight) and oblique (tilted) shock. The element values for the left-hand matrix of the template are set to 0 (black), while for the right-hand matrix to be 1 (white). The framework matches the template pixel by pixel, row by row, and calculates a matching error for each location. The matching error is the sum of the differences between the pixel values of the template and the footprint template created around the considered pixel in the image. For example, in Fig. 6a-c, around the top left of the image where only background is present, the left-side matrix of the template will match perfectly. However, the right-side matrix will not match at this location, thus resulting a huge error. Moving further along the columns, the template will give minimum error around the front where both the matrices seem to match perfectly. To eliminate the loss of result due to boundary pixels in the lower boundary, the last row was replicated for five more rows than the actual data. After the matching process, the pixel in each row where the template matches the best, i.e., the point with the minimum error was chosen. For most of the images, the tracking was accurate; however, for a few images in Exp. 2558, some mistracked points were observed, as seen in Fig. 7a. To overcome this problem, the following a priori information was incorporated: 1. The oblique shock wave is ahead of the normal shock. 2. The position of triple point does not goes above row no. 150.
The median of the normal shock (row 50:150), tracked at the first (green curve in Fig. 7a) was calculated which is shown by the yellow line in Fig. 7b. Later, during the second tracking, the framework tracks the front only till the calculated median, such that the noise present beyond the median will be excluded. The median line with the final front is shown in Fig. 7b.

Estimating triple points based on segmented regression
To estimate the triple point and the shock angle, the tracked front was divided into a normal shock and an oblique shock, by using a segmented regression technique. Segmented regression is a method of fitting multiple lines from a single dataset [15]. In case of fitting two lines from a single dataset, it requires only one break point BP, and the model can be written as, In our case, the x-y points of the tracked front served as a dataset from which two straight lines were fitted: one for the normal shock and one for the oblique shock, as illustrated in Fig. 8. Please note that, the segmented regression process was conducted in the x-y coordinate system. Figure 9 summarizes the segmented regression process of determining the optimum BP by consecutively fitting two lines to the underlying dataset. For the sake of understanding, the process is demonstrated in the images. The blue curve represents the tracked front, and the yellow mark yields the breaking point BP, while the orange line represents the first line fitted as in (1). The white line is then fitted to the remaining points below BP as in (2). The process follows a brute force approach, from top to bottom, starting BP at row no. 10 and for each iteration, the BP moves down by one row. After each iteration, the least square error is calculated for both fitted lines with their respective location and the errors are summed up and stored. For example, in Fig. 9a the orange line will give a small amount of error as it almost coincides with blue curve above the yellow point, while the white line is misplaced and causes large errors. Figure 9b shows the fitted line when the separating point is located near the triple point and thereby the fitted lines representing the normal and the oblique shock are rather accurate. The total error gathered from both lines seems to be at a minimum with this location. Further down, see Fig. 9c, the error increases. For each image, after the process was finished, the BP that gave the minimum error was considered as the separating point between the normal and the oblique shock. A close look around the triple point is presented in Fig. 9d, where the point which gives the minimum error between two fitted lines is visualized with a yellow mark. Now, based on the selected point (yellow), two straight lines are refitted, the red line with respect to points from the tracked front below the yellow point, whereas the green one is fitted with respect to 100 points above it. (By trial and error 100 points were used in this study, but other selections may be used). The white circular point shows the crossover between the two fitted lines, i.e., the triple point, while the slope of the red line gives the shock angle.

Results and discussion
In this work, an image processing framework including frequency filtering, template matching and segmented regression was developed to track the wave front in multiple highspeed videos from shock experiments. The fronts tracked in Exp. 2516 are presented in Fig. 10. Figure 11 shows the comparison between the proposed algorithm with previous work [8,9]. The segmented images in Fig 11a and b demonstrate that the segmentation with FFT filtering performs better than median filtering used in the segmented images demonstrates that the segmentation with FFT filtering performs better than median filtering used in [8,9]. In addition, the average time taken by the proposed method to track the front is about 4 sec per image, which is a second less than by [9] and around 20 sec (45 iterations) less than [8]. Comparison between results from proposed framework with results from [8,9]; a the segmented image and a tracked front by [9]; b the segmented image and the front tracked by proposed method; c the fronts tracked by all three methods (blue- [8], green- [9] and redproposed method; d enlarged around triple point (color figure online) the actual front. However, the differences are small and is a subjective matter. The first few frames are excluded while calculating presenting the heights of triple points and shock angles as the number of points representing the oblique shock was too sparse to have a valid calculation. The heights of the triple points from the lower wall have a likewise behavior in both experiments as can be observed in Fig. 12. The height increases as the front moves away from the reflecting wall, however at the later part beyond 50 mm from reflecting wall, and the height flattens out. A similar behavior has earlier been observed by [16]. The shock angles calculated for both exper-  iments are plotted in Fig. 13. The black diamond shape marks in Fig. 13 represent the manual measurement performed on the images of Exp. 2516. The images were selected such that the position of foot of shock lies around 10,20,30…mm from the reflecting wall. The comparison with the manual measurements shows that the shock angle calculations are in good agreement. Further, to compare the results derived from the image processing with actual pressure recording, a graphic method known as the shock polare [11] was used. Basically, it makes possible to find the shock strength, i.e., pressure ratio across the shocks (both normal and oblique), from knowing the Mach number (speed) of the incident shock and the angle of the oblique shock. The incident shock speed was calculated by using nonlinear regression for line fitting on the position of the incident shock [2]. The shock polare for Exp. 2516 is shown in Fig. 14 (left). The blue curve in Fig. 14 (left), gives the relation between deflection angle (theta) and the pressure ratio across the oblique shock (p/p1). When the shock angle, as shown in Figs. 2 and 13, is known, the deflection angle can be calculated with shock equation [11]. Figure 14 (right) shows the pressure recordings of the transducer located 50 mm from the reflecting wall. We manually select the image from the high-speed video in which the foot position of the wave was located closest to the selected transducer. The corresponding shock strength in the shock polare agrees with the shock strength recorded in the pressure recording as shown by red marks. Another shock polare (red curve) originated at shock strength point in the blue curve is the shock polare for the rear shock (see Fig. 2). Behind the normal shock and the rear shock, deflection angle and pressure will be the same [11]. Therefore, the upper crossing of the blue and the red polare represents this solution for these parameters.
So, from image processing we are able to predict the pressure and the state in the triple point configuration (above the boundary layer). This information will not be available from a pressure recordings, since pressure will be influenced by boundary layer interactions as they were mounted in the wall. To this end, it can be concluded that the high-speed videos in combination with image processing can enrich and give new detailed insight to shock wave boundary layer interactions.

Conclusions
In this paper we have used high-speed video in combination with image processing as a framework for studying reflected shock wave and shock wave boundary layer interactions in detail, within time intervals of 300 µs. Frequency filtering showed an overall accurate performance both with respect to robustness and precision for the application at hand. Template matching is a simple, but still powerful technique for identifying specified features in an image. In case of front tracking, misplacements of the front was seldom observed in this study. Segmented regression proved to be an efficient tool for dividing the tracked wave front into two parts: the normal and the oblique shock. However, the template matching and the segmented regression took almost 80 percent of total processing time. Hence, further work can be done to obtain a speed-up concerning these steps. The pressure ratio across the shock wave was calculated with a shock polare, which then was compared with the measurements from actual pressure transducers. The comparison revealed that our estimates are in good agreement with the pressure recordings. The overall processing of one high-speed video takes less than an two hour to complete which is comparatively lesser than any small computer simulation done for the same study. With the proposed framework, we were able to identify and study shock wave properties that occurred within less than 300 µsec and to track evolvement over a distance of 100 mm.