# Compressing Vision Transformers for Low-Resource Visual Learning

Eric Youn
&Sai Mitheran J
&Sanjana Prabhu
&Siyuan Chen &Carnegie Mellon University
{eyoun, sjagades, sprabhu2<EMAIL_ADDRESS>

##  1 Motivation

Vision transformer (ViT) and its variants have swept through visual learning
leaderboards and offer state-of-the-art accuracy in tasks such as image
classification, object detection, and semantic segmentation by attending to
different parts of the visual input and capturing long-range spatial
dependencies. However, these models are large and computation-heavy. For
instance, the recently proposed ViT-B [1] model has 86M parameters making it
impractical for deployment on resource-constrained devices. As a result, their
deployment on mobile and edge scenarios is limited. In our work, we aim to
take a step toward bringing vision transformers to the edge by utilizing
popular model compression techniques such as distillation, pruning, and
quantization.

Our chosen application environment is an unmanned aerial vehicle (UAV) that is
battery-powered and memory-constrained, carrying a single-board computer on
the scale of an NVIDIA Jetson Nano with 4GB of RAM. On the other hand, the UAV
requires high accuracy close to that of state-of-the-art ViTs to ensure safe
object avoidance in autonomous navigation, or correct localization of humans
in search-and-rescue. Inference latency should also be minimized given the
application requirements. Hence, our target is to enable rapid inference of a
vision transformer on an NVIDIA Jetson Nano (4GB) with minimal accuracy loss.
This allows us to deploy ViTs on resource-constrained devices, opening up new
possibilities in surveillance, environmental monitoring, etc. Our
implementation is made available at https://github.com/chensy7/efficient-vit.

##  2 Related Work

Knowledge distillation is a technique that improves the performance of a
smaller student model by using "soft" labels from a larger teacher model: this
does not compress the model in any dimension of the network [2]. Soft labels
are thought to be more informative than hard labels, which can benefit student
training. Soft labels effectively transfer the knowledge of the teacher to the
student model [3].

Pruning techniques are classified into two types: unstructured [4] and
structured [5]. Unstructured pruning entails removing insignificant weights
based on specific criteria, resulting in irregular sparsity that makes
hardware acceleration difficult. Structured pruning addresses this by
assigning a weight to groups of parameters, such as convolutional channels or
matrix rows [6]. Sensitivity can be determined using a variety of metrics,
including batch normalization layer scaling factors, channel-wise summation
over weights, channel rank, and the geometric median of convolutional filters.
Common structures to prune in transformer-based models include blocks,
attention heads, and fully-connected matrix rows [7]. Several studies have
proposed various methods for pruning redundant heads, entire layers, and
extracting shallow models during inference [8].

Quantization-Aware Training (QAT) [9] and Post-Training Quantization (PTQ)
[10] are the two most common quantization methods. To achieve highly
compressed quantization, QAT relies on training, which frequently necessitates
significant expert knowledge and GPU resources for training or fine-tuning.
PTQ has gained popularity as a way to reduce the costs associated with
quantization since additional training is not required. Many excellent works
have resulted from PTQ, including OMSE [11], which reduces quantization errors
by determining the value range of activation, and AdaRound [12], which adapts
data and task loss via a novel rounding mechanism. Liu et al. [13] proposed a
post-training quantization method for vision transformers that employs
similarity-aware and rank-aware strategies, but it excludes Softmax and
LayerNorm modules, resulting in incomplete quantization.

One promising direction to reduce the computational complexity of the self-
attention mechanism in Vision Transformers (ViT) while maintaining or even
improving their performance is the use of low-rank approximation of attention
matrices. This can significantly reduce the number of computations needed for
each self-attention operation. Several studies have investigated different
variants of low-rank attention mechanisms, such as structured and unstructured
low-rank approximations, tensor decompositions, and dynamic rank selection
[14, 15, 16, 17]. These methods have been shown to be effective in reducing
the computational cost of ViT while achieving competitive performance on
various image classification benchmarks.

Efficient/sublinear-regime attention mechanisms in the Natural Language
Processing domain attempts to reduce the quadratic complexity of standard
Multi-Head Self-Attention (MSA). To accomplish this, various methods have been
developed, including low-rank decomposition [18], kernelization [19], memory,
and sparsity mechanisms [20]. These methods, however, have shown drawbacks in
Computer Vision (CV) tasks, where Spatial Reduction Attention (SRA) [21] and
Local Windowed Attention [22] are used instead. However, these mechanisms only
pay attention to either local or global issues. Some approaches addressed this
limitation by introducing additional global tokens or combining local window
attention with depthwise convolutional layers [23]. Other attention mechanisms
consider both local and global attention [24], but due to inefficient
operations, they may be slow on GPUs.

Figure 1: Our framework for edge-oriented, rapid on-device disaster scene
parsing

##  3 Methodology

###  3.1 Background and Formulation

To meet the requirements of utilizing no more than 4GB of RAM system-wide, we
incorporate most of the above approaches to achieve a compressed model with
Swin Transformer level accuracy as shown in Figure 1. In order to balance both
memory usage and inference latency improvements, we propose the use of
structured pruning to remove redundant heads from the attention mechanism and
neurons from the linear layers. [25] observed that many traditional
distillation techniques from CNNs do not carry over to vision transformers.
They propose a mimicking and generation-based distillation approach. We
utilize a similar approach to our distillation framework, using logit
knowledge transfer at the prediction stage, coupled with mimicking feature
distillation at early layers and generation-based feature distillation in
later layers. Our pruning technique focuses on convolution layers, linear
layers, and entire heads, but the actual attention matrices are not
compressed. Therefore, we can potentially build on this and use low-rank
approximation to further compress the model. Our final avenue for compression
leverages post-training quantization which is quite limited in PyTorch due to
the lack of robust support for quantized CUDA operations. Nevertheless, we
quantize our model to the maximum possible extent to reduce memory
requirements while keeping the performance at level.

Our target is ambitious - a tiny Swin Transformer has at least 40M parameters
while we estimate the maximum parameter size for a model running on the Jetson
Nano 4GB is about 5M. However, by applying the approaches specified above, we
push the accuracy of the model deployed on the edge device as close as
possible to the state-of-the-art level (Swin Transformer). Our task involves
segmentation of search-and-rescue categories and therefore is more specialized
and easier than datasets such as ADE20K or COCO. Therefore, a smaller network
will intuitively be able to capture the information across these classes.

###  3.2 Relation to prior work

As mentioned in the previous sections, much of the knowledge for model
compression in previous architectures and NLP-based transformers does not
translate well to vision transformers. Additionally, much of the prior work in
this area focuses only on certain compression techniques rather than unified
models combining various techniques together. One notable exception is the
work by Yu et al. [26] which proposes a unified framework using a structured
pruning approach to compress the attention module, knowledge distillation, and
using skip connections as a mechanism to drop entire blocks. However, they do
not include quantization and only conduct distillation at the prediction
logits stage. Similarly, certain works [27] focus on a combined structured
pruning and low-rank approximation framework but do not include distillation
for further compression. Our contribution to this line of research is to
combine modern techniques spread throughout the literature into a unified
compression framework and apply this approach to an ambitious real-world
problem: our on-device compression framework is applied for disaster scene
parsing, where the goal is to segment and parse images of UAV-view disaster
scenes using a model that runs efficiently on a Jetson Nano device. We test
our compression framework by performing rapid inference using the transformer
model on an NVIDIA Jetson Nano 4GB, finding the sweet spot between performance
and compute requirements.

###  3.3 Dataset

Figure 2: Samples from the LPCV dataset (CVPR 2023 challenge)

Our work utilizes the LPCV 2023 challenge dataset for the on-device disaster
scene parsing task. This dataset, which is already publicly available online,
is a UAV-view disaster scene dataset consisting of 1,720 images collected from
various disaster scenes.

Out of the total images, we use only 1,120 images along with their ground-
truth labels for training our model and testing its performance, while the
remaining images are used for benchmarking submissions on the challenge
server. The images are densely annotated with 14 categories, including
Background, Avalanche, Building_undamaged, Building_damaged,
Cracks/fissure/subsidence, Debris/mud/rock_flow, Fire/flare,
Flood/water/river/sea, Ice_jam_flow, Lava_flow, Person, Pyroclastic_flow,
Road/railway/bridge, and Vehicle, all resized to a resolution of 512x512
pixels.

Our use of this dataset will aid in the development of more accurate and
efficient disaster response systems and contribute to the existing body of
research in the field of building efficient, performant systems for on-the-
edge use cases, such as in the case of natural disasters and remote area
access.

###  3.4 Evaluation Metrics

Our work is inspired by the 2023 LPCV Challenge), for CVPR. Our target task is
dense prediction, i.e., scene parsing or semantic segmentation. Following the
challenge, we plan to evaluate our framework using the following metrics: for
accuracy, we use mean IoU, as shown below ($\mathrm{n}$ is the total number of
test images; TP, FP, and FN correspond to the True Positive, False Positive,
and False Negative pixels). Then, execution time is measured as the average
inference time over all test images. Our final score is computed as accuracy
per latency, as shown in the following formulae.

mean IoU
$\displaystyle=\frac{\sum_{i=1}^{n}\left(\frac{TP_{i}}{TP_{i}+FN_{i}+FP_{i}}\right)}{n}$
(1)

Final Score $\displaystyle=\frac{\text{ mean IoU }}{\text{ Execution Time }}$
(2)

In addition, we measure model size given general memory constraints on the
edge. Although we have set a clear upper bound for memory (4GB), these metrics
tell us if our solution can be extended to even more resource-constrained
scenarios.

Furthermore, we conduct an ablation study to understand how each compression
method (distillation, pruning, quantization) contributes to the accuracy
versus latency tradeoff.

###  3.5 Dataset and Software

We use PyTorch for training and inference procedures. While faster inference
engines are available, they are mostly orthogonal to our solution. We use the
LPCV dataset for knowledge distillation, quantization-aware training, and
finetuning after pruning. Our implementation is made available at
https://github.com/chensy7/efficient-vit.

###  3.6 Baselines

We use Swin Transformer [22] and MobileViT [28] as our primary encoder
backbones. The Swin Transformer is well-designed with unique features such as
the use of shifted windows to compute the representation, which allows it to
capture more global information while maintaining local details. MobileViT is
tailored to be used in mobile devices that have limited memory and
computational resources. It accomplishes this by replacing fully connected
layers with convolutional layers in the MLP blocks of the standard ViT
architecture. This strategy reduces the number of parameters in the model and
results in a more lightweight model, without significantly compromising
performance.

For the segmentation decoder, our three candidates are FANet [29], UPerNet
[30] and DeepLabv3 [31]. FANet is a segmentation model that can perform real-
time segmentation and it achieves this by using a feature aggregation network
to produce high-level features, which are then used in a feature upsampling
network for pixel-level predictions. UPerNet, on the other hand, enhances the
resolution of the predictions by bringing together different levels of
features using a multi-level feature fusion module and a pyramid pooling
module. Finally, DeepLabv3 is a cutting-edge model that captures multiscale
context using atrous convolution and incorporates low-level features using
skip connections.

Table 1: Experimental results for different backbones on semantic segmentation Backbone | Params | Mean IoU
---|---|---
MobileViT | 13.4M | 37.06
Swin-v2-T | 60M | 44.51

Results of experiments conducted on ADE20K, a common benchmark for semantic
segmentation, using the aforementioned networks as backbones have been
reported in Table 1. Note that for a difficult dataset like ADE20K, a
difference of 7 in mean IoU is considered very large, and we intend to deploy
a MobileViT on the Jetson Nano with similar accuracy as Swin Transformer by
applying our compression methods.

##  4 Results and Analysis: LPCV Challenge Dataset

Given the large volume of existing network architectures for semantic
segmentation, our first step is to identify a segmentation decoder
architecture with the best trade-off between efficiency and accuracy on the
LPCV dataset. In addition, by comparing the accuracy of different backbones,
we utilize (1) a teacher network for distillation, and (2) our target backbone
that we compress (pruning + quantization).

###  4.1 Quantitative Evaluation

The results are shown in Table 2. We include throughput (frames per second,
FPS), accuracy, mean IoU, and model size in the study. For all results, we use
an ImageNet-1k pretrained backbone and a randomly initialized decoder. We then
finetune the models for 20,000 iterations on the LPCV dataset.

As shown in Table 2, out of the three decoders - FANet, DeepLabv3 and UPerNet,
FANet offers the least mean IoU with the fewest parameters, while DeepLabv3
offers the highest mean IoU with the largest model. We choose UPerNet, a
middle ground between the two, to conduct our backbone comparison. Out of the
four backbones we compared, we noticed that: (1) although ResNet18 has
significantly more parameters than MobileViT, it has faster runtime. This is
due to the self-attention matrix multiplication computation not contributing
to parameter count. Also, convolution kernels in PyTorch seem to be more
optimized than attention kernels. (2) Swin-v2-Tiny offers similar, although
slightly higher IoU than Swin-v2-Base at just 40% of the parameters. We
believe that this is due to the LPCV dataset’s small size. Swin-v2-Tiny is
already capable of learning from data given our training configuration (data
augmentation, etc.).

Table 2: Comparison of different models on the LPCV Challenge Dataset

Conventional Method
---
Arch. + Backbone | FPS (GPU: RTX 4090) | FPS (Jetson Nano) | Acc. | MIoU | # Params | Score
FANet (ResNet18) | 120.48 | - | 0.82 | 0.4752 | 13.9M | -
DeepLabv3 (MobileViT) | 64.93 | - | 0.85 | 0.5872 | 6.4M | -
UPerNet (ResNet18) | 76.49 | 6.02 | 0.84 | 0.5200 | 16.4M | 3.130
UPerNet (MobileViT) | 56.75 | 3.11 | 0.85 | 0.5365 | 5.6M | 1.668
UPerNet (Swin-v2-T) | 47.80 | - | 0.88 | 0.6468 | 40M | -
UPerNet (Swin-v2-B) | 33.70 | - | 0.88 | 0.6436 | 108M | -
Distillation-based
UPerNet (MobileViT, KD) | 56.75 | 3.11 | 0.86 | 0.6056 | 5.6M | 1.883
UPerNet (ResNet18, KD) | 76.49 | 6.02 | 0.85 | 0.5493 | 16.4M | 3.307
Distillation and Quantization
UPerNet (MobileViT, KD, fp16) | 57.07 | 4.49 | 0.86 | 0.6056 | 5.6M | 2.719
UPerNet (ResNet18, KD, fp16) | 81.75 | 5.94 | 0.85 | 0.5493 | 16.4M | 3.262

We choose UPerNet with Swin-v2-T as the teacher model and apply knowledge
distillation with simple MSE Loss between the final segmentation logits. As
shown in Table 2, UPerNet with MobileViT has a significant increase in IoU
(0.06) after distillation, while ResNet18 + UPerNet sees a 0.02 increase. We
also experiment with logit-based knowledge transfer using both simple MSE Loss
and a KL-divergence loss of softmax of different temperatures. We found that
consistently, the simple MSE Loss works better. On the other hand, using KL-
divergence loss with temperatured softmax led to unstable training. Our final
network after knowledge distillation: UPerNet with MobileViT, achieves 0.61
mean IoU while using only 5.6M parameters.

In addition, we report the throughput of these models on a Jetson Nano.
Consistent with FPS results on desktop GPUs, ResNet18 achieves much better
runtime than MobileViT although having a bigger model size. Since we are
limited by PyTorch’s support of datatypes, we quantize the model to half-
precision floats and report their FPS. Interestingly, while MobileViT saw a
1.4x increase in throughput, the throughput for ResNet18 at fp16 decreased. We
believe this is still associated with PyTorch’s implementation of convolution
vs. attention kernels in CUDA.

In our final verification experiments, ResNet18 still achieves the best score
(3.307). The best score achieved by a MobileViT is 2.719. However, MobileViT
achieves a much higher mean IoU. Additionally, MobileViT still has a
theoretical advantage in both model size and FLOPs, which implies that with
the continuing optimization of software kernels for more recent computations
like the attention mechanism, it has the potential of achieving better
accuracy vs. runtime trade-off than conventional ConvNets. Our final model
(distilled and quantized UPerNet + MobileViT) consumes 3742 MB of RAM with
1030 MB of swap on the Jetson Nano, while when idle, the OS alone consumes
1424+558 MB of RAM+swap. With optimized baremetal solutions and inference
engines like ONNX or TensorRT, the memory footprint can be further optimized.

###  4.2 Pruning Ablation

Our first pruning attempt is filter pruning for the Conv layers in MobileViT.
Our initial estimation is to get at least 0.5 sparsity since traditional
ConvNets with similar size (MobileNet) are able to achieve that with
negligible loss in accuracy. The results are shown in Figure 3. As observed in
the plot, starting at 0.2 sparsity results in a significant drop in mean IoU.

Our conjecture is that MobileViT is already compressed to a high degree, in
order to be able to tolerate any additional pruning. In order to understand
this hypothesis further, we attempt channel pruning, unstructured pruning on
MobileViT, as well as filter pruning on the much larger Swin-v2-T and Swin-B
models. As shown in Figure 3, channel pruning results in a similar accuracy
degradation compared to filter pruning, while unstructured pruning can achieve
much high sparsity. When pruning Swin models, accuracy only starts to degrade
at around 0.7 sparsity, corroborating our hypothesis about MobileViT’s
inability to further compress. We note that Swin Transformers contain MLP
layers instead of Conv layers, which have more inherent redundancy than
MobileViT.

Figure 3: Effects of Pruning Convolutional and Linear layers: MIoU vs Sparsity
trends for different models Table 3: Effects of head pruning on different
models

Model | # Attention Heads | Mean IoU | # Params | FPS (GPU)
---|---|---|---|---
UPerNet (MobileViT) | 4 | 0.545 | 5.6M | 62
UPerNet (MobileViT) | 2 | 0.533 | 4.6M | 62
DeepLabv3 (MobileViT) | 4 | 0.601 | 6.4M | 64
DeepLabv3 (MobileViT) | 2 | 0.590 | 5.4M | 64

Next, we perform pruning of attention heads based on L2 norm of the weights on
the models with transformer-based backbones. The results are shown in Table 3.
The results show that reducing the number of heads from 4 (the default value
for MobileViT) to 2 had minimal impact on model performance while offering a
significant reduction in the number of parameters. Of note, we pruned the
weights by setting the selected heads to zero but did not actually change the
dimensions of the attention block. Therefore, the parameter reduction is
theoretical, and we do not see an actual improvement in speed as the number of
operations does not change.

###  4.3 Qualitative Evaluation

Using the UPerNet with MobileViT distilled from Swin-v2-T and fp16
quantization as our final framework, the segmentation results on the LPCV
dataset are shown in Figure 4. We can see from the resulting images that our
compressed model actually does a fairly good job at capturing the dominant
features of the image, but fails to accurately capture finer details. For
example, on the second image from the right, our compressed model has
difficulty identifying distant roads. This task is particularly difficult as
many classes look very similar, even to the human eye. The differences between
mudflow and a dirt road are subtle and often not apparent.

Figure 4: Segmentation results on the LPCV dataset

This issue becomes further complicated as we have been evaluating our models
on global performance. Our compression techniques are only preserving the
parameters and interactions that are the most performant at a global scale.
Therefore, it is reasonable to expect that not all classes would be impacted
in a similar manner as classwise distributions are, by the nature of disaster
scenes, not equal. This is evident even in our example images as only a few
classes, including the background, constitute the overwhelming majority of the
pixels. This issue of classwise imbalance becomes problematic because, in
disaster scene parsing, the low-frequency classes are actually the most
important in a real-world scenario.

###  4.4 Unsuccessful Experiments: Iterative Pruning

Figure 5: MIoU vs Sparsity trends for iterative pruning

We perform iterative pruning in order to try and achieve better pruning
results on MobileViT as shown in Figure 5. Our approach involves initially
pruning a small percentage of filters remaining in each convolution layer,
followed by 1k or 2k fine-tuning iterations. We then repeat this cycle and
collect performance metrics for comparison. Our hypothesis was that by not
removing many filters at once, the model may be able to better preserve
performance by redistributing the weights during the fine-tuning iterations.
We find that this process gives marginal improvements over the bulk pruning
described previously, but still does not provide much avenue for compression.

##  5 Summary

In our work, we establish baseline ConvNet and Transformer models for
segmentation on the LPCV dataset. We devise a framework to compress these
models to meet the memory usage requirements of no more than 4GB of RAM
system-wide. Our compression approach builds on structured pruning to remove
redundant heads from the transformer layers, along with linear layer weights.
Furthermore, we also leverage the higher expressive power of a stronger
teacher network using distillation, which guides our smaller, more efficient
student network. In order to strike the optimal balance between memory usage
and inference latency improvements while maintaining segmentation accuracy and
mean IoU, we further explore quantization strategies and employ the most
effective method on the student network at inference time on the Jetson Nano.
We also explore other compression techniques and their effectiveness in
reducing the model size and inference latency. On the whole, our work focuses
on and successfully develops a framework to enable highly efficient
segmentation on the edge that could be deployed on resource-constrained
devices for real-world applications.

##  6 Future Directions

Moving forward, our future work will focus on further compression and
evaluation of our model on the low-cost NVIDIA Jetson Nano 2GB. This
environment provides the right balance of low power, affordability, and robust
support in a remote execution setting, especially in disaster-stricken
regions.

We also plan to experiment and analyze other quantization techniques and their
effectiveness in reducing inference latency while remaining performant.
Additionally, we will investigate the limitations of PyTorch’s support for
quantized CUDA operations and seek alternative solutions if necessary. Beyond
this, we hope to explore low-rank approximation of the attention mechanism for
our use case.

We also aim to explore adaptive pruning methods such as Layer-Adaptive
Magnitude-based Pruning (LAMP) [32] to achieve a better trade-off between
sparsity and performance without requiring extensive hyperparameter tuning or
heavy computation.

## References

  * [1] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale, 2021.
  * [2] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2015.
  * [3] Xilai Yuan, Qiang Liu, Ming Tang, and Duo Zhang.  Revisiting knowledge distillation via label smoothing regularization.  arXiv preprint arXiv:2002.11042, 2020.
  * [4] Jaejun Lee, Sangdoo Kim, Jongseok Lee, Haecheon Shin, and Jung Kwon Kim.  Snip: Single-shot network pruning based on connection sensitivity.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4364–4373, 2018.
  * [5] Jian-Hao Luo and Jianxin Wu.  Thinet: A filter level pruning method for deep neural network compression.  IEEE transactions on pattern analysis and machine intelligence, 40(6):131–144, 2017.
  * [6] Ji Lin, Rongrong Ji, Yimin Li, Baochang Gao, and Yongjian Li.  Accelerating cnns via global & dynamic filter pruning.  arXiv preprint arXiv:1805.06394, 2018.
  * [7] Tianlong Chen, Yu Cheng, Zhe Gan, Lu Yuan, Lei Zhang, and Zhangyang Wang.  Chasing sparsity in vision transformers: An end-to-end exploration, 2021\.
  * [8] Paul Michel, Omer Levy, and Graham Neubig.  Are sixteen heads really better than one?, 2019.
  * [9] Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, and He Wen.  Dorefa-net: Training low bitwidth convolutional neural networks with low bitwidth gradients.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 442–450, 2016.
  * [10] Ron Banner, Yury Nahshan, Elad Hoffer, and Daniel Soudry.  Post-training 4-bit quantization of convolution networks for rapid-deployment, 2018.
  * [11] Daniel Choukroun, Thibaud Bouvier, and Benoit Gauthier.  Low-latency, low-precision inference with integer neural networks.  In Advances in Neural Information Processing Systems, pages 1295–1305, 2019.
  * [12] Till Nagel, Simon Brunner, Pablo Marquez Neila, Martin Wiesner, and Simon Kornblith.  Adaround: A lossy compression for the masses.  arXiv preprint arXiv:2006.08924, 2020.
  * [13] Yanghao Liu, Dong Li, and Jian Cheng.  Post-training quantization for vision transformers with enhanced rank-aware strategies.  arXiv preprint arXiv:2103.15730, 2021.
  * [14] Sangho Lee, Youngjae Yu, Gunhee Kim, Thomas Breuel, Jan Kautz, and Yale Song.  Parameter efficient multimodal transformers for video representation learning.  arXiv preprint arXiv:2012.04124, 2020.
  * [15] Beidi Chen, Tri Dao, Eric Winsor, Zhao Song, Atri Rudra, and Christopher Ré.  Scatterbrain: Unifying sparse and low-rank attention.  Advances in Neural Information Processing Systems, 34:17413–17426, 2021.
  * [16] Jyotikrishna Dass, Shang Wu, Huihong Shi, Chaojian Li, Zhifan Ye, Zhongfeng Wang, and Yingyan Lin.  Vitality: Unifying low-rank and sparse approximation for vision transformer acceleration with a linear taylor attention, 2022.
  * [17] Bosheng Qin, Juncheng Li, Siliang Tang, and Yueting Zhuang.  Dba: Efficient transformer with dynamic bilinear low-rank attention, 2022\.
  * [18] Sinong Wang, Belinda Z. Li, Madian Khabsa, Han Fang, and Hao Ma.  Linformer: Self-attention with linear complexity, 2020.
  * [19] Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret.  Transformers are rnns: Fast autoregressive transformers with linear attention, 2020.
  * [20] Rewon Child, Scott Gray, Alec Radford, and Ilya Sutskever.  Generating long sequences with sparse transformers, 2019.
  * [21] Wenhai Wang, Enze Xie, Xiang Li, Deng-Ping Fan, Kaitao Song, Ding Liang, Tong Lu, Ping Luo, and Ling Shao.  Pyramid vision transformer: A versatile backbone for dense prediction without convolutions, 2021.
  * [22] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows, 2021.
  * [23] Qiang Chen, Qiman Wu, Jian Wang, Qinghao Hu, Tao Hu, Errui Ding, Jian Cheng, and Jingdong Wang.  Mixformer: Mixing features across windows and dimensions, 2022.
  * [24] Shitao Tang, Jiahui Zhang, Siyu Zhu, and Ping Tan.  Quadtree attention for vision transformers, 2022.
  * [25] Zhendong Yang, Zhe Li, Ailing Zeng, Zexian Li, Chun Yuan, and Yu Li.  Vitkd: Practical guidelines for vit feature knowledge distillation, 2022\.
  * [26] Shixing Yu, Tianlong Chen, Jiayi Shen, Huan Yuan, Jianchao Tan, Sen Yang, Ji Liu, and Zhangyang Wang.  Unified visual transformer compression, 2022.
  * [27] Ankur Kumar.  Vision transformer compression with structured pruning and low rank approximation, 2022.
  * [28] Sachin Mehta and Mohammad Rastegari.  Mobilevit: Light-weight, general-purpose, and mobile-friendly vision transformer.  ArXiv, abs/2110.02178, 2021.
  * [29] Ping Hu, Federico Perazzi, Fabian Caba Heilbron, Oliver Wang, Zhe L. Lin, Kate Saenko, and Stan Sclaroff.  Real-time semantic segmentation with fast attention.  IEEE Robotics and Automation Letters, 6:263–270, 2020.
  * [30] Tete Xiao, Yingcheng Liu, Bolei Zhou, Yuning Jiang, and Jian Sun.  Unified perceptual parsing for scene understanding.  ArXiv, abs/1807.10221, 2018.
  * [31] Liang-Chieh Chen, George Papandreou, Florian Schroff, and Hartwig Adam.  Rethinking atrous convolution for semantic image segmentation.  ArXiv, abs/1706.05587, 2017.
  * [32] Jaeho Lee, Sejun Park, Sangwoo Mo, Sungsoo Ahn, and Jinwoo Shin.  Layer-adaptive sparsity for the magnitude-based pruning.  In International Conference on Learning Representations, 2020.


# Fast automatic deforestation detectors and their extensions for other
spatial objects.

Jesper Muren2, Vilhelm Niklasson2,
Dmitry Otryakhin1,2, Maxim Romashin

(1<EMAIL_ADDRESS>
2 Department of Mathematics, Stockholm University, SE-10691 Stockholm, Sweden
)

###### Abstract

This paper is devoted to the problem of detection of forest and non-forest
areas on Earth images. We propose two statistical methods to tackle this
problem: one based on multiple hypothesis testing with parametric distribution
families, another one- on non-parametric tests. The parametric approach is
novel in the literature and relevant to a larger class of problems– detection
of natural objects, as well as anomaly detection. We develop mathematical
background for each of the two methods, build self-sufficient detection
algorithms using them and discuss numerical aspects of their implementation.
We also compare our algorithms with those from standard machine learning using
satellite data.

##  1 Introduction

Within the last decade, there have been some attempts to create remote
deforestation detectors. Contemporary literature on the topic seems to be
geared towards machine learning approaches. Just to name a few examples we
refer to [10] that employed deep learning for deforestation detection in the
Brazilian Amazon and [15] that presented k-nearest neighbour classification to
track deforestation in Peru. Bayesian approach has also seen use, e.g. in [12]
on time series data from Fiji. The present work aims to create methods that
not only show excellent performance, but are also, unlike machine learning
methods, fast, simple and interpretable. To achieve these characteristics we
employ sample means and covariances and empirical characteristic functions
(ECFs) of colour intensities of image samples. The idea behind it is that
forests seen from above, like many other natural objects, show rather simple
homogeneous structure. Moreover, in many cases natural objects are
distinguishable due to their colour features rather than their spatial
patterns. Thus, operating on distributions and their ECFs and completely
erasing all the spatial dependencies seems to be reasonable.

The general statistical setting of the problem is as follows. In visible
spectrum images of geographical objects are represented as 3 matrices. Each of
them consists of intensities of one of the red, green or blue colours. We
employ supervised learning where tags (forest / non-forest) are assigned to
entries of the training data and learn features of the distributions of the
matrices in those entries. During classification, test images are split into
much smaller square sub images for the purpose of obtaining better resolution,
and these sub-images are tested against reference images that are known to
contain forest.

We intend to investigate two different approaches to the problem: non-
parametric multiple hypothesis testing and parametric one. In the case of
parametric multiple hypothesis testing, we assume parametric distributions for
the colour intensities of pixels and use training data to estimate the
parameters of said distributions. New data is then tested against this
distribution. In the case of the non-parametric setting, we do not assume any
known distribution for the colour intensities but instead test if the
multivariate distributions of the colour intensities of the training data and
new data have equal means and covariance matrices.

We build both the parametric and the non-parametric algorithms in such a way
that they work on images in visible spectrum, in particular, RGB composites.
Thus, the algorithms do not require any certain pre-processing technique for
the data and can be applied to a broad range of sources of images of Earth.
There are at least two ways of obtaining large datasets of Earth images:
satellite data obtained through open access Earth observation missions, e.g.
Landsat [13] and Sentinel-2 [4] and aerial photo imaging. In this work,
Sentinel-2 images are used because they provide good spatial and temporal
resolution. Sentinel-2 mission provides data from visible, near infrared and
short-wave infrared sensors in 13 corresponding spectral bands. Although we
use RGB composites obtained via a post-processing procedure of Sentinel-2
original tiles, our methods could work on all of the data in those 13 bands.

###  1.1 A practical example

In order to give more intuition of what task we tackle, we present a real-life
example of deforestation detection. Figure 1a shows a Sentinel-2 image of a
piece of land near the city of Arkhangelsk in Russia suffering from logging.
Figure 1b is a black-and-white image representing the detection result. White
represents forest and black is used for all types of terrains not covered by
it.

(a) Forest image analysed

(b) The result of detection. White- forest, black- non-forest

Figure 1: A plot of land with deforestation and the result of detection
produced by the non-parametric algorithm.

One can see that the detector highlighted the areas with logging traces (black
rectangles) as well as natural objects: lakes, bare ground and pieces of
rivers.

The rest of the paper is organized as follows. Section 2 briefly mentions the
key steps in data preparation and discusses what features forest data possess.
Sections 3 and 4 describe correspondingly novel non-parametric and parametric
methods. Section 5 discusses theoretical properties of the two methods, while
the last section compares the two with well-known methods from machine
learning using Satellite data.

Our main contribution contains the application of the known theory of
Hotelling $T^{2}$ statistic to the problem of deforestation detection, which
occupies Section 3. There, we prove Theorem 1 and using it, build
classification and training algorithms in 3.2 and 3.3. Another major
contribution is Section 4 (excluding 4.1), where we develop a new setting of
hypothesis testing for stable models based on empirical characteristic
functions in 4.2 and use this setting to build algorithms for detection of
forest in 4.3 and 4.4.

##  2 Data and exploratory analysis

To obtain data in visible spectrum we use a simple preprocessing procedure,
collecting red, green and blue bands from Sentinel-2 Level-1C products and
combining them into one GeoTIFF file. It is then handled by GDAL [5], which
creates three numeric matrices (one per each of the red, green and blue
channels) with values between 0 and 1. The procedure is described in appendix
A in more detail.

Initial inspection of the empirical density of the forest data gave us some
additional ideas for potentially suitable distributions. We considered
distributions including, but not limited to, stable, Normal, Log-Normal,
Gamma, Cauchy, Weibull and Levy. We estimated parameters for each of the
colour intensity distributions and calculated the root-mean-square-error
(RMSE)

(a) Forest image analysed

(b) Comparison of densities for red colour intensity

(c) Comparison of densities for green colour intensity  (d) Comparison of
densities for blue colour intensity

Figure 2: Plots of empirical and estimated probability density functions for
distribution of colour intensities in an image of a forest. The means of the
red, green and blue intensities are respectively 0.16, 0.24 and 0.20.

$\text{RMSE}=\sqrt{\frac{\sum_{i=1}^{n}(\hat{f}(x_{i})-f^{*}(x_{i}))^{2}}{n}}$

where $\hat{f}(\cdot)$ is the probability density function using the estimated
parameters, $f^{*}(\cdot)$ is the empirical density and $x_{i}$ is colour
intensity of data point $i=1,\dots,n.$ To choose the best fitting
distributions we then plotted the densities of the three distributions with
the lowest RMSE’s together with the empirical densities and compared. These
distributions appeared to be Normal, Gamma and Stable consistently throughout
the set of images. An example of a forest image with such comparison can be
seen in Figure 2, where we can observe the three best fitting distributions
for a forest image with all parameters estimated. Further from the plots we
see that the Normal distributions do not seem to fit very well while the Gamma
and Stable distributions are quite close, with Stable edging out Gamma
especially for the intensity of the colour red.

Since the focus in this work is to be able to detect if an area has forest or
not, we describe detection of presence or lack of this type of natural
objects. However, our exploratory data analysis shows us that the framework
developed could be used to detect other natural objects from satellite images
as well. In Figures 11 and 12 in appendix B, we can see promising results of
using the stable distribution to describe the colour intensities of images of
mountains and sea. When it comes to man-made objects such as, for example a
city, our approach does not work, and we have to at least assume a multimodal
mixture of some kind to be appropriate since satellite images contain
distinctly different parts, i.e., an orange roof, green grass and grey
pavement.

##  3 Detection with a non-parametric method

In this section, the non-parametric method for forest detection is developed.
Later, in Section 4, the ideas presented here are used to create the theory
and implementation for the parametric one.

###  3.1 Comparing images with Hotelling T-square statistic

The Mahalanobis distance $D$ was originally introduced as a distance between a
given multivariate point $\boldsymbol{p}$ and a distribution. For a point
$\boldsymbol{p}$ and a distribution with mean vector $\boldsymbol{\mu}$ and
covariance matrix $\boldsymbol{\Sigma};$ the distance is given by

$D=\sqrt{(\boldsymbol{p}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{p}-\boldsymbol{\mu}).}$

Where, in practice, the mean and covariance matrix of the distribution would
be estimated by the sample mean and covariance matrix. Further extensions can
be made to calculate the distance between two multivariate populations with

$D=\sqrt{(\boldsymbol{\bar{X}}_{1}-\boldsymbol{\bar{X}}_{2})^{T}\boldsymbol{\hat{\Sigma}}^{-1}(\boldsymbol{\bar{X}}_{1}-\boldsymbol{\bar{X}}_{2}}).$
(1)

Where $\boldsymbol{\bar{X}}_{1}$ and $\boldsymbol{\bar{X}}_{2}$ are the sample
mean vectors of the two populations. The variable
$\boldsymbol{\hat{\Sigma}}^{-1}$ is the inverse of the pooled sample
covariance matrix, given by
$\boldsymbol{\hat{\Sigma}}=(n_{1}\boldsymbol{\hat{\Sigma}_{1}}+n_{2}\boldsymbol{\hat{\Sigma}_{2}})/(n-2)$
for samples with respective sample sizes $n_{1},n_{2}$, total sample size
$n=n_{1}+n_{2}$ and sample covariance matrices $\hat{\Sigma}_{1}$,
$\hat{\Sigma}_{2}$. An interested reader can find more on the topic in [3].

The Mahalanobis distance is often used in cases where there is correlation
present in the considered multivariate distribution. As the reader might have
noticed, the Mahalanobis distance expression is closely related to the
statistical concept of standardization, i.e, subtracting the mean and dividing
by the standard deviation. The idea is to decorrelate and scale each
dimensions variance to be 1 and then calculate the Euclidean distance.

######  Theorem 1.

Let $\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$ be independent samples of
sizes $n_{1}\times p$ and $n_{2}\times p$ composed from i.i.d p-dimensional
random vectors from distributions with mean $\boldsymbol{\mu}$ and covariance
matrix $\boldsymbol{\Sigma}$. Then, for $D$ from (1) we have that
$\frac{n_{1}n_{2}}{n_{1}+n_{2}}D^{2}\xrightarrow{d}\chi^{2}(p)$.

###### Proof.

By the multidimensional central limit theorem

$(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})\xrightarrow{d}\mathcal{N}(\boldsymbol{0},\gamma\boldsymbol{\Sigma}),\quad\text{where
}\gamma=n_{1}^{-1}+n_{2}^{-1}.$ (2)

After introducing $\boldsymbol{A}:=(\gamma\boldsymbol{\Sigma})^{1/2}$ and
$\boldsymbol{\hat{A}}:=(\gamma\boldsymbol{\hat{\Sigma}})^{1/2}$ we can write

$(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})\xrightarrow{d}\boldsymbol{A}\boldsymbol{Z},\quad\text{where
}\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{p}).$ (3)

Taking into account (3) and that
$\boldsymbol{\hat{A}}^{-1}\xrightarrow{P}\boldsymbol{A}^{-1}$, where the later
is a matrix of constants, by a property of weak convergence

$\left[\boldsymbol{\hat{A}}^{-1},(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})\right]\xrightarrow{d}\left[\boldsymbol{A}^{-1},\boldsymbol{AZ}\right].$
(4)

Then, we introduce a new quantity

$\begin{split}\boldsymbol{Q}:=\left[\boldsymbol{\hat{A}}^{-1}(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})\right]^{T}\left[\boldsymbol{\hat{A}}^{-1}(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})\right].\end{split}$
(5)

The multivariate continuous mapping theorem applied to (4) implies
$\boldsymbol{Q}\xrightarrow{d}\boldsymbol{Z}^{T}\boldsymbol{Z}\stackrel{{\scriptstyle
d}}{{=}}\chi^{2}(p)$. On the other hand,

$\boldsymbol{Q}=(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})^{T}(\boldsymbol{\hat{A}}^{-1})^{T}\boldsymbol{\hat{A}}^{-1}(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}}),$
(6)

and using the following chain of transformations

$(\boldsymbol{\hat{A}}^{-1})^{T}\boldsymbol{\hat{A}}^{-1}=(\boldsymbol{\hat{A}}^{T})^{-1}\boldsymbol{\hat{A}}^{-1}=(\boldsymbol{\hat{A}}^{T}\boldsymbol{\hat{A}})^{-1}=(\gamma\boldsymbol{\hat{\Sigma}})^{-1}$
(7)

we obtain

$\boldsymbol{Q}=(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}})^{T}(\gamma\boldsymbol{\hat{\Sigma}})^{-1}(\boldsymbol{\bar{X}_{1}}-\boldsymbol{\bar{X}_{2}}).$
(8)

∎

######  Remark 1.

Note, that $\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$ do not have to be
drawn from the same distribution. If they are drawn from different ones, but
possessing the same mean and covariance matrix, the result of the theorem
still holds.

###  3.2 Hypothesis testing based on Mahalanobis distances between images

Suppose we have an RGB image of forest (Figure 3). It consists of three
matrices corresponding to each of the red, green and blue channels. The
matrices contain values of colour intensities of pixels in the range $[0,1]$.
Then, each matrix is split into column-vectors that are subsequently stack
together to form one column-vector of the length $p^{2}$. The three obtained
large vectors are bound to form a $(p^{2})\times 3$ matrix $\boldsymbol{X}$.
Since the same operation is applied to every matrix, the triplets of pixel
intensities are kept unchanged, so it is possible to analyse dependencies
between colours. On the contrary, matrix $\boldsymbol{X}$ does not contain any
information about spatial locations of the pixels. From $\boldsymbol{X}$, a
three-dimensional mean vector and a $3\times 3$ covariance matrix are
calculated.

Figure 3: Data processing for the non-parametric testing

The assumptions of the theorem hold for data of type of $\boldsymbol{X}$:
since all entries are bounded between $0$ and $1$, the means and covariances
will always be finite. This naturally leads us to a hypothesis testing
framework. Having two images- one picture that is known to depict forest and
one test image, we can use the theorem to accept or reject the null hypothesis
that the test image is forest. This is done by first converting RGB composites
into matrices $\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$, obtaining their
sizes, means and covariance matrices, and second, computing the Mahalanobis
distance between $\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$, and checking
whether or not $T:=\frac{n_{1}n_{2}}{n_{1}+n_{2}}D^{2}$ is lower than a
certain threshold.

Figure 4: Classification in the non-parametric model

A scheme for multiple hypothesis testing framework is depicted on Figure 4.
When it comes to classification of a certain image, rescaled Mahalanobis
distances $T_{1},T_{2},\dots$ between the image and every forest picture in
the training set are computed. Out of all $T$’s, the minimal one is chosen,
and is then compared to a threshold obtained during model training. If the
minimal $T$ is smaller that the threshold, the test image is tagged forest, if
greater, then– non-forest.

######  Remark 2.

This setting is equivalent to multiple hypothesis testing with the thresholds
set to be equal. This restriction can be dropped, in which case the precision
of predictions will generally increase. But the downside of such a
generalization would be that the additional computational complexity during
both training and classification would increase linearly with the number of
forest images.

######  Remark 3.

This algorithm describes dependences between colours by means of the pooled
sample covariance matrix. On the contrary, all the information about spatial
features are erased. Also, distributions of colour intensities in each channel
are described only by their location parameters.

###  3.3 Training and cross-validation of the non-parametric model

Training of the non-parametric model is schematically depicted on Figure 5.
One can observe two sets of images there: reference forest images in the left
column (set 1) and images of both types in the upper row (set 2). For every
pair i-th image from set 1 and j-th image from set 2, statistic
$\hat{T}_{ij}=\frac{n_{1}n_{2}}{n_{1}+n_{2}}D_{ij}^{2}$ is computed, where
$D_{ij}$ is the Mahalanobis distance between the two pictures. Then, for every
image from set 2, the minimal statistic $\min_{i}\hat{T}_{ij}$ is found.
Lastly, the threshold $T$ is chosen to minimize misclassification error when
classifying images from set 2 according to a result of evaluation of Boolean
expression $\min_{i}\hat{T}_{ij}<T$. As previously, True corresponds to
forest, while False– to non-forest.

Aside from simple training on a set of images, we adopt the setting of k-fold
cross-validation for this model. Initially, we shuffle pictures and group them
into k folds. In this procedure, in order to preserve internal features, we
don’t split or merge pictures. We follow a similar procedure to what was
described above. Once a fold is chosen, set 1 consists of all forest pictures
of a training set, while set 2 is a holdout set. Then statistics
$\min_{i}\hat{T}_{ij}$ are computed for each of the $k$ holdout sets and $T$
is chosen to minimize the misclassification error on all the combinations.
Cross-validation is summarized in Algorithm 1.

Figure 5: Training in the non-parametric model Algorithm 1 Non parametric
cross validation procedure

1:Inputs: Data set of $N$ images labelled forest or non-forest, number of
folds $k$ & max threshold range $T_{max}$.

2:Compute means and covariances for every image.

3:Randomly shuffle images and group them into $k$ folds.

4:for $i=1$ to $k$ do

5: for each image in fold $i$ do

6: Compute $\hat{T}=\frac{n_{1}n_{2}}{n_{1}+n_{2}}D^{2}$ from Theorem 1,
between the image and all images labelled forest in the remaining $k-1$ folds.

7: Keep only the smallest $\hat{T}$ for the image.

8: end for

9:end for

10:for $t$ in $0$ to $T_{max}$ do

11: For all images and all $\hat{T}$’s, tag images with $\hat{T}<t$ as forest
and $\hat{T}>t$ as non forest.

12: Compute the accuracy for the given threshold $t,$ $\text{Acc}(t).$

13:end for

14:return $\text{arg}\max_{t}\text{Acc}(t)$ and $n$’s, means and covariances
of forest images.

##  4 Detection with a parametric method

The non-parametric approach described in Section 3 is simple and universal in
terms of underlying data, it doesn’t make any assumptions on the data coming
from RGB composites. On the other hand, it possesses a theoretical drawback-
it doesn’t take advantage of the distribution of the data, only its mean and
covariance matrix. The parametric approach to the detection was designed
basing on the opposite ideas- to capture the actual distributions of the data
and parametrize them in order to make use of parametric theoretical machinery.

###  4.1 Stable distributions for signal processing

The second classification method we develop in this study is based on
univariate stable distributions as they proved to fit well to the empirical
data (Section 2). In this short section we give some context for use of stable
distributions in the field of signal processing, describe how they are fit to
satellite data in the present paper and how we estimate sampled data.

The stable distribution family is a class generalising the Gaussian one. In
signal processing it appears due to two factors. First, stable distribution is
a result of the Generalised Central Limit Theorem which says that if there
exists any limit of a mean of i.i.d. variables with possibly unbounded
variations, then this limit must have a stable distribution. Thus, this class
is a natural choice for modelling signals with potentially large deviations.
Second, stable distributions might be asymmetric, and it is controlled by one
of the parameters. This family of random variables has already seen use in
remote sensing areas such as ship detection [16] and video background
subtraction [2]. While it is unclear to us what role the Generalised Central
Limit Theorem plays in shaping the distributions of forest colour intensities,
it is certain that parametric control of their shape provides accurate fitting
(Figure 2) when the general stable distribution is assumed.

######  Definition 1.

A random variable $X$ is said to have a stable distribution if there are
parameters $0<\alpha\leq 2,$ $\sigma\geq 0$, $-1\leq\beta\leq 1$ and
$\delta\in\mathds{R}$ such that its characteristic function has the following
form:

$\phi(t)=E[e^{itX}]=\begin{cases}\text{exp}\bigg{\\{}-\sigma^{\alpha}|t|^{\alpha}(1-i\beta(\text{sign
}t)\text{tan}\frac{\pi\alpha}{2})+i\delta t\bigg{\\}}&\text{if $\alpha\neq
1$}\\\ \text{exp}\bigg{\\{}-\sigma|t|(1+i\beta\frac{2}{\pi}(\text{sign
}t)\text{ln}|t|)+i\delta t\bigg{\\}}&\text{if $\alpha=1$.}\end{cases}$ (9)

Where sign $t=0$ if $t=0.$ The parameter $\alpha$ is known as the index of
stability and controls the decay of the tails. The parameters $\beta$ is not
the classical statistical skewness but it indicates how the distribution is
skewed, indicating left skewness if $\beta<0$, right skewness if $\beta>0$ and
symmetry if $\beta=0.$ The parameter $\sigma$ is a scale parameter controlling
variability and $\delta$ a location parameter which shifts the distribution.
The location parameter only corresponds to the mean when $\alpha>1$, otherwise
the mean is undefined. In fact, for $\alpha<2$ the $p$-th absolute moment
$E[|X|^{p}]<\infty$ iff $p<\alpha.$

Apart from a few cases, the closed form probability function is unknown for
stable laws. There exists a method for computing stable densities numerically
[9] but it is computationally demanding. That is why we developed a
mathematical foundation of our parametric method based on CF’s and ECF’s. The
most famous special cases where the probability density function can be
expressed explicitly is when the stable distribution follows a Gaussian,
Cauchy or Lévy distribution. A stable distribution is Gaussian when
$\alpha=2$.

####  4.1.1 Parameter estimation

This section will cover estimation of the four parameters of the univariate
stable distribution used in this article. There is a number of techniques
available for estimation of these parameters, including quantile-based
estimators, Maximum-likelihood-based estimators and those based on the method
of moments. For this work, the Koutrouvelis regression-type estimator [6] was
chosen, as it shows a great balance between precision and computational
performance. Estimators which performed well but were very slow were also
studied, including ones based on the generalized method of moments. Similarly,
estimators which were faster than the Koutrouvelis one were assessed, but
these did not perform well enough for our purposes: the McCulloch estimator
[8] based on quantiles and the one developed by Press [11]. A more in-depth
study on the performance of different parameter estimators may be of interest
to individuals interested in improving the algorithm proposed in this thesis.

The Koutrouvelis regression-type estimator is quite fast, but one big downside
is that the estimator of the location parameter $\delta$ is often unreliable.
For our purposes this turns out not to be a big issue though because
experimentally we have seen that the parameter $\alpha>1$ in all our forest
data and when this is the case the location parameter $\delta$ equals the mean
of the distribution [14, Property 1.2.16]. Which implies that the sample mean
is a consistent estimator for the location parameter $\delta$ for our data.

###  4.2 A test for empirical characteristic functions of stable distributions

In this part, we develop a statistical test that rejects a null hypothesis
that a sample comes from a stable distribution defined by (9) with certain
parameters.

First, we provide a limit theorem for vectors composed of parts of ECFs of
stable random variables. For a sequence $\\{X_{k}\\}_{k=1,\dots}$ of iid
stable random variables the ECF is given by

$\phi_{n}(t):=\frac{1}{n}\sum_{k=1}^{n}[\cos(tX_{k})+i\sin(tX_{k})],\text{where
we denote}$ (10) $u_{n}(t):=\cos(tX_{k})\quad v_{n}(t):=\sin(tX_{k}),$ (11)

and define new random vectors

$\boldsymbol{Z_{n}}(t):=\frac{1}{n}\sum_{k=1}^{n}\begin{bmatrix}\cos(tX_{k})\\\
\sin(tX_{k})\end{bmatrix}\quad\boldsymbol{Z_{0}}(t):=\begin{bmatrix}\mathrm{Re}[\phi(t)]\\\
\mathrm{Im}[\phi(t)]\end{bmatrix}$ (12)

######  Proposition 1.

Fix $t\in R$, then
$\sqrt{n}\left[\boldsymbol{Z_{n}}(t)-\boldsymbol{Z_{0}}(t)\right]\xrightarrow{d}\mathcal{N}(0,\boldsymbol{\Sigma_{z}})$,
where

$\boldsymbol{\Sigma_{z}}=\begin{bmatrix}\sigma_{11}&\sigma_{12}\\\
\sigma_{21}&\sigma_{22}\end{bmatrix}$ (13)

and the elements of the matrix are given by

$\begin{split}\sigma_{11}=&\frac{1}{4}\left[\phi(2t)+2+\phi(-2t)-(\phi(t))^{2}-2\phi(t)\phi(-t)-(\phi(-t))^{2}\right]\\\
\sigma_{22}=&\frac{1}{4}\left[[\phi(t)]^{2}-2\phi(t)\phi(-t)+[\phi(-t)]^{2}\right]-\frac{1}{4}[\phi(2t)+\phi(-2t)-2]\\\
\sigma_{12}=\sigma_{21}=&\frac{1}{4i}\left[\phi(2t)-[\phi(t)]^{2}-\phi(-2t)+[\phi(-t)]^{2}\right].\end{split}$
(14)

######  Remark 4.

Note that when $X_{j},j=1,\dots,n$ are also symmetric and centered around the
origin $\boldsymbol{\Sigma_{z}}$ reduces to the form

$\boldsymbol{\Sigma_{z}}=\begin{bmatrix}\frac{1}{2}\left[1+\phi(2t)-2(\phi(t))^{2}\right]&0\\\
0&\frac{1}{2}[1-\phi(2t)]\end{bmatrix}$ (15)

This case was described in Section 3 of [11].

###### Proof of proposition 1.

First, rewrite $u_{j}$ and $v_{j}$:

$\begin{split}u_{j}=\frac{1}{2}(e^{itX_{j}}+e^{-itX_{j}})\\\
v_{j}=\frac{1}{2i}(e^{itX_{j}}-e^{-itX_{j}})\end{split}$ (16)

Since both of them are bounded, then $\boldsymbol{Z_{n}}(t)$ is a mean of
random vectors with components having the second moment. Thus, the
multidimensional central limit theorem immediately applies. Therefore, it is
only left to find the entries of the matrix $\boldsymbol{\Sigma_{z}}$ (13).

$\begin{split}\operatorname{Cov}(u_{j},v_{j})=\operatorname{Cov}\left[\frac{1}{2}(e^{itX_{j}}+e^{-itX_{j}}),\frac{1}{2i}(e^{itX_{j}}-e^{-itX_{j}})\right]&=\\\
\frac{1}{4i}\left[\operatorname{Cov}(e^{itX_{j}},e^{itX_{j}})-\operatorname{Cov}(e^{itX_{j}},e^{-itX_{j}})+\operatorname{Cov}(e^{itX_{j}},e^{-itX_{j}})-\operatorname{Cov}(e^{-itX_{j}},e^{-itX_{j}})\right]&=\\\
\frac{1}{4i}\left[\operatorname{Cov}(e^{itX_{j}},e^{itX_{j}})-\operatorname{Cov}(e^{-itX_{j}},e^{-itX_{j}})\right]&\end{split}$
(17)
$\begin{split}\operatorname{Var}\left(e^{itX_{j}}\right)&=\mathrm{E}\left[e^{2itX_{j}}\right]-\left(\mathrm{E}\left[e^{itX_{j}}\right]\right)^{2}=\phi(2t)-[\phi(t)]^{2}\\\
\operatorname{Var}\left(e^{-itX_{j}}\right)&=\phi(-2t)-[\phi(-t)]^{2}\end{split}$
(18)

Having substituted covariances in (17) using (18) one has

$\begin{split}\operatorname{Cov}(u_{j},v_{j})=\frac{1}{4i}\left[\phi(2t)-[\phi(t)]^{2}-\phi(-2t)+[\phi(-t)]^{2}\right]\\\
\end{split}$ (19)

Other elements of the covariance matrix $\operatorname{Var}(u_{j})$ and
$\operatorname{Var}(v_{j})$ are computed in the same manner:

$\begin{split}\operatorname{Var}(u_{j})=\mathrm{E}[u_{j}^{2}]-\left(\mathrm{E}[u_{j}]\right)^{2}=\frac{1}{4}\mathrm{E}\left(e^{2itX_{j}}+2e^{0}+e^{-2itX_{j}}\right)-\frac{1}{4}\left[\mathrm{E}\left(e^{itX_{j}}+e^{-itX_{j}}\right)\right]^{2}&=\\\
\frac{1}{4}\left[\phi(2t)+2+\phi(-2t)-(\phi(t))^{2}-2\phi(t)\phi(-t)-(\phi(-t))^{2}\right]&\end{split}$
(20)
$\begin{split}\operatorname{Var}(v_{j})=\mathrm{E}[v_{j}^{2}]-\left(\mathrm{E}[v_{j}]\right)^{2}=-\frac{1}{4}\mathrm{E}\left[e^{2itX_{j}}-2e^{0}+e^{-2itX_{j}}\right]-\left(\frac{1}{2i}\mathrm{E}\left[e^{itX_{j}}-e^{-itX_{j}}\right]\right)^{2}&=\\\
-\frac{1}{4}\left[\phi(2t)-2+\phi(-2t)\right]-\frac{-1}{4}\left[\phi(t)-\phi(-t)\right]^{2}&=\\\
\frac{1}{4}\left[[\phi(t)]^{2}-2\phi(t)\phi(-t)+\phi(-t)^{2}\right]-\frac{1}{4}[\phi(2t)+\phi(-2t)-2]&\end{split}$
(21)

∎

Proposition 1 provides a limiting distribution for a two-dimensional vector
$\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}}$, making the result inconvenient for
statistical inference. To address that, we define an analogy of Mahalanobis
distance– a statistic characterising how far $\boldsymbol{Z_{n}}$ is from
$\boldsymbol{Z_{0}}$. The following theorem provides a limiting distribution
for this statistic and plays the key role in the parametric inference.

######  Theorem 2.

Fix $t$. Assume that
$\boldsymbol{\hat{\Sigma}_{z}}\xrightarrow{d}\boldsymbol{\Sigma_{z}}$ and that
$\boldsymbol{Z^{(1)}_{n}}(t)$, $\boldsymbol{Z^{(2)}_{n}}(t)$ come from iid
$\boldsymbol{X_{1}}$ and $\boldsymbol{X_{2}}$ with distribution determined by
$\boldsymbol{Z_{0}}(t)$. Then

$\displaystyle\begin{split}\text{a)
}&\hat{T}_{i,0}:=n\left[\boldsymbol{Z^{(i)}_{n}}(t)-\boldsymbol{Z_{0}}(t)\right]^{{}^{\prime}}\boldsymbol{\hat{\Sigma}_{z}}^{-1}\left[\boldsymbol{Z^{(i)}_{n}}(t)-\boldsymbol{Z_{0}}(t)\right]\xrightarrow{d}\mathcal{\chi}^{2}(2)\quad
i=1,2\\\ \text{b)
}&\hat{T}_{1,2}:=n\left[\boldsymbol{Z^{(1)}_{n}}(t)-\boldsymbol{Z^{(2)}_{n}}(t)\right]^{{}^{\prime}}\boldsymbol{\hat{\Sigma}_{z}}^{-1}\left[\boldsymbol{Z^{(1)}_{n}}(t)-\boldsymbol{Z^{(2)}_{n}}(t)\right]\xrightarrow{d}\mathcal{\chi}^{2}(2)\end{split}$
(22)

###### Proof.

a) According to Proposition 1 and properties of weak convergence

$\left(\frac{1}{n}\boldsymbol{\hat{\Sigma}_{z}}\right)^{-1/2}(\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}})\xrightarrow{d}\mathcal{N}(0,I).$
(23)

Hence, by Continuous mapping theorem,

$\begin{split}\left[\left(\frac{1}{n}\boldsymbol{\hat{\Sigma}_{z}}\right)^{-1/2}(\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}})\right]^{\prime}\left[\left(\frac{1}{n}\boldsymbol{\hat{\Sigma}_{z}}\right)^{-1/2}(\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}})\right]=\\\
n(\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}})^{\prime}\left(\boldsymbol{\hat{\Sigma}_{z}}^{-1/2}\right)^{\prime}\boldsymbol{\hat{\Sigma}_{z}}^{-1/2}(\boldsymbol{Z_{n}}-\boldsymbol{Z_{0}})\xrightarrow{d}\mathcal{\chi}^{2}(2)\end{split}$
(24)

b) follows from independence of $\boldsymbol{Z^{(1)}_{n}}$ and
$\boldsymbol{Z^{(2)}_{n}}$. ∎

###  4.3 Hypothesis testing based on ECFs of colour intensities

Now, let us return to the data extraction (Figure 6). Assume we have an image
of forest and a test image. Having extracted matrix $\boldsymbol{X}$, one can
estimate parameters of the stable distributions of the forest picture channels
and compute vectors $\boldsymbol{Z_{0}}^{r}(t),\boldsymbol{Z_{0}}^{g}(t)$ and
$\boldsymbol{Z_{0}}^{b}(t)$. Similarly, after obtaining $\boldsymbol{X}$ from
the test image, empirical characteristic functions are computed and
subsequently $\boldsymbol{Z_{n}}^{r}(t),\boldsymbol{Z_{n}}^{g}(t)$ and
$\boldsymbol{Z_{n}}^{b}(t)$ are obtained. Note that we use CFs and ECFs of
each channel separately, so dependencies between colours are ignored so as
spatial dependencies.

Figure 6: Data processing for the parametric testing

Hypothesis testing in the parametric setting is depicted on Figure 7. A test
image is compared with each of the forest images in the training set. When
comparing a test image and a forest one,
$[\boldsymbol{Z_{0}}^{r}(t),\boldsymbol{Z_{0}}^{g}(t),\boldsymbol{Z_{0}}^{b}(t)]$
and
$[\boldsymbol{Z_{n}}^{r}(t),\boldsymbol{Z_{n}}^{g}(t),\boldsymbol{Z_{n}}^{b}(t)]$
are used to test separately that the colour intensities of the test image are
those of forest:

$\begin{split}T^{r}:=n(\boldsymbol{Z_{n}}^{r}-\boldsymbol{Z_{0}}^{r})^{{}^{\prime}}\boldsymbol{\hat{\Sigma}_{z;r}}^{-1}(\boldsymbol{Z_{n}}^{r}-\boldsymbol{Z_{0}}^{r})\xrightarrow{d}\mathcal{\chi}^{2}(2)\\\
T^{g}:=n(\boldsymbol{Z_{n}}^{g}-\boldsymbol{Z_{0}}^{g})^{{}^{\prime}}\boldsymbol{\hat{\Sigma}_{z;g}}^{-1}(\boldsymbol{Z_{n}}^{g}-\boldsymbol{Z_{0}}^{g})\xrightarrow{d}\mathcal{\chi}^{2}(2)\\\
T^{b}:=n(\boldsymbol{Z_{n}}^{b}-\boldsymbol{Z_{0}}^{b})^{{}^{\prime}}\boldsymbol{\hat{\Sigma}_{z;b}}^{-1}(\boldsymbol{Z_{n}}^{b}-\boldsymbol{Z_{0}}^{b})\xrightarrow{d}\mathcal{\chi}^{2}(2)\end{split}$
(25)

It implies that in order to be tagged as forest, a test image must have all
three colour samples passing their corresponding test. Since the limiting
distributions in (25) are the same, thresholds for each channel are set to be
equal.

For each pair the test image and a forest image,
$T_{tf}:=\min(T^{r},T^{g},T^{b})$ is computed and then
$T_{min}:=\min_{f}(T_{tf})$ is obtained and compared with the threshold.

Figure 7: Classification in the parametric model

###  4.4 Training and cross-validation of the parametric model

The training procedure starts with finding
$\boldsymbol{\hat{\Sigma}_{z;r}}$’s, $\boldsymbol{Z_{0}}$’s and
$\boldsymbol{Z_{n}}$’s for forest images and $\boldsymbol{Z_{n}}$’s for non-
forest ones in formulas (25). After that, training proceeds as in 3.3: T
statistics for pairs of images are computed and the threshold is set to
minimise the misclassification rate. Throughout classification,
$\boldsymbol{\hat{\Sigma}_{z;r}}$’s and $\boldsymbol{Z_{0}}$’s will be
available in advance, so that only $\boldsymbol{Z_{n}}$’s will be left to
compute to perform hypothesis testing based on (25). For that reason we do not
build the parametric framework on Theorem 2 b) with
$\boldsymbol{\hat{\Sigma}_{z;r}}$’s obtained from both the test image and a
training forest image: estimation of parameters of stable distribution is
quite computationally demanding. Although, this approach would provide better
accuracy.

Cross-validation for the parametric model is summarised in Algorithm 2.

Algorithm 2 Parametric cross validation procedure

1:Inputs: Data set of images labelled forest or non-forest, number of folds
$k$, max threshold $T_{max}$, parameter for CF’s ($t$ in (9)).

2:For every forest image compute $Z_{0}$, parameters and the covariance
matrix.

3:For every non-forest image compute $Z_{n}$.

4:Randomly shuffle images and group them into $k$ folds.

5:for $i=1$ to $k$ do

6: for each image in fold $i$ do

7: Compute 3 $\hat{T}$’s from (25), between the image and all images labelled
forest in the remaining $k-1$ folds.

8: Keep only the smallest $\hat{T}$ for the image.

9: end for

10:end for

11:for $t$ in $0$ to $T_{max}$ do

12: For all images and all $\hat{T}$’s, tag images with $\hat{T}<t$ as forest
and $\hat{T}>t$ as non forest.

13: Compute the accuracy for the given threshold $t,$ $\text{Acc}(t).$

14:end for

15:return $\text{arg}\max_{t}\text{Acc}(t)$ and $Z_{0}$’s, parameters and the
covariance matrices of forest images.

##  5 Comparison of the methods

There are features typical for both of our methods. First, if by accident
there is a non-forest image in the training forest set, test images of its
type start to be tagged as forest because they pass the corresponding test. As
a consequence, additional false positive results will occur. Specially
designed cross-validation may be exploited to eliminate this problem. Second,
the setting is multiple hypothesis testing, so it could naturally be re-set to
classify forests images into multiple types such as rain forest, coniferous
forest, etc. Another thing is that clouds and sunbeam scattering affect the
algorithms because they introduce noise to the distributions of objects. It
could be cured by atmospheric correction [7]. Also, both models increase
accuracy with growth of the numbers of points in the training data and test
data. Improvement of prediction accuracy could be achieved in three ways.
Increasing resolution, e.g. by switching from satellite data to aerial
imagery, would add more pixels to objects under consideration while adding new
tagged images to the training set would provide better training of thresholds
and more tagged classes. In practice, merely increasing the size of an image
with fixed resolution (adding more area) works well only in the case of
training pictures: enlarging test images would imply observing pictures with
objects belonging to different classes and tagging some areas incorrectly
because of that. And lastly, unlike classical machine learning methods, the
models developed in the present article are robust to data transformations
such as rotation and reflection.

There is also a number of differences between the two models. An advantage of
the parametric model is that it can work on incomplete data: e.g. when only 1
channel is available or different channels are available at different times.
The non-parametric method needs at least two channels simultaneously. On the
other hand, the disadvantage of the parametric model is that the algorithm
generally requires more pixels per image in training as it involves parameter
estimation and because we chose to use theorem 2 a). The estimation procedure
also makes the parametric method slower than the other one. A con of the non-
parametric method is that the model does not take into account the actual
distributions of colour intensities, which can contribute to additional false
positives.

To summarise the comparison, we point out that at the current stage of
research it is not possible to claim superiority of one method over the other
as they both have pros and cons.

##  6 Empirical application

Sentinel-2 satellite images are used to illustrate how the new Mahalanobis
distance classifier (MDC) and stable distribution classifier (SDC) perform in
comparison to some other commonly used classification methods. The methods
used for comparison are support vector machine (SVM), logistic regression (LR)
and deep neural network (DNN). Three experiments with real data are shown
here: two comparing accuracy of the methods and one giving visual
representation of performance. As the source of terrain images, the region
Occitanie in France is used, more precisely, vicinity of the city of Albi.

In each of the experiments SVM is implemented with the radial basis function
kernel and the DNN has four hidden layers of sizes 350, 200, 100 and 50
together with ReLu as activation function. The feature space for the
comparison methods consists of all pixels and their RGB colour values. In
practise, one would typically use some dimensionality reduction technique,
e.g. PCA, before implementing the classification methods, but since no such
technique is used for our new methods we want to be consistent with using all
of the available data. Considering that the images are rather small this
should not cause any issues.

In the first two experiments we look at how the methods perform using a
training and test split of the available dataset. In the first one the
training set consists of 510 images of sizes $20\times 20$ pixels with an
equal number forest and non-forest images. The test set is also equally
balanced between forest and non-forest images but it consists of 200 images of
sizes $20\times 20$ pixels. The second experiment is similar, except that each
$20\times 20$ image is split into four equal pieces, making 2040 images of
sizes $10\times 10$ pixels in the training set and 1000 $10\times 10$ images
in the test set. It should be noted that all of the images are collected from
the same region and under similar conditions. This means that the models are
not trained for different seasons or different types of forests. However, it
would be straightforward to make the classification models more general by
including more diverse data.

The performance measure used to evaluate the methods is the error rate, i.e.,
the number of incorrectly classified images in the test set divided by the
total number of test images. Figure 8 illustrates how the methods perform when
the training set varies between 10 % and 100 % of the total training set using
images of sizes $10\times 10$ pixels and $20\times 20$ pixels, respectively.
Even though the training set size varies, an equal number of forest and non-
forest images are always used for training. The test set is kept unchanged.

(a) $10\times 10$ pixels

(b) $20\times 20$ pixels

Figure 8: Forest classification error rates depending on the training set size
using MDC, SDC, SVM, LR and DNN for different image sizes.

It is clear in Figure 8 that the new methods outperform all of the comparison
methods. The error rate drops to zero almost immediately in the $10\times 10$
pixel case when using SDC whereas the error rate corresponding to MDC
gradually decreases from 2.4% and reaches zero when the training set size is
about 1800. All of the standard classification methods perform worse and their
error rates do not fall below $2\%$. The best performing methods among the
comparison ones are DNN and SVM. They show similar results apart from the fact
that SVM is more robust to changes in performance when adding new data to the
training set. Analogous results are presented in the $20\times 20$ pixel case.
Here, MDC generally has a lower error rate, and it drops to zero after the
size of the training set reaches 300. On the contrary, the error rate of SDC
is higher than in the previous case. Moreover, when the size of the training
set is lower than 62, the error rate is higher than those of DNN and SVM. When
the size of the training set reaches 62, the error sharply drops and becomes
lower or comparable to those two.

In the third experiment we visualise how the new methods and DNN (one of the
best performing comparison methods) behave when every sub-image of size
$10\times 10$ pixels is classified of an image of size $2000\times 2000$
pixels. Sub-images that are classified as forest are marked as white and non-
forest sub-images are marked as black. This makes it possible to get a better
feeling for how the methods differ and identify their strengths and
weaknesses. The same (complete) training set as in the first part is used to
calibrate the methods. This training set does not contain images that are part
of the $2000\times 2000$ pixels image but they were collected form a nearby
region under similar conditions. The result of the classification is
illustrated in Figure 9.

(a) Original image

(b) MDC

(c) SDC

(d) DNN

Figure 9: Forest classification using MDC, SDC and DNN of an $2000\times 2000$
pixel image using sub-images of sizes $10\times 10$ pixels. White corresponds
to forest and black corresponds to non-forest.

Generally, all three models performed well, though there are differences in
their outcomes. One can see that unlike MDC and SDC, DNN is not so sensitive
to the intensity of the green colour. It is best illustrated by performances
of the models on forest covered by a shadow of the cloud in the right part of
Figure 9. MDC and SDC incorrectly treat the shadowed vegetation as non-forest,
while DNN tags it as forest. Shades from clouds generally could be accounted
for by cloud detection, which is a solved problem in Earth image processing
(see e.g. [7]). In some other cases DNN tags non-forest vegetation as forest.
In particular, in Figure 10 there is forest rounded by non-forest vegetation.
Unlike MDC and SDC, DNN fails to detect some non-forest parts of the circle.

(a) Original image

(b) MDC

(c) SDC

(d) DNN

Figure 10: Central parts of images in Figure 9.

## References

  * [1] PEPS - Plateforme d’Exploitation des Produits Sentinel.  URL https://peps.cnes.fr.
  * Bhaskar et al. [2010] H. Bhaskar, L. Mihaylova, and A. Achim.  Video foreground detection based on symmetric alpha-stable mixture models.  _IEEE Transactions on Circuits and Systems for Video Technology_ , 20(8):1133–1138, 2010.
  * Bibby et al. [1979] J. Bibby, K. Mardia, and J. Kent.  _Multivariate analysis_.  Academic Press, 1979.
  * Drusch et al. [2012] M. Drusch, U. Del Bello, S. Carlier, O. Colin, V. Fernandez, F. Gascon, B. Hoersch, C. Isola, P. Laberinti, P. Martimort, et al.  Sentinel-2: ESA’s optical high-resolution mission for GMES operational services.  _Remote sensing of Environment_ , 120:25–36, 2012.
  * GDAL/OGR contributors [2021] GDAL/OGR contributors.  _GDAL/OGR Geospatial Data Abstraction software Library_.  Open Source Geospatial Foundation, 2021.  URL https://gdal.org.
  * Koutrouvelis [1980] I. A. Koutrouvelis.  Regression-type estimation of the parameters of stable laws.  _Journal of the American statistical association_ , 75(372):918–928, 1980.
  * Main-Knorn et al. [2017] M. Main-Knorn, B. Pflug, J. Louis, V. Debaecker, U. Müller-Wilm, and F. Gascon.  Sen2Cor for Sentinel-2.  In _Image and Signal Processing for Remote Sensing XXIII_ , volume 10427, page 1042704. International Society for Optics and Photonics, 2017\.
  * McCulloch [1986] J. H. McCulloch.  Simple consistent estimators of stable distribution parameters.  _Communications in Statistics-Simulation and Computation_ , 15(4):1109–1136, 1986.
  * Nolan [1997] J. P. Nolan.  Numerical calculation of stable densities and distribution functions.  _Communications in statistics. Stochastic models_ , 13(4):759–774, 1997.
  * Ortega Adarme et al. [2020] M. Ortega Adarme, R. Queiroz Feitosa, P. Nigri Happ, C. Aparecido De Almeida, and A. Rodrigues Gomes.  Evaluation of deep learning techniques for deforestation detection in the Brazilian Amazon and cerrado biomes from remote sensing imagery.  _Remote Sensing_ , 12(6):910, 2020.
  * Press [1972] S. J. Press.  Estimation in univariate and multivariate stable distributions.  _Journal of the American Statistical Association_ , 67(340):842–846, 1972.
  * Reiche et al. [2015] J. Reiche, S. De Bruin, D. Hoekman, J. Verbesselt, and M. Herold.  A Bayesian approach to combine Landsat and ALOS PALSAR time series for near real-time deforestation detection.  _Remote sensing_ , 7(5):4973–4996, 2015.
  * Roy et al. [2014] D. Roy, M. Wulder, T. Loveland, W. C.E., R. Allen, M. Anderson, D. Helder, J. Irons, D. Johnson, R. Kennedy, T. Scambos, C. Schaaf, J. Schott, Y. Sheng, E. Vermote, A. Belward, R. Bindschadler, W. Cohen, F. Gao, J. Hipple, P. Hostert, J. Huntington, C. Justice, A. Kilic, V. Kovalskyy, Z. Lee, L. Lymburner, J. Masek, J. McCorkel, Y. Shuai, R. Trezza, J. Vogelmann, R. Wynne, and Z. Zhu.  Landsat-8: Science and product vision for terrestrial global change research.  _Remote Sensing of Environment_ , 145:154–172, 2014.  ISSN 0034-4257.  doi: https://doi.org/10.1016/j.rse.2014.02.001.  URL https://www.sciencedirect.com/science/article/pii/S003442571400042X.
  * Samorodnitsky and Taqqu [2017] G. Samorodnitsky and M. S. Taqqu.  _Stable Non-Gaussian Random Processes. Stochastic Models with Infinite Variance. Stochastic Modeling_.  Routledge, 2017.
  * Shermeyer and Haack [2015] J. Shermeyer and B. Haack.  Remote sensing change detection methods to track deforestation and growth in threatened rainforests in Madre de Dios, Peru.  _Journal of Applied Remote Sensing_ , 9(1):096040, 2015.
  * Wang et al. [2008] C. Wang, M. Liao, and X. Li.  Ship detection in SAR image based on the alpha-stable distribution.  _Sensors_ , 8(8):4948–4960, 2008.

##  Appendix A Data pre-processing technology

We have built a preprocessing procedure of creating RGB composites from
Sentinel satellite images using Python programming language.

First, we download Sentinel 2 Level-1C products from PEPS[1] with the means of
EODAG. The PEPS (Plateforme d’Exploitation des Produits Sentinel) platform,
redistributes the products of Sentinel satellites, S1A, S1B, S2A and S2B, S3A
and S3B from COPERNICUS, the European system for Earth monitoring. EODAG
(Earth Observation Data Access Gateway) is a command line tool and a Python
package for searching and downloading remotely sensed images while offering a
unified API for data access regardless of the data provider.

Second, for each product we do a preprocessing procedure of creating RGB
composites from Sentinel 2 Level-1C products. Using rasterio
package111https://rasterio.readthedocs.io/en/latest/ we extract red, green and
blue bands from channels B04, B03 and B02 respectively and combine them in a
single file using GeoTIFF file
format222https://gdal.org/drivers/raster/gtiff.html#raster-gtiff.

Third, we compress RGB-composites with lossless compression algorithm to
decrease the amount of storage space needed for data using gdat_translate
utility [5].

Lastly, in order to obtain images with sufficient brightness, we normalize
values in the channels to range from 0.0 to 1.0 by dividing them by the
constant 2000. We truncate values with normalised intensity more than 1.0 by
setting them to 1.0, which primarily affects bright white clouds and conforms
with typical ”crop and scale” image processing procedure. Deforestation
detection algorithms use this normalised values.

##  Appendix B Distributions of colour intensities of other natural objects

In this section we explore two RGB images: one- of a mountain and another one-
of sea. Each picture is accompanied by empirical density functions plotted for
all channels. For comparison, best fitting stable, gamma and normal
distributions are shown.

(a) Mountain image analysed

(b) Comparison of densities for red color intensity

(c) Comparison of densities for green color intensity

(d) Comparison of densities for blue color intensity

Figure 11: Plots of empirical and estimated probability density functions for
distribution of colour intensities in an image of mountain.

(a) Sea image analysed

(b) Comparison of densities for red color intensity

(c) Comparison of densities for green color intensity

(d) Comparison of densities for blue color intensity

Figure 12: Plots of empirical and estimated probability density functions for
distribution of colour intensities in an image of sea.