\section{Introduction}

The exact relationship between local brain motion and traumatic brain injury (TBI) remains  poorly understood in the large strain and high rate regime, due to challenges of obtaining direct and quantitative measurements of brain motion during injurious impacts. TBI studies typically rely on acceleration measurements of the head \citep{rimel1981disability,greenwald2008head,Campolettano2019,Canrillo2015,CRIPTON20141,Daniel2012,Rowson2013}. However, local brain motion under impact is far more complex than the overall rigid motion of the skull \citep{bayly2005deformation,sabet2008deformation}.

Since injurious experiments on the living brain are infeasible due to ethical reasons and since tissue degrades rapidly post mortem, TBI studies rely extensively on simulations \citep{Miller2019,zhang2004proposed,ji2014head,kleiven2007}. Finite element-based models for instance, use material properties of brain samples where the elastic constants vary by 2-3 orders of magnitude \citep{goriely2015mechanics}. There is therefore a need to directly measure local brain motion during impact to establish its link with injury and in order to validate TBI simulations.

Shear shocks are fundamentally unlike acoustical shocks which have been studied extensively including within the context of TBI. The soft tissue in the brain has nonlinear shear properties that are several orders of magnitude larger than its compressional properties. A typical Mach number (particle-velocity/wave-speed) for compressional waves in soft tissue is on the order of $10^{-4}$, and on the order of 1 for shear waves \citep{pinton2010nonlinear,Espindola2017}. This is due to the very low value of the shear wave velocity ($\sim$2 m/s) which in the case of a violent impact is comparable to the particle velocity (typically 2 m/s
or higher).  Consequently these extremely nonlinear shear waves can generate shock fronts within a single propagation wavelength i.e. areas with smooth shear waves can be adjacent to areas with violent shear shocks.

Impacts typically last $<$50~ms \citep{cloots2013multi} during which, the brain undergoes large, nonlinear strains (0.07) at high rates (40~$s^{-1}$) \citep{Sullivan2015}, which are challenging to sample adequately using conventional imaging techniques.
X-rays can image the human brain at high frame-rates of 10,000 images/s \citep{hardy2007study}, but are limited by poor soft tissue contrast. Although high speed MRI offers high resolution, the frame-rates are too low (167 images/s \citep{sabet2008deformation}) to sample brain motion at rapid time scales $<$50~ms. The high compressional wave speed of soft tissue compared to the shear wave speed enables the use of acoustics to measure shear wave dynamics in the brain \citep{catheline2003observation}.
Recently, using high frame-rate ultrasound (at 6000 images/s) and high sensitivity displacement tracking ($<$ 1~$\rm \mu m$ ), we showed that, planar shear waves propagated in the \textit{ex vivo} porcine brain develop into destructive shear shock waves \citep{Espindola2017}. The 35$g$ excitation wave developed into a 300$g$ shock-front deep inside the brain which we hypothesize, may be responsible for certain types of TBI. In the past studies, the imaging was performed on \textit{ex vivo} porcine brains embedded in gelatin, which presented a) fewer sources of image degradation as compared to the reflective skull surfaces around the \textit{in situ} brain and b) an indirect mechanism of transmitting shear waves into the brain. In this previous work, a flash focus imaging sequence was designed to improve tracking performance as compared to  plane-wave compounding approaches. Multiple focused ultrasound propagation events were used to image partial sections of the \textit{ex vivo} brain during each shear impact, to maintain a high frame-rate as well as high signal-to-noise ratio (SNR) \citep{espindola2017flashfocus}. However, this approach required multiple impacts (i.e. 16) to be generated in the brain, limiting our potential to investigate injury in the live brain in the realistic, single impact regime.

Here we develop an imaging method that a) measures shear shock wave formation in the challenging \textit{in situ} environment presented by the skull b) requires only one impact and c) uses a shear wave source that is directly attached to the brain. A high frame-rate ultrasound imaging sequence (2193 images/s) was designed to interleave multiple transmit-receive events to image different sections of the brain but collectively synchronize with a single mechanical shake i.e., only a single impact is necessary to quantitatively image brain motion rather than multiple impacts. A wide-beam emission was adopted for a long field of view and to optimize the number of ultrasound  events necessary to maintain a high SNR throughout the image width. The wide emission pulse was designed using a finite difference time domain-based ultrasound simulation tool (referred to as ``Fullwave") that we have previously developed \citep{pinton2009heterogeneous,pinton2011erratum,pinton2014spatial}, in conjunction with a time-reversal technique \citep{pinton2011effects,Fink1992,pernot2007vivo,soulioti2019super,pinton2012direct}.

Brain displacements were measured from the beamformed  radio-frequency (RF) data using an adaptive Normalized Cross Correlation (NCC)-based tracking algorithm that has been previously validated for shear shock wave motion \citep{pinton2014adaptive}. The experimentally detected shear shock waves were characterized in terms of particle velocity, acceleration, strain-rate, and amplification ratio and were compared to theoretical predictions of nonlinear shear wave propagation in a viscoelastic brain medium, using a Piecewise-Parabolic-Method based simulation tool \citep{Tripathi2019_PPM2D}.

\begin{figure}[t]
\centering
\includegraphics[width=0.659\textwidth,trim=0 0 0 0,clip]{figure1.pdf}
\end{figure}

\section{Methods}

Porcine heads were extracted 20 minutes post-mortem from 4 month old Yorkshire pigs (weighing $\sim$70~lbs). A craniectomy and durectomy was performed a) to directly attach the \textit{\textit{in situ}} brain to an external excitation and b) to provide access to the brain for ultrasound imaging. The head was secured in place using a stereotactic frame. A 3D-printed plate conforming to the brain surface (circular arc cross section, 28~mm radius) was adhered with cyanoacrylate to one hemisphere of the brain (Fig.1) and vibrated along the $z$ axis by coupling to an electromechanical shaker (Vibration Test Systems, Aurora, Ohio, USA), which propagated shear waves in the brain along $x$.  The ultrasound imaging was performed using a research scanner (Verasonics, Kirkland, Washington, USA) by placing a 7.8 MHz ultrasound probe (ATL L11-5) $\sim$6~mm above the brain with its lateral axis oriented along $x$, over a layer of ultrasound gel for acoustical coupling.

\subsection{Shear wave propagation in brain using a mechanical shaker}

To impart shear impacts to the brain, a 10 cycle, 40~Hz sinusoidal train was enveloped using a -80 dB Chebychev window, amplified and sent to the shaker. The shaker response was recorded using an attached linear accelerometer (PCB Piezotronics, Inc., Depew, NY, USA) (Fig. 2). Although other excitations including non-oscillatory ones are possible, this pulse was chosen to have control over the spectral content and to simplify the analysis of the shock wave physics.  A 40~Hz impact was chosen from the range of frequencies (10-100~Hz) in head impact responses measured using an instrumented mouth-guard in various contact sports \citep{laksari2015resonance,hernandez2015six}. Moreover at this frequency, the shear wavelength is 5.4~cm, allowing the shock formation distance to be well within the brain for accessibility to imaging.

\begin{figure}[h]
\vspace{20mm}
\setlength{\unitlength}{\textwidth}
\begin{picture}(0.2,0.17)(0.27,0.051)
\put(0.50,-0.){\includegraphics[width=0.67\textwidth,trim=0 0 0 5,clip]{figure2.pdf}}
\put(0.77,0.15){\scriptsize \color{black}{{\bf (a)}}}
\put(1.1,0.15){\scriptsize \color{black}{\bf {(b)}}}
\end{picture}
\label{fig:input}
\end{figure}

\subsection{Ultrasound imaging and motion tracking}
Brain motion was captured using a custom imaging sequence with 8 transmit-receive events focused at different lateral locations (Fig.3(a)).  The 8 events were interleaved to obtain one image of the brain in motion, i.e. firing ordering was at locations 1,5,2,6,3,7,4,8. Subsequently, the next set of 8 events generated the next image of the moving brain and in this manner, 4560 transmit events generated 570 images at 2193 images/s (movie included as supplementary material).
Each frame of the RF data was beamformed using a conventional delay-sum algorithm (Fig.3(a)), producing 256 B-mode lines (i.e. 32 lines were beamformed in a parallel-receive fashion per emission) \citep{espindola2017flashfocus,shattuck1984explososcan}. The interframe displacements between the beamformed images were tracked using an adaptive algorithm validated to measure shear shock wave displacements with an accuracy of 1$ \rm \mu m$ \citep{pinton2014adaptive,Espindola2017}. The algorithm maximizes correlation values by adapting the window size based on the displacement and further, rejects poorly correlated displacement estimates using a quality-weighted median filter.

\begin{figure}[t]
\setlength{\unitlength}{\textwidth}
\begin{picture}(0.172,0.57)(0.3,0)
\put(0.4,-0.0){\includegraphics[width=0.67\textwidth,trim=0 0 0 0,clip]{figure3.pdf}}
\put(0.36,0.48){\scriptsize \color{black}{\textcolor{white}{ \textbf{(a)}}}}
\put(0.087,0.29){\scriptsize \color{white}{\textbf{(b)}}}
\end{picture}
\end{figure}

As opposed to using a conventional focused emission, a custom wide-beam emission was designed with a consistent width of 5mm throughout its axial propagation (Fig.3(b)), in order to maintain a high SNR near the beam edges. The beam was designed using Fullwave ultrasound simulations by placing a 5mm wide, thin scatterer in the simulation field. By placing a receiver at the top of the acoustical simulation field, the RF data emitted from the scatterer was collected and time-reversed to produce the emission pulse used in the \textit{in situ} brain imaging sequence (Fig.3). The acoustical pressure field of the wide beam was measured in water at a transmit voltage of 4V, using a needle hydrophone with a frequency range of 1-10~MHz (HNA-0400, ONDA, Sunnyvale, California, USA), over a depth range of 5-30~mm using a sampling step-size of 500~$\rm \mu m$ laterally and axially. Fig.3(b) shows a single "pencil" beam (grid-size interpolated by a factor of 25) uniformly illuminating the area within its lateral extent, providing a long depth of field. Together, the eight beams illuminate the whole field of view (Fig.3(c)).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.69\textwidth,trim=0 0 0 0,clip]{figure4.pdf}
\label{fig:snapshots}
\end{figure}

\subsection{Validation using shear shock wave simulations in brain}
Experimental results were compared to theoretical predictions, by simulating the nonlinear brain deformation due to the propagation of shear shock waves, based on the formulation and numerical methods proposed by Tripathi \textit{et al.} \citep{Tripathi2019_PPM2D}.  Briefly, a system of nine first-order balance laws was proposed to describe the propagation of linearly-polarized nonlinear shear waves in a plane orthogonal to the axis of particle displacement in a homogeneous, isotropic, relaxing soft solid. The nonlinear propagation is governed by a cubic nonlinear term resulting from a fourth-order expansion of the strain energy density function, with a nonlinear parameter $\beta = 3(\mu + A/2 +D)/2\mu$ where $A$ and $D$ are third and fourth order elastic constants, respectively \citep{Landau1986}. Since the brain is highly attenuating and dispersing, a three-body generalized Maxwell body was used to model the relaxing media i.e. the attenuation power law was $\alpha(\omega) = a\omega^b$ and the dispersion was given by the Kramers-Kronig causality conditions (more details in \citep{Tripathi2019_PPM2D}).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.59\textwidth,trim=0 0 0 0,clip]{figure5.pdf}
\label{fig:velocity}
\end{figure}

The simulations were performed in a $95\times75$ mm$^2$ domain, using a circular arc (radius 26~mm, arc length 35~mm) as the source plate. The experimental data recorded near the source plate was used as the excitation pulse for the simulation, with an amplitude of 1~m/s at a frequency of 40~Hz. The spatial domain was discretized with $\Delta x = \Delta y = 130~\mu$m and $\Delta t = 23~\mu$s. The linear shear speed in the brain medium was defined as $c_T=2.10$~m/s at 75 Hz, with a density of 1000 kg/m$^3$, nonlinear parameter $\beta = 44.24$ and attenuation $0.06\omega^{1.05}$ Np/m \citep{Tripathi2020}.

\begin{figure}[t]
\centering
\includegraphics[width=0.59\textwidth,trim=0 0 0 0,clip]{figure6.pdf}
\label{fig:accln}
\end{figure}

\section{Results}
By dividing the tracked brain displacements by the interframe sampling period the local particle velocity was calculated (Fig.4(a)). The rightward propagation of the wave-front into the brain is clearly visible in the cross-section view in Fig.4(a).

Acceleration being a more relevant metric for brain injury, was calculated using a frequency-based derivative \citep{oppenheim1997signals} as, ${a(t)=\mathcal{F}^{-1}(2{\pi}i f \mathcal{F}(v(t))})$, where $a(t)$ is the acceleration time signal, ${i=\sqrt{-1}}$, $f$ is spectral frequency and $v(t)$ is the velocity time signal.
Numerical differentiation methods can be noise sensitive and it is critical that our acceleration estimates are accurate particularly at the shear shock-front. Noise in the experimental velocity was reduced prior to the derivative calculation using a low-pass Butterworth filter, with a cut-off frequency at the 13 \textsuperscript{th} harmonic (520~Hz for the shear wave central frequency of 40~Hz). This spectrum-based differentiation method was previously validated using a simulated known velocity signal with a similar noise level and sampling rate as experimental ultrasound-based velocity \citep{Espindola2017}. It was shown that the method preserved the shock-front acceleration peaks and underestimated the magnitude by $\sim$15\%.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.89\textwidth,trim=0 0 0 0,clip]{figure7.pdf}
\label{fig:sims}
\end{figure}

Snapshots of the particle acceleration movie (included as supplementary material) show that the acceleration amplitude reaches up to $>$110$g$ (Fig.4b) and that the development and decay of shock-front is a rapid event lasting $<$5ms and $<$15mm of propagation distance.

The spatio-temporal map of particle velocity is shown in Fig.5(a), by averaging the velocity estimates over a depth range of 13 to 15mm (Fig.4(a)). The shock-front is evident in Fig.5(a) as the sharp phase-change in the velocity profile at the two central positive cycles, up to $>$10mm from the plate source. From the velocity waveforms shown in Fig.5(b) at different distances from the plate source, the development of the shock-front is evident at 7.1~mm from the plate (orange plot in Fig.5(b) at the location of the orange marker in Fig.3(a)) from the sudden rise in slope at 0.11 s. In contrast, the shear wave at 3.3~mm from the source (blue plot, Fig.5(b)) has distorted much lesser from the smooth sinusoidal excitation imparted at the source. At farther distances of propagation, the velocity amplitude decreases due to the frequency-dependent attenuation within the medium. As a result, the slope at the shock-front wanes beyond 7.1 mm of propagation (yellow and purple plots, Fig.5(b)).

The spatio-temporal map of the particle acceleration (Fig.6(a), shows the finely localized peak in-between the peak and trough of the wave, i.e. at the shock-front. The corresponding  waveforms at different distances from the source are shown in Fig.6(b). Positive acceleration at the shock-front grows and reaches a peak of 102$g$ at 7.1mm from the plate, amplified by a factor of 4.0 with respect to the surface excitation of 25$g$. This shock wave formation distance is much smaller than the shear wavelength (5.4~cm) and the shock-front is highly localized not just in time, but also spatially as seen by the peak acceleration diminishing by a factor of 1.6 between 7.1~mm (orange plot, Fig.6(b)) and 12.5~mm (yellow plot, Fig.6(b)) from the source.

{Theoretical predictions from shear shock wave} simulations were compared to experimental results for four different distances of propagation = 5.2, 7.8, 10.4, 13.0 mm. An excellent match was seen between the experimental and simulated velocity at all the four locations (Fig.7(a)). The shark-fin shape of the cubically nonlinear waveform \citep{catheline2003observation} and the overall amplitude tracked closely. Correspondingly, the comparison between the experimental and simulated acceleration (Fig.7(b)) and strain-rate (Fig.7(c)) waveforms are shown. Both acceleration and strain-rate being rate-dependent metrics, show that the wave shape is far from monochromatic hence, demonstrating the physics of the nonlinear wave model. It is further seen that the simulations preserve the localized peaks in acceleration and strain-rate at the rise and fall of the shock-front (Fig.7(b) and (c)). The peak strain-rate at the shock-front reaches up to 250 $s^{-1}$ (Fig.7(c)). This exceeds by far previously reported brain strain-rates in concussed football players ($<100~  s^{-1}$) \citep{kleiven2007}.

\section{Discussion}
The shock wave formation in the brain is driven by high Mach numbers of $M=c_T/v_0=$~0.47 where, $c_T=$~2.14~m/s \citep{Espindola2017} is the shear wave speed and $v_0=$~1m/s is the amplitude particle shear velocity. The simplest nonlinear viscoelastic representation of brain motion during a linearly-polarized plane wave shear excitation can be described by a cubically nonlinear version of Burgers' equation~\citep{catheline2003observation,Espindola2017}. This type of motion has a characteristic odd harmonic signature wherein only the third, fifth, seventh, etc. multiples of the fundamental frequency are produced by the nonlinear propagation.
As the wave propagates, due to the cubic nonlinearity of the brain, the shear wave speed changes as a function of the local velocity amplitude \citep{Burgers1948}, resulting in an abrupt increase in slope at the peaks and troughs.

The acceleration amplification of the shear wave due to nonlinearity, is almost instantaneous in the brain i.e., at 3.3 mm (blue plot), a peak acceleration of
65$g$ was measured, while the initial condition given at the plate was just 25$g$ (Fig.6(b)). Experimentally the linear acceleration at the surface is not well captured owing to the noise in the displacement estimates at the plate boundary. However, the fully developed shock-front was clearly captured within the brain at 7.1~mm, close to the \textit{ex vivo}-measured shock formation of 6~mm, due to an 75~Hz impact. However, in the \textit{in situ} brain, a lower amplification factor of 4.0 was observed in the peak acceleration (as compared to a factor of 8.1 in the \textit{ex vivo} brain), for the same velocity amplitude of 1m/s. This disparity may be attributed to the difference in the input energy at the plate source between the two experimental scenarios. The attachment area between the plate source and the \textit{in situ} brain was $\rm \sim 6~cm^2$ in contrast to $\rm \sim 150~cm^2$ of embedded area of the plate inside the \textit{ex vivo} brain-gelatin medium. Besides larger contact area, anatomically, the spherical geometry of the skull is likely to have a strong focusing effect on the propagating shear shock waves in brain, which we are investigating using a human head gelatin phantom in a separate study \citep{tripathi2020shear}.

In Fig.6(b) it is seen that there is an asymmetry between the positive and negative cycles of the shear shock waveform. This is likely to be due to the shear wave reflections from the bone boundary underneath the brain. Furthermore, the particle motion is not strictly vertically polarized due to reflections from the boundaries of the small volume of the brain, which are unaccounted for in the shear shock wave simulations.

One of the principal guiding parameters for the design of the custom imaging sequence was measuring displacements at the maximum frequency contained in shear shock waves. The number of transmit events was chosen maintain a high frame-rate of 2193 images/s. Based on our previous experimental results in fresh porcine brain ~\citep{Espindola2017} shear wave propagation beyond frequencies of 1000 Hz are very rapidly attenuated and quickly fall below the sub-micron displacement sensitivity of ultrasound-based tracking algorithms~\citep{pinton2014adaptive}.
For the proposed impact velocity and acceleration the chosen imaging frequency is conservatively high. Nevertheless, at larger impact magnitudes in different intended scenarios, higher sampling frequencies may be necessary to measure shock waves. Injurious impacts in the National Football League, for example, are in the 98 $\pm$ 28$g$ range with impact velocities in the 9.3 $\pm$ 1.9~m/s range \citep{pellman2003concussion}. The average Olympic boxing hit to the head generates accelerations that are 58 $\pm$ 13$g$ \citep{walilko2005biomechanics}. In vehicular impact tests, translational head accelerations of 29 to 120$g$ have been measured \citep{zhang2006brain}. The 25$g$ brain surface accelerations used here are thus comparatively low. However, even in this moderate to mild regime, the accelerations deep inside the {\it in situ} porcine brain grow rapidly. Locally, at the thin shock-front, they exceed 100$g$. The shock-front accelerations are expected to increase for the realistic, stronger impacts. This is driven by the Mach numbers which are orders of magnitude larger than in fluids or in classical solids.  This behavior was also previously observed in the \textit{ex vivo} porcine brain which demonstrates that the skull, which contains and supports the brain, does not significantly alter this type of nonlinear viscoelastodynamic behavior occurring on short and rapid spatio-temporal scales.

\section{Conclusion}

It was shown that shear shock wave formation and propagation can be imaged \textit{in situ} in the porcine brain. Ultrasound imaging sequences were designed to reduce off-axis clutter and to record brain motion from a single impact. Time-reversed acoustical simulations using the Fullwave tool, were used to generate transmit profiles that were optimized to insonify a beamwidth of 1/8 of the imaging field of view.  Compared to previous shock wave imaging sequences, which were designed for a less challenging imaging environment, the proposed technique is better suited for \textit{in vivo} imaging because it only requires a single injurious impact, which is a more realistic model for TBI. Using quantitative motion estimation it was shown that an initially smooth monochromatic impact of 40~Hz frequency and amplitude of 25$g$ was amplified into a 102$g$ shear shock wave within 8~mm of propagation.
In the \textit{in vivo} brain, the availability of biological response markers to axonal injury  due to the occurrence of high accelerations at the shock-front, can establish the link between brain injury and the shear shock waves. These ultrasound displacement imaging methods can also be used to measure shear wave propagation from traumatic impacts in other vulnerable organs in the body.


\section{Introduction}
	In the electronics manufacturing industry, SMD assembly lines produce printed circuit boards (PCB) with high density at a very fast rate, which makes them prone to errors. Detecting solder joint errors in the early stages of the assembly is critical to reduce rework cost in production and scrap rates. However, this task is difficult due to the highly varied appearance of PCBs, given that even the electronic components used for the same purpose differ enormously in appearance. Furthermore, different types of solder joint errors can occur such as solder bridging, excessive/insufficient solder, pseudo solder, solder skips etc. Additionally, quality requirements standardized by the Institute for Printed Circuits should be met for high-quality production. In order to fulfil this need, Automated Optical Inspection (AOI) devices are used in assembly lines. These devices achieve very high recognition rates but require dedicated hardware. There are surveys \cite{moganti1996automatic,janoczki2013automatic} that elaborate on the assembly process and the importance of solder joint inspection in the manufacturing industry.

	With the advancement in deep learning, CNN, a powerful algorithm for feature extraction, is used for the inspection of solder joints. On the other hand, despite achieving high accuracy, supervised CNN based approaches require a lot of labelled defective (abnormal) samples for training, but there is a lack of abnormal samples although plenty of error-free (normal) samples are available in many datasets used for solder joint inspection. Moreover, abnormalities may occur due to various reasons, and take on various shapes and sizes which might be hard to capture in a data collection. Yet, to our knowledge, there is no open dataset that has enough abnormal samples, possibly due to confidential reasons.

	Therefore, in this paper, we propose to treat the solder joint defect detection as an anomaly detection problem, and resort to deep generative models, specifically to $\beta$-VAE \cite{Higgins2017betaVAELB}, which enable learning solely from normal samples.

	$\beta$-VAE is a variation of VAE \cite{kingma2013auto,rezende2014stochastic} that maximizes the probability of generating real data while balancing the likelihood and Kullback–Leibler (KL) divergence terms with an adjustable hyperparameter $\beta$.

	The $\beta$ term puts an emphasis on learning statistically independent latent factors which lead to disentangled representations where one single latent unit is sensitive to changes in a single generative factor. This disentanglement property of the $\beta$-VAE leads to more successful and more interpretable representations \cite{bengio2013representation}.

	Once the model describing the normal data is estimated via $\beta$-VAE, we generate an anomaly score for each test sample and apply a decision rule to detect anomalies. For generating the anomaly scores, we employ gradient-based representations as well as activation-based representations such as the reconstruction error.
	Finally, the model specifically learned to describe normal data is expected to yield a low anomaly score for normal data, and a high anomaly score for abnormal data. Hence, we detect these anomalies by comparing the resultant score against a threshold. With better-disentangled representations obtained using $\beta$-VAE, we show that anomalies on solder joints can be detected with high accuracy via a model trained directly on normal samples. Further, we compare this algorithm to two state-of-the-art approaches, namely the VAE and the Convolutional Autoencoder (CAE), and show the increase in detection rates.

	The main contributions of this study are the following:
	\begin{enumerate}[label=(\roman*)]
		\item  We attempt to solve a unique solder joint inspection problem in unrestricted domains, as such, it does not require specialized lighting, feature engineering, abnormal samples in training nor an error-free reference board.

		\item  We design a $\beta$-VAE architecture tailored for this problem, analyze it with different anomaly scoring techniques, and investigate the effect of $\beta$ for solder joint inspection. Further, we show the benefits of using $\beta$-VAE and show how the latent space representations of normal and abnormal solders become visibly separated.

		\item Our model can inspect both Integrated Circuit (IC) and non-IC component solder joints with a single architecture. To the best of our knowledge, there is no such study within the scope of anomaly detection.

	\end{enumerate}

	In the rest of this paper, we describe the related work in Sec. \ref{Sec:relWork}, present our proposed method in Sec. \ref{sec:Methods}, discuss our data set, parameter settings and experiments in Sec. \ref{sec:Experiment} and provide the results in Sec. \ref{sec:Results}.

	\section{Related Work} \label{Sec:relWork}

	The literature for solder joint inspection can be divided into two main groups, namely, the methods based on feature extraction and the methods based on deep learning.

	Feature extraction based methods can mostly either detect IC or non-IC solder joint errors since there are no discriminative features presented for both types of components. Besides, their successful applications mostly rely on a  circular or hemispherical tiered RGB illumination for inspection of both non-IC component solder joints \cite{kim1999visual,hongwei2011solder,wu2013classification,wu2008aoi,wu2014inspection,song2019smt}
	and IC component solder joints \cite{wu2011feature,ko2000solder}. On the other hand, there are some studies that do not require an RGB illuminator, but mostly an error-free reference image to inspect solder joints \cite{crispin2007automated, acciani2006application, ulger2019standalone, mar2011design}. The tiered illumination is used to generate colour contours on the solder paste, making use of the specular surface characteristics of the solder defects \cite{capson1988tiered}. These color contours are then analyzed via handcrafted features based on both color and shape to classify solder-joint errors. The most prominent color features are the average intensity and the percentage of highlight of each binarized RGB channel \cite{hongwei2011solder,  wu2013classification}.
	Geometric (shape) features in frequent use can be listed as solder and polarity areas \cite{wu2008aoi}, barycenter and distribution features \cite{wu2014inspection} as well as the template matching feature obtained by normalized cross correlation (NCC) \cite{crispin2007automated}. Other proposed features include Wavelet coefficients \cite{acciani2006application} and Gabor features \cite{mar2011design}. Additionally, feature selection mechanisms are employed to reduce the number of features and to improve the classification performance \cite{crispin2007automated,  song2019smt}.
	The classifiers employed in the more traditional literature include the Bayes classifier \cite{kim1999visual,wu2013classification}, multi-layer perceptron (MLP) \cite{ko2000solder, acciani2006application}, learning vector quantization (LVQ) \cite{wu2014inspection,ko2000solder} and Support Vector Machines (SVM) \cite{wu2013classification,dai2020soldering,song2019smt}; followed by the multi-stage classifiers \cite{kim1999visual,wu2013classification} to further improve accuracy. In practice, obtaining such reflections on a solder joint to extract these features is difficult; it requires special lighting and prior or manual input to parse the solder regions as shown in Fig. \ref{fig:nonIC_IC}. Additionally, traditional feature extraction based methods require determining distinctive features for each type of solder-joint errors.

	\begin{figure*}
		\begin{center}
			\begin{tabular}{c}
				\includegraphics[width=3.0cm]{Images/non_IC_parcellation.png} \hspace{1.5cm}
				\includegraphics[width=3.5cm]{Images/IC_parcellation.png}
				\\
				\hspace{-.5cm} (a) \hspace{4.25cm} (b)
			\end{tabular}
		\end{center}
		\caption
		{
			Surface Mount Device (SMD) parcellation: (a) non-IC component image and its solder pad and component body, (b) IC component image and its body, lead and solder pad regions.}
		\label{fig:nonIC_IC}
	\end{figure*}

	In recent years, there have been a few studies that employ deep learning based methods to inspect solder joints \cite{cai2018smt,dai2020soldering,mujeeb2019one}. Dai et al. \cite{dai2020soldering} employed the YOLO detector for solder localization and used semi-supervised learning for classification with SVM. Cai et al. \cite{cai2018smt} used cascaded CNNs: one learns the region of interest where solder defect is probable and the others are used to classify the solder joint of IC component.

	These studies employed both anomalous and normal data, and  still used the tiered RGB illumination. Since anomalous data is rare, it is difficult to gather a representative training set for a supervised CNN to learn to identify all types of solder joint errors. Additionally, commonly inspected PCB defects in the literature are good solder, insufficient solder, excessive solder, pseudo solder, missing solder, missing component and tombstone; however there are more different types of possible PCB defects such as shifted and misplaced component, solder balling, solder flags, pin holes, solder discoloration, billboarding, dewetting and the like. Therefore, an alternative way to handle solder error detection would be the anomaly detection perspective.

	For anomaly detection, Autoencoder \cite{zong2018deep} is a widely employed network that can learn to model complex data distributions. Autoencoders are trained directly on normal (error-free) data to minimize the reconstruction error such that, an encoding network, which strongly reduces the input to a compressed form, and a decoding network, which reconstructs its output to resemble the original input, are simultaneously trained. In the context of solder joint inspection, Mujeeb et al.\cite{mujeeb2019one} used a linear autoencoder on RGB illuminated solder joints and proposed using the $L2$ distance between the features of a reference image and a defective image as a similarity measure. However, the method requires a reference image that is a good (error-free) sample of the defective sample and tiered RGB illumination.

	Although Variational Autoencoders (VAE) are already being used for anomaly detection in different contexts such as dermatology \cite{lu2018anomaly} and medical imaging \cite{zimmerer2019unsupervised}; to the best of our knowledge, we are the first to use it for solder error detection and to offer a customized architecture for this problem. Also, our problem is unique, in that, we consider both IC and non-IC solder joints, we do not require any special lighting and neither do we require an error-free reference image. In doing so, we want to develop a more general solder-defect classifier that is not restricted to laboratory environments.

	\section{Methods} \label{sec:Methods}
	We briefly give an introduction to VAE and its variation, $\beta$-VAE. Latterly, scoring techniques that assign an anomaly score to each test data are presented.

	\subsection{VAE and $\beta$-VAE}
	VAE is a generative model that maximizes the marginal likelihood
	$p_{\theta}(\mathbf{X})$  where $\mathbf{X}$ is a data point (image) in a high dimensional space. It consists of a probabilistic encoder and a decoder where $\phi$ and $\theta$ are the encoder and decoder parameters respectively that include weights and biases, as shown in Fig. \ref{fig:VAEflow}. The encoder approximates the true posterior distribution $p_{\phi}(\mathbf{z}|\mathbf{X})$ whose distribution is assumed to be a Gaussian, where $\mathbf{z}$ is the latent (unobserved) vector of variables. The mean and standard deviation of the approximate posterior $q_{\phi}(\mathbf{z}|\mathbf{X})$ are the output of the encoder.

	We also let the prior distribution have a standard Gaussian distribution $p(\mathbf{z})= N(0, \mathbf{I})$. Then, $\mathbf{z}$ is sampled from the posterior distribution as $\mathbf{z} \raise.17ex\hbox{$\scriptstyle\mathtt{\sim}$} q_{\phi}(\mathbf{z}|\mathbf{X})= N(\boldsymbol{\mu}, \boldsymbol{\sigma})$.  The decoder is a generator that reconstructs $\tilde{\mathbf{X}}$ from $\mathbf{z}$.  However, we can not backpropagate the loss function to train the model since sampling operation is not differentiable. Therefore, a reparametrization trick is introduced. This trick uses a random variable $\boldsymbol{\epsilon}$ that is sampled from $N(0,I)$ and it is transformed to random variable z shifted by the mean and scaled by the standard deviation i.e. $\mathbf{z}= \boldsymbol{\mu} + \boldsymbol{\sigma} \odot \boldsymbol{\epsilon}$ that becomes a sample from the distribution $\mathbf{z} \raise.17ex\hbox{$\scriptstyle\mathtt{\sim}$} N(\boldsymbol{\mu}, \boldsymbol{\sigma})$.
	The illustration of VAE architecture is given in Fig. \ref{fig:VAEflow}.

	\begin{figure*}[ht]
		\centering
		\includegraphics[width=5.5in]{Images/VAE_Recons.png}
		\caption{Variational Autoencoder Architecture with isotropic Gaussian prior distribution over latent vector.}
		\label{fig:VAEflow}
	\end{figure*}

	To train VAE, the objective is finding the parameters $\theta$ and $\phi$ that maximize the evidence lower bound (ELBO) $\mathcal{L}$ of the marginal likelihood of a data point $p_{\theta}(\mathbf{X})$. The first term on the right-hand side of Eq. \ref{objective:VAE} is maximizing the probability of reconstructing the input $\mathbf{X}$ which corresponds to reconstruction loss. The second term is minimizing the KL divergence, which is a distance measure for probability distributions, between approximated posterior (encoder's distribution) and fixed Gaussian distribution with zero mean and unit variance. This leads to the lower bound on the probability density function of our data as is shown in Eq. \ref{objective:VAE}. We want to find the parameters (weights, biases) that maximize his objective function. Note that in the original VAE, $\beta=1$.

	\begin{figure*}
		\begin{equation}
			log p_{\theta}(\mathbf{X})\geq max_{\phi,\theta} \frac{1}{N} \sum_{n=1}^{N} log p_{\theta}(\mathbf{X_{n}} | \mathbf{z}) - \beta KL(q_{\phi}(\mathbf{z}|\mathbf{X}) || N(\mathbf{z}; 0, \mathbf{I}))
			\label{objective:VAE}
		\end{equation}
	\end{figure*}

	$\beta$-VAE introduces a weighting factor $\beta$ that balances reconstruction accuracy and latent channel capacity in Eq. \ref{objective:VAE}. With $\beta>1$ stronger pressure for the posterior $q_{\phi}(\mathbf{z}|\mathbf{X})$ to be close to the standard Gaussian prior $p(\mathbf{z})$ is introduced. This limits the capacity of latent $\mathbf{z}$ to learn a better disentanglement.  $\beta$-VAE with $\beta=1$ corresponds to the original VAE. Higgins et al. \cite{Higgins2017betaVAELB} demonstrated that $\beta$-VAE with the hyperparameter $\beta>1$ outperforms VAE on various datasets since it forces learning a more efficient latent representation of the data. This provides increased disentangling performance. However, there is a trade-off between reconstruction accuracy and learnt disentanglement quality within the learnt latent representations. A high $\beta$ value can lead to poorer reconstruction due to the loss of information while passing through the restricted latent channel. Therefore, it is important to find a balance that ensures a disentangled representation while still being able to reconstruct $\mathbf{X}$.
	On the other hand, $\beta$-VAE requires only tuning $\beta$ hyperparameter that can be estimated heuristically and does not suffer from the training instability of GAN-based approaches.

	\subsection{Anomaly Score}
	In testing, samples $\mathbf{X}$ at the test set are inputted to the encoder of the trained VAE. The encoder outputs a vector of means and standard deviations of each image. Reparametrization trick is applied to get the latent vector $\mathbf{z}$ which is sampled from Gaussian distribution with the mean and standard deviation computed by the encoder. Then the decoder reconstructs the input image from this distribution. The first term in Eq. \ref{objective:VAE} forms the reconstruction loss $(\mathcal{A(\mathbf{X})}_{Recon})$. It is computed as the anomaly score by calculating the mean squared error between the input image and reconstructed image for each test sample, as given in Eq. \ref{sRecons} where N is the mini-batch size. VAE trained on only normal samples is expected to be able to reconstruct normal samples with low reconstruction error and yield low anomaly score, whereas abnormal samples are reconstructed with high reconstruction error and yield high anomaly score.

	\begin{equation}
		\mathcal{A(\mathbf{X})}_{Recon}= log p(\mathbf{X}|\mathbf{z})=\frac{1}{N} \sum_{n=1}^{N} (\tilde{\mathbf{X}}_{n} - \mathbf{X_{n}})^2\\
		\label{sRecons}
	\end{equation}

	The objective is probability of occurrence of $\mathbf{X}$, hence anomaly score can also be calculated by $\mathcal{A(\mathbf{X})}_{ELBO}= -logp(\mathbf{X})$. Since the data is continuous and encoder's distribution is assumed to have a Gaussian distribution, KL divergence can be computed in closed form as the first term of the right hand-side of Eq. \ref{sELBO} as derived by Kingma and Welling \cite{kingma2013auto}:

	\begin{figure*}
		\begin{equation}
			\mathcal{A(\mathbf{X})}_{ELBO}=  \frac{1}{2}\sum_{j=1}^{J}(1+log(\sigma^2_{j})-\mu_{j}^2-\sigma_{j}^2) - \frac{1}{N}\sum_{n=1}^{N} (\tilde{\mathbf{X_{n}}} - \mathbf{X_{n}})^2\\
			\label{sELBO}
		\end{equation}
	\end{figure*}

	Kwon et al. \cite{kwon2020backpropagated} proposed training autoencoders with gradient constraint to model normal data distributions and using gradient-based representations for anomaly detection, motivated to capture information unavailable in the activation-based network representation. Anomalies require more model parameter updates to be represented compared to normal data, hence backpropagated gradients are used to characterize how much model update is required by input to detect anomalies.
	The gradients of each decoder layer with respect to model parameters, $\frac{\partial\mathcal{L}}{\partial\theta_i}$ are calculated by backpropagating the reconstruction loss. The algorithm combines reconstruction and gradient loss as an anomaly score is given by, \[\mathcal{A(\mathbf{X})}_{GradCon}= \mathcal{A(\mathbf{X})}_{Recons} + \gamma \mathcal{L}_{grad} \] where $\gamma$ is a scalar and $\mathcal{L}_{grad}$ is cosine similarity i.e. the cosine of the angle between average training gradients and gradient loss of current ($k_{th}$) iteration as given in Eq. \ref{eq:gradCon}.

	\begin{figure*}
		\begin{equation}
			\mathcal{L}_{grad}= - \mathbb{E}_i\Bigg[cosSIM(\frac{\partial{J^{k-1}}}{\partial\theta_{i_{avg}}}, \frac{\partial \mathcal{L}^k} {\partial \theta_{i}})\Bigg], \hspace{0.5cm}
			\frac{\partial{J^{k-1}}}{\partial\theta_{i_{avg}}}= \frac{1}{k-1} \sum_{t=1}^{k-1} \frac{\partial{J^t} }{\partial{\theta_i}}
			\label{eq:gradCon}
		\end{equation}
	\end{figure*}

	\noindent where $J$ is the loss function that is calculated during training and defined as, \[ J= \mathcal{L} + \alpha \mathcal{L}_{grad} \] where $\mathcal{L}$ is ELBO and $\alpha$ is a weight for the gradient loss. This anomaly detection algorithm using gradient constraint for training autoencoders and using a combination of reconstruction error and gradient loss as an anomaly score is called GradCon.

	\subsection{Decision Rule}
	After calculating the anomaly score, a decision rule is used to determine whether the sample is an anomaly or normal by comparing against a threshold \cite{Boracchi2020}. Anomaly score is expected to be low for the samples that the model is trained on i.e. the normal samples, whereas the model is not able to model the never seen abnormal samples which yield a high anomaly score. The samples with higher anomaly score than the threshold are evaluated as an anomaly and the ones below the threshold are normal samples. The threshold for the decision rule is considered as a hyperparameter and determined on the validation set. It is set to a value close to the anomaly score obtained on the validation set.

	\section{Experiment}\label{sec:Experiment}

	In this section, we describe our dataset and the preprocessing methods we employ. Then, we provide comparisons between VAE, $\beta$-VAE, and the CAE. The CAE is an unsupervised model that learns a low-dimensional representation of the input \cite{hinton2006reducing}, selected due to its characteristic to separate normal from abnormal data.

	\subsection{Dataset and Preprocessing}
	Our dataset has $2513$ normal training, $567$ normal validation and $150$ test samples ($66$ abnormal and $84$ normal). Some abnormal samples are obtained from the open dataset \cite{lu2020fics}. The abnormal (defective solder joints) samples in the test set include solder bridge, insufficient solder, excessive solder, solder flags, solder dewetting, pin holes, flux residues, missing and shifted components. Ground truth labels of these samples are only provided during testing to evaluate anomaly detection performance. The dataset is obtained by segmenting both IC and non-IC solder joints from solder pad regions of several PCBs as shown in red in Fig. \ref{fig:nonIC_IC}. For data preprocessing, our data is cropped from the short side to preserve the integrity of solder joint images and resized to $64$x$64$. Then, the data is normalized to zero mean, unit variance to have the range from $-1$ to $1$, and horizontal flipping is applied as a random transformation to avoid overfitting but without increasing the dataset size.

	\subsection{Architecture and Hyperparameters}

	Both VAE and $\beta$-VAE share the same architecture since the only difference is the added $\beta$ term in the objective function. These architectures were shown in Fig. \ref{fig:VAEflow} and in more detail with all the parameter sizes in Fig. \ref{fig:VAEarch}.  For comparison, the CAE has the same encoder-decoder structure except that there is a single convolutional layer at the bottleneck to compress the input. The encoder, which approximates the mean and standard deviation, has $5$ convolutional layers with an increasing number of filters from $3$ to $256$ and two separate convolutional layers at the bottleneck have $10$ filters. The decoder has transposed convolutional layers as a counterpart of the encoder. The convolutional layers have a filter size of $3$x$3$ and stride $1$. Convolutional layers are followed by batch normalization layers that normalize the feature map to zero mean and unit variance \cite{ioffe2015batch}. Leaky ReLU is applied for non-linearity to all the layers except for the output layer of the decoder, which uses Tanh. Max pooling layers are used to reduce the spatial dimension of feature maps by downsampling, while the upsampling layers increase the spatial size of the image by using bicubic interpolation. The VAE architecture is given in Fig. \ref{fig:VAEarch} with all the hyperparameters. As a convention, Conv$_{a,b,c}$ and Tconv$_{a,b,c}$ are the convolutional and the transposed convolutional layers, respectively, where $a$ is the number of channels in the input image, $b$ is the number of channels produced by the convolution, and $c$ is the size of the $c \times c $ convolution kernel. Max pooling$_{2,2}$ shows the max pooling filter size and stride, while upsample$_{2}$ provides bicubic interpolation sampling with a scale factor of $2$. The latent vector is reshaped to the channel size $10\times8\times8$ to generate the input image from latent vector $\mathbf{z}$.

	The model is trained with the ADAM optimizer \cite{kingma2017adam} for $100$ epochs with an initial learning rate of $1e-2$ that is decayed as the validation loss stops decreasing. To prevent overfitting to the training dataset, we resorted to $L2$ regularization that penalizes the high weights on features. The selected batch size is $64$. $\beta$ (KL weight) is selected as greater than one to learn better disentanglement of data, the effect of which will be more deeply analyzed in the next section. Lastly, the latent $\mathbf{z}$ dimension with the size of $640$ is selected since too small $\mathbf{z}$ can not model the marginal likelihood of a data point and too big $\mathbf{z}$ degrades the overall accuracy since learnt latent representation becomes indistinguishable for abnormal and normal samples. These hyperparameters given in Table \ref{hyperparameter1} are determined on the validation set.

	\begin{figure*}[ht]
		\centering
		\includegraphics[width=5.2in]{VAEarch.png}
		\caption{Sequence of blocks in the first row denote the encoder with $5$ convolutional (conv), $3$ max pooling layers and a convolutional bottleneck to obtain the parameters for the sampling of the latent vector $\mathbf{z}$.
			Blocks in the second row form the decoder that uses $6$ transposed convolution (Tconv) and $3$ upsampling layers to reconstruct the input image from $\mathbf{z}$. The terms bn and LRELU stand for batch normalization and leaky RELU, respectively.}
		\label{fig:VAEarch}
	\end{figure*}

	\begin{table}[ht]
		\caption{Hyperparameter Settings for $\beta$-VAE}
		\begin{center}
			\begin{tabular}{|c|c|}
				\hline
				\textbf{Hyperparameters}     & \textbf{Value} \\
				\hline
				Initial learning rate & $1e-2$   \\
				\hline
				Learning rate decay when loss plateous & $\gamma= 0.1$ \\
				\hline
				Number of Epochs     & $100$       \\
				\hline
				Optimizer     & Adam \\
				\hline
				Weight decay & $\lambda= 1e-4$ \\
				\hline
				Batch Size (N)     & $64$        \\
				\hline
				$\beta$	(KL weight)	   & $3$ \\
				\hline
				Latent space dimension & $640$ \\
				\hline
			\end{tabular}
			\label{hyperparameter1}
		\end{center}
	\end{table}

	\section{Results} \label{sec:Results}
	Original versus reconstructed normal and abnormal samples are given in Fig. \ref{fig:recons}.
	Normal samples yield a low anomaly score, whereas abnormal samples yield a high anomaly score. Precision, recall and F1-score are used as evaluation metrics. Recall rate shows what proportion of actual anomalies are detected correctly and precision is the rate for the correctness of these estimations. F1-score is calculated to show the overall accuracy of the model from these metrics. All the results are at least an average of $50$ runs.

	The generative models, VAE and $\beta$-VAE, are compared based on the evaluation metrics given in Table \ref{methods}; namely the reconstruction loss (Recon), ELBO and the GradCon anomaly score (the combination of reconstruction and gradient loss). The CAE model was also evaluated with reconstruction loss and GradCon as an anomaly score. The highest accuracy is obtained for $\beta$-VAE  for $\beta= 3$ with a score of $0.809$ with GradCon anomaly, outperforming CAE and VAE. The ELBO loss resulted in lower accuracy than reconstruction loss for VAE and $\beta$-VAE. GradCon made improvement over all the models with activation-based representations.

	\begin{table*}[ht]
		\caption{Comparison of reconstruction-based and generative models for different anomaly scores in terms of precision, recall and F1-score (average of 50 runs).}
		\begin{center}
			\begin{adjustbox}{width=\textwidth}
				\begin{tabular}{|c|c|c|c|c|c|c|c|c|}

					\hline
					\multirow{2}{*}{\diagbox[innerwidth=2.7cm, height=1.0cm]{\textbf{Evaluation}}{\textbf{Method}}} & \multicolumn{2}{|c|}{\textbf{CAE}}  & \multicolumn{3}{|c|}{\textbf{VAE}} & \multicolumn{3}{|c|}{\textbf{$\beta$-VAE}}    \\ [0.5ex]
					\cline{2-9}
					& Recon & GradCon & Recon & ELBO & GradCon & Recon & ELBO & GradCon    \\ [0.5ex]
					\cline{2-9}
					\hline
					\multicolumn{1}{|c|}{Precision} & $0.728 \pm 0.01$ & $0.708$ & $0.779 \pm 0.05$  & $0.805 \pm 0.03$ & $0.760 \pm 0.03$ & $0.783 \pm 0.04$ & $0.825 \pm 0.03$ & $0.804 \pm 0.04$ \\
					\hline
					\multicolumn{1}{|c|}{Recall} & $0.573 \pm 0.01$ & $0.780 \pm 0.02$ & $0.802 \pm 0.03$ & $0.773 \pm 0.04$ & $0.836 \pm 0.04$ & $0.827 \pm 0.03$ & $0.777 \pm 0.04$ & $0.816 \pm 0.05$ \\
					\hline
					\multicolumn{1}{|c|}{F1-score} & $0.641$ & $0.742 \pm 0.01$ & $0.789 \pm 0.03$ & $0.787 \pm 0.03$ & $0.796 \pm 0.03$ & $ 0.803 \pm 0.03$ & $0.799 \pm 0.03$ & $\bf{0.809} \pm 0.04$\\  [1ex]

					\hline
				\end{tabular}
				\label{methods}
			\end{adjustbox}
		\end{center}
	\end{table*}

	\begin{figure*}[ht]
		\centering
		\includegraphics[width=5.8in]{Images/GradCon_results4.png}
		\caption{Reconstructed normal and abnormal samples. The first row shows original input images with corresponding ground-truth labels and the bottom row shows the reconstructions along with predicted labels and anomaly scores.
		\label{fig:recons}
	\end{figure*}

	Disentangled representations emerge when the right balance is found between the reconstruction accuracy and latent capacity restriction \cite{Higgins2017betaVAELB}, hence the effect of different $\beta$ term is investigated on $\beta$-VAE that is trained with gradient constraint as shown in Table \ref{beta}. $\beta=3$ is selected for $\beta$-VAE, which yields the highest F1-score. Using larger or smaller $\beta$ values also led to lower accuracy while larger beta values resulted in poorer reconstruction accuracy.

	\begin{table*}[ht]
		\caption{The effect of hyperparameter $\beta$ on accuracy for $\beta$-VAE trained with gradient constraint.}
		\centering
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline

			\diagbox[innerwidth= 2.9cm, height=1.2cm]{\textbf{Evaluation}}{\textbf{Beta} ($\beta$)} & \bf{0.01} & \bf{0.1} & \bf{1 (VAE)} & \bf{3} &  \bf{10} \\ [0.5ex]
			\hline
			Precision & $0.704 \pm 0.01$ & $0.734 \pm 0.01$ & $0.760 \pm 0.03$ & $0.804 \pm 0.04$ & $0.720 \pm 0.05$ \\
			\hline
			Recall  & $0.786 \pm 0.03$ & $0.824 \pm 0.02$ & $0.836 \pm 0.04$ & $0.816 \pm 0.05$ & $0.754 \pm 0.05$ \\
			\hline
			F1-score   & $0.743 \pm 0.02$ & $0.776 \pm 0.01$ & $0.796 \pm 0.03$ & $\bf{0.809} \pm 0.04$ & $0.735 \pm 0.04$ \\ [1ex]
			\hline
		\end{tabular}
		\label{beta}
	\end{table*}

	In order to see the effect of different $\beta$ values on the latent space representations, t-SNE \cite{van2008visualizing}, which is a nonlinear dimensionality reduction algorithm that is useful for visualizing high dimensional data, is employed. Test images are encoded in the beta-VAE that is trained with gradient constraint (GradCon) and ELBO objective function for different $\beta$ values, then the resultant mean vector of the approximated posterior distribution is passed to t-SNE with the perplexity of $5$ since the dataset is small. Then t-SNE maps $640$-dimensional latent space into $2$-dimensional space. The embeddings are scaled to bring all the values to $[0,1]$. Latent space visualization of $\beta$-VAE that is trained with gradient constraint and ELBO objective is given in Fig. \ref{fig:tsne_Backprop} and \ref{fig:tsne_ELBO} respectively. t-SNE is run $100$ times and the embeddings with the lowest KL divergence is selected for each model. As shown in Fig. \ref{fig:tsne_Backprop} and \ref{fig:tsne_ELBO}, samples from the same class are closer to each other and appear in clusters. For $\beta$-VAE that is trained with gradient constraint, clusters of the classes are entangled i.e. not easily separable for the $\beta$ values of $1$ and $10$ as shown in Fig. \ref{fig:tsne_Backprop} (b) and (d). This behaviour is observed for the $\beta$ values of $0.1$ and $10$ for $\beta$-VAE trained with ELBO objective as given in Fig. \ref{fig:tsne_ELBO} (a) and (d). However, we can see that the classes are more separable for $\beta= 3$ although they are split into parts as given in Fig. \ref{fig:tsne_Backprop} and Fig. \ref{fig:tsne_ELBO} (c) for both models.

	For the real-world complex samples, it is difficult to achieve clear disentanglement or visualize it in two dimensions, but we can deduce that a better-disentangled representation is learnt for $\beta= 3$.

	\begin{figure*}
		\begin{center}
			\begin{tabular}{c}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/EmbeddingVAE_GradCon_checkpoint_beta01.pth_KL_loss_0.755.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/EmbeddingVAE_GradCon_checkpoint_beta1.pth_KL_loss_0.627.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/EmbeddingVAE_GradCon_checkpoint_beta3_KL_loss_0.747.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/EmbeddingVAE_GradCon_checkpoint_beta10.pth_KL_loss_0.544.png}
				\\
				\hspace{-1.5cm} (a) \hspace{3.5cm} (b) \hspace{3.5cm} (c) \hspace{3.5cm} (d) \hspace{-1.5cm}
			\end{tabular}
		\end{center}
		\caption
		{ \label{fig:tsne_Backprop}
			$2$-dimensional latent space representation of $\beta$-VAE that is trained with gradient constraint for (a) $\beta= 0.1$, (b) $\beta= 1$, (c) $\beta= 3$ and (d) $\beta= 10$ after applying t-SNE. Blue circles represent normal test data and the red circles represent abnormal test data.}
	\end{figure*}

	\begin{figure*}
		\begin{center}
			\begin{tabular}{c}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/Embeddingmodel_beta01_z10_Adam_ep100_KL_loss_0.886.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/Embeddingmodel_beta_1_z10_Adam_ep100_KL_loss_0.693.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/Embeddingmodel_beta3_z10_Adam_ep100_KL_loss_0.612.png}
				\includegraphics[height=3.0cm]{Images/GradCon_embeddings_150testdata/Embeddingmodel_beta10_z10_Adam_ep100_KL_loss_0.488.png}
				\\
				\hspace{-1.5cm} (a) \hspace{3.5cm} (b) \hspace{3.5cm} (c) \hspace{3.5cm} (d) \hspace{-1.5cm}
			\end{tabular}
		\end{center}
		\caption
		{ \label{fig:tsne_ELBO}
			$2$-dimensional latent space representation of $\beta$-VAE that is trained with ELBO objective function for (a) $\beta= 0.1$, (b) $\beta= 1$, (c) $\beta= 3$ and (d) $\beta= 10$ after applying t-SNE. Blue circles represent normal test data and the red circles represent abnormal test data.}
	\end{figure*}

	\section{Conclusions}
	In this study, solder joint defect detection in PCBs is considered as an anomaly detection problem. We propose using $\beta$-VAE, which is a generative model for disentangled representation learning, for anomaly detection to detect errors in both IC and non-IC component solder joints. beta-VAE disentangles factors of variation in data with that disentangled representations understanding and classifying data become easier. We compared different anomaly scoring techniques based on activation and gradient-based representations and the effect of different parameter beta on accuracy is investigated. The highest accuracy is achieved for the combination of reconstruction and gradient loss as an anomaly score to anomaly detection of solder joints.

	Despite the challenging environment of the problem due to the varying shape and small size of solder joints, we show that high accuracy can be achieved without using special lighting, hand-crafted features or labelled data for supervised learning.

	In the context of solder joint inspection where normal samples are in abundance and there are few abnormal samples, we show that the designed beta-VAE architecture is effective and leads to disentangled latent space representations.

	%

	\ifCLASSOPTIONcaptionsoff
	\newpage
	\fi

	\bibliographystyle{IEEEtran}