Determination of droplet size from wide-angle light scattering image data using convolutional neural networks

The WALS image data can be described as a map $I:W\rightarrow\mathbb{R},$ where the pixel space $W = \{0,1,\dots,999\}\times \{0,1,\dots,999\}$ is a discretized square and $I(v)\in\mathbb{R}$ denotes the grayscale value of pixel $v = (v_1,v_2)\in W$. However, due to the nature of the WALS device, only those pixels within a ring around a pole $(c_1,c_2)\in [0,999]\times[0,999]$ are illuminated, pixels outside of that ring contain no information with respect to the droplet, see figure 2(left). Therefore, several steps were taken to prepare the images for analysis. The most important step involves the application of a polar transformation to the images, which is useful for unwrapping circular objects (such as WALS image data).

More specifically, the polar transformation is applied to the input image I in order to obtain a transformed image $I_\mathrm{p}:W_\mathrm{p}\rightarrow\mathbb{R},$ where the polar pixel space $W_\mathrm{p}$ is a set with equidistant points that spans the entire two-dimensional polar coordinate system $(0,500] \times [0,2\pi)$. For that purpose, Cartesian coordinates $(v_1,v_2)\in W$ are transformed into polar coordinates $(r,\theta) \in (0,500] \times [0, 2\pi)$ with respect to a given center point $(c_1,c_2)\in[0,999]\times[0,999]$. More precisely, the polar coordinates $(r,\theta)$ are calculated from the Cartesian coordinates $(v_1,v_2)$ using the following equations:

$$\begin{align*} r\left(v_1,v_2\right) & = \sqrt{\left(v_1-c_1\right)^2 + \left(v_2-c_2\right)^2} \\ \theta\left(v_1,v_2\right) & = \text{tan}^{-1}\left(\frac{v_2-c_2}{v_1-c_1}\right). \end{align*}$$

In this manner we obtain the pair of polar coordinates $(r(v_1,v_2),\theta(v_1,v_2))$ and the corresponding intensity value $I(v_1,v_2)$ for each pixel $(v_1,v_2) \in W$. The intensity values of the transformed image $I_\mathrm{p}$, which are defined on the set $W_\mathrm{p}$, are then determined by linear interpolation.

In order to maintain a high image quality while also reducing the image size, the size of the polar transformed image is chosen to be $500\times500.$ Furthermore, in order to determine the region of interest in the polar transformation of I (i.e. the ring in which WALS signals are measured), a reference image was additionally transformed. More precisely, this reference image was obtained by using a separate image in which the entire mirror area was well-lit, highlighting the pixels that could potentially contain relevant information. We obtained the range of r and θ values that corresponded to the region of interest by identifying the illuminated region in the transformed reference image. Then, the polar transformation $I_\mathrm{p}$ is cropped to only include this region, thus resulting in a reduced image size of $56\times 496$, while retaining all relevant information. Note that, in addition to identifying the region of interest, the reference image was also used to obtain the center point exploited for the transformation into polar coordinates.

However, before applying the polar transformation and cropping the images, a common type of imaging artifacts is identified and removed. These artifacts manifest as singular bright spots with intensity value equal to the maximum intensity value of 255. By 'singular,' we mean that these bright spots are isolated and not part of a larger pattern or structure in the image, as can be seen by the singular gray spots (corresponding to a grayscale value of 255) in the bottom-left magnified cutout of figure 2(left). It is important to remove these bright spots because, during the training of a CNN, they can strongly influence convolutions and lead to overfitting. This is implemented by determining all pixels with intensity value equal to 255. For each such pixel $v \in W$ with $I(v) = 255$, we compute the average grayscale value of its neighboring pixels, where a pixel is considered a neighbor if it is horizontally, vertically, or diagonally adjacent to the pixel under consideration. If the average grayscale value of these neighboring pixels is below a threshold value of 128, we replace I(v) by the calculated average value.

In summary, we applied several preprocessing steps to the WALS image data in order to reduce the input size, while retaining all relevant information and improving the quality of the images for further analysis. From this point on, all analyses and model training steps are performed exclusively on this preprocessed WALS image data.

It is important to reiterate that two different networks are used in order to predict the class labels $s_1,s_2$ and the quantitative descriptors of droplet size $s_3,s_4$. Nonetheless, both networks possess the same basic architecture, which consists of two main parts: a convolutional part and a dense part. The convolutional part is designed to extract and process features from the input data, while the dense part is used to perform classification on the extracted features. A schematic representation of the architecture is shown in figure 3.

Zoom In Zoom Out Reset image size

The convolutional part of the architecture includes several residual blocks and a dilated residual block, similar to the downsampling part in the U-net architecture used in [4]. Each residual block consists of several convolutional layers with a specified number of filters, followed by batch normalization and rectified linear unit (ReLU) activation layers. The input of each residual block is added to its output in order to form a residual connection, allowing the residual block to learn residual mappings between the input and output, see [6]. The use of such residual blocks has been shown to improve learning in deep neural networks, see [5]. The residual blocks also include a max-pooling layer with a specified pool size. These max-pooling layers are used to reduce the spatial dimensions of the output from each block, as in [18]. The dilated residual block is designed to increase the receptive field of the network without increasing the number of parameters. The dilated residual block consists of several dilated convolutional layers with a specified dilation rate and number of filters, each of which are followed by batch normalization and ReLU activation layers. In this design, each successive residual block doubles the number of filters, starting from an initial hyperparameter value n_f, the model architecture for ${n}_{\text{f}} = 64$ is depicted in figure 3. This approach ensures a progressive increase in the network's capacity to extract increasingly complex features at each residual block.

The dense part of the architecture is used to perform classification and prediction on the extracted features of the convolutional part. The output of the final convolutional layer is flattened and passed through three dense layers with $2\cdot {n}_{\text{f}},{n}_{\text{f}}$ and $\frac{{n}_{\text{f}}}{2}$ units, respectively, and ReLU activation functions. The final layer, i.e. the output layer, of the base architecture consists of two units, where the activation function is either a sigmoid (for the architecture which determines the class labels s₁ and s₂) or a ReLU function (for the architecture which predicts the droplet sizes s₃ and s₄).

We aim to identify an effective choice for hyperparameter n_f that strikes a balance between the network's ability to capture essential features and reduced model complexity. For that purpose, we systematically vary n_f by putting it to values in the set $\{8,16, 32, 64\}$, details of this selection process are described in section 2.3.5. The remaining neural network parameters, denoted by φ , which we will refer to as trainable parameters, are chosen via gradient descent as described in the next section.

In order to accurately predict the presence and size of droplets from preprocessed WALS images, we need to train our neural networks on datasets of labeled examples. This training process involves adjusting the network's trainable parameters to minimize the discrepancy between the predicted and true outputs.

The available data, denoted by $\mathcal{D}_{\text{}}$, consists of n pairs of input images and their corresponding output vectors for some integer n > 1, i.e. $\mathcal{D}_{\text{}} = \{(I_\mathrm{p}^{(i)},s^{(i)})\}_{i = 1}^n$ where $I_\mathrm{p}^{(i)}$ is the ith input image and $s^{(i)} = (s_1^{(i)},s_2^{(i)},s_3^{(i)},s_4^{(i)})$ is the corresponding ground truth vector. We use two neural networks: one for predicting the binary labels $s_1,s_2$, which indicate whether a droplet was detected on the left and right ring-shaped segments, respectively; and another one for predicting the quantitative descriptors $s_3,s_4$ of droplet size, corresponding to the ground truth droplet size measured in the left and right ring-shaped segments, respectively. Thus, the respective training sets $\mathcal{D}_{\text{train,1}}$ and $\mathcal{D}_{\text{train,2}}$ for the two networks consist of the following pairs of input images and their corresponding output vectors: $\mathcal{D}_{\text{train,1}} = \{(I_\mathrm{p},(s_1,s_2)): (I_\mathrm{p},s)\in\mathcal{D}_{\text{train}} \}\}$ and $\mathcal{D}_{\text{train,2}} = \{(I_\mathrm{p},(s_3,s_4)): (I_\mathrm{p},s)\in\mathcal{D}_{\text{train}} \,\text{with } s_1+s_2\gt0\}\}$, where $\mathcal{D}_{\text{train}}\subset\mathcal{D}_{\text{}}$ denotes the data that can be utilized during training. The condition that $s_1+s_2\gt0$ for $\mathcal{D}_{\text{train,2}}$ ensures that the training set for the droplet size prediction network includes only those images for which a droplet size can be detected on at least one ring-shaped segment. In order to increase the variety in the training data, we use training data augmentation [13, 21]. This involves the application of random transformations to the input images in order to generate additional training data. Specifically, we use reflection and translation as options for data augmentation.

Reflection involves flipping the image vertically with a probability of 0.5. If the image is flipped vertically, the ground truth values are adjusted accordingly. This adjustment is necessary due to the polar transformation applied to the original image. In this polar transformation, the y-axis of the transformed image corresponds to the angular coordinate in the polar system, which is defined with respect to a reference axis drawn from the center to the top edge of the original image. Therefore, a vertical flip of the transformed image is equivalent to a horizontal flip of the original image in Cartesian coordinates. Because the ground truth labels are assigned specifically to the left and right ring-shaped segments of the WALS image, this change in orientation necessitates an adjustment of the ground truth values to maintain the correct association of the labels. By reflecting a WALS image, we essentially create a mirror image of the droplet, i.e. as if it were captured in a different spot. Therefore, applying reflections during training improves the model's ability in detecting droplets consistently, no matter where they appear in the measurement area. For translation, we shift the image in both the x- and y-directions by a maximum of two pixels which corresponds to an angular change of approximately 1.5 degrees in the original WALS image. Importantly, all movements within this range are equally likely, ensuring a uniform distribution of translations for the data augmentation process. This method is chosen to train models that remain consistent despite minor variations within the experimental setup. This characteristic is crucial for the dependable and practical analysis of future experimental data.

During the training phase for each neural network, we employ a mini-batch gradient descent approach in order to update the trainable parameters. More specifically, for a given batch size $b \unicode{x2A7D} n$, we repeatedly select batches of b random images from their designated training sets and input them into the corresponding neural network, producing predictions. These predictions are then compared to the ground truth labels and the discrepancy is quantified via a specific loss function. Subsequently, the gradient of this loss function with respect to the trainable parameters is computed to refine the model. Based on this gradient, adjustments to the parameters are made using the Adam optimizer, as referenced in [8]. This optimizer adapts the learning rate for each parameter depending on the first and second moments of the gradients. This iterative process of batch selection, prediction generation, loss computation, and parameter updating is conducted for up to 150 epochs, each containing 100 batches. The maximum number of epochs was chosen due to computational limitations. For more details, please refer to section 2.3.5.

During network training, the quality of our predictions is continuously assessed. This assessment is primarily driven by the loss functions, which highlight the difference between the predictions and corresponding ground truth values. In the next section, the loss functions considered in this paper are described in detail.

To train our neural networks in order to accurately predict the binary labels $s_1,s_2$ and the quantitative descriptors of droplet size $s_3,s_4$, we need to consider suitably chosen loss functions. They measure the discrepancy between the predicted values $\hat{s}_j$ and the true values s_j , for $j = 1,2,3,4$, providing a feedback mechanism for the neural networks to adjust their parameters during training. Given the distinct nature of our two tasks—one of them performing classification and the other one performing quantitative scalar prediction—we employ different loss functions for each case. This allows us to tailor our approach according to the specific requirements of each task.

For training the droplet detection network, we use the mean absolute error (MAE) as loss function, where the MAE calculates the mean absolute difference between the first two components of the predicted labels $\hat s = (\hat{s}_1,\ldots,\hat{s}_4)$ and the corresponding ground truth labels $s = ({s}_1,\ldots,{s}_4)$. Specifically, for a set of b predictions $\hat{s}^{(1)},\hat{s}^{(2)},\ldots,\hat{s}^{(b)}$ and corresponding ground truth labels $s^{(1)},s^{(2)},\ldots,s^{(b)}$, the MAE for the jth component is given by

$$\begin{equation*} \text{MAE}_j = \frac{1}{b} \sum_{i = 1}^{b} |\hat{s}_j^{\left(i\right)} - s_j^{\left(i\right)}| \qquad\text{for } j\in\left\{1,2,3,4\right\}.\end{equation*}$$

Moreover, for the training of the droplet detection network, we consider both $\text{MAE}_1$ and $\text{MAE}_2$ as they relate to the binary labels for droplet detection on the left and right ring-shaped segments, respectively. Therefore, the training loss for the detection network is the average of $\text{MAE}_1$ and $\text{MAE}_2$, i.e. the value of $({\text{MAE}_1+\text{MAE}_2})/2$.

For the droplet size estimation network, using the MAE is not optimal because the relative impact of an absolute error is much more pronounced for smaller droplets. Therefore, we need a loss function that equally emphasizes the importance of errors across all droplet sizes during training. This can be achieved by a loss function which is based on the symmetric mean absolute percentage error (SMAPE). The SMAPE calculates the absolute percentage difference between predicted labels $\hat{s}$ and the ground truth s. Specifically, the SMAPE for the jth component is defined as

$$\begin{equation*} \text{SMAPE}_j = \frac{1}{b} \sum_{i = 1}^{b} \frac{2|\hat{s}_j^{\left(i\right)} - s_j^{\left(i\right)}|}{|\hat{s}_j^{\left(i\right)}| + |s_j^{\left(i\right)}|}\qquad \text{for } j\in\left\{3,4\right\}.\end{equation*}$$

As already mentioned above, the advantage of using the SMAPE for the droplet size estimation network is that it is a scale-invariant measure of accuracy. This means that it assigns equal importance to relative errors, regardless of the magnitude of the true values. This property is particularly useful when dealing with data that spans a wide range, such as droplet sizes.

However, as outlined in section 2.2, our approach for the generation of ground truth labels involves analyzing both ring-shaped segments of the WALS image separately. Consequently, a significant number of images in our dataset only has ground truth values for droplet sizes on one ring-shaped segment, e.g. $s_1^{(i)} = 1$ and $s_2^{(i)} = 0$, in which case we can calculate a meaningful size discrepancy for the left ring-shaped segment, but not for the right one. To effectively utilize this data as well, we employ a weighted SMAPE loss function, where the weights are used to exclude data points without ground truth droplet size values by multiplying the computed discrepancy with a weight of 0. This ensures that the SMAPE is computed only for the ring-shaped segments of the WALS image in which droplets were detected. As a result, this approach enables us to accurately evaluate the performance of our size estimation network on all available data. The weights of 1 and 0 correspond to the cases whether or not a droplet was detected in the ground truth. For example, if $s_1^{(i)} = 1$ then $s_3^{(i)}$ is a valid droplet size. Therefore, the weighted SMAPE ($\text{SMAPE}_{\text{w}}$) of the jth component is given by

$$\begin{equation*}\text{SMAPE}_{\text{w},j} = \frac{1}{b}\sum_{i = 1}^{b} s_{j-2}^{\left(i\right)}\frac{2 |\hat{s}^{\left(i\right)}_j - s^{\left(i\right)}_j|}{|\hat{s}^{\left(i\right)}_j| + |s^{\left(i\right)}_j|} \qquad\text{for } j\in\left\{3,4\right\}.\end{equation*}$$

Again, finally, we consider the average of the losses computed on the left and right ring-shaped segments, i.e. the value of $(\text{SMAPE}_{\text{w},3}+\text{SMAPE}_{\text{w},4})/2$.

This combination of the loss functions described above allows us to accurately measure the performance of the neural networks while taking into account all available data.

The goal of this section is to describe the methodology used to determine suitable model parameters, i.e. the number n_f of filters (see section 2.3) and trainable parameters such that the resulting neural network performs well on previously unseen data.

More specifically, we aim to determine a suitable choice for hyperparameters ${n}_{\text{f}} \in \{8,16, 32, 64\}$ and trainable parameters $\boldsymbol{\phi}\in \mathbb{R}^{d_{{n}_{\text{f}}}}$, where the dimension of φ depends on the choice of n_f, such that

$$\begin{align}\left({n}_{\text{f}}, \boldsymbol{\phi}\right) = \mathop{\text{arg min}}_{\left({n}_{\text{f}},\boldsymbol{\phi} \right)}\ \mathcal{L}_{{n}_{\text{f}}, \boldsymbol{\phi}; \mathcal{D}_{\text{test}}}. \end{align}$$

Here the $\text{arg min}$ is taken over all admissible parameters $({n}_{\text{f}},\boldsymbol{\phi})\in\{8,16, 32, 64\}\times\mathbb{R}^{d_{{n}_{\text{f}}}}$ and the value of

$$\begin{equation*}\mathcal{L}_{{n}_{\text{f}}, \boldsymbol{\phi}; \mathcal{D}_{\text{test}}} = \sum_{\left(I_\mathrm{p},s\right)\in\mathcal{D}_{\text{test}}} L\left(M_{{n}_{\text{f}}, \boldsymbol{\phi}}\left(I_\mathrm{p}\right),s\right)\end{equation*}$$

represents the loss on the test dataset $\mathcal{D}_{\text{test}}\subset\mathcal{D}_{\text{}}\setminus\mathcal{D}_{\text{train}}$ for the loss function $L,$ where $M_{{n}_{\text{f}}, \boldsymbol{\phi}}(I_\mathrm{p})$ denotes the output of the neural network with hyperparameter n_f and trainable parameters φ .

For that purpose, we iterate through hyperparameters ${n}_{\text{f}} \in \{8,16, 32, 64\}$ and obtain a suitable φ by applying the training procedure as described in section 2.3.3. This training proceeds for up to 150 epochs, continually updating the trainable parameters, resulting in different models $M_{{n}_{\text{f}}, \boldsymbol{\phi}e^{({n}_{\text{f}})}}$, where $\boldsymbol{\phi}e^{({n}_{\text{f}})}$ denotes the trainable parameters of the neural network that has been trained for e epochs on the dataset $\mathcal{D}_{\text{train}}$ with hyperparameter n_f. In order to evaluate the goodness-of-fit for each of the considered models at these epochs, we evaluate their performance on a validation dataset $\mathcal{D}_{\text{val}}\subset\mathcal{D}_{\text{}}\setminus(\mathcal{D}_{\text{train}}\cup\mathcal{D}_{\text{test}})$. More precisely, we compute the validation loss at the eth epoch $\mathcal{L}_{e,{n}_{\text{f}}}^{(\text{val})} = \mathcal{L}_{{n}_{\text{f}}, \boldsymbol{\phi}e^{({n}_{\text{f}})}; \mathcal{D}_{\text{val}}}.$ To optimize the training process and prevent overfitting, we employ an early stopping mechanism based on this validation loss as discussed in [17]. More specifically, if $\mathcal{L}_{e,{n}_{\text{f}}}^{(\text{val})}$ does not improve in three consecutive checks performed every fifth epoch, we terminate training. Formally, if $e\in\{20,15,\ldots,150\}$ is the smallest number of epochs such that

$$\begin{equation*}\mathcal{L}_{e-15,{n}_{\text{f}}}^{\left(\text{val}\right)}\unicode{x2A7D}\mathcal{L}_{e-10,{n}_{\text{f}}}^{\left(\text{val}\right)}\unicode{x2A7D}\mathcal{L}_{e-5,{n}_{\text{f}}}^{\left(\text{val}\right)}\unicode{x2A7D}\mathcal{L}_{e,{n}_{\text{f}}}^{\left(\text{val}\right)},\end{equation*}$$

then training is stopped at epoch e, denoted by $e_{\text{end},{n}_{\text{f}}} = e$. If no such epoch e exists, then training terminates at epoch 150 and $e_{\text{end},{n}_{\text{f}}} = 150$.

Finally, we select the hyperparameter n_f and corresponding trainable parameters φ such that

$$\begin{equation*}{n}_{\text{f}} = \mathop{\text{arg min}}_{{n}_{\text{f}}\in\left\{8,16, 32, 64\right\}}\mathcal{L}_{e_{{n}_{\text{f}}},{n}_{\text{f}}}^{\left(\text{val}\right)}\quad\quad\text{and }\quad\quad \boldsymbol{\phi} = \boldsymbol{\phi}{e_{{n}_{\text{f}}}}^{\left({n}_{\text{f}}\right)},\end{equation*}$$

where $e_{{n}_{\text{f}}}$ denotes the epoch for which the lowest validation loss has been obtained when training with the hyperparameter ${n}_{\text{f}},$ i.e. $e_{{n}_{\text{f}}} = \text{arg min}_{e\unicode{x2A7D} e_{\text{end},{n}_{\text{f}}}}\mathcal{L}_{e,{n}_{\text{f}}}^{(\text{val})}.$ Meaning that the model with the lowest validation error is chosen, indicating good performance and generalization capabilities.