Defect identification of bare printed circuit boards based on Bayesian fusion of multi-scale features

mage preprocessing

The proposed model uses image enhancement and Gaussian filtering algorithms to preprocess the original image of the bare PCB to reduce the impact of lighting and noise on the image. Next, the PCB template images and PCB test images are registered using the Förstner corner detection algorithm (Naji et al., 2023). The pose correction and image registration of the image to be measured are completed through projection transformation. As shown in Fig. 1, the rectangular box represents the defect area obtained after calculating the difference between the template image and the test image.

The defect area of the image to be tested is then cropped on the bare PCB board. To reduce the impact of interfering pixels on feature extraction, the pixel value of the image area outside the defect area is set to 0. Some examples of defect images are shown in Fig. 2.

Feature extraction

Image recognition is actually a classification process. To identify the category to which an image belongs, the image must be distinguished from images of different categories. Feature extraction requires the selected features of an image to both describe the image well and distinguish it from different types of images. In order to increase the gap between feature information classes and improve the classification performance of feature vectors for bare PCB defect images, this article proposes a defect extraction algorithm based on multi-scale feature fusion.

Figure 2 shows that the cropped bare PCB defect image makes it difficult to extract rich feature information at a single resolution. Therefore, the proposed model performs down-sampling on the defect images. Based on the Gaussian pyramid, a multi-scale defect image is constructed. The image resolution is gradually reduced with the down-sampling, and the information is gradually enriched. Four down-sampling operations are performed on each image, 1–4 layers of pyramid images are constructed, and the features from each scale image are extracted, separately.

The technology roadmap of this model is shown in Fig. 3.

Multi-scale grey scale co-occurrence matrix feature extraction

The grey scale co-occurrence matrix features are extracted in the 0° direction of each layer pyramid image. This model uses the following statistics as features of the grey scale co-occurrence matrix:

Width represents the width of the co-occurrence matrix; c_ij represents the term in the co-occurrence matrix.

Angular second moment (Asm) refers to the sum of squares of the values of the elements in the grayscale co-occurrence matrix, which reflects the uniformity of the image grayscale distribution and the texture thickness: (1)${\displaystyle \mathrm{Asm}=\sum _{i,j=0}^{\text{width}}{c}_{ij}^2}$

Correlation (Cor) measures the similarity of elements in the spatial gray level co-occurrence matrix in the row or column direction: (2)${\displaystyle \mathrm{Cor}=\frac{\sum _{i,j=0}^{\text{width}}\left(i-u_x\right)\left(j-u_y\right)c_{ij}}{s_x{s}_y}}$

where, (3)${\displaystyle u_x=\sum _{i,j=0}^{\text{width}}i\cdot c_{ij}}$ (4)${\displaystyle u_y=\sum _{i,j=0}^{\text{width}}j\cdot c_{ij}}$ (5)${\displaystyle s_x^2=\sum _{i,j=0}^{\text{width}}{\left(i-u_x\right)}^2{c}_{ij}}$ (6)${\displaystyle s_y^2=\sum _{i,j=0}^{\text{width}}{\left(i-u_y\right)}^2{c}_{ij}}$

Homogeneity (Hom) measures the number of local changes in image texture: (7)${\displaystyle \mathrm{Hom}=\sum _{i,j=0}^{\text{width}}\frac{1}{1+{\left(i-j\right)}^2}{c}_{ij}}$

Contrast (Con) reflects the clarity of the image and the depth of the texture grooves: (8)${\displaystyle \mathrm{Con}=\sum _{i,j=0}^{\text{width}}{\left(i-j\right)}^2{c}_{ij}}$

The grey scale co-occurrence matrix features of the 1–4 layers of pyramid images are calculated separately. The multi-scale grey scale co-occurrence matrix feature vectors are obtained by concatenating the feature vectors.

Multi-scale directional projection feature extraction

Directional projection refers to the horizontal and vertical projection of regional grayscale values. The calculation formula is as follows: (9)${\displaystyle P_{\text{Vert}}\left(c\right)=\frac{1}{n\left(c+c^{\prime }\right)}\sum _{\left(r+r^{\prime },c+c^{\prime }\right)\in R}I\left(r+r^{\prime },c+c^{\prime }\right)}$ (10)${\displaystyle P_{\mathrm{Hor}}\left(r\right)=\frac{1}{n\left(r+r^{\prime }\right)}\sum _{\left(r+r^{\prime },c+c^{\prime }\right)\in R}I\left(r+r^{\prime },c+c^{\prime }\right)}$

where R refers to region and represents the target area; I refers to image, which represents the image where the region is located; $\left(r^{\prime },c^{\prime }\right)$ represents the top left corner of the bounding rectangle parallel to the smallest axis in the input area; $n\left(x\right)$ represents the number of region points in the corresponding row (r + r′) or column (c + c′). Horizontal projection returns a one-dimensional function that reflects changes in vertical grayscale values. Similarly, vertical projection returns a function that reflects changes in horizontal grayscale values. Horizontal and vertical projection functions for 1–4 layers of pyramid images are then drawn separately, using open circuit defects as an example, as shown in Fig. 4.

The directional projection function cannot effectively characterize the shape features of defects. The directional projection function is also directly input into the classifier as a feature vector, which may cause data redundancy. The known one-dimensional curve can be regarded as a one-dimensional signal in time-domain space. Therefore, this model extracts the time-domain features of the multi-scale directional projection functions of each defect image as feature vectors. The time-domain feature representation method is as follows:

Root mean square: (11)${\displaystyle X_{\mathrm{RMS}}=\sqrt{\frac{1}{N}\sum _{i=1}^N{g}_i^2}}$

Root amplitude: (12)${\displaystyle X_r={\left[\frac{1}{N}\sum _{i=1}^N\sqrt{\left|g_i\right|}\right]}^2}$

Skewness: (13)${\displaystyle X_{\mathrm{sk}}=\frac{1}{N}\sum _{i=1}^N{g}_i^3}$

Average amplitude: (14)${\displaystyle X_{\text{mean}}=\frac{1}{N}\sum _{i=1}^N\left|g_i\right|}$

Peak value: (15)${\displaystyle X_{\text{peak}}=max\left|g_i\right|}$

Kurtosis: (16)${\displaystyle K_u=\frac{\frac{1}{N}\sum _{i=1}^N{\left(g_i-\overline{g}\right)}^4}{X_{\mathrm{RMS}}^4}}$

Waveform factor: (17)${\displaystyle S_F=\frac{X_{\mathrm{RMS}}}{X_{\text{mean}}}}$

Pulse factor: (18)${\displaystyle I_F=\frac{X_{\text{peak}}}{X_{\text{mean}}}}$

Margin: (19)${\displaystyle L=\frac{X_{\text{peak}}}{X_r}}$

Crest factor: (20)${\displaystyle I_{\mathrm{CF}}=\frac{X_{\text{peak}}}{X_{\mathrm{RMS}}}}$

where N represents the number of pixel sequences, g_i represents the directional projection function value corresponding to the ith pixel sequence. The time-domain characteristics of the directional projection function of the 1–4 layers of pyramid images are then calculated and the time-domain features are concatenated to obtain the multi-scale directional projection feature vectors of the image.

Multi-scale gradient direction information entropy feature extraction

Image gradient

The image gradient represents the axis and partial derivative of the current pixel point to the axis, and can also be understood as the change speed of the pixel gray value. In mathematics, the Sobel mask operator is used to calculate the approximate gradient of changes in the horizontal and vertical directions of an image. The following common operators can be used, with A defined as the source image, and G_x and G_y defined as the approximate horizontal and vertical gradients of an image, respectively. The calculation method is as follows: (25)${\displaystyle G_x=\left[\begin{array}{ccc}-1\hfill & 0\hfill & 1\hfill \\ -2\hfill & 0\hfill & 2\hfill \\ -1\hfill & 0\hfill & 1\hfill \end{array}\right]\ast A,G_y=\left[\begin{array}{ccc}1\hfill & 2\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill \\ -1\hfill & -2\hfill & -1\hfill \end{array}\right]\ast A}$

where ∗ represents the convolution operation between the operator and the image.

The gradient direction feature is reflected in the image in the form of grayscale values. First, the original image is converted to a grayscale image, and then the pixels under the mask are calculated according to the formula. The obtained results replace the values of pixels at the center of the mask, and after calculating all pixels, the original image gradient direction feature map is obtained. One image was selected for each type of defect for experimentation, and the obtained gradient direction feature map is shown in Fig. 5.

This model calculates the grayscale information entropy of the gradient direction feature map of bare PCB defects based on the two-dimensional grayscale entropy calculation method of the image. Then, based on the multi-scale image construction method studied in the previous section, the gradient direction information entropy is calculated for the 1–4 layers of pyramid images. The obtained results are then serialized as multi-scale gradient direction information entropy features for bare PCB defects.

Deep semantic feature extraction based on transfer learning

To further explore the deep semantic features of bare PCB defect images, the model proposed in this article uses a convolutional neural network to extract features from the datasets. The VGG-16 network, proposed by Simonyan & Zisserman (2014) was selected as the pre-training model. The VGG-16 model consists of a convolutional layer, a pooling layer, and a fully connected layer. The network model structure is shown in Table 1.

Each convolutional layer applies a series of convolutional checks to the input data through convolutional operations for feature extraction, generating a feature mapping set. The maximum pooling layer aggregates input data through sliding windows to reduce feature mapping. Pooling operations can also improve the ability of features to describe images. The fully connected layer is used to generate feature vectors. After each layer, there is a nonlinear layer, ReLU, represented as f(x) = max(0, x), used to achieve rapid convergence of network training.

In the process of transfer learning, there are two models for fine-tuning. The first method is to freeze all convolutional feature extraction layers of the pre-trained model and only perform fine-tuning training operations on the classification layer; the second method is to fine tune all convolutional feature extraction and classification layers of the pre-trained model (Lamsiyah et al., 2023).

This model uses the first fine-tuning method. The known convolutional layer feature maps are output in the form of matrix sets. If directly converted into a one-dimensional vector, the dimensionality of the vector is very large, leading to overfitting problems and affecting classification efficiency (Bai & Tang, 2018). In order to reduce the risk of overfitting, the number of parameters and computational complexity are reduced. This model removes the last five layers of the VGGl6-Net network model and replaces them with a global average pooling layer and a fully connected layer with 512 nodes.

The function of global average pooling is to add all the pixel values of the feature map and find a balance to obtain a numerical value. By using this numerical value to represent the corresponding feature map, the spatial parameters are reduced, and the robustness of the network is improved (Hsiao et al., 2019). The schematic diagram is shown in Fig. 6.

The feature maps output from the global average pooling layer are used as inputs to the fully connected layer, and the 512-dimensional feature vectors outputted from the fully connected layer are used as the deep features of the bare PCB defect images. The structure of the adjusted VGGl6 net model is shown in Fig. 7.

Feature fusion

The different semantic features extracted from the same bare PCB defect image reflect different image information. This model integrates multi-scale grayscale co-occurrence matrix features, multi-scale directional projection features, multi-scale gradient directional information entropy features, and deep semantic features extracted by the VGG-l6 neural network to increase the diversity of feature information. If each group of feature vectors were simply added or subtracted, it would cause problems, such as high dimensionality of feature vectors and high computational complexity, so it is necessary to reduce the dimensionality of the feature vectors. It is also necessary to assign different weights to different feature vectors because different feature vectors have different scale ranges and classification performance, and direct addition and subtraction would ignore the imbalance of the data.

This model uses principal components analysis (PCA) to reduce the dimensionality of the extracted feature vectors (Jiang & Sun, 2023), and the high contribution features are extracted as effective features for bare PCB defect images.

Assuming there are M samples $\left\{X_1,X_{2,},\dots ,X_M\right\}$ andeach sample has an N-dimensional feature vector $X_i={\left(x_1^i,x_2^i,\dots ,x_N^i\right)}^T$, then each feature x_j has its own characteristic value.

The steps for the PCA are as follows:

Step 1: All features are centralized using demeaning.

Each feature is averaged and then its own mean is subtracted. (26)${\displaystyle \overline{x_n}=\frac{1}{M}\sum _{i=1}^M{x}_n^i}$

Step 2: The Covariance matrix C, is calculated taking features A and S as examples. (27)${\displaystyle C=\left[\begin{array}{cc}\mathrm{cov}\left(x_1,x_1\right)\hfill & \mathrm{cov}\left(x_1,x_2\right)\hfill \\ \mathrm{cov}\left(x_2,x_1\right)\hfill & \mathrm{cov}\left(x_2,x_2\right)\hfill \end{array}\right].}$

In the above matrix, the variances of features x₁ and x₂ are on the diagonal, and the covariance is found on the non-diagonal. A covariance greater than 0 indicates that if one of x₁ and x₂ x₂ increases, the other also increases; a covariance less than 0 means when one increases, the other decreases; when the covariance is 0, the two are independent. The larger the absolute value of covariance, the greater the impact the two have on each other. The solution formula for $\mathrm{cov}\left(x_1,x_1\right)$ is shown in Eq. (28), and the same applies to the others. (28)${\displaystyle \mathrm{cov}\left(x_1,x_1\right)=\frac{\sum _{i=1}^M\left(x_1^i-\overline{x_1}\right)\left(x_1^i-\overline{x_1}\right)}{M-1}}$

The covariance matrix C of M samples under N-dimensional features are obtained in this covariance calculation formula.

Step 3: The eigenvalues of the covariance matrix are and the corresponding eigenvectors are identified.

The eigenvalue λ of the Covariance matrix C and the corresponding eigenvector µ(each eigenvalue corresponds to a eigenvector) is calculated with the following formula: (29)${\displaystyle C\mu =\lambda \mu .}$

There are N eigenvalues λ, and each λ_i corresponds to a feature vector μ_i. The feature values are sorted in descending order, with the largest k selected, producing a set of $\left\{\left({\lambda }_1,{\mu }_1\right),\left({\lambda }_2,{\mu }_2\right),\dots ,\left({\lambda }_k,{\mu }_k\right)\right\}$.

Step 4: The original features are projected onto the selected feature vector to obtain the new K-dimensional features after dimensionality reduction.

Dimensionality reduction is performed by selecting and projecting the largest eigenvalues and corresponding eigenvectors. For each sample Xⁱ, the original feature is ${\left(x_1^i,x_2^i,\dots ,x_n^i\right)}^T$, and the projected new feature is ${\left(y_1^i,y_2^i,\dots ,y_k^i\right)}^T$. The new feature is calculated using Eq. (30): (30)${\displaystyle \left[\begin{array}{c}\begin{array}{c}{y}_1^i\hfill \\ y_2^i\hfill \end{array}\hfill \\ \cdot \hfill \\ \cdot \hfill \\ \begin{array}{c}\cdot \hfill \\ y_k^i\hfill \end{array}\hfill \end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{\mu }_1^T{\left(x_1^i,x_2^i,\dots ,x_n^i\right)}^T\hfill \\ {\mu }_2^T{\left(x_1^i,x_2^i,\dots ,x_n^i\right)}^T\hfill \end{array}\hfill \\ \cdot \hfill \end{array}\hfill \\ \cdot \hfill \\ \cdot \hfill \\ {\mu }_k^T{\left(x_1^i,x_2^i,\dots ,x_n^i\right)}^T\hfill \end{array}\right]}$

so that each sample Xⁱ is changed from the original $X^i={\left(x_1^i,x_2^i\right)}^T$ to the current $X^i=y_1^i$.

This model selects appropriate features by determining the amount of information contained in the feature values. According to Eq. (30), the new feature vectors calculated by principal component analysis are arranged in descending order. The size of the eigenvalues corresponds to the amount of information contained in the eigenvalues, which also indicates the contribution rate of the principal components. This model selects the feature values with a cumulative contribution rate of over 99.8% and ranking first as the fusion feature values. After the principal component calculation, the first six principal components of the multi-scale gray level co-occurrence matrix feature vector, the first seven principal components of the multi-scale directional projection feature vector, the first two principal components of the multi-scale gradient directional information entropy feature vector, and the first 68 principal components of the deep feature vector were used as the dimensionality reduced feature vectors for feature fusion.

The naive Bayesian theory is then used to calculate the weights of each group of feature vectors (Shitharth et al., 2022). The posterior probability of each group of eigenvectors is calculated through the dataset, and then their respective weight values are calculated. The dataset is first divided and then the conditional probability and prior probability of each category are calculated according to the training set. The posterior probability of the eigenvector is then calculated using the Bayesian formula. The maximum posterior probability is used as the final posterior probability of the eigenvector (Krishna & Divya, 2022).

The specific steps to calculate the posterior probability of eigenvectors and assign weights to each group of eigenvectors are as follows:

Step 1: Assuming there is a feature vector $x=\left\{x_1,x_2,\dots ,x_n\right\}$ and a category set C = c₁, c₂, …, c_k, then the posterior probability calculation formula of the Naive Bayes classifier is as follows: (31)${\displaystyle P\left(c_i|x\right)=\frac{P\left(c_i\right)P\left(x|c_i\right)}{P\left(x\right)}}$ where P(c_i) represents the prior probability of the category in the training set, P(x|c_i) represents the conditional probability of eigenvector x under category c_i, and P(x) represents the marginal probability of feature vector x. According to the naive Bayesian theory, denominator a is a constant in classification and can be omitted, so the posterior probability can be simplified as P(c_i|x) ∝ P(c_i)P(x|c_i).

Step 2: The prior probability P(c_i) and conditional probability P(x|c_i) are calculated. The prior probability P(c_i) represents the frequency of occurrence of category c_i in the training set, and the calculation method is as follows: (32)${\displaystyle P\left(c_i\right)=\frac{\left|D_{c_i}\right|}{\left|D\right|}}$ where |D_ci| represents the number of samples belonging to category c_i in the training set and |D| represents the total number of samples in the training set. The calculation of conditional probability P(x|c_i) is based on the naive Bayesian classifier assumption, also called the feature conditional independence assumption. According to this assumption, conditional probability can be defined as the product of the conditional probability of each feature, as shown in the following formula: (33)${\displaystyle P\left(x|c_i\right)=\prod _{j=1}^{n}P\left(x_j|c_i\right)}$ where x_j represents the jth feature of feature vector x and P(x_j|c_i) represents the conditional probability of the occurrence of feature x_j under category c_i. For discrete features, P(x_j|c_i) can be directly defined as the frequency of feature x_j appearing under category c_i, and the calculation formula is as follows: (34)${\displaystyle P\left(x_j|c_i\right)=\frac{\left|D_{c_i,x_j}\right|}{\left|D_{c_i}\right|}}$ where |D_ci,xj| represents the number of samples in the training set that belong to category c_i and feature x_j appears, and |D_ci| represents the number of samples belonging to category c_i in the training set. For continuous features, assuming that they conform to the Gaussian distribution, the conditional probability P(x_j|c_i) can be calculated using the probability density function of the Gaussian distribution, as shown in the following formula: (35)${\displaystyle P\left(x_j|c_i\right)=\frac{1}{\sqrt{2\pi {\sigma }_{c_i,j}^2}}exp\left(-\frac{{\left(x_j-{\mu }_{c_i,j}\right)}^2}{2{\sigma }_{c_i,j}^2}\right)}$ where μ_ci,j and σ_ci,j represent the mean and standard deviation of feature x_j under category c_i, respectively, as calculated from the training set.

Step 3: The posterior probability of the feature vector is calculated according to the naive Bayesian classifier, as follows: (36)${\displaystyle P\left(c_i|x\right)\propto P\left(c_i\right)\prod _{j=1}^{n}P\left(x_j|c_i\right)}$ where P(c_i) represents the prior probability of category c_i in the training set and P(x_j|c_i) represents the conditional probability of the occurrence of feature x_j under category c_i. Then the maximum a posteriori estimation is the final posterior probability of the eigenvector, and is calculated as follows: (37)${\displaystyle P_{\max }\left(c_i|x\right)={\mathrm{argmax}}_{c_i\in C}P\left(c_i\right)\prod _{j=1}^{n}P\left(x_j|c_i\right)}$

Step 4: The weight coefficients w_i of each group of feature vectors, defined as the normalization result of the maximum a posteriori estimation, are calculated as follows: (38)${\displaystyle w_i=\frac{P_{\max }\left(c_i|x\right)}{\sum _{i=1}^n{P}_{\max }\left(c_i|x\right)}}$

Due to the presence of four sets of feature vectors in this model, n = 4.

After feature dimensionality reduction and weight allocation, the feature vectors need to be connected together to complete feature fusion. This model uses parallel fusion as the feature vector stitching method. Assuming there is a bare PCB defect image sample I, after dimensionality reduction and weighting, $A=\left\{A_1,A_2,\dots ,A_n\right\}$ and $B=\left\{B_1,B_2,\dots ,B_n\right\}$ are obtained, and the feature vector after parallel fusion is $V=\left\{C_1,C_2,\dots ,C_n\right\}$, where C_i is defined as the following equation: (39)${\displaystyle C_i=\sqrt{{\left(A_i\right)}^2+{\left(B_i\right)}^2}}$

If the number in digits of feature vectors A and B is not equal, then 0 is added to the low dimensional feature vectors.

Before dimensionality reduction, the parallel fusion dimension of feature vectors had a total of 512 dimensions. After dimensionality reduction, the parallel fusion dimension of feature vectors was reduced to 68 dimensions. The feature dimensionality was reduced significantly, and the noise was eliminated. Different weight coefficients were assigned to each group of feature vectors based on different classification performance; the important features were strengthened and the non important features were weakened.

The pseudocode of the proposed algorithm is shown in Table 2.