Distribution-preserving data augmentation

Extraction of features

DPDA method works with color images with three channels (red, green, and blue), namely RGB image. For an input RGB image in size W × H there are n pixels where n = WH. If we flatten the pixels of this RGB image, then our feature matrix X has a size of n × d where d = 3. Although n features are sufficient, we further enriched the feature space using the image pyramid approach (Adelson et al., 1984). During the image pyramid generation, we halved the original image in width and height three times. This creates an image pyramid with four levels where each level contains four times fewer pixels than the higher level in the pyramid. We used Lanczos interpolation over 8 × 8 neighborhood for the down-sampling operation (Turkowski, 1990). Using an image pyramid with four levels increases the number of features in X by 32.8%. In Fig. 5A, there are structural missing regions in feature space. This is due to quantization error since decimal parts of colors are quantized in 8 bits RGB images. However, refined feature space in Fig. 5B is denser, and the effect of image quantization errors is reduced, which demonstrates another benefit of the employed image pyramid approach.

Creation of density decreasing paths

Let x be the color of a pixel in the image as a starting point of a path we aim to find. Then probability density function (PDF) on the color feature space X constructed from the image is given in Eq. (1).

(1) $$\begin{array}{c}P(x)={\displaystyle \frac{1}{n}\sum _{j=1}^{n}K(x-x_j)}\hfill \end{array}$$

where K(.) is a kernel function and x_j are data points in X where we used Epanechnikov kernel.

We define a density decreasing path $T=\{x^{(0)},x^{(1)},\dots ,x^{(i)},\dots \}$ (Fig. 6) where x^{(i + 1)} = x⁽ⁱ⁾ + s⁽ⁱ⁾ for i ≥ 0. Here, x⁽⁰⁾ is the starting point of the path and s⁽ⁱ⁾ is a density decreasing direction at point x⁽ⁱ⁾. Also, x⁽ⁱ⁾ is only defined in color space domain where 0 ≤ x⁽ⁱ⁾_r, x⁽ⁱ⁾_g, x⁽ⁱ⁾_b ≤ 255.

Now, we can define the pdf for the point x^{(i + 1)} as below:

(2) $$\begin{array}{c}P(x^{(i+1)})={\displaystyle \frac{1}{n}\sum _{j=1}^{n}K(x^{(i+1)}-x_j)={\displaystyle \frac{1}{n}\sum _{j=1}^{n}K(x^{(i)}+s^{(i)}-x_j)}}\hfill \end{array}$$

Since the points x⁽ⁱ⁾ and x_j are constant, we can rewrite above equations as below:

(3) $$\begin{array}{c}P(x^{(i+1)})=P(x^{(i)}+s^{(i)})={\displaystyle \frac{1}{n}\sum _{j=1}^{n}K(s^{(i)}-{\hat{x}}_j)where{\hat{x}}_j=x_j-x^{(i)}}\hfill \end{array}$$

So, ${\hat{x}}_j$ points are centered to x⁽ⁱ⁾ where x⁽ⁱ⁾ is shifted to the origin; thus, s⁽ⁱ⁾ becomes a direction vector. Finally, we define a gradient descent direction as below that will lead to a density decreasing path:

(4) $$x^{(i+1)}=x^{(i)}-\mathrm{\nabla }J(x^{(i)}+s^{(i)})$$

where J(x⁽ⁱ⁾ + s⁽ⁱ⁾) is the cost function to minimize which is defined as below:

(5) $$\begin{array}{cc}J(x^{(i)}+s^{(i)})\to P(x^{(i)}+s^{(i)})\hfill & subjectto{s}^{(i)}\in \mathrm{\Omega }\hfill \\ \mathrm{\Omega }=\{\parallel s^{(i)}\parallel \le S_{length}\hfill & and1-⟨{\displaystyle \frac{s^{(i)}}{S_{length}},{\hat{s}}^{(prior)}⟩\le S_{angle}\}}\hfill \\ \hfill & with{\hat{s}}^{(i)}={\displaystyle \frac{s^{(i)}}{S_{length}}and\hat{s}(prior)={\displaystyle \frac{s^{(i-1)}}{\parallel s^{(i-1)}\parallel }}}\hfill \end{array}$$

In Eq. (5), S_length is the maximum length and S_angle is the maximum angle for the direction vector s⁽ⁱ⁾. Here, second constraint is only defined for i > 0 where ${\hat{s}}^{(prior)}$ is the prior direction in unit length and considered as constant. Owing to Eq. (5), density decreasing direction s⁽ⁱ⁾ is regularized in length and orientation to avoid chaotic shifts in sparse data regions.

The gradient descent method is not practical to minimize cost functions with constraints. In such cases, one can use the projected gradient descent (PGD) method, which can minimize a cost function subject to a constraint where this constraint defines a domain (Boyd & Vandenberghe, 2004). Although PGD works fine for cost functions with a single constraint, we can still use it efficiently since our cost function’s length, and orientation constraints only form a single domain, Ω. Therefore, we use PGD to obtain a gradient descent path on the cost function J as defined in Eq. (6).

(6) $$\begin{array}{cc}\underset{s^{(i)}}{min}P(x^{(i)}+s^{(i)})\hfill & \text{subject to}{s}^{(i)}\in \mathrm{\Omega }\hfill \\ \hfill & x^{(i+1)}={\mathcal{P}}_{\mathrm{\Omega }}(x^{(i)}-\mathrm{\nabla }P(x^{(i)}+s^{(i)}))\hfill \\ \hfill & {\mathcal{P}}_{\mathrm{\Omega }}(x_{new})=\underset{s^{(i)}\in \mathrm{\Omega }}{argmin}\parallel (x^{(i)}+s^{(i)})-x_{new}\parallel \hfill \\ \hfill & \hfill \end{array}$$

After doing some algebraic manipulations one can see that s⁽ⁱ⁾ equals to the opposite of the mean-shift direction m⁽ⁱ⁾ such that s⁽ⁱ⁾ = −m⁽ⁱ⁾. First, we will limit the number of iterations in the PGD to L/2 since we aim to find a density decreasing path with a limited number of points. Next, we will stop the PGD iteration if (a) the norm of mean-shift direction is becoming smaller than a tolerance value (C_tolerance) or (b) next point x^{(i + 1)} exiting from the image color space domain. Also, we set default value for C_tolerance as 10⁻²d. The final density decreasing path generation method is presented in Algorithm 1.

Calculation of the mean-shift direction m⁽ⁱ⁾ at point x is as given as below:

(7) $$m^{(i)}={\displaystyle \frac{\sum _{x_j\in \mathcal{N}(x)}K(x_j-x)x_j}{\sum _{x_j\in \mathcal{N}(x)}K(x_j-x)}-x}$$

where K(.) is the kernel function and x_j are k nearest neighbours of x. We used FLANN proposed by Muja & Lowe (2014) to have fast k nearest neighbor search operations for efficiency. In this study, we used 256 as the default value of k. Note that, one need to put value of x⁽ⁱ⁾ into the point x in the Eq. (7) given in the Algorithm 1. However, kernel functions require the selection of the bandwidth parameter h. Since each image’s characteristic is different, we estimate bandwidth parameter h from the image to balance differences between images as an approximation to median pair-wise distances to closest points. First, we find the Euclid distances of each pixel with its 4 neighbors. Then, we use Quick Select (Cormen et al., 2009) algorithm to find the median value of these distances as our bandwidth h.

We used the PGD method, which first does a gradient descent step, then back-projection of gradient descent result to the domain Ω. In Fig. 7, blue vectors are prior directions; green vectors are new directions in the domain, and red vectors are new directions out of the domain. Domain Ω is the region between dotted gray lines, determined by the constraints in each gradient descent step.

Note that prior direction and new direction form a plane where its normal is the cross product of these two vectors. Therefore, a rotation matrix can be formed, which aligns this normal vector to the canonical z-axis where the prior direction and new direction vectors transform onto xy-plane. Once prior and new directions are rotated, all the back-projection operations can be done in 2D easily then back-projected direction can be rotated back to the original space. We used the method proposed by Möller & Hughes (1999) to construct a rotation matrix that aligns normal vector to z-axis as given in Eq. (8).

(8) $$\begin{array}{cc}v=f\times t,u\hfill & ={\displaystyle \frac{v}{\Vert v\Vert }}\hfill \\ c=f\cdot t,r\hfill & =(1-c)/(1-c^2)\hfill \end{array}\to \begin{array}{c}\left[\begin{array}{ccc}c+r{v}_x^2& r{v}_x{v}_y-v_z& r{v}_x{v}_z+v_y\\ r{v}_x{v}_y+v_z& c+r{v}_y^2& r{v}_y{v}_z-v_x\\ r{v}_x{v}_z-v_y& r{v}_y{v}_z+v_x& c+r{v}_z^2\end{array}\right]\hfill \end{array}$$

where f is plane normal calculated by cross product of prior direction and new direction, and t is z-axis.

Creation of augmented images

Each pixel has its corresponding density decreasing path with L colors. 0^th color (first path node) has the largest color deviation from the original pixel color towards the tail of image color distribution. (L − 1)^th color (last path node) equals to original image pixel color. So, we can take a different node (color) from the corresponding path for each pixel of an augmented image. Here all 0 indices will yield to augmented image with the most perturbation, while L − 1 indices will yield to the original image. For each pixel, we can randomly choose an index number between 0 and L − 1, which will lead to different augmented images that allow the generation of any number of augmented images. However, utterly random selection will result in unnatural results. Thus, we want to sample from path nodes in a random but spatially smooth manner. We modified the Perlin noise generator, which is proposed by Perlin (1985) to obtain a smooth but random index map as seen in Fig. 8. Here, we choose parameter values randomly from a predefined range where parameters are roughness (N_roughness), noise scale (N_scale), and noise center (N_center). Finally, modified Perlin noise is generated using C_x,y = 0.5 (tanh(N_scale * (N_x,y − N_center)) + 1) where N_x,y = Perlin.generate(xN_roughness, yN_roughness, 1) is original Perlin noise generation function.

Datasets

We used Pxfuel¹ for qualitative experiments, and three different datasets for quantitative experiments, namely the UC Merced Land-use Yang & Newsam (2010), the Intel Image Classification (https://www.kaggle.com/puneet6060/intel-image-classification/), and the Oxford-IIIT Pet datasets (Parkhi et al., 2012).

UC Merced Land-use dataset consists of satellite images of size 256 × 256 and 0.3-m resolution that are open to the public. There are a total of 21 classes and 100 images in each class (see Fig. 9).

The Intel Image Classification dataset contains about 25,000 images of size 150 × 150 pixels, classified under six categories (buildings, forest, glacier, mountain, sea, and street) of natural scenery worldwide (see Fig. 10). The training set is around 17,000 images and the test set is the rest.

The Oxford-IIIT Pet dataset includes 37 dog and cat classes with 25 dog and 12 cat categories. There are approximately 200 images for each class with significant variations in image size, exposure, and lighting. The total number of images in the dataset is just over 7,000. The dataset also includes labels as bounding boxes and segmentation masks as trimaps. In trimap, absolute background is shown as white, absolute foreground is shown as black, and mixed region is shown as gray (see Fig. 11).

Qualitative results

To demonstrate our data augmentation results qualitatively, we first downloaded sample images from Pxfuel, which provides high-quality royalty-free stock photos. Note that augmented images’ brightness is slightly increased to emphasize the difference between original images and augmented images.

In Fig. 12, augmentation results are shown for man, forest, food, car, and urban images. In the first row, in the augmented male image, the male’s skin color becomes lighter, and the eye color shifts to green. Also, there are some color changes in the background and t-shirt of the man. In the second row, the augmented forest image contains red trees, although there are no red trees in the original image. Here, red trees occur in the augmented image because the original image’s color distribution contains colors towards the red tones in the distribution’s tail. In the third row, in the augmented image, each olive type’s colors in the original image are changed differently while the background is not changed. In the fourth row, the old car’s rust tones in the original image are changed naturally in the augmented image. In the fifth row, the trees and the building roof’s colors become greener with slight color changes in the buildings and the road in the augmented urban image.

DPDA results shown in Fig. 12 are all plausible image augmentations. In all these visible results, some image pixel colors are shifted to the tail of the image’s color data distribution. Thus, the image is transformed into a less occurring version of itself. This is quite useful to increase data variability of the training dataset since the proposed data augmentation approach generates fewer occurring images, and thus original dataset is enriched. Therefore, the over-fitting problem is reduced while increasing the training accuracy. Since the image color data guides data augmentation, the algorithm does not require different parameter selections for different images, i.e., images with different content, resolution, or camera characteristics. Accordingly, default DPDA parameters are used for the data augmentations in Fig. 12 (as qualitative experiments) and also for all the quantitative experiments.

Quantitative results

Training a DL network from scratch requires a considerable amount of data and computational power. Therefore, researchers and practitioners with limited data and computational resources prefer to reuse existing DL architectures, which are trained with millions of data and using server farms. This reuse methodology employs a transfer learning approach where a well-proven DL model is fine-tuned with a new dataset (Shao, Zhu & Li, 2015). Pre-training a DL network with transfer learning yields successful results, even with a small train dataset. However, transfer learning provides excellent results if the data and pre-trained model are on a similar domain (Yosinski et al., 2014). A model pre-trained with the Imagenet dataset gives better outcomes for the datasets in the same domain, such as CIFAR-10 or Caltech-101. On the other hand, if the model is tuned using a small amount of training data that is not in a similar domain, the performance benefits of transferring features decrease. So, data augmentation helps increase dataset size and variety to remedy such problems (Shao, Zhu & Li, 2015).

Note that our aim is not to give an extensive study of the architecture of CNNs as done by Szegedy et al. (2015), or He et al. (2016) but to briefly use them for evaluating the performance of the proposed DPDA method in transfer learning settings. Resnet50 (He et al., 2016) and DenseNet201 (Huang, Liu & Weinberger, 2016) network weights trained on the ImageNet are used as starting weights in the classification task since they are widely used in the current studies (Khan et al., 2020). Then the models are fine-tuned during training (Vrbančič & Podgorelec, 2020) since initial layers of CNNs preserve more abstract, generic features. We just copy the weights in convolutional layers rather than the entire network, excluding fully connected layers. MobileNetV2 (Sandler et al., 2018) network weights trained on the ImageNet are used as starting weights in segmentation task as a base model and trained with CNN architecture based on U-Net (Silburt et al., 2019). For all the experiments, we used an Stochastic Gradient Descent (SGD) solver with a momentum of 0.9. Weights are initialized from a Gaussian distribution $\mathcal{N}(\mu ,\sigma )$ for μ = 0 and σ = 10⁻². We found 20 epochs and a batch size of 32 typically sufficient for convergence.

The following methodology was utilized to create train and validation sets for all datasets used in this study. First, we randomly selected 20 images from each class as a validation set and used the same validation set in all tests. Then, we created different train sets in various sizes (N = 20, 30, 40, 50, 60, 70, 80) to investigate the effect of training dataset size on classification performance using the data augmentation approaches. For segmentation tests, we only used the train set size of 80. To avoid sample imbalance in the training datasets, we randomly selected the training datasets in equal numbers from each class. We evaluated the final classification performance of each dataset with the average accuracy over 10 runs. We increased the original training dataset size 5-fold, utilizing random erase (RE), flip image (FI), gamma correction (GC), histogram equalization combined with gamma correction (HE+GC), the proposed DPDA method, and the DPDA method combined with the flip image (DPDA+FI) separately. We implemented color-based augmentation methods as done in CLoDSA (Casado-Garca et al., 2019) library.

UC merced land-use dataset

First, we compare the performance increase obtained with various data augmentation methods and DPDA on the UC Merced Land-use dataset using DenseNet201 architecture. As seen in Table 1, average accuracy improvement ranges from 1.51% to 4.43%. All data augmentation methods provide performance increase compared to baseline performance for all the train set sizes. However, the DPDA method and the DPDA combined with the flip image consistently provide the best performances in every test. The results also show that data augmentation in data sets with fewer elements contributes more to accuracy. For example, the highest accuracy increase is 6.98% in the training set consisting of 20 images per class, which is obtained with DPDA+FI augmentation.

Next, we compared the performance increase with various data augmentation methods, including the DPDA, using ResNet50 architecture. As seen in Table 2, average accuracy improvement ranges from 2.35% to 5.84%. All data augmentation methods provide performance increase compared to baseline performance for all the train set sizes. However, the DPDA method itself and in combination with the flip image method, consistently provide the best performances in every single test. The results also show that data augmentation in data sets with fewer elements contributes more to accuracy. For example, the highest accuracy increase is 8.49% in the training set consisting of 20 images per class, which is obtained with DPDA+FI augmentation.

Intel image classification dataset

Like the UC Merced Land-use dataset, first, we compare the performance increase obtained with various data augmentation methods, including DPDA, on the Intel Image Classification dataset using DenseNet201 architecture. As seen in Table 3, average accuracy improvement ranges from 2.25% to 4.94%. All data augmentation methods provide performance increase compared to baseline performance for all the train set sizes. However, the DPDA method provides the best performances in every single test. The results also reveal that data augmentation in data sets with fewer elements contributes more to accuracy. For instance, the highest accuracy increase is 7.23% in the training set consisting of 30 images per class, which is obtained with DPDA augmentation.

Next, we compare the performance increase with various data augmentation methods, including the DPDA, using ResNet50 architecture. As seen in Table 4, average accuracy improvement ranges from 1.82% to 4.80%. Every data augmentation methods provide a performance increase compared to baseline performance for all the train set sizes. However, the DPDA method provides the best performances in every single test. The results also indicate that data augmentation in data sets with fewer elements contributes more to accuracy. For instance, the highest accuracy increase is 7.27% in the training set consisting of 20 images per class, which is obtained with DPDA augmentation. This result is also in compliance with the DenseNet comparison study.

Oxford-IIIT pet dataset

First, we compare the performance increase with various data augmentation methods, including DPDA, on the Oxford-IIIT Pet dataset using DenseNet201 architecture. As seen in Table 5, average accuracy improvement ranges from 2.34% to 3.34%. All data augmentation methods provide performance increase compared to baseline performance for all the train set sizes while DPDA being superior. The results also reveal that data augmentation in data sets with fewer elements contributes more to accuracy. For example, the highest accuracy increase is 6.68% in the training set consisting of 20 images per class, which is obtained with DPDA+FI augmentation.

Next, we compare the performance increase with various data augmentation methods, including the DPDA, using ResNet50 architecture. As seen in Table 6, average accuracy improvement ranges from 3.03% to 6.33%. Every data augmentation methods provide a performance increase compared to baseline performance for all the train set sizes. However, the DPDA+FI method provides the best performances in every single test. The results also show that data augmentation in data sets with fewer elements contributes more to accuracy. For instance, the highest accuracy increase is 12.66% in the training set consisting of 20 images per class, which is again obtained with DPDA+FI augmentation.

We used the U-Net architecture on top of the MobileNetV2 architecture for the Oxford-IIIT Pet Dataset segmentation experiments. DPDA provides by far the best performance in these experiments (see Table 7). The accuracy improvement obtained with DPDA is 8.49%. The second highest accuracy improvement achieved with the FI augmentation is 0.81%, which is much less than the DPDA accuracy improvement. On the other hand, every data augmentation method does not provide a performance increase compared to baseline performance. For instance, the GC decreases the accuracy by −2.20%.

The above results indicate that models trained with the augmented UC Merced Land-use, Intel Image Classification, Oxford-IIIT Pet datasets, with the DPDA method, significantly improve classification performance. Besides, DPDA also provides superior performance in an image segmentation task. Thus, we can infer that the proposed DPDA method can improve DL performance for different datasets, different DL architectures (ResNet, DenseNet, and MobileNetV2), and different image analysis tasks.

Execution time analysis

There are n pixels in an image. During the execution of DPDA, for each pixel, we find a path with a length of up to L/2. For each point in a path, the nearest neighbor search using FLANN method (Muja & Lowe, 2014) is done to retrieve k neighbor points. Our data is 3 channel so dimension d is 3. For each image, the FLANN tree is constructed for once where tree construction has a computational complexity of $\mathcal{O}(ndKI(logn/logK))$, where I is a maximum number of iterations, and K is the branching factor. We used exact search in FLANN, which leads to $\mathcal{O}(Md(logn/logK))$ for single nearest neighbor search where M is a maximum number of points to examine. However, we need to do a separate neighbor search for L/2 times for n pixels, which leads to nL/2 neighbor search operations. Thus, computational complexity of the all neighbour search operations is $\mathcal{O}(nLMd(logn/logK))$. Therefore, computational complexity of the DPDA is $\mathcal{O}(nd(KI+LM)(logn/logK))$ including tree construction and neighbor search operations.

The average execution time for 10 augmentations of DPDA, RE, GC, and FI methods concerning image size (# of pixels) is shown in Fig. 13. Although FLANN provides efficient nearest neighbor search operations, as shown in Fig. 13, the execution times of the DPDA method are longer compared to RE, GC, and FI methods. Fortunately, in DL training, images are generally in small sizes, i.e. 435 × 387 for Oxford dataset, 250 × 250 for UCMerced dataset, and 150 × 150 for Intel dataset (see Fig. 13).