1 Introduction

Infrared (IR) thermography has been commonly used for global detection of laminar-turbulent transition fronts. This transition-detection technique exploits the fact that shear stress correlates directly with convection rate. Given nonzero heat transfer through the skin of the model under study, a temperature differential will develop between the laminar and turbulent regions. This differential can be effectively imaged with the use of an IR camera. Once the surface temperature field is measured, these data must then be interpreted in some manner to extract the transition location.

Traditionally, this has been done visually by a human. Unfortunately, humans are notorious for inducing bias error when subjectively evaluating images. Additionally, in a crossflow-dominated environment, the transition front is typically quite jagged. In effect, the human interpreting the images must read some single representative transition location from a very complex signal. In general, measurements taken in this fashion can only be trusted to accuracy on the order of \(\pm 5\,\%\) chord. Furthermore, IR images are not solely functions of shear stress. Variations in heating or cooling, surface emissivity, internal conduction due to structural members, reflections from nearby hot bodies, and numerous other phenomena introduce significant noise into the IR image. Finally, since each image must be analyzed by hand, a significant amount of time is required to adequately reduce the data from a substantial number of images.

As a result of these difficulties, increasing efforts have been made in processing IR thermography in an automated fashion. Raffel and Merz (2014) have developed a technique using differential IR thermography that provides high-contrast images of transition that are well suited for automated processing. In order to generate these high-contrast images, the motion of the transition front is leveraged. However, this is not suitable for analyzing a stationary crossflow-induced transition front where the transition location is quite stable. Additionally, Richter and Schülein (2014) have developed a technique for measuring Tollmien–Schlichting-induced transition on a moving helicopter rotor using IR thermography. This demonstrates one of many scenarios where IR thermography can provide excellent quantitative data, even in a dynamic or otherwise difficult to measure environment. One of the main remaining challenges is the ability to measure transition in a crossflow-dominated environment with the corresponding characteristic sawtooth transition pattern. In a crossflow-dominated transition front, the transitional region is a much smaller area and varies strongly as a function of span. This introduces a different set of challenges.

The technique described in this article was developed during testing on the SWIFTER (Swept Wing In-Flight Excrescence Research) and SWIFTEST (Swept Wing In-Flight Testing Excrescence Stability Theory) research models at Texas A&M University in both the wind tunnel and flight environment. The sample reduced data presented were acquired using these models. Additional details about SWIFTER and SWIFTEST can be found in Duncan et al. (2014a, b) and Crawford et al. (2015), respectively. In a typical 1.5 h SWIFTEST flight, approximately 300 test points are measured. If a skilled person is reducing the data, 6 s per image is a reasonable pace. At that rate, reducing 300 images takes about half an hour. In practice, bursts of approximately 20 images are used to get more accurate measurements. This would take approximately 10 h to reduce just one flight worth of data. Due to all of these factors, a more-rigorous computer-automated analysis technique was developed.

2 Experimental setup

IR thermography has two key experimental requirements: there must be heat transfer through the surface being studied, and the surface must be capable of holding a crisp temperature gradient.

2.1 Temperature differential generation

In order to generate a measurable temperature gradient for IR thermography, the model under study must be either heated or cooled to drive the laminar-turbulent surface temperature differential. Generally, this temperature differential must be induced by external means. One option is to heat or cool the surface via means such as heating wire or water cooling. The other option is to heat or cool the freestream fluid, whether by rapid altitude changes in flight, or via cycling of the cooling system in a wind tunnel. Only a few degrees of heating or cooling is required to obtain a very clear image. This can be thought of as a \(T_{\mathrm{wall}}/T_{\mathrm{wall,adiabatic}}\) variation on the order of one percent for most conditions. Such minimal temperature variation does not significantly affect the boundary-layer stability for most phenomena. Verification of this fact can be readily achieved by disabling the heat/cooling source, and carefully watching for motion of the transition front as the model approaches equilibrium.

Fig. 1
figure 1

Heating-sheet diagram

Fig. 2
figure 2

Notional temperature map

For the SWIFTER/SWIFTEST model, internal, electric heating wire was used, along with an internal insulating blanket, to direct the majority of the heat transfer through the test surface of the model. Figure 1 shows a notional graphic of the configuration used. The power density required to achieve sufficient temperature differential at a flight \(Re' = 6.0\,\times \,10^6/\mathrm {m}\) was approximately 300 \(\hbox {W/m}^2\). Off-the-shelf heating wire designed for heating hardwood and tile floors is ideal for this purpose. It is both cheap and readily available. It also achieves these power densities with 25-mm wire spacing, which results in fairly uniform heating through a 3-mm aluminum skin. Furthermore, the wire is already shielded and sheathed, which serves as both a safety precaution against nicking and shorting, as well as electrical shielding which helps reduce any interference with neighboring electronics when the wire is switched on and off to control the temperature. Closed-loop temperature control is easily achieved using an off-the-shelf proportional-integral-derivative (PID) controller and a pair of resistance temperature detectors (RTDs); one is on the inside of the heated surface of the model, and the other is on the unheated side of the model. The second RTD serves as a reference temperature, so that temperature differential relative to an unheated wall can be maintained, rather than absolute temperature. This allows for operation in widely varying freestream temperatures without user intervention.

2.2 Surface coating

The other primary requirement for IR thermography is a surface coating with high emissivity in the IR band being imaged, as well as a low thermal conductivity. Additionally, the coating must be at least on the order of 300 \(\upmu \hbox {m}\) thick to support a crisp temperature gradient. If the coating is too thin, the conductive substrate washes out the temperature gradient. However, it is also advantageous for the underlying substrate to be highly thermally conductive when using internal heating/cooling in order to help distribute the heat for a more uniform internal temperature. A notional drawing of the temperature gradients involved is shown in Fig. 2. Most of the thermal gradient is in the paint itself, while the underlying aluminum substrate is much more uniform in temperature and serves to spread out the heat from the discrete wires. The high emissivity of the coating is required in order to accurately read the surface temperature using an IR camera, as well as to reduce reflections on the surface due to neighboring hot or cold bodies.

To fulfill the coating requirements on the SWIFTER/SWIFTEST model, six coats of Sherwin Williams F93 lusterless black aircraft paint, conforming to MIL-PRF-85285 Type 1, Class H, were applied to the model over the corresponding recommended primer. This particular paint has the added benefit of being representative of what would be applied to a typical military transport aircraft. Furthermore, it is both relatively inexpensive and readily available through most aircraft paint suppliers. It is important to note that the surface must not be sanded. Even the lightest sanding destroys the high-emissivity properties of the paint.

2.3 Camera

For these tests, a FLIR SC8100 IR camera was utilized. The SC8100 has a spectral range of 3–5 \(\upmu \hbox {m}\) and a sensor resolution of 1024 \(\times\) 1024. The camera is capable of detecting temperature differences as small as 25 \(\mathrm{m} \hbox {K}\). The spatial and temperature resolution assist in detecting smaller scale features such as crossflow and the fine sawtooth pattern. However, lower resolution cameras can still successfully resolve the general transition front as the difference between the fully laminar and fully turbulent transition regions are on the order of a couple degrees and span the entire model. For the SWIFTER/SWIFTEST wind tunnel testing, a 17-mm lens was utilized, while a 50-mm lens was utilized for the flight testing.

3 OpenCV

The OpenCV library is utilized to dramatically simplify many of the operations used in the following image processing steps. OpenCV provides functions for many generic filtering operations, as well as a relatively thorough model for handling camera and lens distortion. Several functions for estimating position and point (PnP) are also provided, utilizing the aforementioned camera model. The image filtering is fairly straightforward; however, the camera model and PnP estimation deserve some brief explanation for clarity since they are used heavily throughout the following algorithm. A much more thorough treatment of these algorithms can be found in the OpenCV manual.

The camera model utilized is nonlinear and accounts for variables including image-sensor geometry, lens distortion, and perspective transformation. The full model is as follows:

$$\begin{aligned} \begin{bmatrix} x_c \\ y_c \\ z_c \end{bmatrix}&= \mathbf R \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} + \mathbf {t} \end{aligned}$$
(1)
$$\begin{aligned} \begin{bmatrix} x_c' \\ y_c' \end{bmatrix}&= \frac{1}{z_c} \begin{bmatrix} x_c \\ y_c \end{bmatrix} \end{aligned}$$
(2)
$$\begin{aligned} r&= \sqrt{x_c'^2 + y_c'^2} \end{aligned}$$
(3)
$$\begin{aligned} x_c''&= x_c' \frac{1 + k_1 r^2 + k_2 r^4 + k_3 r^6}{1 + k_4 r^2 + k_5 r^4 + k_6 r^6} \nonumber \\&\qquad + 2 p_1 x_c' y_c' + p_2 \left( r^2 + 2 x_c'^2 \right) \end{aligned}$$
(4)
$$\begin{aligned} y_c''&= y_c' \frac{1 + k_1 r^2 + k_2 r^4 + k_3 r^6}{1 + k_4 r^2 + k_5 r^4 + k_6 r^6} \nonumber \\&\qquad + p_1 \left( r^2 + 2 y_c'^2 \right) + 2 p_2 x_c' y_c' \end{aligned}$$
(5)
$$\begin{aligned} \begin{bmatrix} u \\ v \end{bmatrix}&= \begin{bmatrix} f_x&0 \\ 0& \quad f_y \end{bmatrix} \begin{bmatrix} x_c'' \\ y_c'' \end{bmatrix} + \begin{bmatrix} c_x \\ c_y \end{bmatrix} \end{aligned}$$
(6)

Equation 1 defines the standard rotation–translation transformation between the model coordinates (XYZ) and the camera coordinates \((x_c,y_c,z_c)\). Equation 2 transforms the camera coordinates into homogeneous camera coordinates \((x_c',y_c')\). Equations 3, 4, and 5 remove the distortion due to camera lens effects. The variables \(k_{1-6}\) are the radial distortion coefficients, while \(p_{1-2}\) are the tangential distortion coefficients. These coefficients are determined a priori via imaging multiple fiducial grids of known geometry at different orientations. Solving for the coefficients then becomes a nonlinear least-squares optimization problem, which is readily solved using the Levenberg–Marquardt algorithm. Equation 6 is a standard pinhole-camera model, which takes the undistorted homogeneous coordinates, \((x_c'',y_c'')\), and transforms them into their image-sensor coordinates, (uv). The \(f_x\) and \(f_y\) coefficients are the focal depth of the camera, and the \(c_x\) and \(c_y\) coefficients are the image-sensor origin, in pixel coordinates. These values are readily available in the spec sheet of the FLIR SC8100 camera and its lenses. However, to help account for manufacturing tolerances, these values were solved for using the same technique as the lens distortion coefficients. These coordinate systems are also indicated in Fig. 3. The homogeneous and undistorted systems are not indicated in the figure as they are nominally identical to the \((x_c,y_c,z_c)\) frame with the exception of scale and distortion effects.

Once the camera is fully calibrated, the unknown rotation matrix, \(\mathbf R\), and translation vector, \(\mathbf {t}\), can be solved for using nonlinear least-squares given at least six corresponding image-sensor-coordinate to model-coordinate pairs. The OpenCV function that performs this calculation is called solvePnP.

4 Image filtering

The data directly from the SC8100 camera are in the form of raw integrated sensor counts that relate nominally to the intensity of the IR signal received. Before every test, a uniform cold plate is imaged, followed by a uniform warm plate. These calibration images are used to normalize the gain and offset of the entire sensor on a per-pixel basis. This removes much of the static noise from the image sensor due to pixel-to-pixel variation in sensitivity. Additionally, these images are used to identify any dead, stuck, or flickering pixels for later removal. Periodically, a more in-depth temperature calibration is also performed. This calibration consists of imaging several plates of different known uniform temperatures. These known temperature plates are used to develop a set of calibration curves transforming the raw sensor counts into radiance values, and then from radiance into temperature.

When processing these data, the span and offset correction from the two-plate test is first applied to the raw data. These corrected data are then converted into temperature readings using the multi-plate calibration curves. Next, any damaged pixels are reconstructed by interpolating from neighboring good pixels. Up to this point, the process is identical to that performed by the off-the-shelf FLIR ExamineIR software provided with the camera. The freestream temperature from a freestream temperature probe is then subtracted from the image. A copy of the image is then spatially low-passed using a Gaussian filter. The original image is then divided by the low-passed copy. The result is the local percent variation in heat transfer. In dividing by the low-passed temperature, any spatially slow-varying phenomena such as heating-sheet variation and varying heat-transfer rate due to structural members are suppressed. All that remains is the local variation in heat transfer due to shear stress. Since shear stress is directly proportional to convection rate, the percent variation in convection rate is also the percent variation in shear stress. As such, the image is now in terms of local percent variation in shear stress. It should also be noted that this conversion to percent variation in shear stress is completely agnostic of whether the model is heated or cooled. It is also largely independent of the degree of heating or cooling. The only constraint is that if the degree of heating or cooling is too low, the denominator becomes small, and the signal becomes small, causing the remaining sensor noise to dominate, thereby washing out the signal; however, a few tenths of a degree differential is sufficient for a clear image.

5 Image coordinate mapping

The first step in rigorously detecting a quantitative transition location from an IR image is to map each pixel on the camera’s image sensor to the corresponding model coordinate that it is measuring. In general terms, one must first track some identifying features on the model, such as fiducials. These fiducials are then used to determine the PnP of the camera. The image is finally undistorted using a combination of the measured PnP of the camera relative to the model in addition to the a-priori knowledge about the geometry of the model, the geometry of the camera sensor, and the distortion caused by the non-ideal behavior of the camera lens.

Fig. 3
figure 3

SWIFTEST in flight (top), rotated to match IR camera coordinates; IR camera FOV in solid blue; Region of interest in dashed white. Sample IR image (bottom)

Figure 3 shows a visible spectrum image of the SWIFTEST model in flight on top, with the corresponding IR image on bottom. The image is rotated such that flow is from lower left to upper right in order to avoid confusion as the images are transformed into their final coordinate system with flow from left to right. This keeps the images consistent in orientation between the flight, wind tunnel, and computational fluid dynamics solutions. The approximate field of view of the IR camera is outlined in solid blue, and the region of interest for these tests is outlined in dashed white. Once the undistortion process is complete, the output image is bounded by the dashed white line and referenced to the xz coordinate system indicated. The fiducials to be tracked are visible in the image on the left as three silver squares, and as three white squares in the IR image.

5.1 Fiducial tracking

The fiducials used for SWIFTER/SWIFTEST consist of three 51 mm squares of silver mylar tape. These were chosen to have a very low emissivity to provide a sharp contrast with the high-emissivity surface of the model. They were placed in a triangle configuration just outside the region of interest on the model. Their locations were determined relative to the known pressure-tap locations on the model to within 100 \(\upmu \hbox {m}\) using calipers. This provides four vertices of known location per fiducial for a total of 12 vertices to be tracked with known locations in the model-coordinate space.

In order to facilitate tracking of the fiducials, a copy of the image of interest is first spatially high-passed on the scale of the fiducials themselves. This removes any features in the image larger than the fiducial which makes the fiducials themselves significantly easier to track. Additionally, the image is spatially median-filtered to remove much of the high-frequency noise present in the IR sensor itself, while preserving the edges of the fiducials. Finally, a Gaussian filter is used to further remove any high-frequency noise. The result is shown in Fig. 4. After filtering, the magnitude of the spatial gradient is taken using a Sobel filter and stored for later use.

Fig. 4
figure 4

Filtered fiducial tracking image

At this point, a coarse set of fiducial locations is generated by thresholding the image into black and white at a low threshold. Any closed contours are stored, and the process is repeated for a slightly higher threshold. This process continues until the image is completely black. The resulting set of closed contours is then reduced to quadrilaterals using the Ramier–Douglas–Peuker algorithm. Any contours that exhibit large error when reduced to quadrilaterals are discarded. If a resulting quadrilateral is convex, it is also discarded, as there is no orientation of a concave fiducial that would result in a convex projection in the image space. Finally, any highly skewed quadrilaterals are discarded, as the model is not expected to be at a highly oblique viewing angle.

The resulting quadrilaterals are then further refined using the previously computed gradient data. The first step of the refinement consists of marching along each side of the quadrilateral on a per-pixel basis and taking a slice from the image normal to the side. The coordinate of the maxima of this slice is then stored. Once all of the local maxima locations are computed along the side, a linear least-squares fit is computed for that side. This results in four lines per quadrilateral. The four intersections are computed from these lines, and those are stored in floating point as the exact vertex locations. Since the entire side of the fiducial is examined in aggregate, sub-pixel accuracy can be obtained. Additionally, this method is robust to peeling and/or damaged corners on the fiducials (Fig. 5).

Fig. 5
figure 5

Filtered fiducial tracking image with fiducials in arbitrary order

5.2 Model tracking

At this point, the locations of the fiducial vertices are known in both model space and image-sensor space. However, each vertex in model space must still be matched to its corresponding vertex in image-sensor space. Since the fiducials are not labeled, they must be identified by their relative location. This is complicated by the fact that the model is currently in a completely unknown orientation. This problem can be broken into two separate sub-problems. First, the three fiducials must be identified. Second, each vertex of each fiducial must be identified relative to the fiducial itself.

In general, there are six permutations mapping three fiducials between the two systems. If the convention of numbering the fiducials counter-clockwise about their collective centroid is assumed, the number of possible mappings is reduced from six to three since the fiducials cannot be viewed through the back of the model. Next, the affine transform is computed that maps the three fiducials from the image coordinates of their centroids, to the known span/chord coordinates for each of the three remaining permutations. If the sum of the products of the square portion of the affine matrix is taken for each of the three transforms, the value will only be positive for the case where the spanwise leg of the triangle formed by the three fiducials is within \(\pm 45^\circ\) of the spanwise direction. This value will only be positive for one of the remaining permutations, provided that the first two fiducials are approximately vertical with respect to each other within the camera frame. Given that one constraint on camera placement, the fiducials are now identified. The result is shown overlaid on the image in Fig. 6.

Fig. 6
figure 6

Image with fiducials ordered

Each individual fiducial has 24 possible permutations mapping the four vertices between the two systems. If the vertices are defined to be numbered counterclockwise about the individual fiducial’s centroid, similar to the fiducials themselves. The number of permutations is thereby reduced from 24 to 4. Next, an arbitrary permutation is chosen for each of the fiducials. The normal vector from the camera to each fiducial’s centroid is then computed using the camera model. Then, the solvePnP function is used on each fiducial individually. Since a random permutation of the four viable options was chosen, the rotation matrix is likely incorrect, and the translation vector is likely in the wrong coordinate system. However, the magnitude of the translation vector will be correct due to the fact that the fiducial is only incorrect by a rotation. The magnitude of each fiducial’s translation vector is then applied to the respective normal vectors. This corrects the rotation and gives a good estimate of the location of the three fiducials in three-dimensional camera coordinates. The Triad method is then applied to the three fiducials to obtain an estimate rotation matrix. Next, the rotation matrix is calculated for each of the four permutations for each of the fiducials individually. Whichever of the four permutations that results in the closest match to the rotation matrix calculated by the Triad method is the correct permutation. Once this step is complete, each vertex in the image is matched to its corresponding coordinate in model space, as shown in Fig. 7.

Fig. 7
figure 7

Image with fiducials and vertices ordered

At this point, solvePnP is used one final time to solve for all 12 coordinate pairs simultaneously. This gives the final, refined translation and rotation values. The model location and thereby the full nonlinear transform from image-sensor coordinates to model coordinates are now known.

5.3 Undistortion

Once the full transform is computed, the next step is to use it to undistort the image from the sensor. First, a uniform grid of desired span/chord points is chosen. For SWIFTER/SWIFTEST a grid at 0.5 mm pitch between 205 and 1029 mm in the streamwise direction, and \(\pm 203\,\hbox {mm}\) along the span (\(\pm 234\,\hbox {mm}\) measured parallel to the attachment line) is chosen. This grid pitch was chosen to be a round number slightly smaller than the finest physical pixel pitch on the model (\(\sim\)0.9 mm) in order to interpolate, rather than decimate, the sensor data. This grid of coordinates is then transformed using the camera model in order to obtain a grid of coordinates in image-sensor coordinates. The output image is then generated by retrieving the interpolated sensor value at each of the grid points. Lanczos resampling is used to perform the interpolation.

Figure 8 shows a representative uncropped, undistorted image from flight on the left, and the wind tunnel on the right. Even though the flight configuration has a relatively straight view of the model at optimal range, there is still a fairly visible amount of distortion. The wind tunnel corrections show a significantly more extreme case. Due to tunnel geometry constraints, the camera must be placed close to the model, at an oblique angle, with a wide-angle lens. This results in images that are difficult to orient by eye due to the heavy distortion, but makes little difference to the automated undistortion algorithm presented. In either case, the result is an accurate, distortion-free image.

Fig. 8
figure 8

Sample flight (top) and wind tunnel (bottom) undistorted images

6 Transition detection

Now that a clean image with uniform pixel spacing in span/chord coordinates has been obtained, a copy is made, which will be operated on for transition detection. First, the image is histogram-equalized. This operation linearizes the cumulative distribution function (CDF) of the histogram, which in effect provides a nonlinear, but still monotonic, mapping from percent shear stress into a normalized zero to one intensity value. The main benefit of this operation is that it strongly enhances small gradients in the image, making the transition front much easier to detect. Next, the result is median-filtered and low-passed at an approximately 2-mm cutoff wavelength. This is done to even further remove any residual sensor noise since the next step (Sobel gradient) is highly sensitive to high-frequency noise. Additionally, there is no worry of destroying useful small features in the image at this point due to the fact that the operations are currently being performed on a copy, and the only feature currently of interest is the transition front, which occurs on a much longer scale. The Sobel gradient of the image is then taken in both the x (streamwise) and z (spanwise) directions, giving an approximation of the vector gradient of the image. The result is then projected onto the characteristic transition propagation direction(s). For crossflow in the favorable pressure gradient present on SWIFTER, this is approximately \(\pm 10^\circ\). For Tollmien–Schlichting (T–S) dominated transition, \(90^\circ\) would be a pragmatic choice. Alternately, a polar function selecting for both \(90^\circ\) and \(\sim \pm 7^\circ\) (depending on pressure gradient) may be necessary in order to select for both T–S and for breakdown due to bypass or other phenomena. The easiest way to find the desired value is to retrieve a raw image, and measure the propagation angle of the phenomena of interest directly. This technique will function well as long as the measurement is within several degrees of the true characteristic angle(s) for the test conditions of interest. The purpose of this projection is to heavily select for any gradients related to the transition front and suppress all other gradients. In practice, for cases such as crossflow where there are two characteristics (positive and negative), the absolute value of the z gradient is taken, and the positive projection is then used, thereby selecting for both characteristics. The top image of Fig. 9 shows an example output of this process.

Fig. 9
figure 9

Transition detection at five selected iterations

The high-intensity values of this result are likely to be locations of transition. However, there is also a fair degree of noise and interference from features such as exhaust flare and sun reflecting off the model. Since these effects have been reduced by choosing an appropriate model coating, the highest intensity points will correspond to transition locations with high certainty. These points are marked as likely transition locations; Any point between the characteristic lines cast by turbulent wedges, in both the upstream and downstream directions, is reduced in intensity as transition is unlikely in those areas. The process is repeated for successively lower intensity values, until the transition front is fully defined. Figure 9 shows this process from top to bottom at five uniformly spaced iterations. For clarity, the chosen high-intensity points are shown highlighted in green. As the images progress from top to bottom, the penalizing of regions inside the characteristic lines of the previously selected points is visible as a significant reduction in intensity of those regions. It should be noted that if T–S is of interest, the characteristic lines cast by turbulent wedges should still be chosen in case there are some unintended trips in the flow. Once a full transition front is obtained, a mild median filter is applied to remove outliers where some feature other than the transition front was selected.

7 Human-readable image generation

At this point, processing of the images into a quantitatively usable state is complete. The images are spatially co-located to sub-pixel precision, allowing multiple images at similar conditions to be averaged together to remove sensor noise without compromising the clarity of the physical phenomena. Additionally, the magnitude of the images is equivalent to the local variation in shear stress. This allows for multiple useful measurements. One example is spatial FFTs with meaningful relative amplitudes that can be used to measure growth rates of crossflow. High- and low-speed streaks from isolated roughness are also clearly visible and can be measured as well.

Visually, the resulting images no longer exhibit a distinction between the laminar and turbulent regions due to the effective low-passing of the image during the conversion to local percent variation in shear stress. Instead, there is a sharp rise in the local variation in shear stress as the transition location is approached, followed by a sharp drop after transition occurs. While this is useful when processing the images quantitatively, an image with the DC component added back in is more visually intuitive. In order to generate an approximation of this DC component, the transition front is used to generate an image where the laminar region is zero and the turbulent region is a value chosen to correspond with the approximate difference in the shear stress between the laminar and turbulent regions. This image is then low-passed with the parameters previously used to high-pass the original image. The result is then added to the original images depicting the local variation in shear stress. The result is an image that still depicts an accurate variation in local shear stress, but also includes an approximation of the global shear-stress variation, making the image visually intuitive. These images have a dynamic range significantly higher than the human eye can perceive, let alone be represented by a 16-bit image. To map these values from their floating point representations to a 16-bit grayscale value, while visually preserving the high dynamic range, histogram equalization is utilized once again. This serves to enhance small variations in the image, such as crossflow streaking, while keeping the already strong gradients near transition acceptably strong.

Fig. 10
figure 10

Example single processed IR image (top) compared to averaged burst of 20 images (bottom)

Bursts of 20 processed images are then taken at each test point at constant conditions and averaged together. Since the camera and model are vibrating, the camera view will vary. Due to the undistortion process, any real features inherent to the model or flow remain spatially fixed. However, any features inherent to the camera sensor itself, such as static sensor noise, will move with the camera. This average therefore has the desirable effect of washing out the noise, but leaving the physical phenomena intact. Furthermore, since the images are co-located to sub-sensor-pixel-size precision, when the images are interpolated onto the finer grid, and then averaged together, features smaller than individual sensor pixels can be resolved. This technique is called superresolution and is also suitable for use on the non-human-readable images. Figure 10 compares a single processed IR image to the corresponding averaged burst of 20 images. In the unaveraged (top) image, there is a significant amount of high-frequency sensor noise apparent, which results in the image appearing much more muted. The high-frequency noise also somewhat obscures the crossflow streaking. In the averaged image (bottom), much of the noise is eliminated by the averaging, resulting in a much clearer image, including finer details of the crossflow streaking.

8 Statistical quantification

Once a transition front has been located as a function of span, one of the most basic measurements desired is a single representative transition location. This allows different conditions to be compared by a single number. The most obvious metric is a mean transition location. However, real test data frequently have imperfections such as bugstrikes that can contaminate part of the span. Since the mean considers the whole data set blindly, a single bugstrike or other similar contaminant will distort the measurement. Median is the next logical metric. It is more robust to these sorts of contaminants. However, it is still not immune.

A more robust metric is required. The probability density function (PDF) fills this requirement. The PDF can be considered analogous to the probability that transition has occurred near a streamwise location. The abscissa of the maxima of the PDF is therefore the most probable location of transition. In order to compute the PDF, the transition data are first sorted from most forward to most aft. The data are then uniformly numbered from zero to one in fractional increments. This value as a function of the transition location is called the cumulative distribution function (CDF). The derivative of the CDF is the PDF. Figure 11 shows the previously presented example IR image on bottom, with the corresponding PDF aligned in the top graph. The image is an averaged burst of 20 images. Transition detection is performed on each of the 20 images individually, and the individual PDFs are plotted as colored lines. Additionally, all of the transition data are combined into a single PDF, which is plotted in black. The abscissa of the maxima of the PDF is the dominant transition location. The presented values for transition location are taken from the aggregate PDF. The SEM (standard error of the mean) of the dominant transition location from the 20 individual PDFs is used, in conjunction with the estimated uncertainty in fiducial alignment, to compute the displayed uncertainty. It should be noted that this is simply a detection uncertainty; it does not take into account any transition location uncertainty due to variation in environmental factors such as angle of attack or Reynolds number.

In the event of a contaminated data set, the resulting wedge will either stretch one of the tails of the function, or even induce an additional peak. However, the dominant peak remains unmoved. As such, the PDF inherently rejects contamination with no user input or bias. Figure 12 shows a typical IR image with a bugstrike. Visually, the center of the wedge structure is at approximately 0.34 x/c. The PDF indicates that transition occurs at 0.348 x/c, whereas the mean transition location is 0.304 x/c and the median transition location is 0.323 x/c. As to be expected, the average is heavily moved forward by the bugstrike and the median is also dragged forward to a lesser extent. However, if the clean region of \(z = 0\) to \(-150\,\hbox {mm}\) is analyzed, the mean, median, and PDF all indicate a transition location of 0.346 x/c, which agrees well with the PDF of the entire front. As such, it is highly desirable to utilize PDFs for these comparisons due to their robustness to these environmental disturbances.

Fig. 11
figure 11

Processed IR image with PDF

Fig. 12
figure 12

PDF with bugstrike at \(z=75\,\hbox {mm}\)

9 Conclusions

IR thermography provides global transition data without disturbing the flow under study. It has two primary experimental requirements: a relatively high-emissivity surface that can hold a crisp temperature gradient, and a source of heating or cooling for the surface under study. Judicial choice of surface coating greatly improves the quality of data. However, IR thermography’s principle weakness has been the necessity of reducing the resulting images by hand. Reducing these data by hand is time intensive, requires a skilled eye, and is fairly subjective. The presented data reduction technique solves these issues. The result is a method that can not only quantitatively determine the global transition front over an entire model, but do so in a reliable and robust fashion; even when the model is subject to significant motion and real-world environmental contamination such as bugstrikes, this technique continues to provide a reliable transition location.

Furthermore, this technique is capable of examining finer details such as crossflow streaking. This opens up the possibility of quantitatively measuring the wavelengths present in the crossflow growth, including environments where more traditional techniques such as traversing hotwire anemometry are infeasible. Additionally, IR thermography provides a significantly faster global measurement than traversing a hotwire probe, and a more global picture than techniques such as hotfilm anemometry.