A Discrete Optimal Transport Based Melding Defect Detection Method for PCB in Gas Meter

1. Introduction

Gas meters are indispensable measurement devices in modern urban infrastructure, widely used for quantifying gas consumption in residential and industrial environments [1,2]. With the rapid advancement of smart cities (where every single natural person in the city is managed by the government in order to provide intelligent service to all residents as a whole [3]), the role and significance of gas meters in energy management and intelligent monitoring have become increasingly prominent. The quality of gas meters directly affects the accuracy of gas measurement and the safety of users [4,5,6]. The control printed circuit board (PCB) of a gas meter, serving as its core component, is responsible for critical functions such as data processing, signal transmission, and control execution. The manufacturing quality of this PCB decisively influences the reliable operation of the gas meter. Therefore, ensuring the manufacturing quality of the gas meter’s control PCB is central and crucial to quality assurance in gas meter manufacturing [7,8].

The quality issues of gas meter control PCBs primarily arise from various defects during the manufacturing process, such as poor welding, component misplacement, and circuit breaks [9,10,11]. Among these, the detection precision of dynamic defects during manufacturing is the key factor determining the final product quality. Consequently, research focusing on enhancing the detection precision of dynamic defects in PCBs has become a hot spot in the field of defect detection.

Current research methodologies for PCB defect detection mainly include machine vision-based methods and deep learning algorithms [12,13,14]. Machine vision-based methods typically employ image processing techniques to identify surface defects on PCBs by comparing captured images with reference images. While these methods are straightforward and computationally efficient, they often lack the precision required to detect subtle defects, especially in complex or noisy backgrounds [15]. Moreover, they generally do not provide information that can guide the tracing of defects back to specific stages in the manufacturing process or inform effective repair strategies.

Deep learning approaches, such as those based on convolutional neural networks (CNNs), have been proposed to improve detection precision by learning complex patterns from large datasets [16,17]. These methods can achieve high accuracy in defect detection and classification. However, they require extensive training on large, labeled datasets, which can be time-consuming and costly to obtain. Additionally, deep learning models often act as black boxes, providing limited interpretability, which hinders their ability to guide defect tracing or repair in the manufacturing process.

In summary, existing research inadequately addresses the challenges of achieving high detection precision without the need for large training datasets and providing actionable insights for defect tracing and repair. Most methods focus solely on defect detection without offering the ability to analyze the discrepancies between the detected image and the desired image in a way that can inform corrective actions in the manufacturing process. This limitation hampers the ability to improve manufacturing processes based on the detection results.

To address these challenges, this paper proposes a melding defect detection method for PCBs in gas meters based on discrete optimal transport (as introduced in Section 2.3, the solution to a convex linear program introduced by Kantarovich [18]). The proposed method is an image comparison approach that contrasts the gathered images of the detecting areas with standard or desired images. Since solving discrete optimal transport problems can become computationally expensive, particularly when dealing with high-resolution images, our method mitigates this issue by focusing solely on the welding areas of the PCB. This localized analysis significantly reduces computational cost while maintaining precision. By utilizing discrete optimal transport theory, the method calculates the optimal mass flow between the detected image and the reference image, enabling precise identification of discrepancies in the welding areas. This high detection precision does not rely on large datasets for training, making it more practical for industrial applications where labeled data may be scarce. Furthermore, the calculated mass flow provides valuable information about the nature and extent of the defects, which can be used as a guide to trace back to potential issues in the manufacturing process and to inform subsequent repairs. This approach overcomes the limitations of existing methods by combining superior detection precision, minimal data requirements, relative low computational cost, and actionable insights for manufacturing process improvement and product repair, thereby offering a robust solution for quality assurance in gas meter manufacturing.

The specific structure of this paper is as follows: The literature review and related work are shown in Section 2. The problem definition and mathematical formulation are established in Section 3. The solver design to solve the optimal transport problem is presented in Section 4. The introduction for the welding defect detection method based on discrete optimal transport is described in Section 5. The experimental results and analysis are provided in Section 6. The conclusion to this paper is given in Section 7.

2. Literature Review and Related Work

2.1. Machine Vision-Based Methods

Machine vision has been widely used for PCB defect detection due to its ability to offer real-time monitoring and non-contact inspection at relatively low cost. Baygin et al. proposed a machine vision-based algorithm for detecting missing components and solder joint issues, demonstrating reliable detection in well-lit environments [14]. Similarly, Chang et al. applied machine vision techniques to detect soldering defects in PCBs, contributing to the faster and more cost-effective detection of welding problems [15].

Despite these benefits, machine vision approaches often struggle with precision, particularly when handling smaller or subtle defects. Liu and Qu found that in noisy environments, machine vision systems often fail to detect minor welding defects, leading to inaccuracies [19]. Ren et al. also noted that traditional machine vision systems face limitations in real-world conditions, such as variable lighting and background noise, which can compromise their detection accuracy [20]. Additionally, Baygin et al. highlighted that these systems provide little insight into the root causes of defects, making it difficult for manufacturers to trace issues back to their source [14].

Furthermore, Liao et al. emphasized that traditional machine vision methods often require manual tuning for different PCB designs, limiting their scalability in modern high-throughput production environments [21]. Despite their speed and cost advantages, these limitations make machine vision less suitable for environments requiring high precision in detecting subtle welding defects.

2.2. Deep Learning-Based Methods

Deep learning approaches, particularly convolutional neural networks (CNNs), have now become increasingly popular for PCB welding defect detection due to their ability to recognize complex patterns and achieve higher precision. Chen et al. demonstrated the effectiveness of deep learning-based models in detecting subtle defects in PCBs, achieving more than 95% accuracy [12]. Similarly, Legon et al. used CNN-based models to classify PCB defects, improving both detection accuracy and speed [22]. Liu et al. proposed a deep learning model that leveraged generative adversarial networks (GANs) to improve the resolution of PCB defect images, further enhancing detection performance [23].

However, deep learning models come with significant drawbacks. Ling et al. pointed out that these models are highly data-dependent, requiring large labeled datasets to perform effectively, which is a challenge in industrial environments where labeled defect data is scarce [24]. Pham et al. also emphasized the difficulties of acquiring enough training data, noting that deep learning models struggle when the dataset is small or unbalanced [25].

Moreover, the black-box nature of deep learning algorithms poses interpretability issues, as described by Sezer and Altan. This makes it difficult for manufacturers to understand the root causes of detected defects and to implement targeted repairs [10]. Although deep learning models significantly improve detection precision, their dependence on large datasets and lack of transparency make them challenging to deploy in many real-world manufacturing environments.

2.3. Optimal Transport

2.3.1. Monge’s Problem

The concept of optimal transport can be traced back to the period right before the French revolution. At that time, Monge proposed an intuitive idea of the optimal transport problem in a very practical reason: how to move a terrain from a given current landscape into another desired landscape with the least amount of total work done [26].

The initial formulation of the problem is as follows:

$$\begin{array}{ccccccc}\hfill \underset{T:X\to X}{inf}{\int }_{X}c(x,T\left(x\right))f\left(x\right)dx& \hfill & \hfill s.t.& \hfill & \hfill {\int }_{T^{-1}\left(B\right)}f\left(x\right)dx={\int }_{B}g\left(x\right)dx& \hfill & \hfill \forall B\subset X\end{array}$$

(1)

where $X\in {\mathbb{R}}^2$ and $f,g$ represents the current and desired landscape (as shown in Figure 1). Since the terrain does not change in size, $f,g$ must satisfy ${\int }_{X}f\left(x\right)dx={\int }_{X}g\left(x\right)dx$. $c(a,b)$ represent the work of transporting a to b is convex.

However, since the function $f,g$ are zero everywhere except in the landscape, the integral over X will lead to zero, thus the mass conservation constraint cannot be properly expressed. To solve this problem, replacing with measure instead is helpful. By using sum of weighted dirac masses to represent terrain, the problem can be solved. Then the Monge’s problem with measure is formulated as follows:

$$\begin{array}{ccccccc}\hfill \underset{T:X\to Y}{inf}{\int }_{X}c(x,T\left(x\right))df& \hfill & \hfill s.t.& \hfill & \hfill g\left(B\right)=f\left(T^{-1}\left(B\right)\right)& \hfill & \hfill \forall B\subset X\end{array}$$

(2)

where $X,Y$ are measurable, and $f,g$ are two measures on $X,Y$ respectively, and the cost function c is convex. Note that the constraint $g\left(B\right)=f\left(T^{-1}\left(B\right)\right)\forall B\subset X$ is the same as saying g is the push forward of f by T, $g=T\#f$

2.3.2. Kantorovitch’s Relaxed Problem

Nevertheless, the Monge’s problem with measure has a problem that the mass cannot be splited or merged as shown in Figure 2. To solve it, Kantarovich introduced the idea of marginal and relaxed the problem into the following form [27]:

$$\begin{array}{c}\hfill \underset{\gamma }{inf}\left\{{\int }_{X\times Y}c(x,y)d\gamma (x,y)|\gamma \in \Pi (f,g)\right\}\end{array}$$

(3)

where

$$\begin{array}{c}\hfill \Pi (f,g)=\{\gamma \in \mathcal{P}(X\times Y)|{\pi }_X\#\gamma =f,{\pi }_Y\#\gamma =g\}\end{array}$$

(4)

Here $\gamma \in \Pi$, a measure on the product of space $X,Y$ with marginal $f,g$, represents a transport plan between two measures $f,g$ on $X,Y$. And ${\pi }_X:X\times Y\to X$ is the projection map, so as ${\pi }_Y:X\times Y\to Y$. Note that the space of transport plans $\Pi (f,g)$ is convex.

In this case, in contrast with Figure 2, the Kantorovitch’s relaxed problem allows mass to be splitted, and we can have optimal transport transport as shown in Figure 3.

2.3.3. Discrete Optimal Transport in 1D with Finite Element

The discrete optimal transport problem in 1D with finite element can be easily driven from the Kantarovich’s relaxed problem as follows:

$$\begin{array}{ccccc}\hfill \underset{\Phi }{min}\{Flat\left(C\right)·Flat\left(\Phi \right)\}& \hfill & \hfill s.t.& \hfill & \hfill \left\{\begin{array}{c}{M}_X\Phi =F\hfill \\ M_Y\Phi =G\hfill \\ {\varphi }_{i,j}\ge 0\forall i,j\hfill \end{array}\right.\end{array}$$

(5)

where (i) $F={\left(f_i\right)}_{i=1,...,m}\in {\mathbb{R}}^m,G={\left(g_j\right)}_{j=1,...,n}\in {\mathbb{R}}^n$ are the vectors representing measures on $X,Y\subset S\subseteq \mathbb{Z}$. (ii) $\Phi ={\left(\varphi \right)}_{i,j}\in {\mathbb{R}}^{m\times n}$ represents the transport plan which transport ${\left(\varphi \right)}_{i,j}$ amount of mass from $f_i$ to $g_j$. (iii) $C={\left(c\right)}_{i,j}\in {\mathbb{R}}^{m\times n}$ represents the unit cost of transporting from $f_i$ to $g_j$. (iv) $M_X\in {\mathbb{R}}^{m\times mn},M_Y\in {\mathbb{R}}^{n\times nm}$ represents the constraints to satisfy marginal. (v) $Flat\left(\alpha \right)=\beta$ is some function that flatten matrix $\alpha \in {\mathbb{R}}^{m\times n}$ into vector $\beta \in {\mathbb{R}}^{k=m\times n}$. Note that flatten function is one-to-one.

For example, as shown in Figure 4, the sum of colum of $\varphi$ must equal to the mass above, and the sum of row of $\varphi$ must equal to the mass at the left. The C is the unit cost of each possible transport pair. The problem is to minimize the cost while satisfying the (1) equality of $\varphi$ and (2) mass above and on the left.

3. Problem Formulation

3.1. Core Problem

The core problem this research seeks to address is the accurate detection of welding defects in PCBs in gas meter while overcoming the limitations of traditional methods. To tackle this problem, we reformulate the PCB welding defect detection challenge as an image comparison task using optimal transport (OT) theory. Specifically:

• Input: Two local images—one representing the reference (defect-free) welding area on PCB and the other representing the welding area on PCB under inspection.
• Objective: Identify any discrepancies between the reference image and the test image, which may correspond to welding defects.
• Constraints:

  (a)

  Achieving high precision in detecting small, subtle welding defects.

  (b)

  Reducing dependency on large, labeled datasets for training detection algorithms

  (c)

  Providing interpretability of defect to enable actionable insights for defect tracing, process improvement, and product repair.

3.2. Formulation of Discrete Optimal Transport Problem

We suppose the 2 dimensional image of the detecting image and the reference image are both in the size of $d_1\times d_2$ pixels for some $d_1,d_2\in \mathbb{N}$, say

$$X,Y\subseteq \mathbb{Z}\times \mathbb{Z}$$

respectively. Since $X,Y$ are in same size, we can easily transform X to Y by elementary arithmetic, thus we assume $X=Y$. Similarly, we assume

$$X=Y=\{(x,y):x\in [0,d_1-1],y\in [0,d_2-1]\}$$

Consider the gray scale

$$S_{gray}=\{0,1,2,...,\kappa \}$$

for some $\kappa \in \mathbb{N}$ represents the brightness of the image, where 0 and $\kappa$ represent the darkest and brightest respectively.

We consider $f^{*},g^{*}$ to be densities on $X,Y$ such that for all $x\in X,y\in Y$, we have $f^{*}\left(x\right),g^{*}\left(y\right)\in S_{gray}$. By normalizing $f^{*},g^{*}$, we have the normalized image (density/measure) of the detecting image f and reference image g where

$$\begin{array}{ccc}\hfill f\left(x\right)={\displaystyle \frac{f^{*}\left(x\right)}{{\sum }_{\alpha \in X}{f}^{*}\left(\alpha \right)}}& \hfill & \hfill g\left(y\right)={\displaystyle \frac{g^{*}\left(y\right)}{{\sum }_{\beta \in Y}{g}^{*}\left(\beta \right)}}\end{array}$$

Then the optimal transport problem is as follows:

$$\begin{array}{ccc}\hfill C_{optimal}=min\left\{{\int }_{X\times Y}c(x,y)d\gamma (x,y)\right|\gamma \in \Pi (f,g)\}& \hfill & \hfill c={\parallel x-y\parallel }_2\end{array}$$

(6)

where the cost function c is the Euclidean distance. If $C_{optimal}$ exist, it is the optimal transport cost and the map $\gamma$ is the optimal transport map.

Table 1. Notations in Problem Formulation.

Notation	Description
$d_1\times d_2$	The size in pixel of the input image
X	The 2 dimensional coordinates for each pixel in detecting image
Y	The 2 dimensional coordinates for each pixel in reference image
$S_{gray}$	The set of gray scale
$\kappa$	Maximum of gray scale
$f^{*}$	The measure on X representing the gray sacle detecting image
$g^{*}$	The measure on Y representing the gray sacle reference image
f	The normalized measure of $f^{*}$ on X
g	The normalized measure of $g^{*}$ on Y
$C_{optimal}$	The optimal cost of transporting f to g

Table 1. Notations in Problem Formulation.

Notation	Description
$d_1\times d_2$	The size in pixel of the input image
X	The 2 dimensional coordinates for each pixel in detecting image
Y	The 2 dimensional coordinates for each pixel in reference image
$S_{gray}$	The set of gray scale
$\kappa$	Maximum of gray scale
$f^{*}$	The measure on X representing the gray sacle detecting image
$g^{*}$	The measure on Y representing the gray sacle reference image
f	The normalized measure of $f^{*}$ on X
g	The normalized measure of $g^{*}$ on Y
$C_{optimal}$	The optimal cost of transporting f to g

3.3. Solvability of the Assignment Problem

Clearly, any $f\left(x\right),g\left(y\right)$ can be considered as

$${\displaystyle \frac{n}{{\sigma }_X{\sigma }_Y}}={\displaystyle \frac{n}{\left({\sum }_{\alpha \in X}{f}^{*}\left(\alpha \right)\right)\left({\sum }_{\beta \in Y}{g}^{*}\left(\beta \right)\right)}}\in \mathbb{Q}$$

So both $f,g$ could be considered as some sum of Dirac masses in same scale ($1/{\sigma }_X{\sigma }_Y$). Note that

$$\sum _{x\in X}f\left(x\right)=\sum _{y\in Y}g\left(y\right)=1$$

which implies that the discrete optimal transport problem could be solved [18]. This could be cone by algorithm such as transportation simplex [28], Hungarian [29], auction [30], etc.

3.4. Properties of the Assignment Problem

As shown above, the assignment problem has some properties as follows:

• To solve the assignment problem does not need any outside information beside the detecting image and reference image.
• By solving the assignment, any slight difference between two image would be included in the transport plan and also represented by the optimal transport cost.
• The optimal transport plan is two-way since the cost is Euclidean, thus the plan could show that how can the detecting image be transported to the reference image and vice versa.

4. Solver Design

4.1. 2 Dimensional to 1 Dimensional

For further simplicity, we could transform the assignment problem in Section 3 to 1 dimensional form. We first consider a transform function $\mathcal{T}:X\to \mathbb{N}$ such that

$$\begin{array}{c}\hfill \mathcal{T}(x,y)=x+d_{2}y\end{array}$$

(7)

Clearly, $\mathcal{T}$ is bijective. Denote $k=d_1{d}_2-1$, and we let $\mathcal{D}=\mathcal{T}\left(X\right)=[0,k]\in \mathbb{N}$ to be the 1 dimensional frame for images of detecting image and reference image, and we let $\widehat{f},\widehat{g}$ to be the 1 dimensional representation of the images where $\widehat{f}\left(x\right)=f\left({\mathcal{T}}^{-1}\left(x\right)\right)$ and $\widehat{g}\left(y\right)=g\left({\mathcal{T}}^{-1}\left(y\right)\right)$. Then the assignment problem becomes

$$\begin{array}{ccccc}\hfill min\sum _{i,j\in \mathcal{D}}{x}_{i,j}{c}_{i,j}& \hfill & \hfill s.t.& \hfill & \hfill \left\{\begin{array}{cc}\sum _{j\in \mathcal{D}}{x}_{i,j}=\widehat{f}\left(i\right),\hfill & \forall i\in \mathcal{D}\hfill \\ \sum _{i\in \mathcal{D}}{x}_{i,j}=\widehat{g}\left(j\right),\hfill & \forall j\in \mathcal{D}\hfill \end{array}\right.\end{array}$$

(8)

where $x_{i,j}$ represents the mass transports from ${\mathcal{T}}^{-1}\left(i\right)$ to ${\mathcal{T}}^{-1}\left(j\right)$

$$\begin{array}{c}\hfill c_{i,j}={\parallel {\mathcal{T}}^{-1}\left(i\right)-{\mathcal{T}}^{-1}\left(j\right)\parallel }_2\end{array}$$

(9)

Table 2. Notations in Encoding and Decoding.

Notation	Description
$\mathcal{T}$	The bijective (2D to 1D) flatten function
k	The number of pixels minus one
$\mathcal{D}$	The 1 dimensional coordinates for the gray scale of both input images
$\widehat{f}$	The measure on $\mathcal{D}$ representing the gray scale detecting image
$\widehat{g}$	The measure on $\mathcal{D}$ representing the gray scale reference image
$c_{i,j}$	The unit cost for transporting from ith coordinate to jth coordinate in $\mathcal{D}$
$x_{i,j}$	The mass to transport from ith coordinate to jth coordinate in $\mathcal{D}$
${\sum }_{i,j\in \mathcal{D}}$	The cost of a transport
${\sum }_{j\in \mathcal{D}}{x}_{i,j}$	Mass conservation constraints for detecting image at ith coordinate in $\mathcal{D}$
$=\widehat{f}\left(i\right)$
${\sum }_{i\in \mathcal{D}}{x}_{i,j}$	Mass conservation constraints for reference image at jth coordinate in $\mathcal{D}$
$=\widehat{g}\left(j\right)$
$\mathcal{C}$	The cost map
$\mathcal{X}$	The transport plan
$\mathcal{A}$	The constraint matrix which project $\mathcal{X}$ to vector in ${\left[\begin{array}{ccc}\mathcal{D}& & \mathcal{D}\end{array}\right]}^T$
$\mathcal{B}$	The constraint vector represents the gray scale at each point for both image

Table 2. Notations in Encoding and Decoding.

Notation	Description
$\mathcal{T}$	The bijective (2D to 1D) flatten function
k	The number of pixels minus one
$\mathcal{D}$	The 1 dimensional coordinates for the gray scale of both input images
$\widehat{f}$	The measure on $\mathcal{D}$ representing the gray scale detecting image
$\widehat{g}$	The measure on $\mathcal{D}$ representing the gray scale reference image
$c_{i,j}$	The unit cost for transporting from ith coordinate to jth coordinate in $\mathcal{D}$
$x_{i,j}$	The mass to transport from ith coordinate to jth coordinate in $\mathcal{D}$
${\sum }_{i,j\in \mathcal{D}}$	The cost of a transport
${\sum }_{j\in \mathcal{D}}{x}_{i,j}$	Mass conservation constraints for detecting image at ith coordinate in $\mathcal{D}$
$=\widehat{f}\left(i\right)$
${\sum }_{i\in \mathcal{D}}{x}_{i,j}$	Mass conservation constraints for reference image at jth coordinate in $\mathcal{D}$
$=\widehat{g}\left(j\right)$
$\mathcal{C}$	The cost map
$\mathcal{X}$	The transport plan
$\mathcal{A}$	The constraint matrix which project $\mathcal{X}$ to vector in ${\left[\begin{array}{ccc}\mathcal{D}& & \mathcal{D}\end{array}\right]}^T$
$\mathcal{B}$	The constraint vector represents the gray scale at each point for both image

4.2. Assignment Problem to Linear Programming Problem

Since for any $i,j$, $c_{i,j}$ is fixed, we can consider

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathcal{C}=[c_{0,0},c_{0,1},...,c_{0,k},c_{1,0},...,c_{k,k}]\hfill \\ \mathcal{X}=[x_{0,0},x_{0,1},...,x_{0,k},x_{1,0},...,x_{k,k}]\hfill \end{array}\right.\end{array}$$

(10)

and

$$\begin{array}{ccc}\hfill \mathcal{A}=\left[\begin{array}{cccc}\mathcal{I}& 0& \cdots & 0\\ 0& \mathcal{I}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \mathcal{I}\\ {\mathbf{I}}_0& {\mathbf{I}}_0& \cdots & {\mathbf{I}}_0\\ {\mathbf{I}}_1& {\mathbf{I}}_1& \cdots & {\mathbf{I}}_1\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{I}}_k& {\mathbf{I}}_k& \cdots & {\mathbf{I}}_k\end{array}\right]& \hfill & \hfill \mathcal{B}=\left[\begin{array}{c}\widehat{f}\left(1\right)\\ \widehat{f}\left(2\right)\\ ⋮\\ \widehat{f}\left(k\right)\\ \widehat{g}\left(1\right)\\ \widehat{g}\left(2\right)\\ ⋮\\ \widehat{g}\left(k\right)\end{array}\right]\end{array}$$

(11)

such that

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathcal{I}=[1,1,...,1,1]\in {\mathbb{R}}^k\hfill \\ {\mathbf{I}}_0=[1,0,...,0,0]\in {\mathbb{R}}^k\hfill \\ {\mathbf{I}}_1=[0,1,...,0,0]\in {\mathbb{R}}^k\hfill \\ ⋮\hfill \\ {\mathbf{I}}_k=[0,0,...,0,1]\in {\mathbb{R}}^k\hfill \end{array}\right.\end{array}$$

(12)

Then we can reformulate the assignment problem as follows

$$\begin{array}{ccccc}\hfill \underset{\mathcal{X}\ge 0}{min}\mathcal{C}·\mathcal{X}& \hfill & \hfill s.t.& \hfill & \hfill \mathcal{A}\mathcal{X}=\mathcal{B}\end{array}$$

(13)

Note that the problem here is a standard linear programming (LP) problem [31], thus it can be solved by standard method such as dual simplex method [32], interior point method [33], etc.

5. Discrete Optimal Transport Based Welding Defect Detection Method

5.1. Method Overview

As shown in Figure 5, we gather two images from various ways and by various equipment. One of them is the detecting image and another is the reference image. Then we find the welding area for the component to detect, and cut both images into the same welding area.

To lower the time consumption, we then resize both images into some convenient size without losing too much precision, and we transform the image into gray scale image. Next, we use matrix to represent every pixel in gray scale image, and we normalize the matrix.

The matrix could be considered as measures over some domain $Dom\subset {\mathbb{Z}}^2$, and we could let Euclidean distance to be the transport cost. By doing so, we form an optimal transport problem. To solve the problem, we project the domain $Dom$ to one dimension, say $Do{m}_2\subset \mathbb{Z}$, by some bijection. Then we formulate the transport map $T_{map}$ in size of $Do{m}_2\times Do{m}_2$ representing how much mass is transport from detecting image to the reference image at each pixel. By such a variable matrix, we then generate a corresponding cost matrix C representing the unit cost for every transport. Flattening the map $T_{map}$ and cost C, taking dot product would result the cost of transport. Adding constraints of mass conservation, then the optimal transport problem becomes a simple minimizing problem subject to constraints.

By solving the minimization problem, we get the optimal transport map. Using the map, we can have the optimal transport cost, and optimal transport path from detecting image to reference image.

5.2. Image to Matrix

We assume that the input image is in same size as it is assumed in Section 3.2. We also assume that the input image can be formulated as pixels of some ratio (0-255) of red, green, and blue (RGB). Therefore, an image in RGB can be considered as a matrix where each element is a vector in ${\mathbb{N}}^3$. For example:

5.3. Find Welding Area

Since the components on PCB have distinct color in comparison to the color of PCB, we consider matrix $R,G,B$ having same size as the RGB matrix such that $R,G,B$ contain only R,G,B value in RGB matrix respectively. We define

$$\begin{array}{ccccc}\hfill \mathcal{N}(M,i,j,k)={\displaystyle \frac{1}{{(2k+1)}^2}}\sum _{q\in [-k,k]}{m}_{i+q,j+q}& \hfill & \hfill s.t.& \hfill & \hfill \left\{\begin{array}{c}i\pm k\in [0,m]\hfill \\ j\pm k\in [0,n]\hfill \end{array}\right.\end{array}$$

(14)

for $i,j,k\in \mathbb{N},M\in {\mathbb{R}}^{m\times n}$, which represents the average value for elements near the element at ith row jth column of matrix M. Then for the $R,G,B$, say they all have size $m\times n$, we consider $c=(c_1,c_2)=(\lfloor {\displaystyle \frac{m}{2}}\rfloor ,\lfloor {\displaystyle \frac{n}{2}}\rfloor )$ represents the relative center element of the matrix. Let

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{l}_1=\{(0,c_2),(1,c_2),...,(m,c_2)\}\hfill \\ l_2=\{(c_1,0),(c_1,1),...,(c_1,n)\}\hfill \\ l_3=\{(0,0),(1,\lfloor {\displaystyle \frac{1\times n}{m}}\rfloor ),(2,\lfloor {\displaystyle \frac{2\times n}{m}}\rfloor ),...,(m,\lfloor {\displaystyle \frac{m\times n}{m}}\rfloor )\}\hfill \\ l_4=\{(m,0),(m-1,n-\lfloor {\displaystyle \frac{1\times n}{m}}\rfloor ),(2,\lfloor {\displaystyle \frac{2\times n}{m}}\rfloor ),...,(0,n-\lfloor {\displaystyle \frac{m\times n}{m}}\rfloor )\}\hfill \end{array}\right.\end{array}$$

(15)

be four sets of elements in R (horizontal, vertical, and diagonal element w.r.t. c in R). Let some $k\in \mathbb{N}$ such that $k\le min\{{\displaystyle \frac{m}{8}},{\displaystyle \frac{n}{8}}\}$, and consider

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{v}_1\left(x\right)=(x+k,c_2),\hfill & x\in [0,m-2k]\hfill \\ v_2\left(x\right)=(c_1,x+k),\hfill & x\in [0,n-2k]\hfill \\ v_3\left(x\right)=(x+k,\lfloor {\displaystyle \frac{(x+k)\times n}{m}}\rfloor ),\hfill & x\in [0,m-2k]\hfill \\ v_4\left(x\right)=(m-(x+k),n-\lfloor {\displaystyle \frac{(m-(x+k\left)\right)\times n}{m}}\rfloor ),\hfill & x\in [0,m-2k]\hfill \end{array}\right.\end{array}$$

(16)

Define $V_1,V_2,V_3,V_4$ with same domain as $v_1,v_2,v_3,v_4$ respectively such that

$$\begin{array}{ccc}\hfill V_{\alpha }\left(x\right)=\mathcal{N}(R,i,j,k)& \hfill & \hfill (i,j)=v_{\alpha }\left(x\right)\end{array}$$

(17)

Let some $\delta ,\gamma ,\omega \in \mathbb{R}$ to be the threshold that represents the red value range for the component body and the red value for text on it or other color on it if exist. Consider $V_1^{\prime },V_2^{\prime },V_3^{\prime },V_4^{\prime }$ such that

$$\begin{array}{c}\hfill V_{\alpha }^{\prime }=\left\{\begin{array}{cc}0,\hfill & \delta \le V_{\alpha }\left(x\right)<\gamma \hfill \\ 0,\hfill & \omega \le V_{\alpha }\left(x\right)\hfill \\ 1,\hfill & otherwise\hfill \end{array}\right.\end{array}$$

(18)

Then for $V_1^{\prime },V_2^{\prime },V_3^{\prime },V_4^{\prime }$, let ${\lambda }_{11},{\lambda }_{12},...,{\lambda }_{41},{\lambda }_{42}$ to be the values farthest (but connected) from the center of domain which maps to 0 respectively. For rectangle component, we could have the welding area $A_1,...,A_8$ as follows:

(19)

And for component in other shape such as round, we could have the welding area $A_1,...,A_8$ as follows:

(20)

The welding area matrix $A_i$ is defined to be in same shape with R, but every element outside of the welding area is valued by $-1$.

Similarly, we take some $k,\delta ,\gamma ,\omega$, and then we could find two other sets of $A_1,...,A_8$ for green matrix and blue matrix $G,B$. Denoting the welding area of $R,G,B$ to be $A_1^R,...,A_8^R$, $A_1^G,...,A_8^G$, $A_1^B,...,A_8^B$ respectively. They we say $A_1^{*},...,A_8^{*}$ defined by

$$\begin{array}{cc}\hfill & {\left(a_{\alpha }^{*}\right)}_{i,j}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & =\left\{\begin{array}{cc}(R=-1,G=-1,B=-1),\hfill & ifmin\{{\left(a_{\alpha }^R\right)}_{i,j},{\left(a_{\alpha }^G\right)}_{i,j},{\left(a_{\alpha }^B\right)}_{i,j}\}=-1\hfill \\ (R={\left(a_{\alpha }^R\right)}_{i,j},G={\left(a_{\alpha }^G\right)}_{i,j},B={\left(a_{\alpha }^B\right)}_{i,j}),\hfill & otherwise\hfill \end{array}\right.\hfill \end{array}$$

(22)

to be the welding area.

5.4. Refine Welding Area

Now suppose we have $W_1,W_2,...$ to be the welding areas where $W_{\beta }$ is some combination of element in $A_1^{*},...,A_8^{*}$. Since the PCB is in green, we consider function $S_r(M,j),S_c(M,i)$ which take the average of elements with non-negative value in jth row and ith column of matrix M, such that the average is defined by $avg(x_1,...,x_n)=(R={\displaystyle \frac{{\sum }_{i\le n}{x}_{iR}}{n}},G={\displaystyle \frac{{\sum }_{i\le n}{x}_{iG}}{n}},B={\displaystyle \frac{{\sum }_{i\le n}{x}_{iB}}{n}})$ for $x_i=(R=x_{iR},G=x_{iG},B=x_{iB})$

For each welding area, we could change the value of every element in ith column to $(R=-1,G=-1,B=-1)$ if

$$\begin{array}{ccccc}\hfill S_c(R,i)=(R=S_1,G=S_2,B=S_3)& \hfill & \hfill satisfy& \hfill & \hfill \left\{\begin{array}{c}{S}_1\in Gree{n}_R\hfill \\ S_2\in Gree{n}_G\hfill \\ S_3\in Gree{n}_B\hfill \end{array}\right.\end{array}$$

for $Gree{n}_R,Gree{n}_G,Gree{n}_B$ is the color range of the PCB. Similarly, we change value at jth column to $(R=-1,G=-1,B=-1)$ if

$$\begin{array}{ccccc}\hfill S_r(R,i)=(R=S_1,G=S_2,B=S_3)& \hfill & \hfill satisfy& \hfill & \hfill \left\{\begin{array}{c}{S}_1\in Gree{n}_R\hfill \\ S_2\in Gree{n}_G\hfill \\ S_3\in Gree{n}_B\hfill \end{array}\right.\end{array}$$

Taking (5.6) as an example, the process would be as follows:

(23)

5.5. RGB to Normalized Gray Scale

We use the formula $Grayscale=0.299R+0.587G+0.114B$ as discussed in [34] to transform welding area matrix $W_1,...$ into gray scale matrix. For example, if we now let the example discussed in Section 5.2 to be the non-negative part of welding area matrix, we will have the following:

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}R=219,G=48,B=122& & R=103,G=12,B=77\\ R=66,G=200,B=187& & R=0,G=255,B=0\end{array}\right]\end{array}$$

where the element with $(R=-1,G=-1,B=-1)$ is transformed to $-1$

Next, we normalize the gray scale matrix where we let the sum of all non-negative element in the matrix to be 100 by letting the element at ith row jth column in normalized gray scale matrix to be

$$\begin{array}{c}\hfill m_{i,j}^{*}=\left\{\begin{array}{cc}-1,\hfill & m_{i,j}=-1\hfill \\ {\displaystyle \frac{100\times m_{i,j}}{{\sum }_{r,s}max\{m_{r,s},0\}}},\hfill & otherwise\hfill \end{array}\right.\end{array}$$

(24)

where $m_{i,j}$ is the element at ith row jth column in gray scale matrix. In this example, we will have the

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}23.2671& & 10.0838\\ 34.2724& & 32.3767\end{array}\right]\end{array}$$

(25)

5.6. Matrix to 2D Measure

While we have a normalized gray scale welding matrix, we can consider the matrix as a space in ${\mathbb{N}}^2$ where the gray scale value to be the measure on the space. Namely, let

$$\begin{array}{c}\hfill Z=[0,m-1]\times [0,n-1]\subseteq {\mathbb{N}}^2\end{array}$$

(26)

and

$$\begin{array}{ccccccc}\hfill h:Z\to {\mathbb{R}}_{+}& \hfill & \hfill s.t.& \hfill & \hfill h(i,j)=\left\{\begin{array}{cc}{m}_{i,j}^{*},\hfill & m_{i,j}^{*}\ge 0\hfill \\ 0,\hfill & otherwise\hfill \end{array}\right.& \hfill & \hfill \forall (i,j)\in X\end{array}$$

(27)

To eliminate the unnecessary part and keep only the welding area, we consider

$$\begin{array}{c}\hfill X^{*}=[a,b]\times [c,d]\subseteq Z\end{array}$$

(28)

to be the smallest rectangle to contain every element $z\in Z$ satisfy $h\left(z\right)>0$. Consider mapping $m_{origin}(x,y)=(x-a,y-b)$ that maps $X^{*}$ to $[0,b-a]\times [0,d-c]$. We could then use

$$\begin{array}{ccc}\hfill \Phi \left(x\right)=h\left(m_{origin}\left(x\right)\right)& \hfill & \hfill \forall x\in [0,b-a]\times [0,d-c]\end{array}$$

(29)

to represent the welding area.

With the same example as above subsection, we could have

$$\begin{array}{ccc}\hfill [0,b-a]\times [0,d-c]=[0,1]\times [0,1]\subset {\mathbb{N}}^2& \hfill & \hfill \left\{\begin{array}{c}\Phi (0,1)=23.2671\hfill \\ \Phi (1,1)=10.0838\hfill \\ \Phi (0,0)=34.2724\hfill \\ \Phi (1,0)=32.3767\hfill \end{array}\right.\end{array}$$

(30)

5.7. Resizing 2D Measure

In order to avoid high computational cost, we could resize the detecting welding area by considering

$$\begin{array}{c}\hfill X=[0,{\sigma }_1-1]\times [0,{\sigma }_2-1]\end{array}$$

(31)

and

$$\begin{array}{ccc}\hfill \left\{\begin{array}{c}{\mu }_1=\lfloor {\displaystyle \frac{b-a}{{\sigma }_1}}\rfloor \hfill \\ {\nu }_1={\displaystyle \frac{b-a}{{\sigma }_1}}-{\mu }_1\hfill \end{array}\right.& \hfill & \hfill \left\{\begin{array}{c}{\mu }_2=\lfloor {\displaystyle \frac{d-c}{{\sigma }_2}}\rfloor \hfill \\ {\nu }_2={\displaystyle \frac{d-c}{{\sigma }_2}}-{\mu }_2\hfill \end{array}\right.\end{array}$$

(32)

with

$$\begin{array}{c}\hfill f(x,y)=\left\{\begin{array}{cc}\Phi (x,y),\hfill & (x,y)=(0,0)\hfill \\ \Phi (x,y),\hfill & (x,y)=(0,{\sigma }_2-1)\hfill \\ \Phi (x,y),\hfill & (x,y)=({\sigma }_1-1,0)\hfill \\ \Phi (x,y),\hfill & (x,y)=({\sigma }_1-1,{\sigma }_2-1)\hfill \\ \phi (x,y),\hfill & otherwise\hfill \end{array}\right.\end{array}$$

(33)

for

$$\begin{array}{cc}\hfill \phi (x,y)=& {\displaystyle \frac{{\sum }_{i\le {\mu }_1,j\le {\mu }_2}\Phi (\lfloor ({\mu }_1+{\nu }_1)x\rfloor +i,\lfloor ({\mu }_2+{\nu }_2)y\rfloor +j)}{({\mu }_1+1)({\mu }_2+1)+{\mu }_1+{\mu }_2+1}}\hfill \\ \hfill & +{\displaystyle \frac{{\nu }_2{\sum }_{i\le {\mu }_1}\Phi (\lfloor ({\mu }_1+{\nu }_1)x\rfloor +i,\lfloor ({\mu }_2+{\nu }_2)y\rfloor +{\mu }_2+1)}{({\mu }_1+1)({\mu }_2+1)+{\mu }_1+{\mu }_2+1}}\hfill \\ \hfill & +{\displaystyle \frac{{\nu }_1{\sum }_{j\le {\mu }_2}\Phi (\lfloor ({\mu }_1+{\nu }_1)x\rfloor +{\mu }_1+1,\lfloor ({\mu }_2+{\nu }_2)y\rfloor +j)}{({\mu }_1+1)({\mu }_2+1)+{\mu }_1+{\mu }_2+1}}\hfill \\ \hfill & +{\displaystyle \frac{{\nu }_1{\nu }_2\Phi (\lfloor ({\mu }_1+{\nu }_1)x\rfloor +{\mu }_1+1,\lfloor ({\mu }_2+{\nu }_2)y\rfloor +{\mu }_2+1)}{({\mu }_1+1)({\mu }_2+1)+{\mu }_1+{\mu }_2+1}}\hfill \end{array}$$

(34)

which takes averages for neighbor measure.

Then we could represent the shrank or resized welding area’s measure by $(f,X)$.

5.8. 2D Measure to 1D Measure

As discussed in Section 4.1, since the element is finite, we can easily construct the flatten function $\mathcal{T}$. Applying $\mathcal{T}$ to the space X, we can get the 1D space $\mathcal{D}$. Assign element $d\in \mathcal{D}$ with measure f on ${\mathcal{T}}^{-1}\left(d\right)$, denoting as $\widehat{f}$. Then we have the 1D flat version ($\mathcal{D}$,$\widehat{f}$) from 2D version (X,f).

With the same example as above, we have $d_2=d_2=2$, which leads to

$$\begin{array}{ccccc}\hfill \mathcal{D}=[0,1,2,3]& \hfill & \hfill where& \hfill & \hfill \left\{\begin{array}{c}\widehat{f}\left(1\right)=23.2671\hfill \\ \widehat{f}\left(3\right)=10.0838\hfill \\ \widehat{f}\left(0\right)=34.2724\hfill \\ \widehat{f}\left(2\right)=32.3767\hfill \end{array}\right.\end{array}$$

(35)

5.9. 1D Measure to Discrete OT Problem

Followed by the previous steps, two input images would be transformed into two 1D measures on some $\mathcal{D}$. As discussed in Section 4.1, we then construct the 1D discrete optimal transport problem as formulated in (8) and (9).

If we inherit the example used in previous subsections to be the original image, and we arbitrarily let

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\widehat{g}\left(1\right)=10\hfill \\ \widehat{g}\left(3\right)=30\hfill \\ \widehat{g}\left(0\right)=20\hfill \\ \widehat{g}\left(2\right)=40\hfill \end{array}\right.\end{array}$$

(36)

to be the 1D measure representing the objective image. Then we will have some transport plans and transport cost as shown in Table 3 and Table 4.

Then the discrete optimal transport problem is to find the transport plan that minimizes ${\sum }_{i,j\le 3}{x}_{i,j}{c}_{i,j}$ satisfying the marginal, namely the sum of each column equals to the value above and the sum of each row equals to the value on the left.

5.10. Discrete OT Problem to LP Problem

As discussed in Section 4.2, we could construct the minimizing problem as in (13). Followed by the example, we will have

$$\begin{array}{ccc}\hfill \mathcal{A}=\left[\begin{array}{cccccccccc}1& 1& 1& 1& 0& 0& 0& 0& \cdots & 0\\ 0& 0& 0& 0& 1& 1& 1& 1& \cdots & 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \cdots & 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \cdots & 1\\ 1& 0& 0& 0& 1& 0& 0& 0& \cdots & 0\\ 0& 1& 0& 0& 0& 1& 0& 0& \cdots & 0\\ 0& 0& 1& 0& 0& 0& 1& 0& \cdots & 1\\ 0& 0& 0& 1& 0& 0& 0& 1& \cdots & 1\end{array}\right]\in {\mathbb{R}}^{8\times 16}& \hfill & \hfill \mathcal{B}=\left[\begin{array}{c}34.2724\\ 23.2671\\ 32.3767\\ 10.0838\\ 20\\ 10\\ 30\\ 40\end{array}\right]\in {\mathbb{R}}^8\end{array}$$

(37)

and

$$\begin{array}{cc}\hfill \mathcal{C}& =(0,1,1,\sqrt{2},1,0,\sqrt{2},1,1,\sqrt{2},0,1,\sqrt{2},1,1,0)\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill \mathcal{X}& =(x_{0,0},x_{0,1},x_{0,2},x_{0,3},...,x_{3,0},x_{3,1},x_{3,2},x_{3,3})\hfill \end{array}$$

(39)

5.11. LP Problem Solution to Optimal Cost and Optimal Plan

Solving the LP problem gives argmin, which corresponds to the optimal transport map. Suppose ${\mathcal{X}}^{*}$ is the solution to the LP problem, then ${\mathcal{X}}^{*}\mathcal{C}$ is the optimal transport cost. For $x_{i,j}\in {\mathcal{X}}^{*}$, moving $x_{i,j}$ amount of mass from ${\mathcal{T}}^{-1}\left(i\right)$ to ${\mathcal{T}}^{-1}\left(j\right)$ is the optimal transport plan.

For the same example, the solution to the constrained minimization problem is given by

$$\begin{array}{c}\hfill {\mathcal{X}}^{*}=(20,0,0,14.2724,0,10,0,13.2671,0,0,30,2.3767,0,0,0,10.0838)\end{array}$$

(40)

which can be reformulated as in Table 5.

Applying ${\mathcal{T}}^{-1}$ to the plan, we finally get the transport plan from measure f to measure g.

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}(0,0)\to (0,0)\hfill & 20\hfill \\ (0,0)\to (1,1)\hfill & 14.2724\hfill \\ (0,1)\to (0,1)\hfill & 10\hfill \\ (0,1)\to (1,1)\hfill & 13.2671\hfill \\ (1,0)\to (1,0)\hfill & 30\hfill \\ (1,0)\to (1,1)\hfill & 2.3767\hfill \\ (1,1)\to (1,1)\hfill & 10.0838\hfill \end{array}\right.\end{array}$$

(41)

5.12. Intermediate Images

For each assignment, we can consider moving path $p_{i,j}\left(t\right)$ where $p_{i,j}\left(0\right)={\mathcal{T}}^{-1}\left(i\right)$ is the initial position and $p_{i,j}\left(1\right)={\mathcal{T}}^{-1}\left(j\right)$ is the final position. Then for intermediate status $t^{\prime }\in (0,1)$, the mass is located at $p_{i,j}\left(t^{\prime }\right)$. For each assignment, finding the closest coordinate add the transporting mass $x_{i,j}$ to the position, then subtracting from the initial position would give a intermediate image representing the moving status at time $t=t^{\prime }$ according to the optimal transport plan.

6. Experimental Result and Discussion

The proposed Discrete Optimal Transport Based Welding Defect Detection Method is implemented in Python 3.11.4. The experimental configuration is Intel Core i9-13900K CPU @5.00 GHz with 128 GB memory and RTX 4070 Ti GPU. The test data are gathered from real PCB manufacturing line of intelligent gas meter from Qinchuan IoT Technology Co.,Ltd. In this section, we will show the following advantages of the proposed method from real welding status for components on PCB.

• (Precision) Sensitivity to slight change in image
• (Efficiency) Computational cost comparison for entire component and only welding part
• (Independency) Do not need large training data sets
• (Directivity) Can provide actionable insights for manufacturing process improvement and product repair

6.1. Image Gathering

For each PCB manufactured, we take a picture from the above. Since the position of PCB in image stay in some range, and also the component to detect stay in some range in PCB, we can easily get a relative small area containing the component (as shown in Figure 6).

6.2. Sample And Standard Image

Since for each component, there is a white rectangle-like box containing the component, we can then get the cleaner image of the detecting component by methods such as finding contours. Then we can compare the image (Figure 7) with the standard image (Figure 8). Note that the standard image could come from various origin such as previous qualified product, simulation, etc.

6.3. Find Welding Area

Once we obtained the detecting component, we first transform the RGB image into RGB matrix as discussed in Section 5.2. Since we know the detecting component is in rectangle shape, we then use the Welding Area Finding method discussed in Section 5.3 with

$$\begin{array}{ccc}\hfill k=2& \hfill & \hfill \left\{\begin{array}{c}{\delta }_R={\delta }_G={\delta }_B=30\hfill \\ {\gamma }_R={\gamma }_G={\gamma }_B=60\hfill \\ {\omega }_R={\omega }_G={\omega }_B=100\hfill \end{array}\right.\end{array}$$

(42)

where taking ${\lambda }_{11},{\lambda }_{12},{\lambda }_{21},{\lambda }_{22}$ to be the farthest value mapped to 0 by $V_1^{\prime },V_2^{\prime }$. Then we take $A_1^{*}\cup A_8^{*}\cup A_7^{*}$ and $A_3^{*}\cup A_4^{*}\cup A_5^{*}$ to be detecting welding areas (as shown in Figure 9, the welding areas are L and R). Then, with the same $\lambda$ and $A^{*}$, we also find the detecting welding area for the reference image as follows in Figure 10.

6.4. Defect Detection

Following the method discussed from Section 5.4 to 5.12, we could compare the right detecting welding area with the right welding area of the standard image, and for the left welding area, we compare with the left welding area of the standard image.

As shown in the Table 6, Table 7 and Table 8, the practical time complexity of solving optimal transport problem in the proposed method is approximately $\mathcal{O}\left(n^3\right)$. Thus by localizing the welding area, we significantly lower the computational cost with no harm in precision.

Also, the proposed method is extremely sensitive to slight defect. As we can see that the similar welding images in Figure 7 results in very different optimal transport cost. Any small difference in measure will reflect on the optimal transport plan since the mass change need to be transported out or filled in. This can be illustrated very clearly by Figure 11, Figure 12, Figure 13 and Figure 14. The difference in detecting welding area will be represented by different measure (or we could consider it as different mass distribution) and to find the optimal transport plan will need to deal with these differences. A small difference could result in different transport plan for a large neighborhood of such a difference.

This can be supported by the Table 9, where the result in optimal transport cost change is based on changing the gray scale value for some randomly chosen point(s) in the detecting welding area. Note that if the change exceed the min and max of 0 and 255, then the changed value will be limited in the range

Therefore, another advantage of only detecting the welding area is that it could avoid the error due to the pattern or text on detecting component. For example, in Figure 15 and Figure 16, the optimal transport plan need to consider the mass for the text printed on the component, even the text is unimportant.

Looking at the optimal cost for same set of samples with different resized size in Table 6 and Table 7, we could see that while the resolution of detecting welding area increases 4 times, the optimal cost approximately doubles. This is due to the Euclidean transport cost we chose. It also shows that reasonably resizing the size does not harm precision very much.

Note that the only input data we need is two images: the detecting image and reference image. During the detecting process, there is no outside data used. The welding area finding process, image resizing process, optimal transport problem generating process, etc. are all done locally. The independency on large training data set allow this method to omit the training process and to be applied in data scarce welding defect detecting scene.

Finally, the method not only compare two imagess’ difference, but also gives a path for one image to transform to another. As shown in Figure 11, Figure 12, Figure 13 and Figure 14, the method unveil the possible way for the detecting welding area to become ideal based on 2D image. It may also help to trace back to potential issues in manufacturing process. For example, if a batch of PCB all has problem on the same component, and the optimal transport path shows the similar mass moving trends, then it is very possible that the welding strategy of the welding machine has issue when dealing with this component. By considering the optimal transport path, the welding machine may add deviation to its strategy based on the path.

7. Conclusions

This study addressed the lack of welding defect detecting method that could maintain precision, be independent to training data, relatively cost low in computation, and guide the potential issue tracing in manufacturing process and the afterwards repair. An welding defect detecting method based on discrete optimal transport was proposed. First, the discrete optimal transport problem based on image comparison was established, and a solver transforming discrete optimal transport problem to linear programming problem was designed to solve the discrete optimal transport problem. Second, a method to process the input detecting image to matrix and then finding the welding area was proposed. In this method, images were transformed into RGB matrix, and by detecting the edges of component, the welding areas were derived from the matrix. Third, a way to transform matrix into 2D measure on ${\mathbb{N}}^2$ was proposed. The welding area is defined in measure and the non-welding area is eliminated from the domain. Fourth, a resizing model is proposed to resize the measure on ${\mathbb{N}}^2$ into smaller domain was proposed. It takes the weighted average over the neighborhood in order to summarize the information. Finally, with the resized 2D measure, the discrete optimal transport problem can be formulated and solved by the solver. According to the optimal transport plan, the transport path from the image of detecting welding area to the image of standard welding area could be generated, as well as the intermediate status.

However, there also exist some limitations for this method. Although the time complexity in practice is not bad, the time consumption remains high when the size of image is large since the variable is the square of total pixels of the image. Considering continuous transport or semi-continuous transport may benefits in some cases. Another limitation is that the method proposed in this paper only analyzes a 2 dimensional image. It is possible that the images of same object are different due to the different position of camera. Using multiple images to form a 3 dimensional density and finding the optimal transport may be the way to solve it.