Data augmentation for deep ensembles in polyp segmentation

The transformation consists of stretching the original image. The stretching direction is randomly chosen from the four possible options: left/up, right/up, left/down and right/down. The stretching is applied to the columns first, then the empty spaces are interpolated with a weighted nearest neighbor method, and the same procedure is repeated for the rows. The stretching of the original image is obtained by changing the position of the rows/columns according to the following equations:

$$\begin{eqnarray*}&&\lambda =U(2,17)\end{eqnarray*}$$

$$\begin{eqnarray*}&&k\left[i\right]=\frac{\frac{1}{\lambda }-\lambda }{s-1}t\left[i\right]+\lambda ,\,i\in \left[1,s\right]\end{eqnarray*}$$

$$\begin{eqnarray*}&&k\left[i\right]=\frac{\lambda -\frac{1}{\lambda }}{s-1}t\left[i\right]+\frac{1}{\lambda },\,i\in [1,s]\end{eqnarray*}$$

$$\begin{eqnarray*}&&a=a+k(\,\,j),\,j\,{\rm{\epsilon }}\,(2,s)\end{eqnarray*}$$

$$\begin{eqnarray*}y\left[\,j\right]=\left\{\begin{array}{ll}1 & \,j=1\\ \mathrm{round}\left(a\ast 5\right) & j\in [2,s]\,\end{array},\,\right.\end{eqnarray*}$$

where s is the number of columns and rows (in this work both equal 224), t[i] is a vector containing the numbers [0, s−1], k[i] is the distance between the row/column of the original image and the final image, while y[j] is the final row/column which varies according to j, the original position. The function U(a, b) returns a random number in the range a, b. An example of the application of the spatial stretch technique to an image of the training set and its mask is shown in figure 8.1.

Zoom In Zoom Out Reset image size

The final image is obtained by applying a shadow [36] to the left or the right of the original image, as shown in figure 8.2. In particular, the intensities of the columns are multiplied by the following equation:

$$\begin{eqnarray*}&&y=\left\{\begin{array}{c}{\rm{\min }}\left\{\mathrm{0.2\; +\; 0.8}\sqrt{\frac{x}{0.5}},\,1\right\}{\rm{direction}}=1\,\\ \,\\ {\rm{\min }}\left\{\mathrm{0.2\; +\; 0.8}\sqrt{\frac{1-x}{0.5}},\,1\right\}{\rm{direction}}=0\,\end{array}\right..\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

This transformation [36] is the composite of two alterations: first, it is necessary to modify the contrast of the original image, increasing or decreasing it, then a filter that simulates the movement of the camera is applied. Two functions to modify the contrast are implemented, but only one randomly chosen between the two is applied to the image.

The first contrast function (figure 8.3) is based on the following equation:

$$\begin{eqnarray*}&&y=\frac{\left(x-\frac{1}{2}\right)\sqrt{1-\frac{k}{4}}}{\sqrt{1-k{(x-\frac{1}{2})}^{2}}}+0.5,\,k\leqslant 4.\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

The parameter k controls the contrast, there is an increase in contrast if k < 0, a decrease in contrast if 0 < k ⩽ 4, the image remains the same when k = 0.

In the code k is chosen randomly from a specific range, among the following four:

• U(2.8, 3.8) → Hard decrease in contrast.
• U(1.5, 2.5) → Soft decrease in contrast.
• U(−2, −1) → Soft increase in contrast.
• U(−5, −3) → Hard increase in contrast.

The second contrast function (figure 8.4) is based on the following equation:

$$\begin{eqnarray*}y=\left\{\begin{array}{ll}\frac{1}{2}\left(\frac{x}{0.5}\right)\alpha & 0\leqslant x\lt \frac{1}{2}\\ \,\\ 1-\frac{1}{2}\left(\frac{1-x}{0.5}\right)\alpha & \frac{1}{2}\leqslant x\leqslant 1\end{array}\,\right..\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

The parameter α controls the contrast, there is an increase in contrast if α > 1, a decrease in contrast if 0 < α < 1 and the image remains the same when α = 1. This parameter is chosen randomly from four possible ranges:

• U(0.25, 0.5) → Hard decrease in contrast.
• U(0.6, 0.9) → Soft decrease in contrast.
• U(1.2, 1.7) → Soft increase in contrast.
• U(1.8, 2.3) → Hard increase in contrast.

This technique consists of three color operations (color adjusting, blurring and adding Gaussian noise) and a spatial operator (image rotation). The jitterColorHSV (I, Name, Value) function is used for color adjusting by randomly changing saturation, contrast and brightness levels within a random range:

• Saturation → [−0.3–0.1].
• Contrast → [−0.3–0.1].
• Brightness → [1.2 1.4].

The imgaussfilt function is used for blurring and the imnoise function is used for adding Gaussian noise. The image is finally rotated by an angle in the interval [−90° 90°] (figure 8.5).

Zoom In Zoom Out Reset image size

This technique consists of the segmentation of the image based on three different colors. The image obtained in this way is divided into three images, each of these containing a different color. The two images with the highest number of black pixels are selected, then the third image with the highest brightness is added to these. The images are combined and rotated by an angle selected from the range [−90°, 90°] (figure 8.6).

Zoom In Zoom Out Reset image size

This approach consists of the random selection of two transformations (one color-based and one shape-based) among a set of 21 [37].

Thirteen transformations belong to the color category:

• Create a composite image from two images A and B, by using the function imfuse.
• Add Gaussian noise.
• Adjust the saturation.
• Adjust the brightness.
• Control the contrast.
• Adjust the sharpness by using the function imsharpen.
• Application of the motion filter to blur the image.
• Histogram equalization.
• Conversion from RGB to YUV and histogram equalization.
• Application of the disk filter to blur the image.
• Add salt-and-pepper noise by using the function imnoise.
• Convert from RGB to HSV and adjust the hue by using the function jitterColorHSV(I,Name,Value).
• Adjust the local contrast by using the function localcontrast.

Eight transformations belong to the shape category:

• Rotation of an angle selected from the range [−40°, 40°].
• Vertical flip.
• Horizontal shear by using the function randomAffine2d.
• Vertical shear by using the function randomAffine2d.
• Horizontal translation by using the function randomAffine2d.
• Vertical translation by using the function randomAffine2d.
• Uniform scaling.
• Cutout by using the function imcrop.

The Rand augment technique (figure 8.7) consists of the application of a transformation from the color category with probability 1 and the subsequent application of a technique from the shape category with a probability of 0.3.

Zoom In Zoom Out Reset image size

The RICAP transformation [38] consists of three steps:

• Four images are randomly selected from the training set.
• The selected images are cropped.
• The cropped images are patched together to construct a new image (figure 8.8).

Zoom In Zoom Out Reset image size

Ten artificial images are created (figure 8.9) for each image in the dataset [39]. The operations performed are:

• The image is displaced to the right or the left.
• The image is displaced up or down.
• The image is rotated by an angle randomly selected from the range [0°, 180°].
• Horizontal or vertical shear is applied by using the function randomAffine2d.
• A horizontal or vertical flip is applied.
• Change the brightness levels by adding the same value to each RGB channel.
• Change the brightness levels by adding different values to each RGB channel.
• Add speckle noise by using the function imnoise.
• Application of the technique 'contrast and motion blur', described previously.
• Application of the technique 'shadows', described previously.

Zoom In Zoom Out Reset image size

Three occlusion methods are used. The first technique consists in selecting some rows (or columns) and replacing them with black lines (figure 8.10). The distance between the lines in terms of pixels is chosen randomly from the interval U(15, 30).

Zoom In Zoom Out Reset image size

In the second occlusion technique the images are transformed by replacing some parts of them with black rectangles (figure 8.11). The size of the rectangles, their number and their position are chosen randomly from the following intervals:

• Length of the rectangle n → U(4, 14).
• Height of the rectangle m → U(4, 14).
• Number of rectangles → U(7, 14).
• Coordinate x → U(1, 224 – n).
• Coordinate y → U(1, 224 – m).

Zoom In Zoom Out Reset image size

Given an input image, some of its pixels are removed, as shown in figure 8.12. The equation used is y = x × M, where x is the input image, M is a binary mask that stores pixels to be removed, while y is the final image.

Zoom In Zoom Out Reset image size

For each image x₁ of the training set an image x₂ from the original dataset is chosen randomly. The image x₂ is divided into a 7 × 7 grid, and N elements are cut from this grid. The elements cut from the second image are pasted on top of the first in their original position, as in figure 8.13. The equation used [40] is

$$\begin{eqnarray*}&&y=B\times x2+(1-B)\times x1,\end{eqnarray*}$$

where y is the result of the algorithm,while B is a binary mask (figure 8.13).

Zoom In Zoom Out Reset image size

For each image A of the training set, a second image B is selected randomly, resized and pasted over image A [41].Then the technique 'color change and rotation' is applied (figure 8.14).

Zoom In Zoom Out Reset image size

This technique consists of randomly selecting an image B from the training set for each original image A. Three methods of color normalization of image A versus image B are applied (figure 8.15), using the Stain Normalization toolbox by Nicholas Trahearn and Adnan Khan [42]:

• RGB histogram specification.
• Reinhard.
• Macenko.

Zoom In Zoom Out Reset image size

We have performed the experiments following two different protocols, which are used widely in the literature.

The first testing protocol [13, 43] uses the Kvasir-SEG dataset [44] partitioned into 880 images for training and the remaining 120 for testing. The Kvasir-SEG dataset includes 1000 polyp images acquired using a high-resolution electromagnetic imaging system (the image sizes vary between 332 × 487 and 1920 × 1072 pixels), with a ground-truth consisting of bounding boxes and segmentation masks.

The second testing protocol exploits five polyp datasets, namely Kvasir-SEG, CVC-ColonDB, EndoScene, ETIS-Larib Polyp DB and CVCClinic DB. The CVC-ClinicDB [45] and ETIS-Larib [46] datasets consist of frames extracted from colonoscopy videos, annotated by expert video endoscopists, and include 612 images (384 × 288) and 196 high-resolution images (1225 × 966), respectively. The CVC-ColonDB [47] dataset contains 380 images (574 × 500) representing 15 different polyps. EndoScene is a combination of CVC-ClinicDB and CVC300. The training set for the second testing protocol includes 900 images from Kvasir-SEG and 550 images from CVC-ClinicDB. The testing set is made up of the remaining images from the above-cited datasets: 100 images from Kvasir-SEG, 62 images from CVC-ClinicalDB, 380 from CVC-ColonDB, 196 from ETIS-Larib and 60 from CVC-T, which is the testing set for EndoScene. Note that only a small set of images are taken from EndoScene to avoid the inclusion of images already seen in the training stage.

All the datasets for both protocols can be downloaded from the GitHub of [43].

In all experiments, even when the images are resized to the input size of a CNN model, the predicted masks are resized back to the original dimensions before performance evaluation.We do not include other approaches in the comparison that evaluated the performance on the resized version of the images.

In accordance with most of the works on polyp segmentation we use pixel-wise metrics as performance indicators: accuracy, precision, recall, F₁-score, F₂-score, intersection over union (IoU) and Dice. A mathematical definition of each indicator for a bi-class problem (foreground/background), starting from the confusion matrix (TP, TN, FP and FN refer to the true positives, true negatives, false positives and false negatives, respectively), is the following:

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{accuracy}\,=\,\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}+\mathrm{TN}}\end{array},\end{eqnarray*}$$

which is the number of pixels correctly classified over the total number of pixels in the image,

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{precision}\,=\,\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}},\end{array}\end{eqnarray*}$$

which is the fraction of the polyps that are correctly classified,

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{recall}\,=\,\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\end{array},\end{eqnarray*}$$

which is the fraction of the predicted mask that is actually polyp pixels,

$$\begin{eqnarray*}&&\begin{array}{l}{{\rm{F}}}_{1}{\unicode{x02010}}\mathrm{score}\,=\frac{2\times \mathrm{TP}}{2\times \mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\end{array},\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}{{\rm{F}}}_{2}{\unicode{x02010}}\mathrm{score}\,=\frac{5\times \mathrm{precision}\times \,\mathrm{recall}}{4\times \mathrm{precision}+\,\mathrm{recall}}\end{array},\end{eqnarray*}$$

which are two measures that try to average precision and recall,

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{IoU}\,=\,\frac{\,\left|A\,\cap \,B\right|}{\left|A\,\cup \,B\right|}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\end{array},\end{eqnarray*}$$

which is defined as the area of intersection between the predicted mask A and the ground-truth map B, divided by the area of the union between the two maps, and

$$\begin{eqnarray*}&&\begin{array}{l}\text{Dice}\,=\,\frac{\,\left|A\,\cap \,B\right|}{\left|A{\rm{| }}+{\rm{| }}B\right|}=\frac{2\times \mathrm{TP}}{2\times \mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\end{array},\end{eqnarray*}$$

which is defined as twice the overlap area of the predicted and ground-truth masks divided by the total number of pixels (it coincides with the F₁-score for binary masks)

All the above metrics range in [0,1] and must be maximized. The final performance is obtained by averaging on the test set the performance obtained for each test image.

The first experiment in table 8.1 is aimed at comparing the two different backbone networks: ResNet18 and ResNet50. Since the size of images in the Kvasir dataset is quite large we also evaluate versions of the ResNet with larger input size, i.e. 299 × 299 (ResNet18–299/ResNet50–299) and 352 × 352 (ResNet18–352/ResNet50–352). Clearly, a larger input size improves the performance.

The training of all the models has been performed with the SGD optimizer for 20 epochs, with an initial learning rate of 10e-2, a learning rate drop period of 5 epochs and a drop factor of 0.2. In this experiment we use the same 'base' data augmentation of our previous paper [48], i.e. horizontal and vertical flip, 90° rotation.

The second experiment (table 8.2) is aimed at comparing the different data augmentation approaches proposed in this work. For the sake of computation time all the tests are performed using ResNet18 with an input size of 224 and compared to the 'base' data augmentation reported in the first rows of tables 8.1 and 8.2. The other rows of table 8.2 include the performance obtained using the different data augmentation approaches described in section 8.4 and the last three rows report the performance obtained by evaluating a combination of several data augmentation approaches.

Comb_1 → 'color and shape change' and 'color mapping';

Comb_2 → 'color and shape change' and 'color change and rotation';

Comb_3 → 'color and shape change' and 'color mapping' and 'color change and rotation'.

Clearly, the combination of different data augmentation approaches boosts the performance.

The third experiment (table 8.3) is aimed at designing effective ensembles by varying the activation functions. Each ensemble is the fusion by the sum rule of 14 models (since we use, in [48], 14 activation functions). For sake of computation time, this test is performed only in the Kvasir-SEG dataset augmented by the same method as table 8.1 (base data augmentation, first testing protocol). The ensemble name is the concatenation of the name of the backbone network and a string to identify the creation approach:

• act: each of the 14 networks of the ensemble is obtained by substituting its activation layers by one of the activation functions used in [48] (the same function for all the layers, but a different function for each network).
• sto: 14 stochastic models are generated, repricing the activation layers by randomly selected activation function used in [48] (which may be different for each layer).
• sel: ensembles of 'selected' stochastic models. The network selection is performed using three cross-validations on the training set among 100 stochastic models. The selection procedure is aimed at picking the most performing/independent classifiers to be added to the ensemble. For a fair comparison with other ensembles, we selected a set of 14 networks, which are finally fine-tuned on the whole augmented training set at a larger resolution.
• relu: 14 networks with standard architecture. All the starting models in the ensemble are the same, except for the initialization.

According to the results reported in table 8.3, it is clear that stochastic variation of activation functions (sto) allows a performance improvement with respect to a simple fusion of networks sharing the same architecture (relu) or a set of networks differing by the activation function (act). Such improvement is also noticeable if a selection procedure (sel) is used to include only the most performing architecture. Unfortunately, the selection process is computationally heavy, because it is performed among a pool of 100 networks. The best performance among the ensembles is obtained by ResNet50–352_sel, showing that the input size is critical in this problem.

Finally, in tables 8.4 and 8.5 a comparison with some state-of-the-art results is reported according to the above-cited first and second testing protocols, respectively.

In these final experiments, for sake of computational time, ResNet50 is trained using only the base data augmentation consisting of horizontal and vertical flip, and 90° rotation, while ResNet18 is trained using the here-proposed data augmentation 'color and shape change' and 'color mapping'.

The methods named HarDNet_SGD and HarDNet_Adam are our experiments using HardNet [43] trained by SGD and Adam optimizers (using the code shared by the authors). The reason for this choice is that the authors of [43] proposed a different optimizer for the two protocols.

In table 8.4, in addition to the methods presented in the above tables, we also report the performance of the following fusions by average rule (symbol ⊕):

• Res18&50 = ResNet18_352_sel ⊕ ResNet50_352_sel.
• HardNet_SGD/Adam ⊕ Res18&50, which is the fusion of HardNet (SGD or ADAM version) with the ensemble above.
• HardNet_SGD⊕HardNet_Adam⊕ Res18&50, which is the fusion of both HardNet (SGD and ADAM version) with our ensemble.

As far the latter protocol is concerned, since the training set is larger, to reduce computation time we have run the following methods without supervised selection:

• ResNet18_352, a stand-alone ResNet18 trained using the base data augmentation approach (only horizontal/vertical flip and 90° rotation).
• ResNet18_352_a, a stand-alone ResNet18 trained using the data augmentation 'color and shape change' and 'color mapping'.
• ResNet50_352, a stand-alone ResNet50 trained using base data augmentation.
• ResNet18_352_relu_a, an ensemble of ten standard ResNet18 networks trained using the data augmentation 'color and shape change' and 'color mapping'.
• ResNet50_352_relu, an ensemble of ten standard ResNet50 networks trained using base data augmentation.
• ResNet18_352_sto_a, an ensemble of ten stochastic ResNet18 networks trained using the data augmentation 'color and shape change' and 'color mapping'.
• ResNet50_352_sto, an ensemble of ten stochastic ResNet50 networks trained using base data augmentation.

Moreover, in table 8.5 the performance of the following fusions by the average rule are reported:

• Res18&50sto = ResNet18_352_sto_a ⊕ ResNet50_352_sto.
• HN1&R = HardNet_Adam⊕ Res18&50sto, which is the fusion of HardNet_Adam with the ensemble above.
• HN2&R = HardNet_SGD⊕HardNet_Adam⊕Res18&50sto, which is the fusion of both HardNet (SGD and ADAM version) with our ensemble.

Finally, for the sake of comparison in table 8.5, we report the performance of several state-of-the-art approaches: only methods following the second protocol exactly are included.

The following conclusions can be obtained from tables 8.4 and 8.5:

• From the comparison of ResNet18_352 (or ResNet18_352_a) and the ensembles ResNet18_352_relu_a (or ResNet18_352_sto_a) we can observe the advantage of fusing networks together.
• While in table 8.2 data augmentation granted a performance improvement for ResNet18, from the results in table 8.5 it seems that data augmentation is not helpful. The reason may be the different input sizes of the networks—by using larger input layers the performance improves reducing the need for data augmentation.
• The fusion of ensembles created from different models (ResNet18 and ResNet50) allows a performance improvement over each component.
• The best method for training the HarDNet model is Adam (considering all datasets from the second protocol).
• HarDNet_Adam has a slight advantage over our proposed ensembles (Res18&50 or Res18&50sto), nonetheless, the fusion of both approaches improves both performances.

This work shows the usefulness of ensembles. Our best method is HN2&R which is given by the fusion of HardNet_SGD, HardNet_Adam and the ensemble Res18&50sto. This method closes the gap with transformers, gaining performance almost comparable to that of transformers.