## 2 Details of the Proposed Approach
### Contour Detection Using the Oriented Gradient of Histograms  
The oriented gradient of histograms has been used in Pb contour detectors and has been a successful method in image processing literature [32; 33]. For most classical gradient-based contour operators, such as the Canny operator or the Sobel operator, the gradient at an investigated pixel is calculated based on the gray levels of its neighbor pixels. Different from these methods, the oriented gradient of histograms is calculated based on the intensity histogram of a local region, the center of which is the investigated pixel.  
The calculation of the oriented gradient for a pixel (\(m\), \(n\)) is demonstrated in Figure 2. A circular disc with a specified radius of \(r\) is placed at the investigated pixel and the disc is split into two half-discs by a diameter at angle \(\theta\), as shown in Figure 2a. According to the findings in the previous work [32], when calculating the intensity histogram of each half-disc region, the gray level range of [0, 255] is divided into 64 bins in our work; consequently, each bin covers an intensity range of about 4. Intensity histograms of the two half-disc regions were calculated and are shown in Figure 2b,c, respectively.  
Figure 1: Flowchart of the proposed approach.  
The horizontal ordinate in the intensity histogram describes the intensity bins, and values at the vertical ordinate are the count of pixels falling into each bin. The oriented gradient of histograms at pixel \((m,n)\) along direction \(\theta\) is then described by the \(\chi^{2}\) distance between the histogram \(g(i)\) of the lower half-disc and the histogram \(h(i)\) of the upper half-disc:  
\[G(m,n,\theta)=\frac{1}{2}\sum_{i=1}^{64}\frac{\left(g(i)-h(i)\right)^{2}}{g(i)+ h(i)} \tag{1}\]  
For each pixel, the oriented gradient is calculated along different directions by changing the value of \(\theta\). In our work, oriented gradients are calculated in nine directions by changing \(\theta\) from \(0^{\circ}\) to \(160^{\circ}\) with an interval of \(20^{\circ}\). Consequently, each pixel has nine oriented gradients corresponding to nine directions. Only the maximum gradient of each pixel remains as the prominent gradient, and a gradient map is then constructed with the prominent gradient of each pixel. In the constructed gradient map, pixels around which gray levels change greatly will have large gradient magnitudes. These pixels with large gradient magnitudes constitute the contours of the defects.  
Figure 2: Calculation of the oriented gradient of histograms: (**a**) a patch of a thermogram; (**b**) histogram of the upper half-disc; (**c**) histogram of the lower half-disc; (**d**) oriented gradient map after non-maximum suppression (NMS) processing.  
The contour of a defect presents a wide border in the constructed gradient map. In order to locate the defect contour more accurately, a narrow border with a pixel width should be determined. Non-maximum suppression (NMS) is implemented for this object in this work. NMS is a commonly used method in the image processing literature for obtaining a one-pixel width edge. Details of NMS are given in [35]. Through NMS operation, only the local maxima pixels remain in the gradient map. Consequently, the defect contour is indicated as a one-pixel width curve. The gradient map of Figure 2a after NMS processing is shown in Figure 2d.  
Once the contour of a subsurface defect is detected, the size of the defect can be easily calculated based on the ratio between the physical size and the pixel number in the thermogram of the tested specimen.