Source Design Optimization for Depth Image Reconstruction in X-ray Imaging

1. Introduction

Non-destructive testing plays a key role in industrial quality control [1]. X-ray imaging in particular is a highly effective non-destructive testing technique that can deliver detailed 3D images of the internal structure of products. It is therefore widely used in various industries for quality control [2,3,4]. Conventional tomographic acquisition requires a full rotation around the object, limiting its applicability to off-line inspection of small to medium-sized objects. To overcome this limitation, several approaches to tomographic reconstruction from limited-angle data have been developed. These methods typically assume additional prior knowledge of the object to fill in information that is missing in the measurements. A prime example is discrete tomography, where objects consisting of a few distinct materials can be reconstructed from severely under-sampled data. In some applications, objects are wide and flat, and tomographic data can only be collected for a very limited angular range. This pushes the limits of conventional limited-angle tomography and we have to resort to another class of methods [5]. Laminography offers a method to reconstruct an image of a specific slice of a flat layered object. As shown in Figure 1, it needs a flexible source and detector that can move and collect data from several directions such that during scanning the X-ray tube is located on the same side of the flat object and the detector is on the opposite side of it [6].

In this paper we consider a particular setup for laminographic imaging which requires the acquisition of X-ray projection images from N source positions on a fixed grid, while the detector is fixed on the opposite side of the object. Our aim is to be able to optimally design acquisition protocols with $K\le N$ sources for a given set of objects. By optimally designing the setup we can reduce imaging artifacts as well as potentially speed up the data-acqisition process in cases where moving the source is time consuming (e.g., robotic arm imaging). We pose the resulting optimal experimental design (OED) problem [7,8] as a bi-level optimization problem that employs derivative-based methods to facilitate implementation in high-dimensional settings [9].

Our main contributions are the following.

• We analyze the corresponding 3D and 2D (slice-based) image reconstruction process, present an analog of the Fourier Slice Theorem for this setup and analyze the image artifacts present in slice-images.
• We pose an experimental design problem for selecting the K most informative source positions for laminography of a class of objects, and formulate a bi-level optimization in the Bayes risk minimization framework to solve it.

The remainder of the paper is organized as follows. In Section 2, we present the forward and backward operators in the continuous domain for the depth imaging procedure. We study depth reconstruction in the discrete domain in Section 3. In Section 4, we propose a bi-level optimization problem to find a source design with K sources chosen from a fixed array of $N\ge K$ sources. In Section 5, we propose an algorithm to solve the proposed bi-level optimization problem. In Section 6, we show the results obtained by our method and compare them to optimal source designs, obtained by exhaustive search. In Section 8, we provide a discussion on the limitations and future work of the proposed method.

1.1. Notation

Throughout the paper, scalars and vector elements are denoted by lower-case letters, vectors by lower-case boldface letters, and matrices by upper-case boldface letters. The notation ${·}^{\top }$ stands for transpose of either vector or matrix. Moreover, # denotes cardinality of a vector and $\parallel \mathit{x}\parallel =\sqrt{{\sum }_{i=1}^n{x}_i^2}$ denotes the ${\ell }_2$-norm of a vector.

2. Forward and Backward Operators in the Continuous Domain

In this section, we will first review the linear X-ray projection model, followed by introducing the forward and backward operators for the setup shown in Figure 2.

2.1. Linear Projection Model

During the X-ray imaging process, an object is exposed to a beam of X-ray photons emitted by an X-ray tube. Some of these photons are absorbed by the object, while others reach the detector. The captured measurements can be modeled using the Beer-Lambert law, which describes the exponential relationship between the intensity of the X-ray beam and the distance it travels through the object

$$\begin{array}{c}\hfill I_i=I_{0}exp\left(-\int x\left({\mathit{t}}_i+{\mathit{\eta }}_{i}r\right)\mathrm{d}r\right).\end{array}$$

(1)

Here, $I_i$ denotes the number of photons that reach the detector from a starting point ${\mathit{t}}_i$ in the direction of ${\mathit{\eta }}_i$, where $i=1,\cdots ,m$. $I_0$ is the initial intensity of the beam (i.e., the number of photons emitted by the source), and x is the function representing the attenuation coefficient of the object. This relation is usually written as:

$$\begin{array}{c}\hfill y_i=\int x\left({\mathit{t}}_i+{\mathit{\eta }}_{i}r\right)\mathrm{d}r,\end{array}$$

(2)

with $y_i=-log\frac{I_i}{I_0}$ representing the measured data for a single ray. In the following, we present the forward and backward operators for the X-ray setup shown in Figure 2. To study the depth imaging procedure within this specific setup, we need to understand forward and backward operators for both 2D and 3D objects.

2.2. Forward and Backward Operators for the Slice Imaging Setup

The proposed imaging setup is schematically depicted in Figure 2. X-rays are emitted by a source at lateral position $\mathit{s}\in {\mathbb{R}}^2$ and height h. Without loss of generality, we parametrize the detector position by $\mathit{r}\in {\mathbb{R}}^2$. The forward operator, $\mathcal{A}$, is now given by

$$y(\mathit{r},\mathit{s})=\left(\mathcal{A}x\right)(\mathit{r},\mathit{s})=\int x\left(t\mathit{r}+(1-t)\mathit{s},(1-t)h\right)\mathrm{d}t,$$

(3)

as shown in Figure 2. Here, $x(\mathit{u},z)$ denotes the attenuation coefficient of the object in terms of lateral ($\mathit{u}\in {\mathbb{R}}^2$) and vertical ($z\in \mathbb{R}$) coordinates. In fact, it integrates along the ray from the source at $(\mathit{s},h)$ to the detector located at $(\mathit{r},0)$.

2.2.1. A Fourier Slice Theorem

Defining the measurements in receiver-offset coordinates as $g(\mathit{r},\mathit{o})=y(\mathit{r},\mathit{r}-\mathit{o})$, we find the following relation between the 2D Fourier transform of g and the 3D Fourier transform of x. Let $\mathit{\xi }$ and $\zeta$ represent the spatial frequency coordinates, then it reads: (see the derivation in Appendix A.1)

$$\begin{array}{cc}\hfill {\mathcal{F}}_{2D}\left(g\right)(\mathit{\xi },\mathit{o})& ={\mathcal{F}}_{3D}\left(x\right)(\mathit{\xi },\zeta ),\hfill \\ \hfill \zeta & =h^{-1}{\mathit{\xi }}^{\top }\mathit{o},\hfill \end{array}$$

(4)

where ${\mathcal{F}}_{2D}$ and ${\mathcal{F}}_{3D}$ denote the 2D and 3D Fourier transforms, respectively, and $\mathit{o}=\mathit{r}-\mathit{s}$ signifies the vector difference between the detector position $\mathit{r}$ and the source position $\mathit{s}$. Figure 3 depicts sampled elements ${\mathcal{F}}_{2D}\left(g\right)$ (black circles) with respect to the Fourier object samples ${\mathcal{F}}_{3D}\left(x\right)$ (shown by green squares). Obviously, the sampling procedure is missing a portion of the object spectrum due to the limited values of $\mathit{o}$, i.e., the forward operator samples from $\zeta =h^{-1}{\mathit{\xi }}^{\top }\mathit{o}$. Since the range of $\mathit{o}$ is finite, this leads to a missing wedge in the frequency domain. This missing wedge leads to an ill-posed inverse problem of recovering x from g (or equivalently, y).

2.2.2. The 2D Slice at Depth z

To study the depth-imaging procedure, we introduce the corresponding forward operator ${\mathcal{A}}_z$ which maps a 2D image at depth z, denoted by $x_z\left(\mathit{u}\right)$ to the corresponding projection data $y(\mathit{r},\mathit{s})$ (See Appendix A.2),

$$\begin{array}{c}\hfill y\left(\mathit{r},\mathit{s}\right)=\left({\mathcal{A}}_z{x}_z\right)(\mathit{s},\mathit{r})=x_z\left(\left(1-h^{-1}z\right)\mathit{r}+h^{-1}z\mathit{s}\right).\end{array}$$

(5)

The corresponding adjoint operator ${\mathcal{A}}_z^{*}$ is given by (See Appendix A.2),

$$\begin{array}{c}\hfill \left({\mathcal{A}}_z^{*}y\right)\left(\mathit{u}\right)=\frac{1}{1-h^{-1}z}\int y\left(\frac{\mathit{u}-h^{-1}z\mathit{s}}{1-h^{-1}z},\mathit{s}\right)\mathrm{d}\mathit{s}.\end{array}$$

(6)

We then find that ${\mathcal{A}}_z^{*}{\mathcal{A}}_z$ is given by

$$\begin{array}{c}\hfill \left({\mathcal{A}}_z^{*}{\mathcal{A}}_z{x}_z\right)\left(\mathit{u}\right)=\frac{\int \mathrm{d}\mathit{s}}{1-h^{-1}z}{x}_z\left(\mathit{u}\right),\end{array}$$

(7)

so that for a limited range of $\mathit{s}$, ${\mathcal{A}}_z^{*}{\mathcal{A}}_z\propto \mathcal{I}$ and hence is trivially inverted. Depth imaging of an object consisting of a single thin layer can thus be achieved by a simple back-projection. When multiple layers are present we expect artifacts, which we study next.

2.3. Depth Image Reconstruction

Here, we are interested in studying the artifacts that arise when imaging a slice from measurements of a 3D layered object. In the remainder of this section, we let $h=1$ for ease of notation. As argued above, we can form a depth image by back projection, so it suffices to study the operator ${\mathcal{A}}_z^{*}\mathcal{A}$ (See Appendix A.3). In particular, we expect to image layered objects, which we can model as

$$x(\mathit{u},z)=\sum _{i=1}^M{x}_i\left(\mathit{u}\right)\delta (z-z_i).$$

Using a finite number of sources ${\left\{{\mathit{s}}_j\right\}}_{j=1}^N$, the depth image of the $k^{\mathrm{th}}$ layer is then given by

$$\begin{array}{cc}\hfill {\tilde{x}}_k\left(\mathit{u}\right)& =\sum _{i=1}^M\sum _{j=1}^N{x}_i\left(\frac{1-z_i}{1-z_k}\mathit{u}+\frac{z_i-z_k}{1-z_k}{\mathit{s}}_j\right)\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & =x_k\left(\mathit{u}\right)+N^{-1}\sum _{i\ne k}\sum _{j=1}^N{x}_i\left(\frac{1-z_i}{1-z_k}\mathit{u}+\frac{z_i-z_k}{1-z_k}{\mathit{s}}_j\right)\hfill \end{array}$$

(9)

Hence, the error in the depth image of the $k^{\mathrm{th}}$ layer is a superposition of spatially averaged versions of the other layers. An illustration of the resulting spatial sampling patterns is shown in Figure 4.

The source design problem is now to find a set of source positions ${\mathit{s}}_j$ such that unwanted parts in the depth image are minimized. As we don’t expect to be able to find a closed-form solution to the source design problem in general, we discretize the problem and resort to numerical optimization.

3. Discrete Problem

3.1. Forward and Back Projection

We discretize the forward model in (3) and express it as

$$\begin{array}{c}\hfill \mathit{y}=\mathit{A}\mathit{x},\end{array}$$

(10)

where $\mathit{y}\in {\mathbb{R}}^m$, $\mathit{x}\in {\mathbb{R}}^n$ and $\mathit{A}\in {\mathbb{R}}^{m\times n}$. Note that the number of rows in $\mathit{A}$ is the product of the number of detector pixels and the number of sources. Every block of rows corresponds to a particular X-ray source.

Let ${\mathit{A}}_z$, ${\mathit{x}}_z$ denote the forward matrix and the slice object at depth z, respectively. In general, ${\mathit{A}}_z^{\top }{\mathit{A}}_z$ will not be a multiple of the identity as in the continuous setting (7), so we consider solving ${\mathit{A}}_z^{\top }{\mathit{A}}_z{\mathit{x}}_z={\mathit{A}}_z^{\top }\mathit{y}$. We can view this as an approximation of 3D reconstruction through the Schur complement [10]. To see this, let ${\mathit{A}}_c$ and ${\mathit{x}}_c$ represent the forward operator and the object for all pixels except the ones contained in ${\mathit{x}}_z$. The linear equations for the forward projection (10) can now be expressed as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathit{A}}_z& {\mathit{A}}_c\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_z\\ {\mathit{x}}_c\end{array}\right]=\mathit{y}.\end{array}$$

(11)

By multiplying ${\mathit{A}}^{\top }$, it reads

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathit{A}}_z^{\top }{\mathit{A}}_z& {\mathit{A}}_z^{\top }{\mathit{A}}_c\\ {\mathit{A}}_c^{\top }{\mathit{A}}_z& {\mathit{A}}_c^{\top }{\mathit{A}}_c\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_z\\ {\mathit{x}}_c\end{array}\right]=\left[\begin{array}{c}{\mathit{A}}_z^{\top }\mathit{y}\\ {\mathit{A}}_c^{\top }\mathit{y}\end{array}\right].\end{array}$$

(12)

By utilizing Schur complement, one can re-write this in terms of ${\mathit{x}}_z$ as

$$\begin{array}{cc}\hfill & {\mathit{A}}^{\prime }{\mathit{x}}_z={\mathit{y}}^{\prime },\hfill \end{array}$$

(13)

where

$$\begin{array}{cc}\hfill & {\mathit{A}}^{\prime }=\left({\mathit{A}}_z^{\top }{\mathit{A}}_z-{\mathit{A}}_z^{\top }{\mathit{A}}_c{\left({\mathit{A}}_c^{\top }{\mathit{A}}_c\right)}^{-1}{\mathit{A}}_c^{\top }{\mathit{A}}_z\right),\hfill \\ \hfill & {\mathit{y}}^{\prime }=\left({\mathit{A}}_z^{\top }\mathit{y}-{\mathit{A}}_z^{\top }{\mathit{A}}_c{\left({\mathit{A}}_c^{\top }{\mathit{A}}_c\right)}^{-1}{\mathit{A}}_c^{\top }\mathit{y}\right).\hfill \end{array}$$

However, it is computationally not feasible to compute the corresponding Schur complements, and we can interpret depth imaging as neglecting the additional terms involving ${\mathit{A}}_c$.

3.2. Algebraic Reconstruction

Although the operator ${\mathcal{A}}_z$ is unitary (up to a scaling factor), the corresponding discrete operator ${\mathit{A}}_z$ is generally not. Moreover, it is desirable to stabilize the reconstruction with additional regularization. We therefore pose depth imaging as a regularized least-squares problem:

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_z=\underset{{\mathit{x}}_z}{\mathrm{argmin}}\frac{1}{2}\parallel {\mathit{A}}_z{\mathit{x}}_z{-\mathit{y}\parallel }^2+\frac{\lambda }{2}{\parallel {\mathit{x}}_z\parallel }^2,\end{array}$$

(14)

where $\lambda$ controls the trade-off between the two terms. While the problem is not extremely ill-posed, we utilize $\lambda$ as a safety factor to ensure the uniqueness of the reconstructed depth image. The above optimization can be solved using iterative algorithms such as the LSQR [11]. Various other regularization techniques, such as Total Variation regularization, or sparse regularization, introduce different penalty functions that promote desired properties in the reconstructed image.

4. Efficient Source Design for K Sources

Let $\mathit{\beta }\in {\mathbb{R}}^N$ be the design parameter, where the components of $\mathit{\beta }$ determine the influence of each source in the reconstruction process. The problem of finding the optimal design parameter $\mathit{\beta }$ for a fixed number of sources in the reconstruction of a specific depth plane can be formulated as a bi-level optimization problem as below,

$$\begin{array}{cc}\hfill & \underset{\mathit{\beta }}{min}\underset{\alpha \left(\mathit{\beta }\right)}{min}\frac{1}{2}{\parallel \alpha \left(\mathit{\beta }\right){\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\parallel }^2\mathrm{s}.\mathrm{t}\#\mathit{\beta }=K,\mathit{\beta }\in {\{0,1\}}^N,\hfill \end{array}$$

(15)

where ${\mathit{x}}_{\mathrm{gt}}\in {\mathbb{R}}^n$ represents the ground truth and the hyper-parameter $\alpha \left(\mathit{\beta }\right)\in \mathbb{R}$ serves as scaling factor, ensuring that $\alpha \left(\mathit{\beta }\right){\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)$ is on the same scale as ${\mathit{x}}_{\mathrm{gt}}$. Also, ${\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)$ is obtained by solving

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)=\underset{{\mathit{x}}_z}{\mathrm{argmin}}\frac{1}{2}\parallel {\mathit{T}}_{\mathit{\beta }}({\mathit{A}}_z{\mathit{x}}_z-\mathit{y}){\parallel }^2+\frac{\lambda }{2}{\parallel {\mathit{x}}_z\parallel }^2,\end{array}$$

(16)

where $\mathit{y}$ denotes the data corresponding to ${\mathit{x}}_{\mathrm{gt}}$ and ${\mathit{T}}_{\mathit{\beta }}$ is a linear operator which applies the weights $\mathit{\beta }$ to each source’s corresponding measurements. The optimization problem (15) minimizes the Euclidean distance between the reconstructed image (scaled by $\alpha \left(\mathit{\beta }\right)$) and the ground truth subject to certain constraints on $\mathit{\beta }$ to enforce its binary nature with K elements.

The binary constraint on $\mathit{\beta }$ in (15) can be relaxed using bound constraints [12], and the bi-level optimization is then formulated as

$$\begin{array}{c}\hfill {min}_{\mathit{\beta }}\frac{1}{2}{\parallel \widehat{\alpha }\left(\mathit{\beta }\right){\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\parallel }^2\mathrm{s}.\mathrm{t}{\mathit{\beta }}^{\top }\mathbf{1}=K,\mathbf{0}\le \mathit{\beta }\le \mathbf{1},\end{array}$$

(17)

while the optimal value of the scalar parameter $\alpha \left(\mathit{\beta }\right)$ in (15) is substituted in above optimization with $\widehat{\alpha }\left(\mathit{\beta }\right)=\frac{{\widehat{\mathit{x}}}_z{\left(\mathit{\beta }\right)}^{\top }{\mathit{x}}_{\mathrm{gt}}}{\parallel {\widehat{\mathit{x}}}_z{\left(\mathit{\beta }\right)\parallel }^2}$. Mathematically speaking, this problem seeks a $\mathit{\beta }$ that minimizes the Euclidean distance of $\widehat{\mathit{x}}\left(\mathit{\beta }\right)$ and ${\mathit{x}}_{\mathrm{gt}}$ and lies on the intersection of a hyperplane ${\mathit{\beta }}^{\top }\mathbf{1}=K$ and the hyperbox $0\le \mathit{\beta }\le 1$. Note that we cannot guarantee that this will produce binary solutions. In such cases, we convert the solution to a binary form by setting the K largest elements to one and assigning zeros to the remaining elements.

In practice the ground truth of ${\mathit{x}}_{\mathrm{gt}}$ is unknown, however, a set of calibration data ${\left\{{\mathit{x}}_{\mathrm{gt}}^i\right\}}_{i=1}^M$ of representative images may be available. This can then be utilized to compute an efficient source design such that minimizes the empirical Bayes risk. More precisely, for a given set of calibration data ${\left\{{\mathit{x}}_{\mathrm{gt}}^i\right\}}_{i=1}^M$ with corresponding measurements ${\left\{{\mathit{y}}^i\right\}}_{i=1}^M$, we can generalize (17) to

$$\begin{array}{c}\hfill \underset{\mathit{\beta }}{min}\frac{1}{2M}{\sum }_{i=1}^M{\parallel {\widehat{\alpha }}^i\left(\mathit{\beta }\right){\widehat{\mathit{x}}}_z^i\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}^i\parallel }^2\mathrm{s}.\mathrm{t}{\mathit{\beta }}^{\top }\mathbf{1}=K,0\le \mathit{\beta }\le 1,\end{array}$$

(18)

where ${\widehat{\mathit{x}}}_z^i\left(\mathit{\beta }\right)$ is the solution of (16) for the i-th calibration data ${\mathit{y}}^i$.

5. Implementation

To address the implementation of the Efficient K-Source Design (EKSD) method, it is worth mentioning that the non-linear optimization problem of (17) requires numerical optimization algorithms to start from an initial guess and refine the solution iteratively. This section discusses the solution of (17) using a projected gradient method.

5.1. Projected Gradient Method

In general, the projected gradient (PG) method solves,

$$\begin{array}{c}\hfill \underset{\mathit{\beta }}{min}f\left(\mathit{\beta }\right)s.t.\mathit{\beta }\in \mathcal{B},\end{array}$$

(19)

where $f:{\mathbb{R}}^N\to \mathbb{R}$ is a smooth function and $\mathcal{B}$ is a nonempty closed convex set. The projection onto $\mathcal{B}$ is the mapping $\mathcal{P}:{\mathbb{R}}^N\to {\mathbb{R}}^N$ defined by:

$$\begin{array}{c}\hfill \mathcal{P}\left(\mathit{u}\right)=\underset{\mathit{x}}{argmin}{\parallel \mathit{x}-\mathit{u}\parallel }^{2}s.t.\mathit{x}\in \mathcal{B}.\end{array}$$

(20)

The PG algorithm is defined by,

$$\begin{array}{c}\hfill {\mathit{\beta }}_{k+1}=\mathcal{P}\left({\mathit{\beta }}_k-\gamma \nabla f\left({\mathit{\beta }}_k\right)\right)\end{array}$$

(21)

where $\gamma \in \left(0,2/L\right)$ is the step size, L is the Lipschitz constant of $\nabla f$ , and $\nabla f\left({\mathit{\beta }}_k\right)$ is the gradient of f at ${\mathit{\beta }}_k$. The stopping criteria is defined as the condition where the Euclidean norm of the difference between consecutive design parameters, denoted as $\parallel {\mathit{\beta }}_{k+1}-{\mathit{\beta }}_k\parallel$, is less than or equal to a specified threshold value, denoted as $ϵ$.

5.2. Gradient of the Objective Function

To implement the non-linear optimization of (17) using the PG algorithm, the gradient calculation of $\mathit{F}(\widehat{\mathit{x}}\left(\mathit{\beta }\right),\widehat{\alpha }\left(\mathit{\beta }\right)):=\frac{1}{2}{\parallel \widehat{\alpha }\left(\mathit{\beta }\right)\widehat{\mathit{x}}\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\parallel }^2$ is required. Here, it reads,

$$\begin{array}{cc}\hfill \nabla F\left(\mathit{\beta }\right)& =\widehat{\alpha }{\frac{\partial {\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}^{\top }\left(\widehat{\alpha }{\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\right)\hfill \\ \hfill & =2{\widehat{\alpha }\mathit{C}}^{\top }{\left({\mathit{A}}^{\top }{\mathit{T}}_{\mathit{\beta }}^{\top }\right)}_{\lambda }^{†}\left(\widehat{\alpha }{\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\right),\hfill \end{array}$$

(22)

where ${\left({\mathit{A}}^{\top }{\mathit{T}}_{\mathit{\beta }}^{\top }\right)}_{\lambda }^{†}={\left({\mathit{A}}^{\top }{\mathit{T}}_{\mathit{\beta }}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{A}+\lambda \mathit{I}\right)}^{-1}{\mathit{A}}^{\top }{\mathit{T}}_{\mathit{\beta }}^{\top }$, and $\mathit{C}=\left[{\mathit{c}}_1,\cdots ,{\mathit{c}}_N\right]$

with ${\mathit{c}}_i=\frac{\partial {\mathit{T}}_{\mathit{\beta }}}{\partial {\beta }_i}\left(\mathit{y}-\mathit{A}{\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)\right)$. (See the derivation in Appendix A.4)

5.3. Projection Onto the Simplex Constraint

Here, we investigate the projection operator $\mathcal{P}$ presented in (21) which maps each point onto the simplex set of $\mathcal{B}=\{\mathit{u}\in {\mathbb{R}}^m|{\mathit{u}}^{\top }\mathbf{1}=K,0\le u_i\le 1\}$. (See the derivation in Appendix A.5)

$$\begin{array}{c}\hfill \mathcal{P}\left(\mathit{p}\right)=min\{max\{\mathit{p}-\mu \mathbf{1},\mathbf{0}\},\mathbf{1}\},\end{array}$$

(23)

where $\mathit{p},\mathbf{1}\in {\mathbb{R}}^m$ and $\mu$ can be obtained by solving the following optimization problem:

$$\begin{array}{c}\hfill \sum _{i}min\{max\{p_i-\mu ,0\},1\}=K,\end{array}$$

(24)

which can be solved using the Newton method.

6. Numerical Experiments

Here, we present some results related to our proposed method (EKSD). We evaluate the performance of EKSD using the following criteria:

i: Normalized mean square error (NMSE) can be defined as the error between the reconstructed image and ground truth

$$\begin{array}{c}\hfill \mathrm{NMSE}:=\frac{\parallel \widehat{\alpha }{\widehat{\mathit{x}}}_z-{\mathit{x}}_{\mathrm{gt}}{\parallel }^2}{\parallel {\mathit{x}}_{\mathrm{gt}}{\parallel }^2}.\end{array}$$

(25)

ii: Structural similarity index measure (SSIM): Unlike NMSE, SSIM takes into account the structural information and perceived quality of the images. It compares three aspects of the images: luminance, contrast, and structure, and computes a similarity index between 0 and 1, where 1 indicates perfect similarity [13].

6.1. Data Set

As a benchmark example we consider a two-layer PCB of size $80\times 250\times 250$ pixels. It consists of two distinct wire designs (as illustrated in Figure 5) with a thickness of 25 pixels each. These designs are separated by a gap of 20 pixels. The layer of interest is the one depicted in Figure 5a.

We create 200 defect shapes in various locations, forming a data set of defective wire design for the bottom layer. A few samples of this data set are depicted in Figure 6. We use half of this data set for source design (calibration) and half for validation.

6.2. Implementation Details

We define the X-ray setup as shown in Figure 2. The distance h, between the center of the array and the detector, is equal to $6L$. Here, L represents the width of the phantom. Moreover, we consider $\mathit{r}\in {[-L/2,L/2]}^2$ (discretized uniformly with 250 × 250 pixels) and $\mathit{s}\in {[-D/2,D/2]}^2$ (discretized uniformly with $\sqrt{N}\times \sqrt{N}$ source positions), where D is the fixed array length. The projection data are generated using the ASTRA toolbox [14]. For illustration, we present the projection images from all sources in a $4\times 4$ array of sources in Figure 7. It is challenging to inspect individual layers of a the phantom with a single projection, as overlapping them hinders a proper investigation.

We choose a small regularization parameter $\lambda =0.1$ to ensure the solution’s stability and uniqueness without significantly biasing the reconstruction. We carefully tuned the step size, $\gamma ={10}^{-9}$, to be sufficiently small, reducing the risk of algorithmic overshooting during optimization. Additionally, we manually adjusted the stopping threshold, $ϵ=0.001$, to ensure that the algorithm stops at approximately consistent convergence points. The initial value of the parameter vector $\mathit{\beta }$ is initialized either with a vector ones or with i.i.d. uniform random values in $[0,1]$. We refer to these methods as EKSD(K) and EKSDR(K), respectively. In some experiments, we will compare the results to completely random designs of K out of N source positions. We refer to this as the RANDOM method.

6.3. Convergence

We evaluate the convergence behavior of the EKSD algorithm across various iterations. Here, we let $K=10$, $N=16$, $D=4L$. In Figure 8, we depict the NMSE at each iteration as evaluated on the calibration (left) and validation (right) data. This figures shows that the design is not overfitted on the calibration data.

6.4. Exploring the Effects of Array Length

Here, we study the influence of the array length D on the reconstruction error for $N=16,36,64$. Hence, we compare the average NMSE over the validation data across various array lengths, as depicted in Figure 9. The heightened errors observed for both smaller and larger array lengths indicate that certain source positions within the array contribute to artifacts in the depth image, leading to an elevated NMSE.

6.5. Exhaustive Search

We carried out an exhaustive search to assess the performance of EKSD. We consider $N=16$ source positions with $D=4L$. Figure 10 shows the average NMSE over the validation data for all possible K-source designs for each K, as well as the average NMSE achieved by the EKSD-design. We see the the designs found by EKSD are nearly optimal.

6.6. Source Design 1

We consider three configurations with $N=16,36,64$ at $D=4L$, and we show the NMSE and SSIM results in Table 1. Initially, we display results for depth-reconstructed images utilizing all source positions within the array. Subsequently, we derive $K=10$ efficient source designs with EKSD and EKSDR and present the corresponding NMSE and SSIM results. For each value of N, we observe that EKSD and EKSDR outperform than the case using all source positions. This improvement can be attributed to the reduction of artifacts in the images, which is likely due to the selective nature of measurements in the EKSD and EKSDR approaches. As a comparison we show the average NMSE and SSIM of random designs with $K=10$ sources, and these are significantly worse than all other methods.

Figure 11 shows the corresponding depth images of the bottom layer reconstructed with source positions suggested by EKSD. As demonstrated in Section 2.3, artifacts in these images attributed to the superposition of features from different layers. The lack of rotational acquisition in laminography leads to this artifacts. Clearly, the designs proposed by EKSD minimize artifacts in depth images compared to random designs, facilitating accurate defect detection within the depth images.

6.7. Source Design 2

We repeat the previous experiment with the array length optimally according to Figure 9; $D=2.4L$ for different N values. The results are summarized in Table 2. As the full designs are already optimal in a sense, there is less dramatic improvement when optimizing the design with fewer source positions. However, the results for $N=64$ suggest that increasing the number of positions while maintaining a fixed array length, thereby reducing the distance between sources, may result in a non-binary solution. Then, converting the solution into a sparse vector with $K=10$ elements, results in a less optimal outcome. The proposed source designs and corresponding depth images are shown in Figure 12.

7. Discussion

In this research, we explore the methodology of experimental design for a laminographic setup, acknowledging potential challenges that may arise in real-world applications. Here, we outline several significant challenges with our proposed method: (i) The non-convex nature of the objective function may lead to entrapment in local minima. Addressing this, diversifying initial conditions enhances the likelihood of identifying more optimal solutions or the global minimum. (ii) Our method optimizes the positions of sources within a specified grid. However, this grid may not represent the most efficient set of source positions for optimization. An alternative approach could involve optimizing within a continuous domain of positions, leveraging the closed-form formulation for depth image reconstruction presented in this paper. (iii) While our investigation of the experimental design problem assumes noise-free measurements, this assumption may not hold in real-world scenarios characterized by noisy data. This can be addressed by incorporating noise directly into the optimization process to enhance real-world applicability.

8. Conclusion

In real world scenarios, imaging a particular depth of a planar object can yield valuable information with less effort than a full 3D reconstruction. In a depth image, the information corresponding to the chosen layer appears sharp, while the details of other layers appear blurred. We argue that source design plays a role in mitigating these unwanted artifacts and pose the experimental design problem as a bi-level optimization problem. To solve this problem, we employed a bilvel optimization and implemented it by gradient-based methods. The suggested framework can be tailored more specifically to a particular application by incorporating different regularization terms in lower level optimization or by adopting alternate loss function in the upper level optimization to more effectively capture the characteristics of the image.

Appendix A.1. Forward Operator for the 3D Object

Figure 2 presents the 2D view of acquisition setup for a planar object $x\left(\mathit{u},z\right)$ where $\mathit{u}\in {\mathbb{R}}^2$ is the lateral coordinate and $z\in \mathbb{R}$ the depth coordinate. The line integral of $x\left(\mathit{u},z\right)$ from $\left(\mathit{s},h\right)$ to $\left(\mathit{r},0\right)$ for a fixed h can be expressed as follows:

$$\begin{array}{c}\hfill y\left(\mathit{r},\mathit{s}\right)=\int x\left(\mathit{s}+\left(\mathit{r}-\mathit{s}\right)t,h-ht\right)\mathrm{d}t\end{array}$$

(A1)

where $\mathit{s}\in {\mathbb{R}}^2$, $\mathit{r}\in {\mathbb{R}}^2$ denote the source and detector positions, respectively. Next, we substitute $x(\mathit{u},z)$ by its representation in the Fourier domain ${\mathcal{F}}_{3D}\left(x\right)(\mathit{\xi },\zeta )$:

$$\begin{array}{c}\hfill x(\mathit{u},z)=\int \int {\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },\zeta \right)e^{j2\pi \left({\mathit{\xi }}^{\top }\mathit{u}+\zeta z\right)}\mathrm{d}\mathit{\xi }\mathrm{d}\zeta ,\end{array}$$

(A2)

using (A2) in (A1), it yields:

$$\begin{array}{c}\hfill y(\mathit{r},\mathit{s})=\int {\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },\zeta \right)e^{j2\pi \left({\mathit{\xi }}^{\top }\left(\mathit{s}+\left(\mathit{r}-\mathit{s}\right)t\right)+\zeta \left(h-ht\right)\right)}\mathrm{d}t\mathrm{d}\mathit{\xi }\mathrm{d}\zeta ,\end{array}$$

(A3)

using $\int e^{j2\pi \left({\mathit{\xi }}^{\top }\left(\mathit{r}-\mathit{s}\right)-\zeta h\right)t}\mathrm{d}t=\delta \left({\mathit{\xi }}^{\top }\left(\mathit{r}-\mathit{s}\right)-\zeta h\right)$, it follows,

$$\begin{array}{c}\hfill y(\mathit{r},\mathit{s})=\int {\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },\zeta \right)e^{j2\pi \left({\mathit{\xi }}^{\top }\mathit{s}+\zeta h\right)}\delta \left({\mathit{\xi }}^{\top }\left(\mathit{r}-\mathit{s}\right)-\zeta h\right)\mathrm{d}\zeta \mathrm{d}\mathit{\xi },\end{array}$$

(A4)

and then,

$$\begin{array}{c}\hfill y(\mathit{r},\mathit{s})=\int {\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },h^{-1}\mathit{\xi }(\mathit{r}-\mathit{s})\right)e^{j2\pi {\mathit{\xi }}^{\top }\mathit{r}}\mathrm{d}\mathit{\xi }\end{array}$$

(A5)

using $\mathit{o}=\mathit{r}-\mathit{s}$, it can be written,

$$\begin{array}{c}\hfill g(\mathit{r},\mathit{o})=\int {\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },h^{-1}{\mathit{\xi }}^{\top }\mathit{o}\right)e^{j2\pi {\mathit{\xi }}^{\top }\mathit{r}}\mathrm{d}\mathit{\xi }\end{array}$$

(A6)

which means we can link the 3D Fourier transform of x to the 2D Fourier transform of the data as,

$$\begin{array}{c}\hfill {\mathcal{F}}_{2D}\left(g\right)(\mathit{\xi },\mathit{o})={\mathcal{F}}_{3D}\left(x\right)\left(\mathit{\xi },h^{-1}{\mathit{\xi }}^{\top }\mathit{o}\right).\end{array}$$

(A7)

Appendix A.2. Forward Operator for a 2D Slice on the Depth z

We introduce $x_z$ as a slice of the object x at depth z, i.e.,

$$x(\mathit{u},z^{\prime })=x_z\left(\mathit{u}\right)\delta (z^{\prime }-z).$$

We then define the forward operator ${\mathcal{A}}_z$ as projecting a 2D object to the detector at $z=0$ via,

$$\begin{array}{c}\hfill y\left(\mathit{r},\mathit{s}\right)=x_z\left(\mathit{s}+\left(1-h^{-1}z\right)\left(\mathit{r}-\mathit{s}\right)\right)\end{array}$$

(A8)

or equivalently as,

$$\begin{array}{c}\hfill y\left(\mathit{r},\mathit{s}\right)=x_z\left(\mathit{r}\left(1-h^{-1}z\right)+h^{-1}z\mathit{s}\right).\end{array}$$

(A9)

The adjoint operator ${\mathcal{A}}_z^{*}$ can be obtained by considering the definition ${\langle {\mathcal{A}}_z{x}_z,y\rangle }_{(\mathit{s},\mathit{r})}={\langle x_z,{\mathcal{A}}_z^{*}y\rangle }_{\mathit{v}}$. First, we have

$$\begin{array}{ccc}\hfill & {\langle {\mathcal{A}}_z{x}_z,y\rangle }_{(\mathit{s},\mathit{r})}=\hfill & \hfill \int \int x_z\left(\mathit{r}\left(1-h^{-1}z\right)+h^{-1}z\mathit{s}\right)y\left(\mathit{r},\mathit{s}\right)\mathrm{d}\mathit{r}\mathrm{d}\mathit{s}.\end{array}$$

Now let $\mathit{v}=\mathit{r}\left(1-h^{-1}z\right)+h^{-1}z\mathit{s}$, then $\mathit{r}=\frac{\mathit{v}-h^{-1}z\mathit{s}}{1-h^{-1}z}$ and $\mathrm{d}\mathit{r}=\frac{1}{1-h^{-1}z}\mathrm{d}\mathit{v}$, and it follows,

$$\begin{array}{cc}& {\langle {\mathcal{A}}_z{x}_z,y\rangle }_{(\mathit{s},\mathit{r})}=\hfill \\ \hfill & \frac{1}{1-h^{-1}z}\int \int x_z\left(\mathit{v}\right)y\left(\frac{\mathit{v}-h^{-1}z\mathit{s}}{1-h^{-1}z},\mathit{s}\right)\mathrm{d}\mathit{v}\mathrm{d}\mathit{s}\hfill \\ \hfill & =\langle x_z,{\mathcal{A}}_z^{*}y\rangle .\hfill \end{array}$$

(A10)

We thus find that ${\mathcal{A}}_z^{*}f$ is defined as

$$\begin{array}{c}\hfill \left({\mathcal{A}}_z^{*}y\right)\left(\mathit{u}\right)=\frac{1}{1-h^{-1}z}\int y\left(\frac{\mathit{u}-h^{-1}z\mathit{s}}{1-h^{-1}z},\mathit{s}\right)\mathrm{d}\mathit{s}.\end{array}$$

(A11)

By using (A9) in (A11), the composition ${\mathcal{A}}_z^{*}{\mathcal{A}}_z$ can be written as below

$$\begin{array}{c}\hfill \left({\mathcal{A}}_z^{*}{\mathcal{A}}_z{x}_z\right)\left(\mathit{u}\right)=\frac{x_z\left(\mathit{u}\right)}{1-h^{-1}z}\int \mathrm{d}\mathit{s}.\end{array}$$

(A12)

Appendix A.3. Back Projection for a Depth Image

We use $y\left(\mathit{r},\mathit{s}\right)$ in (A1) to express ${\mathcal{A}}_z^{*}\mathcal{A}$ using (A11),

$$\begin{array}{cc}\hfill & \left({\mathcal{A}}_z^{*}\mathcal{A}x\right)\left(\mathit{u}\right)=\hfill \\ \hfill & \frac{1}{1-h^{-1}z}\int \int x\left(\frac{\mathit{u}-h^{-1}z\mathit{s}}{1-h^{-1}z}t+\left(1-t\right)\mathit{s},h\left(1-t\right)\right)\mathrm{d}t\mathrm{d}\mathit{s}.\hfill \end{array}$$

(A13)

Now let $z^{\prime }=h(1-t)$, then it can be rewritten as

$$\begin{array}{cc}\hfill & \left({\mathcal{A}}_z^{*}\mathcal{A}x\right)\left(\mathit{u}\right)=\hfill \\ \hfill & \frac{1}{h-z}\int \int x\left(\frac{h-z^{\prime }}{h-z}\mathit{u}+\frac{z^{\prime }-z}{h-z}\mathit{s},z^{\prime }\right)\mathrm{d}{z}^{\prime }\mathrm{d}\mathit{s}\hfill \end{array}$$

(A14)

Appendix A.4. Gradient Calculation of the Objective Function

To find the gradient of $\mathit{F}(\widehat{\mathit{x}}\left(\mathit{\beta }\right),\widehat{\alpha }\left(\mathit{\beta }\right)):=\frac{1}{2}{\parallel \widehat{\alpha }\left(\mathit{\beta }\right)\widehat{\mathit{x}}\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\parallel }^2$ with respect to $\mathit{\beta }$, by the chain roll it reads:

$$\begin{array}{c}\hfill \frac{\partial \mathit{F}}{\partial \mathit{\beta }}({\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right),\widehat{\alpha }\left(\mathit{\beta }\right))=\frac{\partial \mathit{F}}{\partial {\widehat{\mathit{x}}}_z}·\frac{\partial {\widehat{\mathit{x}}}_z}{\partial \mathit{\beta }}+\frac{\partial \mathit{F}}{\partial \widehat{\alpha }}·\frac{\partial \widehat{\alpha }}{\partial \mathit{\beta }}\end{array}$$

(A15)

where the second term is equal to zero because $\widehat{\alpha }$ is the optimal point of $\mathit{F}$. Then, the first term can be written as:

$$\begin{array}{c}\hfill \frac{\partial \mathit{F}}{\partial \mathit{\beta }}={\left(\frac{\partial {\widehat{\mathit{x}}}_z}{\partial \mathit{\beta }}\right)}^{\top }\frac{\partial \mathit{F}}{\partial {\widehat{\mathit{x}}}_z}=\widehat{\alpha }{\left(\frac{\partial {\widehat{\mathit{x}}}_z}{\partial \mathit{\beta }}\right)}^{\top }\left(\widehat{\alpha }{\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\right).\end{array}$$

(A16)

To obtain an expression for the derivative of $\widehat{\mathit{x}}$ with respect to $\mathit{\beta }$, we exploit the fact that the lower-level optimization problem (16) has a closed-form solution:

$$\begin{array}{cc}\hfill {\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)=& {\left({\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{A}+\lambda \mathit{I}\right)}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{y}.\hfill \end{array}$$

(A17)

We define ${\mathit{H}}_{\mathit{\beta }}:={\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{A}+\lambda \mathit{I}$, and then the derivative follows,

$$\begin{array}{c}\hfill \frac{\partial \widehat{\mathit{x}}\left(\mathit{\beta }\right)}{\partial {\beta }_i}=\frac{\partial {{\mathit{H}}_{\mathit{\beta }}}^{-1}}{\partial {\beta }_i}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{y}+{{\mathit{H}}_{\mathit{\beta }}}^{-1}\frac{\partial \left({\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{y}\right)}{\partial {\beta }_i},\end{array}$$

(A18)

recall that $\frac{\partial {{\mathit{H}}_{\mathit{\beta }}}^{-1}}{\partial {\beta }_i}=-{{\mathit{H}}_{\mathit{\beta }}}^{-1}\frac{\partial {\mathit{H}}_{\mathit{\beta }}}{\partial {\beta }_i}{{\mathit{H}}_{\mathit{\beta }}}^{-1}$, therefore,

$$\begin{array}{c}\hfill \frac{\partial {{\mathit{H}}_{\mathit{\beta }}}^{-1}}{\partial {\beta }_i}=-2{{\mathit{H}}_{\mathit{\beta }}}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }\frac{\partial T_{\mathit{\beta }}}{\partial {\beta }_i}\mathit{A}{{\mathit{H}}_{\mathit{\beta }}}^{-1}\end{array}$$

(A19)

then (A18) can be rewritten as,

$$\begin{array}{cc}\hfill \frac{\partial \widehat{\mathit{x}}\left(\mathit{\beta }\right)}{\partial {\beta }_i}& =-2{{\mathit{H}}_{\mathit{\beta }}}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }\frac{\partial T_{\mathit{\beta }}}{\partial {\beta }_i}\left(\mathit{A}{{\mathit{H}}_{\mathit{\beta }}}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{y}-\mathit{y}\right),\hfill \end{array}$$

(A20)

where ${{\mathit{H}}_{\mathit{\beta }}}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }{\mathit{T}}_{\mathit{\beta }}\mathit{y}$ is equal to $\widehat{\mathit{x}}\left(\mathit{\beta }\right)$ in (A17) and ${\mathit{H}}_{\mathit{\beta }}^{-1}{\mathit{A}}^{\top }{{\mathit{T}}_{\mathit{\beta }}}^{\top }$ is denoted by ${\left({\mathit{T}}_{\mathit{\beta }}\mathit{A}\right)}_{\lambda }^{†}$, then it can be rewritten,

$$\begin{array}{c}\hfill \frac{\partial \widehat{\mathit{x}}\left(\mathit{\beta }\right)}{\partial {\beta }_i}=2{\left({\mathit{T}}_{\mathit{\beta }}\mathit{A}\right)}_{\lambda }^{†}\left(\frac{\partial {\mathit{T}}_{\mathit{\beta }}}{\partial {\beta }_i}\left(\mathit{y}-\mathit{A}\widehat{\mathit{x}}\left(\mathit{\beta }\right)\right)\right)=2{\left({\mathit{T}}_{\mathit{\beta }}\mathit{A}\right)}_{\lambda }^{†}{\mathit{c}}_i,\end{array}$$

(A21)

for simplicity, we denote $\frac{\partial {\mathit{T}}_{\mathit{\beta }}}{\partial {\beta }_i}\left(\mathit{y}-\mathit{A}\widehat{\mathit{x}}\left(\mathit{\beta }\right)\right)$ as ${\mathit{c}}_i\in {\mathbb{R}}^m$ and it follows,

$$\begin{array}{c}\hfill \frac{\partial \widehat{\mathit{x}}\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}=2{\left({\mathit{T}}_{\mathit{\beta }}\mathit{A}\right)}_{\lambda }^{†}\mathit{C}\end{array}$$

(A22)

where $\mathit{C}=[{\mathit{c}}_1,\cdots ,{\mathit{c}}_N]$. By replacement of (A22) into (22), yields

$$\begin{array}{c}\hfill \nabla F=2\widehat{\alpha }{\mathit{C}}^{\top }{\left({\mathit{A}}^{\top }{\mathit{T}}_{\mathit{\beta }}^{\top }\right)}_{\lambda }^{†}\left(\widehat{\alpha }{\widehat{\mathit{x}}}_z\left(\mathit{\beta }\right)-{\mathit{x}}_{\mathrm{gt}}\right).\end{array}$$

(A23)

Appendix A.5. Projection Onto the Simplex Constraint

Let $\mathcal{P}$ be an operator that maps each point onto the simplex set of $\mathcal{B}=\{\mathit{u}\in {\mathbb{R}}^m|{\mathit{u}}^{\top }\mathbf{1}=K,0\le u_i\le 1\}$. In fact, the orthogonal projection of $\mathit{p}$ is the unique solution of:

$$\begin{array}{c}\hfill \underset{\mathit{t}}{argmin}\frac{1}{2}{\parallel \mathit{t}-\mathit{p}\parallel }^2\mathrm{s}.\mathrm{t}{\mathit{t}}^{\top }\mathbf{1}=K,0\le t_i\le 1,\end{array}$$

(A24)

to solve the above optimization via its dual problem, we introduce the Lagrangian [15],

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{t},\mu ):=\frac{1}{2}{\parallel \mathit{t}-\mathit{p}\parallel }^2+\mu \left({\mathit{t}}^{\top }\mathbf{1}-K\right),\end{array}$$

(A25)

where $\mu \in \mathbb{R}$ is a Lagrangian multiplier. The dual problem can now be written as

$$\begin{array}{c}\hfill \mathcal{J}\left(\mu \right)=\underset{0\le t_i\le 1}{inf}\mathcal{L}\left(\mathit{t},\mu \right),\end{array}$$

(A26)

where $\mathcal{J}\left(\mu \right)$ is the dual function as the minimum value of Lagrangian over $\mathit{t}$. Regarding strong duality the optimal solution of (A24) and (A26) are equal if there exists a dual variable of $\mu \in \mathbb{R}$ that makes (A26) as a tight lower bound for (A24). The solution of the dual problem of (A26) can be written as [16],

$${\mathit{t}}^{★}\left(\mu \right)=min\{max\{\mathit{p}-\mu \mathbf{1},\mathbf{0}\},\mathbf{1}\}$$

where $\mu \in \mathbb{R}$ is obtained by solving ${\mathit{t}}^{★}{\left(\mu \right)}^{\top }\mathbf{1}=K$. Let’s define the function $w\left(\mu \right)={\sum }_{i}min\{max\{p_i-\mu ,0\},1\}-K$ that is non-increasing and can be solved using the Newton method. The derivative is given by,

$$w^{\prime }\left(\mu \right)=\left|\mathcal{I}\right|,\mathcal{I}=\{i\in \{0,\cdots ,n\}:0\le p_i-\mu \le 1\}.$$