DNA Storage: A Promising Large Scale Archival Storage?

Deoxyribonucleic Acid (DNA), with its high density and long durability, is a promising storage medium for long-term archival storage and has attracted much attention. Several studies have verified the feasibility of using DNA for archival storage with a small amount of data. However, the achievable storage capacity of DNA as archival storage has not been comprehensively investigated yet. Theoretically, the DNA storage density is about 1 exabyte/mm3 (109 GB/mm3). However, according to our investigation, DNA storage tube capacity based on the current synthesizing and sequencing technologies is only at hundreds of Gigabytes due to the limitation of multiple bio and technology constraints. This paper identifies and investigates the critical factors affecting the single DNA tube capacity for archival storage. Finally, we suggest several promising directions to overcome the limitations and enhance DNA storage capacity.

INTRODUCTION An indisputable fact is that we are surrounded by more and more available data. Meanwhile, various applications such as IoT, autonomous vehicles, video streaming, scientific discovery, etc. need to preserve the increasingly large available datasets for future data analysis/mining in a longer duration. Unfortunately, the development of storage technologies/devices is not matching with the booming storage demand. Moreover, traditional storage media are not durable enough (e.g., disks for 3-5 years, tapes for 10-15 years) for the increasingly longer data retention duration. Therefore, data must be migrated from the existing obsolete devices to new devices when the current devices wear out every few years. The data migration and device replacement cost can be over half of the total cost of ownership even in cloud archival storage [1]. To meet the burgeoning storage demand, we are looking for a new storage medium with high density and long durability.
DNA has emerged with its ultra-high storage density (i.e., a theoretical storage density of about 1 exabyte/mm 3 [2]). When keeping DNA in a proper environment, it can last for centuries [1] and thus save an enormous amount of maintenance cost for data transfer and device replacement. As DNA synthesis/sequencing (i.e., write/read) prices keep decreasing, DNA storage is identified as a promising archival storage medium.
Several researchers have demonstrated the feasibility of DNA archival storage. Their investigations focused on error-tolerant encoding schemes [3] [4][5] [6] and random accessibility with wet-lab experiments of a small amount data [2] [7]. Encoding schemes are responsible for the conversion between digital data and DNA sequences. Random access allows retrieval of target data without sequencing the entire DNA pool which greatly saves the number of slow and expensive sequencings needed. However, the achievable storage capacity of using DNA for archive based on current biotechnologies remains unknown. In DNA storage, Data is in the form of powder for a longer duration. When sequencing, people liquidize the powders and use a drop of the liquid for sequencing. To avoid long diffusion time and ensure sequencing effectiveness, people usually separate a large powder mixture into multiple physical tubes. A tube is a basic DNA storage device, and DNA storage capacity can be measured by DNA tube capacity. A larger single tube capacity means larger overall DNA storage capacity or higher storage density. In this paper, we study the storage capacity of a single tube.
One important aspect of limiting the real DNA storage capacity is that DNA storage is errorprone. Due to the mutation nature of DNA, errors like insertion, deletion, and substitution can happen during synthesis and sequencing processes. Some special DNA patterns may increase the probabilities of these errors such as long homopolymer, secondary structure, etc. Therefore, it is important to avoid certain biological constraints to decrease the error rate of DNA storage using special DNA encoding schemes. We also convert digital data to error-tolerant DNA sequences with a shorter DNA strand length (the current practical strand length for DNA storage is between 200 to 300 nt) since the longer the DNA strand length, the higher the error rate incurred.
For DNA storage used for archive, we consider random access is a must. The most prevalent sequencing technology, Polymerase Chain Reaction (PCR), may also introduce extra errors. PCRbased random access is a process that amplifies only target DNA strands to enough concentration to be sequenced out. It uses a pair of short synthetic nucleotide sequences called primers to tag a set of DNA strands (the number of DNA strands in a set is considered as a parallel factor since these DNA strands can be read out in one sequencing process). Primers should be sufficiently different from each other such that random access can accurately retrieve the target strands. However, the requirement of sufficient difference between primers reduces the number of usable primers. Similarly, the encoded DNA strands must not contain certain DNA short sequences that are also included in primers. This will further limit the DNA storage tube capacity. We will discuss these capacity limitations in the next section.
In this paper, we investigate the single tube capacity of DNA archival storage. We first list possible factors that influence DNA tube capacity. We then generate a large enough primer library to measure the practical tube capacity with different encoding schemes and different types of data. Due to the extremely high synthesis/sequencing cost for a large amount of data at the current time, the measurements are done without wet labs. We follow a similar setup in existing studies including the value of the parallel factor, bioconstraints, collision definitions, etc., and use NCBI BLAST [8], one of the most widely used genomic sequence analysis tools, to check the feasibility of DNA storage. According to our measurements, DNA tube capacity can only reach hundreds of gigabytes. Finally, based on the thorough analysis of the factors restricting tube capacity, we suggest several potential ways to improve the DNA tube capacity.
The rest of this paper is organized as follows. We first discuss the necessary background. Then we discuss several capacity limitations of DNA archival storage. After that, we present our measurements of the practical DNA tube capacity. Finally, we propose several potential enhancement approaches to enlarge DNA tube capacity.

BACKGROUND
This section explains basic background of DNA storage synthesis and sequencing workflows and PCR procedure for random accesses. Figure 1 summarizes the workflows of synthesis and sequencing of DNA storage. A DNA storage system takes digital data as input. It first adds Error Correction Codes (ECC) [9] to digital data to provide recoverability in case of existing synthesis/sequencing errors. It then coverts the digital data into a sequence of nucleotides based on an encoding scheme and cuts the sequence After that, it flanks a primer pair and an internal index onto each payload sequence to form a DNA strand. The primer pair indicates and enables the beginning and end of amplification during the PCR process. Since one primer pair can be shared among multiple DNA strands, an internal index is used to identify the order of DNA strands that share one primer pair. Each DNA strand will be chemically synthesized nucleotide by nucleotide and stored in physical tubes/pools. To recover digital data back, target DNA strands with a specific primer pair can be selectively amplified via PCR and sequenced via next-generation sequencing (NGS) or other sequencing techniques.

DNA Storage Synthesis/Sequencing Workflows
Finally decoded back to digital information.

PCR Procedure
PCR is a method that can amplify target DNA strands to enough concentration for DNA sequencing. It enables random access of target data without sequencing the whole storage pool. A PCR process usually consists of 25 to 35 cycles. Each cycle doubles the number of target strands with the specific primer pair. As shown in Figure  2(a), a standard PCR cycle follows three steps: denaturation, annealing, and elongation. In the denaturation step, DNA strands are denatured into two single-stranded DNA molecules. In the annealing step, complementary primers will be added to the reaction solution and bind to the target denatured DNA molecules. In the elongation step, with the help of dNTPs, the Taq DNA polymerase adds nucleotides to the growing DNA strand. Typically, the single tube capacity of PCR-based DNA storage is affected by the parallel factor, payload encoding density, payload length, and the number of usable primers.
Parallel factor is the number of strands that share a primer pair. A larger parallel factor produces a larger capacity but requires more PCR cycles to amplify all strands to enough concentration. An overlarge parallel factor can cause data loss if some strands fail to be amplified to enough concentration. Existing studies have adopted 1.55 × 10 6 as a practical parallel factor [3] [10]. Parallel factor together with the number of useable primers can determine how many DNA strands per tube.
Payload encoding density is the number of bits each nucleotide can represent. To reduce the potential errors during synthesis/sequencing, an encoding scheme needs to avoid specific patterns such as long homopolymer, low/high GC content ratio, etc. As a result, an ideal encoding den-sity, 2 bits/nt, cannot be achieved/used. Encoding schemes are required to convert digital data without generating long homopolymers and unbalanced GC content. A theoretical upper bound of encoding density is proofed around 1.98 [4].
Payload length is the length of a DNA strand subtracting the length of primer pair, internal index, and logical redundancy (e.g., ECC). The overall DNA strand length is currently restricted by synthesis technology. A practical DNA strand length is usually no more than 300 nt. Otherwise, it is likely to have more synthesis errors [9]. 15% redundancy is usually enough to recover data from errors [9]. Therefore, excluding redundancy, internal index and primer pair, payload length in a DNA strand can be varied around between 100 nt -250 nt.
The number of usable primers refers to how many primers can be safely used for PCR-based random access in a single tube. Theoretically, with the primer length of 20 bases, the number of primers can reach up to 4 20 . However, due to bio-constraints, the number of primers used in a DNA storage tube is much less than expected. The bio-constraints include homopolymer length shorter than 4, balanced GC content, and the absence of secondary structures. Moreover, since primers serve as unique tags in PCR, they must be sufficiently different from each other to avoid potential crosstalk. The details are shown in the next section when we discuss building a primer library with enough primers.
Even if primers are designed orthogonal to each other, they can still interact with payloads, which are referred to as primer-payload collisions [7]. A primer-payload collision happens when a primer and a payload have a pair of almost identical subsequences. The pair of almost identical subsequences are usually longer than 12 bases and have at most two mismatches or gaps. In each PCR cycle, wrong annealing can happen if any payload has a collision with the target primer. The wrong annealing will amplify some irrelevant payloads and, more importantly, consume the limited PCR reagents (e.g., complementary primers used to bind with the target primer and free DNA bases waiting for complementing new double-helix strands). Figure 3 shows the distribution of primers with different numbers of collisions. The input data is a 135MB video (https://www.youtube.com/watch?v=qybUFnY7Y8w), and the encoding scheme is Blawat code [3] (see the encoding scheme, primer generating, and collision checking details in next section). Although input data is only 135MB, 27,658 primers (i.e., 98.77% primers in our primer library) have collisions with those payloads, and the average number of collisions per collided primer is 155.45. If the input data keeps scaling up, most collided primers likely have hundreds or thousands of collisions. In this case, PCR amplification is very vulnerable to primer-payload collisions. The target primer must compete for the limited PCR reagents with thousands of collided payloads in each PCR cycle. Even if a primer has only a few collisions at the beginning, the collisions can still significantly impact the final result. That is because PCR is an exponential augmentation process. Slight amplification variations can exacerbate cycle by cycle and consume ever more PCR reagents. Eventually, the resultant solution will be saturated with plenty of irrelevant sequences but very few target sequences. As a result, the sequencing for target sequences is inhibited.
Currently, the collision problem has not drawn enough attention for two reasons. Firstly, existing works are mostly based on a small scale data (i.e., several Kilobytes to several Megabytes) [12][5] [2][3] [4]. Secondly, the data in DNA pool is not diverse enough. Existing works build a DNA pool by mixing a few synthetic target data with background non-target data. The non-target data are usually homogenous DNA strands made from PCR amplification because PCR amplification is much cheaper than DNA synthesis. Because non-target data is almost identical, it has less chance of colliding with primers during PCR. However, DNA payloads in a tube should be more diversified when DNA storage is put to practical use and continues to scale up in capacity. In this case, the impact of primerpayload collisions will be very severe. For example, a previous study [7] reported a PCR failure because of primer-payload collisions. They expanded a DNA tube from a single file to multiple files and failed to sequence target data when there were nine files in the tube (the target data is 17.4% of the whole data in the tube). To ensure a precise and effective PCR, it is essential only to use primers that have no collisions with payloads.
Given parallel factor, payload length, and payload encoding density, we can calculate the amount of data each primer pair can accommodate. Since the above three factors are all restricted by the current biological technologies, the number of usable primers becomes a major factor in deciding practical DNA tube capacity. Since the primer-payload collision is input data dependent, we plan to use real-world data as input to check how many primers are usable and what is the practical tube capacity.

Investigation of Practical Tube Capacity
This section illustrates the impact of primerpayload collisions and discovers the practical tube capacity.

Build A Large Primer Library
We first build a large enough primer library since the existing primary libraries are too small for our purpose. We adopt a primer design framework designated for DNA storage to generate our primer library [7]. The framework keeps generating random 20 base DNA sequences and evaluating the sequences with a set of criteria. Sequences can be added to the primer library only if they satisfy all the following criteria: • GC content ∈ [0.45, 0.55] • melting temperature ∈ [55 • C, 60 • C] • homopolymer length ≤ 3 • Hamming distance between any two primers ≥ 6 • no more than 10 bases inter-complementarity • no more than 4 bases intra-complementarity Because GC contributes more bonding power than AT, proper GC content is required to ensure stable binding of primers. Similarly, melting temperature, the temperature at which DNA duplex dissociates to be single-stranded, is also required for reliable primers binding. Hamming distance is required to ensure primers are sufficiently different. It helps reduce the likelihood of non-specific primer binding.
Intra-complementarity refers to two subsequences within a primer complementing each other in reverse order. Similarly, intercomplementarity refers to two complementary subsequences in two primers. More than 4 bases intra-complementarity or 10 bases intercomplementarity can cause a primer hybridizes with itself (i.e., hairpin) or with other primers (i.e., primer dimer) instead of annealing to the binding site of the target payloads. For example, primer "...ATGG...CCAT..." has intra-complementary subsequences with length 4.
This primer design framework has been verified by Organick et al. [7] with a pool of 3,240 small synthetic files. They successfully retrieved up to 48 files in a one-pot multiplex PCR experiment. They also predicted that the number of potentially generated primers is roughly 28,000. We follow their design framework in this study to generate a primer library of 28,000 primers. More primers may be available with different primer design rules, but 28,000 primers are enough for us to investigate the capacity loss due to primerpayload collisions. Hereafter, the practical tube capacity refers to the capacity with the primer library of 28,000 primers and is based on the current biotechnologies.

Encoded Digital Data
We collect 1,500 GB of data with five data types, including image, video, audio, eBook, and software archives (300 GB of each type). Image data is a subset of ImageNet [13]. Audio data contains 187GB of LibriSpeech [14] and some small collections from InternetArchive (https://www.youtube.com/watch?v=qybUFnY7Y8w).

May/June 2022
For other data types, we do not find any single dataset that is large enough and thus collect multiple small datasets from InternetArchive.
We implement the following five highly cited DNA encoding schemes. For some encoding schemes that have no legit name, we name them using the first author's name.
• Church [12] is one of the earliest encoding scheme. It encodes one bit to one base, and every bit has two optional bases to choose from (i.e., 0 to AC and 1 to TG The third base will be selected to ensure the first three bases will not be identical, and the fifth base will be selected to ensure the last two bases will not be identical. With these constraints, the longest possible homopolymer length is 3, and the encoding density is 1.6 bits/base. • Grass [5] first converts binary bits to Galois Field of size 47 (GF(47)). It then bases on a mapping table to encode GF(47) digits to base triplets in which the second base differs from the third base. Since no base triplet has identical second and third bases, the longest homopolymer length in the resultant DNA sequences is three. The encoding density of Grass encoding is 1.77 bits/base. • Fountain code is famous for its robustness [4].
It iterates over Luby transform to select data segments and XOR them bit-wisely. The resultant data will be directly translated to DNA sequences from {00,01,10,11} to {A,C,G,T}.
Then, it screens out the sequences with undesired GC content and long homopolymers. It keeps the Luby transform until generates a predefined number of DNA sequences. For reliability, the number of resultant DNA sequences is recommended to be 5%-10% more than the total input segments (in our case: 10%). Its encoding density is 1.81 bits/base if we count the overhead of the 10% redundancy.
In the encoding process, we first cut the collected five types of digital data into multiple small files of 10 MB. We will add files incrementally for collision checking and capacity identification later. We further partition each file into chunks of 239 bytes and apply Reed-Solomon code RS(255,239) to append 16 parity bytes to each chunk (i.e., each chunk has a size of 255). After that, we use each of the five encoding schemes to encode data chunks to DNA sequences and further partition these sequences into multiple DNA payloads with a length of 198 1 . We chose 198 because it is in the middle of the payload length range and is dividable by three, which is required by Grass code. The DNA payloads will further be appended with an internal index of 12 bases to indicate their ordering associated with a primer pair. We choose 12 bases for internal index length because it is long enough to generate more than 1.55×10 6 base sequences that will not collide with any primers in the primer library. We will choose 1.55 × 10 6 of them in a lexicographical order to form a mapping table to represent the ordering of DNA payloads associated with a primer pair. The payloads with internal indexes will be sent to BLAST to check possible collisions. When decoding DNA sequences back to digital data, the mapping table can be used to identify the order of payloads. Figure 4: Practical tube capacity as increasing input data size using Rotation code. The input data is ImageNet.

Check Collisions and Discover Tube Capacity
After building the primer library and encoding digital data to DNA payloads, the next step is filtering out primers with primer-payload collisions. Unfortunately, synthesizing such a large scale of data for collision checking is currently cost prohibitive (predicted at $21,000/MB by 2027 [16]). Instead, we follow Organick's work [7] and use BLAST [8] to check collisions between primers and payloads. Only primers without any collisions will be considered as usable primers.
We check collisions between all kinds of DNA payloads gathered in the last subsection and all primers from our primer library. Payloads are in multiple partitions and each partition represents 10.66MB input data (i.e., 10MB digital data + 0.66MB ECC). We sequentially check collision on each payload partition (i.e., increase of input data) and accumulate the collided primers. Meanwhile, more usable primers are needed to accommodate data as input data increases. The number of usable primers required to store a data set is the total data size divided by the encoding density, payload length, parallel factor, and two (i.e., two primers to form a primer pair). We keep increasing the input data until the number of remaining usable primers decreases to a number that is just capable of accommodating all input data. We consider the size of digital data in the input data (i.e., 10MB × number of input files) as the practical tube capacity of a specific data type for a particular encoding scheme. Figure 4 shows the achievable tube capacity when use Rotation code and ImageNet. The orange line in Figure 4 shows the increased number of primers needed to accommodate the input data. The blue line shows the decreased number of remaining usable primers. When input digital data grows to 215.42GB, the remaining usable primers are just enough to accommodate the input data. Due to space limitations, we do not show the figures of all five encoding schemes with five data types. The number of usable primers for other encoding schemes drops quickly, even with small amount of input data. We summarize all the results in the next subsection. Table 1 shows the practical, achievable tube capacity of the five encoding schemes with five types of data. Suppose all primers in the library are usable without primer-payload collisions. In that case, the five encoding schemes will achieve a tube capacity varying from 468GB to 848GB, depending on their encoding densities. However, if we remove the primers that have collisions, Blawat/Fountain/Church/Grass codes will suffer more than 99% practical tube capacity reduction and can have only hundreds of Megabytes to a few Gigabytes capacity. Rotation code shows the most collision resistance but still suffers 70% capacity reduction and achieves a capacity of around 215GB.

Summary of Tube Capacity
Blawat/Fountain/Church codes are close to direct conversion between binary data and DNA sequences (from 00, 01, 10, 11 to A, C, G, T). Since no special pattern is applied, the resultant DNA sequences are possible to collide with most primers. An example of the distribution of primerpayload collisions is shown in Figure 3. The encoding scheme is Blawat code, and the input data is the 135MB video. Despite the small data set, Blawat code has already collided with 27,658 primers, and the average number of collisions per collided primer is 155.45. However, Blawat code is still a little better than Fountain and Church codes because Blawat code only directly converts the first 6 bits of every 8 bits. Besides, Fountain code performs slightly better than Church code because of its higher encoding density.
By contrast, Rotation code has the second lowest encoding density but achieves the highest tube capacity among the five encoding schemes. That is because Rotation code sacrifices its encoding density to avoid consecutive identical bases (i.e., no homopolymers). Its DNA payloads are less likely to collide with primers that have consecutive identical bases. As for Grass code, it encodes digital data with base triplets in which the second base is different from the third base. This makes Grass's DNA payloads have less number of three consecutive identical bases and thus less likely to collide with primers that have three consecutive identical bases. However, Grass's DNA payloads still collide with primers that have shorter consecutive identical bases. Therefore, Grass code is worse than Rotation code but better than Blawat/Fountain/Church codes. Besides, the collision problem exists in all five data types, but no significant capacity differences were observed among different data types. That may be because different data types will eventually cover most of the possible digital sequences as the data size increases.

Potential Capacity Enhancement
Although several factors like strand length and parallel factor are temporarily limited by the current technologies, several potential approaches can enhance DNA tube capacity for archive. First, we can investigate new primer design frameworks to generate more primers. The existing primer frameworks are designed for biological constraints rather than the need for more usable primers. They apply strict rules on primer design to ensure PCR specificity. A primer design framework can more aggressively limit some primer design space to ensure both primer reliability and produce more usable primers.
Besides generating more usable primers, we can also work on the payload side. Payload sequences in DNA storage are synthetically generated based on an encoding scheme. Given a primer library, it is possible to design special encoding schemes to encode DNA sequences by avoiding primer-payload collisions. The previous section has shown that the Rotation scheme can avoid some collisions via tightening the longest homopolymer length from 3 to 1. Besides homopolymer length, we can further create other DNA sequence patterns to be different from those of primers. For example, we can tighten DNA sequence to have GC content ∈ [0.45, 0.50] to avoid collisions with primers that have GC content ∈ (0.50, 0.55]. That is, we can co-design primers and encoding schemes together. However, special encoding schemes that create different patterns may lose some encoding density. That is because a smaller number of DNA sequences are qualified to represent digital data.
Another cost-efficient way is to use some payload post-processing. If we can have variablelength payloads and make them break up at the collision locations, we can remove many collisions. If we can change payload content with a simple permutation of nucleotides, we can also remove a considerable number of collisions. These potential approaches work in different aspects and can be applied together to enhance DNA storage capacity.

Conclusion
In this paper, we identify the limitations of DNA tube capacity based on the-state-of-the-art synthesis and sequencing technologies. We then build our own large enough primer library and investigate the tube capacity of different encoding schemes with different types of data. The results show that the tube capacities of the current methods are reduced to hundreds of gigabytes, which is far less than what people expected. By analyzing the DNA storage capacity limitation and bio-constraints, we point out some potential ways to enhance DNA storage tube capacity. We hope with further investigations of these potential approaches and the development of new biotechnologies, DNA tube capacity will be sufficiently large for practical archival storage.


Project to Adapt: Domain Adaptation for Depth Completion from Noisy and Sparse Sensor Data

Depth completion aims to predict a dense depth map from a sparse depth input. The acquisition of dense ground truth annotations for depth completion settings can be difficult and, at the same time, a significant domain gap between real LiDAR measurements and synthetic data has prevented from successful training of models in virtual settings. We propose a domain adaptation approach for sparse-to-dense depth completion that is trained from synthetic data, without annotations in the real domain or additional sensors. Our approach simulates the real sensor noise in an RGB+LiDAR set-up, and consists of three modules: simulating the real LiDAR input in the synthetic domain via projections, filtering the real noisy LiDAR for supervision and adapting the synthetic RGB image using a CycleGAN approach. We extensively evaluate these modules against the state-of-the-art in the KITTI depth completion benchmark, showing significant improvements.

Introduction
Motivation. Active sensors such as LiDAR determine the distance of objects within a specified range via a sparse sampling of the environment whose density decreases quadratically with the distance. RGB cameras densely capture their field of view, however, monocular depth estimation from RGB is an ill-posed problem that can be solved only up to a geometric scale. The combination of RGB and depth modalities form a rich source for mutual improvements where each sensor can benefit from the advantage of the other.
Many pipelines have been proposed for a fusion of these two inputs [2][3][4][5][6][7]. Ground truth annotations for this task, however, require elaborate techniques, manual adjustments and are subject to hardware noise or costly and timeconsuming labeling. The most prominent publicly available data for this task [8] creates a ground truth by aligning consecutive raw LiDAR scans that are cleaned from measurement errors, occlusions, and motion artifacts in a post-processing step involving classical stereo reconstruction. Even after the use of this additional data and tedious processing, the signal is not noise-free as discussed in [8]. To avoid such annotations, some methods perform self-supervision [5,6,9], where arXiv:2008.01034v2 [cs.CV] 5 Aug 2020 (a) Overview of the approach (b) Multicamera set-up in CARLA Fig. 1: We investigate the depth completion problem without ground truth annotations in the real domain, which contains paired noisy and sparse depth measurements and RGB images. We highlight some noise present in the real data: see-through on the tree trunk and bicycle, self-occlusion on the bicycle and missing points on the van. In our approach, we leverage synthetic data with multicamera dense depth and RGB images. An overview of the multicamera set-up in CARLA used to simulate the real projection LiDAR artifacts is included. The Depth Camera acts as a virtual LiDAR and collects a dense depth, which is sparsified using real LiDAR binary masks and projected to either the Left Camera or Right Camera reference frame. Both the Left Camera and the Right Camera collect RGB information, used as part of the input data, and a dense depth map, used for supervision a photometric loss is employed with stereo or video data. The dependence on additional data such as stereo or temporal sequences brings other problems such as line-of-sight issues and motion artifacts from incoherently moving objects. Modern 3D engines are capable of rendering highly realistic virtual environments [10][11][12] with perfect ground truth. However, a significant domain gap between real and virtual scenes prevents from successful training on synthetic data only. Contributions and Outline. In contrast to the self-supervised methods [5,6,9], we propose to use a domain adaptation approach to address the depth completion problem without real data ground truth as shown in Figure 1. We train our method from the synthetic data generated with the driving simulator CARLA [12] and evaluate it on the real KITTI depth completion benchmark [8]. The real LiDAR data is noisy with the main source of noise being the see-through artifacts that occur after projecting the LiDAR's point cloud to the RGB cameras. We propose an approach to simulate the see-through artifacts by generating data in CARLA multicamera set-up, employing random masks from the real LiDAR to sparsify the virtual LiDAR sensor, and projecting from the virtual LiDAR to the RGB reference frame. We further improve the model by filtering the noisy input in the real domain, thus obtaining a set of reliable points that are used as supervision. Finally, to reduce the domain gap between the RGB images, we use a CycleGAN [1] to transfer the image style from the real domain to the synthetic one. We compare our approach to other state-of-the-art depth completion methods and provide a detailed analysis of the proposed components. The proposed domain adaptation for RGB-guided sparse-to-dense depth completion is a novel approach for the task of depth completion, which leads to significant improvements as demonstrated by the results. To this end, our main contributions are: 1. A novel domain adaptation method for depth completion that includes geometric and data-driven sensor mimicking, noise filtering and image style adaptation. We demonstrate that adapting the synthetic sparse depth is crucial for improving the performance, whereas RGB adaptation is secondary. 2. The improved state-of-the-art results for the KITTI depth completion benchmark, amongst ground truth free methods, by 6.4% RMSE and 6.3% MAE, and by 9.2% RMSE and 10.4% MAE when combining our pipeline with video self-supervision.

Related Work
We first review related works on depth estimation using either RGB or LiDAR, and then discuss depth completion methods using Convolutional Neural Networks.

Unimodal Approaches
RGB Images. RGB based depth estimation has a long history [13][14][15] reaching from temporal Structure from Motion (SfM) [16,17] and SLAM [18][19][20] to recent approaches that estimate depth from a static image [21][22][23][24]. Networks are either trained with full supervision [21,25] or use additional cameras to exploit photometric consistency during training [22,26]. Some monocular depth estimators leverage a pre-computation stage with an SfM pipeline to provide supervision for both camera pose and depth [27,28] or incorporate hints from stereo algorthms [29]. These approaches are in general tailored for a specific use case and suffer from domain shift errors, which has been addressed with stereo proxies [23] or various publicly available pre-training sources [24]. The estimated depth often suffers from over-smoothing [30] with wrongly inferred "flying pixels" in the free space close to depth discontinuities. Sparse Depth. While recent advantages in depth super-resolution [31,32] show good performance, they are not directly applicable to LiDAR data which is sparsely and irregularly distributed within the image. Similar to super-resolution, a rectangular grid for the sampling was assumed in [33]. The sampling grid of the sparse depth signal is crucial for the depth completion task [8], which can be provided as a mask to the network, thus helping to densify the input. While classical image processing techniques are used in in [34], an encoder-decoder architecture is applied for this task in [35]. Other approaches [36,37] design more efficient architectures to improve the runtime performance.

Depth Completion from RGB and LiDAR
Most recent solutions to depth completion leverage deep neural networks. These can be divided into (a) supervised and (b) self-supervised approaches. Supervision and Ground Truth. Usually, an encoder-decoder network is used to encode the different input signals into a common latent space where feature fusion is possible and a decoder reconstructs an output depth map [2-4, 7, 38]. Different additional random sampling strategies can increase the density of the input signal [2] while fusing 2D and 3D representations [39] can improve depth boundaries. The noise problem has been targeted with local and global information in [3]. Other methods [4,7] leverage different input modalities such as surface normals to increase the amount of diversity in the input data. The publicly available dataset KITTI [8,40] includes real driving scenes where a stereo RGB camera system is fixed on the roof of a car along with a LiDAR scanner that acquires data while the car is driving. A post-processing stage fuses several LiDAR scans and filters outliers with the help of stereo vision to provide labeled ground truth. While this process is intricate and time-consuming, further error is accumulated from calibration and alignment [8].
Self-Supervised Approaches. Another view either from a second camera or a video sequence can be used for self-supervision. Temporal information and mutually predicted poses between RGB frames were used in [5] for self-supervision with a photometric loss on the reprojected image. A probabilistic formulation was proposed in [6] with a conditional prior within a MAP estimation, which also leverages stereo information. A non-learning method was used in [9] to form a spatially dense but coarse depth approximation from the sparse points, where the coarse approximation was then refined using another network. A photometric loss was also used in [9], where a separate network predicted the poses between RGB frames obtained from a video sequence. Synthetic Data. For monocular depth estimation, two domain adaptation approaches used style-transfer methods [41,42]. Sparse-to-dense methods, however, have used synthetic data without any adaptation [4,6,43] so far. Training on synthetic data requires a high rendering quality [44]. To this end, synthesizing driving scenarios has also been researched: SYNTHIA [10] provides synthetic urban images together with semantic annotations, while Virtual KITTI [11] and its point cloud variant [45] constitutes a dataset with synthetic renderings and point clouds that closely match the videos of the KITTI dataset [40] including pixel-perfect semantic ground truth and depth labels. The CARLA simulator [12] was introduced allowing for photo-realistic simulations of driving scenarios, which we utilize to generate realistic RGB images. A LiDAR simulator using ray-casting and a learning process to drop points was proposed in [46], which was tested in detection and segmentation tasks, but is not publicly available. While LiDAR scans can also be simulated with CARLA via ray-casting, the car shapes are approximated with cuboids losing much detail. We leverage the simulator z-buffer to estimate fine-granular depth and then sparsify the signal to simulate LiDAR scans, thus closely matching the real domain. The main components of our method. We use a simulator with a multicamera set-up and real LiDAR binary masks to transform the synthetic dense depth map into a noisy and sparse depth map. We train in a two-step manner: the green blocks are used in the 1st and the 2nd step of training while the blue blocks in the 2nd step only. In the filtering block, green points are the reliable points S p and red points are dropped. The Image Translation network is pretrained using a CycleGAN approach [1] 3 Method Our method, shown in Figure 2, consists of two main components that include an adaptation of the synthetic data to make it similar to the real data, as well as a retrieval of reliable supervision from the real but noisy LiDAR signal.

Data Generation via Projections
Supervised depth completion methods strongly rely on the sparse depth input, achieving good performance without RGB information [3,5]. To train a completion model from synthetic data that works well in the real domain we need to generate a synthetic sparse input that reflects the real domain distribution. Instead of simulating a LiDAR via ray-casting, which is computationally expensive and hard to implement [47], we leverage the z-buffer of our synthetic rendering engine to provide a dense depth ground truth at first. Previous approaches used synthetic sparse data to evaluate a model in indoor scenes or synthetic outdoor scenes [6,35,48]. To sparsify the data a Bernoulli distribution per pixel is used in some works [2,35,48] which, given a probability p B and a dense depth image x D , samples each of the pixels x D,k by either keeping the value x D,k with probability p B or setting its value to 0 with probability (1 − p B ), thus generating the sparse depth x s B D . We argue that using x s B D does not simulate well a real LiDAR input, thus a model trained with x s B D does not perform well in the real domain. Our results in Section 4 support this observation. There are two reasons for the drop of performance in the real LiDAR data. Firstly, the distribution of the points x s B D does not follow the LiDAR sparse distribution. Secondly, there is no noise in the sampled points, as we directly sample from the ground truth. We now propose an approach to address these two issues. from the real domain similarly to [43]. We use x s R,D to generate a binary mask We then apply the masks to the dense synthetic depth data by This approach adapts the synthetic data directly to the sparsity level in the real domain without the need to tune it depending on the LiDAR used. Generating Projection Artifacts. Previous works use noise-free sparse data to pre-train [4] or evaluate a model [35] with synthetic data. However, simulating the noise of real sparse data can reduce the domain gap and improve the adaptation result. Real LiDAR depth contains noise from several sources including the asynchronous acquisition due to the rotation of lasers, dropping of points due to low surface reflectance and projection errors. Simulating a LiDAR sampling process by modelling all of these noise sources can be costly and technically difficult as a physics-based rendering engine with additional material properties is necessary to simulate the photon reflections individually. We propose a more pragmatic solution and use the z-buffer of a simulator by assuming that the dominating noise is a consequence of the point cloud projection to the RGB camera reference frame. For such a simulation, the error becomes twofold. Firstly, the 3D points are not exactly projected on the pixel center which produces a minor quantization error. Secondly, as we are projecting a sparse point cloud arising from another viewpoint, we do not have a way to filter the overlapping points by depth. This creates the see-through patterns that do not respect occlusions as shown in Figure 3 which is also observed in the real domain [49]. Therefore, a simple point drawing from a depth map at the RGB reference cannot recreate this effect and such method does not perform well in the real domain.
To recreate this pattern, we use the CARLA simulator [12], which allows us to capture multicamera synchronized synthetic data. Our CARLA set-up mimics the camera distances in KITTI [40], as our benchmark is the KITTI depth completion dataset [8]. Instead of a LiDAR, we use a virtual depth camera.
The set-up is illustrated in Figure 1. As the data is synthetic, the intrinsic and extrinsic parameters needed for the projections are known. After obtaining the depth from the virtual LiDAR camera, we sparsify it using the LiDAR masks resulting in x s M D , which is then projected onto the RGB reference with where K L , K RGB are the LiDAR and RGB camera intrinsics and P L RGB is the rigid transformation between the LiDAR and RGB reference frame. The resulting x s P D is the projected sparse input to either left or right camera.

RGB Adaptation
Similarly to domain adaptation for depth estimation methods [41,42,50], we address the domain gap in the RGB modality with style translation from synthetic to real images. Due to the added computational complexity of adapting highresolution images, we first train a model to translate from synthetic to real using a CycleGAN [1] approach. The generator is not further trained and is used to translate the synthetic images to the style of real images, thus reducing the domain gap as shown in Figure 3.

Filtering Projection Artifacts for Supervision
In a depth completion setting, the sparse depth input can be used as supervision data, similarly to [5]. However, the approach from [5] did not take into account the noise present in the data. The given LiDAR input is precise in most points with an error of only a few centimeters. However, due to the noise present, some points cannot be used for supervision, such as the see-through points, which have errors in the order of meters. Another method [49] also used the sparse input as guidance for LiDAR-stereo fusion while filtering the noisy points using stereo information. We propose to filter the noisy input without using additional sensor data such as a stereo pair as this may not always be available.
Our goal is to find a set of reliable sparse points S p , likely to be correct, for supervision based on the assumption used in Section 3.1, i.e., the main source of error are the see-through points after projection. We assume that in any given local window there are two modes of depth distribution, approximated by a closer and a further plane. We show an overview of the idea in Figure 2. The points from the closer plane are more likely to be correct as part of the occluding objects. To retrieve S p we apply a minimum pooling with window size w p yielding a minimum depth value d m per window. Then, we include in S p the points s ∈ [d m , d m + θ] where θ is a local thickness parameter of an object. The number of noisy points not filtered out depends on the window w p and object thickness θ, e.g. larger windows remove more points but the remaining points are more reliable. We use the noise rate η, which is the fraction of noisy points as introduced in noisy labels literature [51][52][53], to select w p and θ in the synthetic validation set, thus not requiring any ground truth in the real domain.
After the filtering step, a certain number of false positives remains. The noisy points in S p are more likely to be further away from the dense depth predictionŷ, hence BerHu will give more weight to those outliers. To provide extra robustness against these false positives we use in the real domain a Mean Absolute Error (MAE) loss instead of the Reverse Huber (BerHu) loss used in the synthetic domain, as MAE weights all values equally, showing more robustness to the noise.

Summary of Losses
Our proposed loss is where L S is the loss used for the synthetic data, L R the loss used for the real data and λ S and λ R are hyperparameters. We use a two-step training approach similar to past domain adaptation works using pseudo-labels [54,55], aiming first for good performance in the synthetic data before introducing noise in the labels. First, we set λ S = 1.0 and λ R = 0.0, to train only from the synthetic data. For L S we use a Reverse Huber loss, which works well for depth estimation problems [21]. Hence, we define L S as where b S is the synthetic batch size, n i the number of ground truth points in image i,ŷ is the predicted dense depth, y is the ground truth depth and L bh is the Reverse Huber loss [56].
In the second step we set λ S = 1.0 and λ R = 1.0 as we introduce real domain data into the training process using S p for supervision. We define L R as where b R is the real domain batch size and #(S p,i ) is the cardinality of the set of reliable points S p for an image i.

Experiments
We use PyTorch 1.3.1 [57] and an NVIDIA 1080 Ti GPU as well as the official implementation of FusionNet [3] as our sparse-to-dense architecture. The batch size is set to 4 and we use Adam [58] with a learning rate of 0.001. For the synthetic data, we train using x s P D by randomly projecting to the left or right camera with the same probability. In the first step of training, we use only synthetic data (i.e., λ S = 1.0, λ R = 0.0, b S = 4 and b R = 0) until performance plateaus in the synthetic validation set. In the second step, we mix real and synthetic images setting λ S = 1.0, λ R = 1.0, b S = 2, b R = 2, w p = 16 pixels and θ = 0.5 m, and train for 40,000 iterations. To test our approach, data from a real LiDAR+RGB set-up is needed as we address the artifacts arising from projecting the LiDAR to the RGB camera. There are no standard real LiDAR+RGB indoor depth completion datasets available. In NYUv2 [59] the dense ground-truth is synthetically sparsified using Bernoulli sampling, while VOID [9] provides sparse depth from visual inertial odometry that contains no projection artifacts. Thus, the KITTI depth completion benchmark [8] is our real domain dataset, as it provides paired real noisy LiDAR depth with RGB images, along with denser depth ground truth for testing. We evaluate our method in the selected validation set and test set, each containing 1,000 images. Following [3], we train using images of 1216x256 by cropping their top part. We evaluate on the full resolution images of 1216x356. The metrics used are Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE), reported in mm, and inverse RMSE (iRMSE) and inverse MAE (iMAE), in 1/km. Synthetic Data. We employ CARLA 0.84 [12] to generate synthetic data using the camera set-up in Figure 1. We collect images from 154 episodes resulting in 18,022 multicamera images for training and 3,800 for validation. An episode is defined as an expert agent placed at random in the map and driving around while collecting left and right depth+RGB images, as well as the virtual LiDAR depth. We use for the virtual LiDAR camera a regular dense depth camera instead of the provided LiDAR sensor in CARLA because the objects in the LiDAR view are simplified (e.g., CARLA approximates the cars using cuboids). The resolution of the images is 1392x1392 with a Field Of View of 90 • . To match the view and image resolution in KITTI, we first crop the center 1216x356 of the image and then the upper part of 1216x256. To adapt the synthetic RGB images, we train the original implementation of CycleGAN [1] for 180,000 iterations.

Ablation Study
We include an ablation study in Table 1 using the validation set. For the result of the whole pipeline, we average the results of three different runs to account for training variability. All of the proposed modules provide an increase in accuracy.
CARLA Adaptation. Table 1 shows that reprojecting the sparse depth is as important as matching the LiDAR sampling pattern, decreasing the RMSE by 32.4% when used jointly, i.e training with x s P D instead of x s B D . Table 1 also shows that training with the reprojected x s B D results in worse performance compared to training with x s B D , showing that it is the combination of using a LiDAR distribution of points and projection to another camera which reduces the domain gap. Even though CycleGAN mostly adapts the brightness, contrast and colors of the images as shown in Figure 3, using CycleGAN adaptation to the style of the real domain further reduces the RMSE by 6.6% when training with x s P D examples. Figure 4 includes some predictions when training using x s B D , x s M D and x s P D for examples with projection artifacts, showing that training using x s P D in the synthetic images is crucial to deal with the noisy input in the real domain. Introducing Real Domain Data. Introducing the reliable points S p as supervision in the real domain alongside the MAE loss function increases the performance as Table 1 shows. If we use BerHu along with S p supervision, the method deteriorates as the noisy points are likely to dominate the loss even if the noise rate η in S p is low. Using MAE without filtering also drops the performance, in this case due to the high noise rate η. These results show that using the noisy LiDAR points for supervision as in [5,9] is detrimental to the depth completion performance. To understand the noise level in the data, we define a point to be noisy if it is at least 0.3 meters away from the ground truth. Using this metric, the noise rate η for the points with available ground truth is 5.8%, with our filtering method decreasing η to 1.7% while dropping 45.8% of input points. The results suggest that the value of η in S p is more important than the total amount of points used for supervision. Our method shows a clear improvement compared to both synthetic baselines (Syn. Baseline) defined in Table 1. Impact of RGB Modality. Contrary to self-supervised methods, which use RGB information to compute a photometric loss, we do not require the RGB image for good performance as shown in Table 2. Including RGB information reduces the error by 0.7% in RMSE, and by using the CycleGAN RGB images the RMSE is reduced by 2.1%. In a fully supervised manner the difference is 16.3% for FusionNet [3], showing that methods aiming to further reduce the RGB domain gap may increase the overall performance. Due to computational constraints, we train the CycleGAN model in a separate step. To test an end-to-end approach, we use the method in [42], which does not use cycle-consistency, however we obtained lower-quality translated images and reduced accuracy as shown in Table 2.

Method Evaluation
Comparison to State-of-the-Art. In Table 3 we compare our method, Ours, with the real domain GT-free state-of-the-art. In the test set 4 our method decreases the RMSE by 6.4%, the MAE by 6.3% and obtains better results for iRMSE and iMAE compared to VOICED [9]. Note that these improvements upon previous methods are obtained by using an architecture with fewer parameters. Table 1 and Table 3 show that we achieve similar results to [9] by training only with synthetic data, i.e., in the first training step, which validates the observation that the main source of error to simulate are the see-through points. DDP [28] uses synthetic ground truth from Virtual KITTI [11] for training, however no  adaptation is performed on the synthetic data, resulting in worse results compared to our method even when using stereo pairs (DDP+Stereo). Both VOICED [9] and SS-S2D [5] use, besides video self-supervision, the noisy sparse input as supervision with no filtering, reducing the achievable performance as shown in Table 1 in No Filter. Non-learning methods in Table 3 perform worse than [5,6,9]. Domain Adaptation Baselines. Following synthetic-to-real depth estimation methods [41,42], we use as a domain adaptation baseline a CycleGAN [1] to adapt the images. To sparsify the synthetic depth, we use the real LiDAR masks [43], shown in Table 1 to perform better than Bernoulli sampling. The performance of this domain adaptation baseline is presented in Table 3 in DA Base. We explore the use of adversarial approaches to match synthetic and real distributions on top of the DA Base. DA Base + Disc. Out. in Table 3 uses an output discriminator using the architecture in [62], with an adversarial loss weight of 0.001 similarly to [63]. Following [42], we also tested a feature discriminator in the model bottleneck in DA Base + Disc. Feat. with weight 0.01. Table 3 shows that the use of discriminators has a small performance impact and that standard domain adaptation pipelines are not capable of bridging the domain gap. Semi-Supervised Learning. In some settings, a subset of the real data may be annotated. Our full pipeline mimics the noise in the real sparse depth and takes advantage of the unannotated data by using the filtered sparse depth S p for supervision. This provides a good initialization for further finetuning with any available annotations as Table 4 shows. Compared to pretraining using the DA Baseline, our method achieves in all cases a better performance after finetuning.  Hyper-Parameter Selection. We do not tune λ S and λ R as we assume we do not have real-domain annotated data. The projected points x s P D in the synthetic validation set are used to choose w p and θ by employing η in S p as the indicator for the filtering process performance. We define a point to be noisy if the distance to the GT is bigger than 0.3 m, but other values tested yield similar results for the optimal w p and θ. Figure 5 shows the noise percentage depending on w p and θ, where we see that curves for η follow a similar pattern in both the synthetic and real domain. We first select w p and then θ as the gain in performance is lower for θ. The optimal values found are w p = 16 pixels and θ = 0.5 m. Figure 5 also shows the MAE depending on the number of iterations in the second step, where we empirically found that any chosen value yields better performance compared to using only synthetic data. After 40,000 training iterations, we did not see a noticeable improvement. Adding Self-Supervision. When real domain video data is available, our approach can be combined with self-supervised methods [5,9]. Ours+SS-S2D in Table 3 adds the photometric loss λ ph L ph from [5] to our pipeline during the second step of training for the real data, with λ ph = 10 to have similar loss values Fig. 7: Failure cases of our method in KITTI, which cannot correct all types of noise. The left side example shows a set of noisy inputs on the wall. The right side example shows dropping of points due to low-reflectance black surfaces as L S + L R . Ours+SS-S2D further reduces the error in the test set and achieves, compared to VOICED [9], a lower RMSE by 9.2% and a lower MAE by 10.4%. Model Agnosticism. We chose FusionNet [3] as our main architecture, but we test our approach with the 18-layers architecture from [5] to show our method is robust to changes of architecture. Due to memory constraints we use the 18-layers architecture instead of the 34-layers model from [5], which accounts for the different parameter count in Table 3 between Ours w/S2D arch and SS-S2D. We set the batch size to 2, increase the number of iterations in the second step to 90,000 (the last 20,000 iterations use a lower learning rate of 10 −4 ), and freeze the batch normalization statistics in the second step. The result is given in Table 3 in Ours w/S2D arch., which achieves state-of-the-art RMSE and MAE. Qualitative Results in PandaSet [64] are shown in Figure 6 for our full method compared to the DA Baseline (DA Base in Table 3) trained for CARLA and KITTI without further tuning. PandaSet contains a different camera set-up with different physical distances compared to the one used in training, e.g., top row in Figure 6 corresponds to a back camera not present in KITTI. Our method is still capable of better correcting projection artifacts (top row and middle row) and completing the missing data (bottom row) compared to the DA Baseline. Limitations. While we addressed see-through artifacts, other types of noise can be present in the real sparse depth as Figure 7 shows. The left side example shows a set of noisy inputs on the wall that is not corrected. The right side example shows missing points in the prediction due to the lack of data in the black hood surface. The fully supervised model deals properly with these cases, suggesting that approaches focused on other types of noise could further decrease the error.

Conclusions
We proposed a domain adaptation method for sparse depth completion using data-driven masking and projections to imitate real noisy and sparse depth in synthetic data. The main source of noise in a joint RGB + LiDAR set-up was assumed to be the see-through artifacts due to projection from the LiDAR to the RGB reference frame. We also found a set of reliable points in the real data that are used for additional supervision, which helped to reduce the domain gap and to improve the performance of our model. A promising direction is to investigate the use of orthogonal domain adaptation techniques capable of leveraging the RGB inputs even more to correct also other types of error in the LiDAR co-modality.


extfieldsofstudy:
Computer Science
Engineering

text:
Hierarchical Sampling based Particle Filter for Visual-inertial Gimbal in the Wild

The gimbal platform has been widely used in photogrammetry and robot perceptual module to stabilize the camera pose, thereby improving the captured video quality. Usually a gimbal is mainly composed of sensors and actuator parts. The orientation measurements from sensor can be inputted directly to actuator to steer camera towards proper pose. But the off-the-shelf custom product is either quite expensive, or depending on highly precise IMU and Brushless DC motor with hall sensor to estimate angles, which is prone to suffer from accumulative drift over long-term operation. In this paper, a CV based new tracking and fusion algorithm dedicated for gimbal system on drones operating in nature is proposed, main contributions are listed as below: a) a light-weight Resnet -18 backbone based network model was trained from scratch, and deployed onto Jetson Nano platform to segment the image into binary parts (ground and sky). b) geometric primitives tracking of the skyline and ground plane in 3D as cues, along with orientation estimation from IMU can provide multiple guesses for orientation. c) spherical surface based adaptive particle sampling can fuse orientation from aforementioned sensor sources efficiently. The final prototyping algorithm is tested on the real-time embedded system, and with both simulation on ground and real functional tests in the air.

II. RELATED WORKS
The stabilization of video problem can be traced back to the Electronic Image Stabilization (EIS) technique [29], which relies on the handcrafted feature points in consecutive image frames, to search for correspondences for image alignment. This renders the cheap hand-held device possibility to stabilize video sequence without too much overhead in hardware or mechanics. The feature points based tracking, [26], or optical flow tracking [31], requires always good association of feature points or raw pixel intensities across video frames [22]. thus, it tends to be brittle to observation noise induced by motion blur. For that reason, humanoid robot vision system [34], or complex video [35] task compensates the rotation and shift of video images by using costly IMU only.
*This work was supported by Flemish government "Vlir-UOS" project under grant VLS-D7314-EC2020SIN278A101. 1  The review work about traditional motion modelling work applied to video stabilization, can be referred to in [30]. But these passive video stabilization work is always constrained on image frame, cannot fully utilizes more freedoms in 3D. The video stabilization problem can also be considered as image deblurring task. The main goal is to infer the pixel intensity in blurry area from transferring or interpolation of neighbouring pixels information on single or multiple images, incorporating both temporal or spatial constraints from the video. The recent deep learning technique boosts video processing performance by introducing prior from dataset at specific domains. [32] warps the selfie image patch locally by decoding the scene representation from multiple encoders, such as foreground contours, feature points on the background, even 3D facial mesh. [33] trains a sensor-driver model by using raw IMU data, to predict the visual odometry in a sliding time window, similar to the filter approach. [36] integrated the video stabilization into large mapping and localization framework via graph in self-supervised manner. To tackle the problem of video stabilization in dynamic scenes, [38] learns a dense warping fields from consecutive video frames, then the warped parts are blended to synthesize the stabilized image. [37] framework even makes use of the motion prediction model, a spatial filter for optical flow tracking, followed by inpainting process to synthesize outputs. All these models involve a lot of image processing overhead, and are quite heavy to deploy on edge device for real-time use.
Video stabilization is also widely applied in aerial surveillance [39]. By Using natural cue to estimate the attitude of UAVs or fixed-wing airplane is literally nothing new, like the work in [2] [6][1] relying all on skyline tracking. But it requires delicate tuning tricks to search for boundary, and input image with ideal condition. Other work depends on the tracking of geometric primitives on 2D image or in 3D world, either by handcrafted feature points, like [23][25] [26], or curved boundary [5], to predict the ego-pose of camera, which can further guide the stabilisation of attitude of UAVs, as work done in [3], [40] even detects the moving objects in the view through tracking filter over the "SIFT" [7] feature points. However, in natural setting, there are a lot of spurious and similar appearance features, all these tracking algorithm relying on traditional feature points cannot work at all in such setting.

III. OVERVIEW
The ultimate goal of our project is to deploy a gimbal platform carried by UAVs for volcanic eruption surveillance in a depopulated region full of mountains, thence, our work differs from prior work in that, a quite robust platform in operation over long-term and with endurance of harsh weather condition should be devised. As shown in Figure 1, the gimbal system is placed underneath airplane outfitted with camera, IMU, and barometer. To overcome the aforementioned correspondence issues of feature points, the gimbal system can take advantage of the natural cues such as skyline and ground for pose tracking. At bottom of Figure 1, there is a binary mask with skyline and points in ground region. The ray direction passing through the red dots position along with height from barometer, can determine the 3D point through trigonometry. More details about perception can be referred to at section IV.

A. Open source Hardware
The gimbal platform design is formed on open-sourced hardware completely. The main processing unit is a Jetson Nano, featured with 2GB GPU memory and Quad-core ARM A57, connected to IMU, camera, and barometer sensors. OpenCR driver board maps the driving command to control commands for two servo motors, in order to compensate roll and pitch angle accordingly.

B. ROS nodes
The software is composed mainly of three parts, The preprocessing part including network model and geometric primitive extraction, followed by tracking module to align the skyline and normal of ground plane in current frame with those at reference frame. The compensation angles from various pipelines are then fed into the proposed particle filter in section V to have fusion output, finally used as input for controller to compensate motion.

IV. PERCEPTION
The perception part is structured into preprocessing and rotation transform estimation, where he the perception part can be further separated into two pipelines, roll and pitch prediction from skyline tracking, and full Euler 3D rotation from ground plane tracking. The following clarification follows this structure.

A. Preprocessing
The general computer vision processing pipeline, like "OpenCV" library is investigated for our purpose at first. a) The feature points are extracted from two image frames to estimate relative pose, b) The boundary between sky and ground is isolated from image by Canny edge detection for alignment. As illustrated in Figure 4, regarding the feature points, the always alike appearance in grassland region and cloudy sky with dynamics all pose a challenge to correspondence match. In Figure 4a, many false inconsistent correspondences are found. As for boundary isolation, the color distribution of sky, clouds, and grass land are significantly contrasting, on top of the HSV image in Figure 4b, the Canny [8] edge detection can be applied to find the boundary between ground and sky region, but some brightness of bluish sky is cast onto grass ground. As a consequence, the final fitting curve are shifted downwards compared to reality as depicted in 4c. The Canny edge detection on raw RGB is quite noisy over whole image, unsuitable for boundary search. To tackle this challenging issue, we referred to the neural network model, and trained from scratch on Skyline dataset [4] by using Resnet-18 [19] as backbone for binary mask segmentation, the model is fine-tuned on self-collected dataset in nature with a few hundred images. Some train data samples along with ground truth masks are shared in Figure  5. The model is exported into "ONNX" and optimized by "TensorRT" to convert to "FP16" precision.

B. Estimation from Skyline Tracking
On top of binary mask results from network model, the skyline can be found along the pixel gradient change direction. The extracted skyline can be further considered as a cue to estimate the roll and pitch of camera with respect to the line in reference frame, as shown in Figure 6. The search of boundary points, as yellow dots in Figure 6c, can be implemented firstly over the whole reference image. Straight line fitting in red as illustrated in Figure 6c, including two parameters, slope m and intercept b in Equation 1. For the subsequent images, we use the constant angular velocity model from skyline tracking results of last two frames, to predict the skyline line position in current frame, like the green line presented in 6d. The assumption is that image never flips upside down, because of airplane maneuverability constraints. The boundary can always be searched along the vertical direction on predicted straight line pixels, for example, the pixel labelled as sky should be searched in downwards direction. In practice, the found points at boundary is further down-sampled to a few hundreds of points to speed up the search. Then the skyline parameters m and k in Equation 2 are calculated by least square algorithm from the sample points. The estimated skyline of current frame is the red line in Figure 6d. To guarantee the skyline has enough length for down-streaming processing, a minimal half length of image width is set as check condition. The roll α can be derived from arctan(m), and pitch β can be derived from subtraction of arctan((h1 − c)/ f y ) − arctan((h2 − c)/( f y )). c is the half size of image height, h1 and h2 is the height of center points at straight line in reference and current image frame. f y is the focus length of camera. To avoid unnecessary operation, roll and pitch have a certain tolerant range, so that processing is only triggered when the movement is out of this range. For the coupling issue of roll and pitch processing, the logic as below is utilized. There are three cases in total. From slope rate m of straight skyline, roll angle compensation can be recovered first. On top of corrected results by roll angle, pitch movement checks further by skyline height shift. In such way, pure roll, pure pitch, or both happening at the same time, are covered totally.

C. Estimation from Ground Plane Tracking
Ground plane tracking, is relying on the normal vector of ground plane, as demonstrated in Figure 1, a set of points in the ground region of binary mask are sampled evenly, followed by a back-projection to 3D camera frame corresponding to Equations 3-9. The height h is measured from barometer, and ray direction passing through the pixel position is P i at left of Equation 5. K is the intrinsic matrix obtained by calibration approach [15]. θ is derived from the dot product of gravitational vector and ray direction. Length scale l i can be obtained by trigonometry in Equation 8. Finally the current normal vector m of ground plane is shaped by cross product of non-Collinear triple points. The ground plane tracing mode is only triggered when the level height of camera is over a 300 meters level in our setting, so that the variance of uneven grass land can be approximated by a flat plane compared to height. The ground region also uses threshold condition check.
cos(θ ) = P i · N G || P i || · || N G || (7) Next, the rotation matrix to align the normal vector, m in current frame and reference normal n at start can be derived following Equations 11-14.
In the end, the 3D Euler angles are retrieved from the Rotation matrix according to  in an order of "yaw ¡-pitch ¡-roll". Only roll and pitch are used for fusion.
V. HIERARCHICAL TEMPLATE BASED PARTICLE FILTER The particle filter is easy to implement, applicable in non-linear problems, in particular for positioning like [20]. Finally, a variant of vanilla particle filter [28], sampling on the spherical surface is devised and presented in pseudo code below. In real configuration, the roll and pitch is virtually constrained by a rotation range because of mechanical kinematics, e.g., roll is -45 • to 45 • .   The general idea behind this adaptive particle filter is straightforward. It comprises of three steps, sampling from orientation measurements of IMU s I (Line 3), sampling according to observation from CV pipeline (Line 9), and resampling proportional to updated weight of each particle(Line 20). This weight is the square root distance of angle between one particle and CV observation. s C 1 , and s C 2 indicate the orientation from skyline and ground plane separately.

Algorithm 2: Particle Filter based on Manifold Template
The manifold surface is discretized into cells at multiscale resolution, where the longitudinal direction is pitch while the latitudinal direction is roll. When both observation are available, more samples will be created around those close to CV observation. In Line 3, particles are sampled from normal Gaussian distribution, with mean centered at IMU measurements plus shift by constant angular velocity propagating through interval and certain bias σ b .
Line 10 of Table 21 manifests the metric for neighbouring check of particles close to CV observation. µ f at Line 9 is coming from fusion Equation 20, a sum of mean of IMU and CV wieghted by inverse of variance.
IMU and CV variance are set according to our practical tests. µ C and δ C are merged results from two CV pipelines in same form of Equation 20 and 21. The fusion of three pipelines are recursive fusion by two times. Each particle life cycle can be represented into a four states' machine, in Figure 7, arrows between them stand for transition condition. The rasterized cells spread over the discretized manifold space of roll and pitch angles. In our setting, the CV approach has smaller variance compared to IMU according to accuracy. Three levels of resolution (Ω 1,2,3 in Table 21) are utilized. The top left is initial state, sampled from normal distribution centered at IMU measurement at time t. Bottom left is prediction based on propagation of previous state by extrapolation in time. When the orientation from single CV pipeline, either skyline or ground plane is available, the new samples in red are generated on finer resolution neighbouring the green dots corresponding to µ f in Equation 20. Each dot stands for the center of cell. If both CV pipelines are at hand, the highest resolution is employed for new samples generation. Each particle's timestamp of creation is recorded. If elapsed time is over certain interval δt, the particles will be moved back to coarser resolution, like the arrow direction from bottom right to top right. The particle set is a mixing states of this machine, each particle with different resolution over time.

VI. EXPERIMENT & RESULTS
In practice, the video is scaled down to 640×480 resolution, to achieve around 20 FPS real-time performance. Intrinsics of camera follow the calibration guide of [15]. Extrinsic calibration between IMU and camera is estbalished on an open source tool "Kalibr" [13] [14]. In Figure 8, the simulation setup at lab can be referred to. There are three parts, motor driver board, Jetson Nano, and gimbal part, the former two will be positioned inside airplane body for outdoor test, only gimbal part outfitted with sensors is exposed in the air. All 3D printed pieces have enhanced connections, also having taken the aerial dynamics into account for flying efficiency. A big natural picture from real scenario is chosen for test. The camera on gimbal is placed to be forward facing the wall. Then gimbal cases along with sensors, are attached onto a pole end, which is not in the view of picture 8, the other end is controlled by hand to simulate the random rotation and vibration. The ground truth data of roll and pitch angles are read from the fixed protractor near the each servo motor accordingly, like the one fixed behind the gimbal for roll reading. In demo video test, we use the fused estimation from our particle filter to compensate the motor rotation conversely for camera stabilization. The closed-loop PID control is leveraged for actuation. Motor is with 50-80ms latency. All the following test sequences are starting from horizontal direction on image. This configuration holds for all tests of our algorithms, we always start from the aligned ground plane orthogonal to gravitational vector, and horizontal skyline. In outdoor test, we always manually align the gimbal platform with proper pose. For comparison, the SOTA visual-inertial frameworks, "ORBSLAM3", "VINS-mono" and "MSCKF" are selected as benchmarks. The open source datasets for SLAM research were captured indoor or within the city or urban region, rarely including the depopulated wild zone, we thus recorded several custom sequences by using the aforementioned setup at lab. The outdoor test can be a bit tough to generate ground truth angle, thus we only performed the video test on real airplane outdoor for verification only. Disturbance to gimbal is simulated by hand-held pole as mentioned before. During recording, each pose configuration of gimbal is kept for a few seconds on a par with protractor readings as ground truth. In the end, the recorded rosbag of raw video and IMU sequences are utilized to run the SOTA SLAM approaches. All results of 6D pose estimation are exported, and only pitch and roll angles are taken out for comparison. There is no displacement estimation in our approach, so all postend loop closure detection and post-optimization of SLAM cannot be triggered.  It can be inferred from Table I, our approach outperforms the other algorithms consistently by a considerable margin, nearly across all the test sequence. The sequences are named after the simulation movement, sequence with pure roll or pitch rotation and both axis rotations. The roll test has bigger movement range, so we recorded more dense poses in Rosbag. The fairly good performance with least RMSE in most tests should be attributed to the good assumption, skyline and ground plane in the wild, unlike the traditional feature point based SLAM, which is struggling or failing to establish the correspondences between video frames, due to wrong match in feature points. Furthermore, our particle based filter even shows quite robust performance over long-term test, compared to other approaches on "full*" sequences. the particle filter models the non-linearity of this problem rather well, fusing the ego-body IMU measurement and CV observation r obustly to get rid of accumulated drift. "full-03" sequence lasts for more than 20 minutes, it is remarkable that drift error on other approaches in last row of Table I has one order of magnitude increase compared to other rows. Figure 9 further justifies our claim. The test keeps the same setting above by repeating 5 times per each sequence, our approach has mean value at 4.30 degrees, and variance at 0.16, which is smallest amongst the all. Thanks to the natural cue assumption we employed, and adaptive particle filter fusion. The approach can always achieve best accuracy and robustness in any configuration, even invariant over the long period operation. It is noteworthy that the ORBs [9] performs worst, which may be explained by the front-end ORB feature extraction, it looks like the optical flow is more robust compared to handcrafted ORB feature points in such scenario.

VII. CONCLUSION & OUTLOOK
A brand new gimbal system is proposed in this paper, which includes natural geometry primitives tracking at frondend, by using skyline and ground plane. The current frame rotation with respect to the manually set reference frame can be derived from the two lines and normal vectors alignment. Next, a specific particle filter with adaptive resolution based sampling strategy can fuse orientation from both CV and IMU pipelines, even jump between different states according to various lifetime phase of one particle. The final experimental results are implemented on the 3D printed gimbal platform. All simulation tests at lab are performed in realtime on a Jetson Nano. To evaluate the performance of our system, several self-collected sequences with static picture as target were recorded. The rotation disturbance to gimbal is simulated by turning a pole at the other end by hand. Three types of popular SOTA visual-inertial based SLAM methods are chosen to compare, revealing the best accuracy and robustness of our approach amongst all comparisons. Our method can even get rid of accumulative drift issue in traditional visual odometry methods, thanks to a global reference frame in use for all frames' tracking. limitations, Our approach depends on the geometric primitive assumption, for challenging wild, the flatness assumption of ground cannot hold, even the abrupt illumination change and low-level cloud occlusion, all pose great challenge to tracking of those visual cues. The skyline now is simplified into a straight line for easy match and real-time performance, on affordable edge device. In reality, the curved mountains in the view make this assumption quite brittle. The traditional iterative closest point (ICP) can be further applied on image level for match. The 360 fisheye camera with full view should be beneficial to robustness of gimbal platform. Even the gimbal motor and processing unit can be upgraded to fulfill more real-time potential for more high-level tasks.


extfieldsofstudy:
Computer Science
Mathematics

text:
An Efﬁcient Piggybacking Design Framework with Sub-packetization l (cid:54) r for All-Node Repair

—Piggybacking design has been widely applied in distributed storage systems since it can greatly reduce the repair bandwidth with small sub-packetization. Compared with other existing erasure codes, piggybacking is more convenient to operate and the I/O cost is lower. In this paper, we propose a new efﬁcient design which can further reduce the repair bandwidth with the sub-packetization l (cid:54) r where r = n − k . Generally, we let l (cid:54) 8 . Compared with other analogous designs, our design has lower l and the value of l is more ﬂexible. Moreover, our design can repair all nodes with small repair bandwidth. Therefore our piggybacking design is more feasible.

I. INTRODUCTION
I N the era of big data, countless amounts of data are generated every day. Therefore distributed storage systems are widely deployed in data centers. However the node failures occasionally happen. Erasure codes are widely used to retain the stability of the systems and MDS code is the most popular. A (k + r, k) MDS storage code, composed of (k + r) nodes, can tolerate the failure of any r nodes. If a node fails, it can be recovered by downloading data from some surviving nodes. The total data need to be downloaded to recover the failed node is termed as repair bandwidth. A trivial method is downloading the total data of any k surviving nodes, but it costs a lot. In 2010, Dimakis et al. [1] proposed the lower bound of the repair bandwidth of MDS codes and introduced a notion of regenerating codes that the failed nodes can be recovered by downloading same amount of data from each of the surviving nodes. Particularly, the ones with the optimal repair bandeidth are termed as MSR codes. MSR codes were studied extensively in [2]- [6].
The number of the sub-packets stored in each node that is termed as the sub-packetization is also an important indicator for regenerating codes. It has been proved that exponential sub-packetization is necessary for MSR codes [7]. Hence MSR codes are extremely difficult to apply into distributed storage systems. In 2013, Rashmi et al. [8], [9] proposed the piggybacking framework based on MDS array codes to reduce the repair bandwidth with samll sub-packetization. The failed nodes can be recovered by solving piggybacking functions which are downloaded from surviving nodes. And the piggybacked codes retain the MDS property. Since then, the piggybacking design has been further developed. Yang et al. [10] utilized the piggybacking design to obtian new MSR codes which are almost optimal repair bandwidth of parity nodes while retaining the optimal repair bandwidth of systematic nodes. Yuan et al. [11] and Shangguan et al. [12] proposed a generalized piggybacking codes (REPB&SPB) with the average repair bandwidth ratio 2 √ r+1 and √ 2r−1 r of systematic nodes respectively. And the sub-packetization l = r in [12]. For repairing all nodes, an efficient one-to-one piggybacking design (OOP) [13] was proposed by Li et al while the sub-packetization l = r − 1 + √ r − 1. Sun et al. [14] reduced the repair bandwidth of parity nodes in [12] by making full use of each parity node while retaining the repair bandwidth of systematic nodes. Jiang et al. [15] reduced the average repair bandwidth of all nodes by sacrificing a part of the systematic nodes repair bandwidth.
In this paper, we propose a new effcient piggybacking design (C). Compared with the designs in [12]- [14], our design has lower sub-packetization l r and can further reduces the repair bandwidth of systematic nodes in [13], [14] when l = r. Compared with the codes in [11], [15], the average repair bandwidth of all nodes and systematic nodes can be further reduced under the same parameters in this paper.
The remaining of the paper is organized as follows. Section II describes the framework of the new piggybacking and gives a toy example to illustrate the repair process of the design. The characteristics of the design are discussed in Section III. Section IV compares our design with some existing ones. The conclusion is given in section v.

II. THE NEW PIGGBACKING DESIGN
In this section, we propose a new explicit piggybacking design and show the repair process of our codes through an example. For any two integers i < j, denote by [i, j] = {i, i + 1, . . . , j} and [i] = {1, 2, . . . , i}.

A. The new piggybacking design
For an (n, k, l) MDS array code over F q , the original data symbols a 1 , . . . , a l are stored in the k systematic nodes, where a i = (a 1,i , . . . , a k,i ) . The corresponding parity symbols are P T 1 a i , . . . , P T r a i T ∈ F r q . In the piggybacking model, each substripe a 1,i , . . . , a k,i , P T 1 a i , . . . , P T r a i T ∈ F n q is an MDS codeword.  Let the sub-packetization l = s + t + u r, and denote k = α 1 u+β. Divide the systematic nodes into u groups where each of the first u−β groups contains α 1 nodes and each of the last β groups contains α 2 nodes, where α 1 = k u , α 2 = k u . As in figure 1, we obtian 1) the piggybacking design of data symbols: • For the ith group, we treat the first l − i symbols of each node as an • Partition the symbols in G i into r − 1 almost equal parts column by column and the number of symbols in each part is Then add each part of G i to P T 2 a l−i+1 , . . . , P T r a l−i+1 by piggybacking function in turn, where l = s + t + u, i ∈ [u].
2) the piggybacking design of parity symbols: • Partition the first s + t symbols of each parity node into two parts, s and t symbols respectively. We treat the first s symbols of the r parity nodes as an r × s matrix H 1 and use permutation g j = (1, 2, . . . , r) j act on the jth and i denotes the position of the ith element of the jth column. Then we have a new matrixH 1 . • Partition the symbols of the mth arrow ofH 1 into t almost equal parts and the number of symbols in each part is required to differ by no more than 1. Then add each part to P T m a s+1 , . . . , P T m a s+t in turn by piggybacking functions, where m ∈ [r]. It is easy to check that the MDS property of the code and the symbols added to each parity symbol will not come from the same node.

III. THE GENERAL REPAIR PROPERTIES
In this section, we will discuss the average repair bandwidth of our piggybacking design and the values of s, t, u to get the optimal repair bandwidth where l = s + t + u r.
Remark 1. Although we consider the total repair bandwidth of systematic nodes in one group, the repair process is one by one, i.e. we still consider the condition of single node failures.

A. The average repair bandwidth of systematic nodes
Define three variables as in section II: 1) The repair bandwidth of systematic nodes in the first u − β groups: Define γ i as the repair bandwidth of the ith group of systematic nodes, where i ∈ [u − β]. Firstly, consider the repair bandwidth of the first α 1 nodes. Let where z 1 is the remainder of α 1 (l − 1) mod r − 1. Therefore, the first r − 1 − z 1 positions of the parity symbols used to piggyback the first group of data symbols each piggyback ζ 1 symbols, and each position piggybacks ζ 1 symbols in the last z 1 positions.
Since only the last one symbol in each node of the first group has not been piggybacked, these α 1 symbols can be recovered only through the MDS property, i.e. download kα 1 symbols totaly. And other symbols are all piggybacked in some parity symbols, they can be recovered by some linear operations. Since kα 1 symbols have been downloaded, the original parity symbols in the last substripe can be recovered. Then in this process, ζ 1 or ζ 1 need to be downloaded for each data symbol recovered. Therefore we need download ((r − 1) − z 1 )·ζ 2 1 +z 1 · ζ 1 2 symbols to recover the remaining α 1 (l − 1) data symbols and γ 1 = kα 1 + ((r − 1) − z 1 ) · ζ 2 1 + z 1 · ζ 1 2 . Define ϕ 1 = · · · = ϕ (r−1)−z1 = ζ 1 , and ϕ (r−1)−z1+1 = · · · = ϕ r−1 = ζ 1 . Then Similarly, we can calculate the total repair bandwidth of the nodes in the ith group, . Then 2) The repair bandwidth of the last β groups: Since the only difference from III-A1 is that the last β groups each contain α 2 nodes, in a similar way, we have 3) The total repair bandwidth of systematic nodes: By (3) and (4), we obtain Theorem 1. The total repair bandwidth of systematic nodes in the new piggybacking design is 4) The lower bound of the average repair bandwidth of systematic nodes: Without loss of generality, we may assume u|k, then α 1 = α 2 = k u α. If u k, only need α 2 = α 1 + 1, and the rest of the process is all the same.
Define the average repair bandwidth and repair ratio of systematic nodes asγ sys = γsys k andγ sys =γ sys kl respectively. Proposition 2.γ sys andγ sys achieve their lower bounds . And Γ is defined as in (8). (6) meets the quality when z i = 0. Therefore the average repair bandwidth and define Then combine (7) with (8),γ sys andγ sys achieve the lower bounds when u = by some simple calculation.

B. The average repair bandwidth of parity nodes
Consider the repair bandwidth of parity nodes on the premise that the average repair bandwidth of the systematic nodes has been optimal in our scheme, i.e. we may assume u = 6l 2 −6l−1 3r−1 . As in Remark 1, consider the total repair bandwidth of parity nodes first. Let ω = s t , ω = s t , y = s − tω, clearly t s. Then the t positions used to piggyback parity symbols in each parity node, each of the first t−y positions piggybacks ω parity symbols, the rest of them each piggyback ω symbols.
1) The total repair bandwidth of parity nodes: Since only the first s symbols of each parity nodes have been piggybacked, the last t + u symbols can be recovered only through the MDS property. Then in this process, (t + u)kr symbols need to be downloaded to recover the original (t + u)r parity symbols first. On the other hand, there are rs+k(s+t)+ k u ((u−1)+· · ·+1) = rs+k(s+t)+ k(u − 1) 2 symbols that have been piggybacked to these (t+u)r symbols, for simplicity, we suppose α 1 = α 2 = k u , i.e. u |k . However only the rs + ks symbols of them haven't been downloaded yet, the rest belong to the rk(t + u) symbols that have been already downloaded in the first step. Hence only the rs + ks more symbols need to be downloaded.
Consider the repair of the first s symbols in each parity nodes next. To recover one symbol in this process, ω or ω symbols need to be downloaded. As in the section A, r (t − y) ω 2 + y (ω ) 2 symbols need to be downloaded totally. Therefore we can obtian Theorem 3. The total repair bandwidth γ par of parity nodes is 2) The lower bound of the average repair bandwidth of parity nodes: Define the average repair bandwidth of parity nodes asγ par = and t|s, ω = s t , thenγ par andγ par achieve their lower bounds (t + u)k + s(k+r) r + tω 2 and t+u l + s(k+r) rkl + tω 2 kl respectively when t = (l − u) r k(r−1) .
By s + t + u = l, we have s = l − u − t, since l and u are known, (10) can be viewed as a univariate polynomial with respect to t. Sinceγ par andγ par are are both constant multiples of γ sys , they achieve the lower bounds at the same point t = (l − u) r k(r−1) .

Remark 2.
Although the two extreme points we consider above do not necessarily make the average repair bandwidth of all nodes optimal, we believe that the repair of systematic nodes is the top priority and the gap is also negligible. The most relevant works to ours are [11], [12], [14], [15]. Define the code in this paper as C and the average repair bandwidth ratio of all nodes asγ all . Clearly,γ all =γ sys k k+r + γ par r k+r . Remark 4. In Yuan et al. [11], in order to reduce the repair bandwidth of parity nodes, the original sub-packetization l needs to be enlarged by m times, where m 2. There are two piggybacking designs that we call C 1 and C 2 here in Jiang et al. [15].
To some extent, the method of [12] is utilized to improve our work, but the repair bandwidth of systematic nodes in [12] and [14] is reduced by using more locations to piggyback the data information symbols in this paper. Furthermore, the designs in [12] and [14] are only suitable for the sub-paketization l = r. The average repair bandwidth ratioγ sys of systematic nodes of our design and SPB in [12] are plotted in Fig. 3. Since the repair bandwidth of systematic nodes is identical in [12] and [14], it is not necessary to make further comparison with [14]. On the other hand, it is not feasible for the subpacketization l 10 in industry, hence it is sufficient to consider l ∈ [5,10] in Fig. 3. And the repair bandwidth ratio can be reduced by about 5%. Sinceγ sys is independent with k, for the sake of high code rate, k in taken as 50 in Fig. 3.
In the following, we focus on comparing our repair bandwidth ratio with [11] and [15] under the same parameters. On the other hand, [15] reduced the average repair bandwidth ratio of all nodes in [11] by sacrificing a part of the systematic nodes repair bandwidth, hence the average repair bandwidth ratio of all nodes is compared with [15] in Fig. 4. And Fig.  5 shows the average repair bandwidth ratio of the systematic nodes in [11] and ours. In Fig. 4 and 5, k = 100 and l = 8. The repair bandwidth ratio can be reduced by about 8% and 6% respectively.
Actually, it is inevitable that our results outperform [11] and [15] , since during the repair process at each node, the fewer  [11] with k = 80, 10 r 20 and l 8 symbols are recovered only through the MDS property in this paper.

V. CONCLUSION
In this paper, we propose a new piggybacking design framework with sub-packetization l r for all nodes. The proposed design further reduces the repair bandwidth of all nodes in [11] and [15]. Moreover, it expands the scope of application of the design in [12] and [14] while reducing the repair bandwidth of systematic nodes. Compared to other piggybacking designs, the proposed design has lower sub-packetization and repair bandwidth, therefore it is more feasible in industry.