Thanh Dat

Maurice Nguyen

Student Member, IEEE Giuseppe Quach

Senior Member, IEEE Pierre Valenzise

Life Fellow, IEEE Duhamel

Lossless Coding of Point Cloud Geometry using a Deep Generative Model

Keywords: Point Cloud Coding, Deep Learning, G-PCC, context model, arithmetic coding

This paper proposes a lossless point cloud (PC) geometry compression method that uses neural networks to estimate the probability distribution of voxel occupancy. First, to take into account the PC sparsity, our method adaptively partitions a point cloud into multiple voxel block sizes. This partitioning is signalled via an octree. Second, we employ a deep auto-regressive generative model to estimate the occupancy probability of each voxel given the previously encoded ones. We then employ the estimated probabilities to code efficiently a block using a context-based arithmetic coder. Our context has variable size and can expand beyond the current block to learn more accurate probabilities. We also consider using data augmentation techniques to increase the generalization capability of the learned probability models, in particular in the presence of noise and lower-density point clouds. Experimental evaluation, performed on a variety of point clouds from four different datasets and with diverse characteristics, demonstrates that our method reduces significantly (by up to 37%) the rate for lossless coding compared to the state-of-the-art MPEG codec.

I. INTRODUCTION

P OINT clouds (PC) are becoming the most popular data structure for many 3D applications such as augmented, mixed or virtual reality, as they enable six degrees of freedom (6DoF) interaction. Typical PCs contain millions of points, each point being represented by x, y, z coordinates, and attributes (e.g. color, normal, etc.). This entails a high transmission and storage cost. As a result, there is a massive demand for efficient Point Cloud Compression (PCC) methods to enable the practical use of this content.

The Moving Picture Expert Group (MPEG) has studied coding solution for various categories of point clouds, including static point clouds (category 1), dynamic point clouds (category 2), and LiDAR sequences (category 3 -dynamically acquired point clouds). As a result, two PCC standards have been developed [START_REF] Schwarz | Emerging MPEG Standards for Point Cloud Compression[END_REF]- [START_REF] Graziosi | An overview of ongoing point cloud compression standardization activities: video-based (V-PCC) and geometry-based (G-PCC)[END_REF]: Video-based PCC (V-PCC) and Geometry-based PCC (G-PCC). V-PCC focuses on dynamic point clouds, and projects the volumetric video onto 2D planes before encoding. The generated 2D videos are then compressed using 2D video coding standards. This approach benefits from efficient 2D video coding solutions which have D. T. Nguyen, M. Quach, G. Valenzise and P. Duhamel are with the Université Paris-Saclay, CNRS, CentraleSupelec, Laboratoire des Signaux et Systèmes (UMR 8506), 91190 Gif-sur-Yvette, France (email: thanh-dat.nguyen@centralesupelec.fr; maurice.quach@l2s.centralesupelec.fr; giuseppe.valenzise@l2s.centralesupelec.fr; pierre.duhamel@l2s.centralesupelec.fr).

been optimized over several decades. On the other hand, G-PCC targets static content, and the geometry and attribute information are independently encoded. Color attributes can be encoded using methods based on the Region Adaptive Hierarchical Transform (RAHT) [START_REF] De Queiroz | Compression of 3d point clouds using a region-adaptive hierarchical transform[END_REF], Predicting Transform or Lifting Transform [START_REF] Graziosi | An overview of ongoing point cloud compression standardization activities: video-based (V-PCC) and geometry-based (G-PCC)[END_REF]. Coding the PC geometry is particularly important to convey the 3D structure of the PC, but is also challenging, as the non-regular sampling of point clouds makes it difficult to use conventional signal processing and compression tools. In this paper, we focus on lossless coding of point cloud geometry.

In particular, we consider the case of voxelized point clouds. Voxelization is the process that quantizes the coordinates of a point cloud to integer precision prior to the coding process. This process is common in many coding scenarios, e.g., when dealing with dense point clouds such as those produced by camera arrays. After voxelization, the point cloud geometry can be represented either directly in the voxel domain or using an octree spatial decomposition. PCs are divided into a fixed number of cubes, which defines the resolution (e.g., 10 bit = 1024 cubes per dimension). Each cube is called a voxel. If a voxel contains at least one point, it is called an occupied voxel. Usually, very few voxels are occupied and a large part of the volume is empty. An octree representation can be obtained by recursively splitting the volume into eight subcubes until the desired precision is achieved. Then, occupied blocks are marked by bit 1 and empty blocks are marked by bit 0. Consequently, at each level, the generated 8 bits represent the occupancy state of an octree node (octant). Our method operates in both the voxel and octree domain. On the one hand, the octree representation can naturally adapt to the sparsity of the point cloud, as empty octants do not need to be further split; on the other hand, in the voxel domain convolutions can be naturally expressed, and geometric information (i.e., planes, surfaces, etc.) can be explicitly processed by a neural network.

In this work, we propose a deep-learning-based method (named VoxelDNN) for lossless compression of static voxelized point cloud geometry. Our main contributions are:

• We employ for the first time a deep generative model in the voxel domain to estimate the occupancy probabilities sequentially using a masked 3D convolutional network. The conditional distribution is then used to model the context of a context-based arithmetic coder. • We propose an optimal rate-driven partitioning and context selection algorithm. The partitioning algorithm adapts to the point cloud sparsity by employing a hybrid octree/voxel representation while the context to encode each block is expanded to the neighboring blocks and the expansion size is optimally selected. • We propose specific data augmentation techniques for 3D point clouds coding, to increase its generalization capability. We demonstrate experimentally that the proposed solution outperforms the state-of-the-art MPEG G-PCC lossless codec in terms of bits per occupied voxel over a set of point clouds with varying density and content type. The rest of the paper is structured as follows: Section II reviews the related work; the proposed method is described in Section III; Section IV presents the experimental results; and finally Section V concludes the paper.

II. RELATED WORK

Relevant work related to this paper includes state-of-the-art PC geometry coding and learning-based methods in image and point cloud compression.

A. MPEG G-PCC and Conventional Lossless Codecs

Most existing methods that compress point cloud geometry, including MPEG G-PCC, use octree coding [START_REF] Schnabel | Octree-based point-cloud compression[END_REF]- [START_REF] Biswas | Muscle: Multi sweep compression of lidar using deep entropy models[END_REF] and local approximations called "triangle soups" (trisoup) [START_REF] Schnabel | Octree-based point-cloud compression[END_REF], [START_REF] Dricot | Adaptive multi-level triangle soup for geometry-based point cloud coding[END_REF].

In the G-PCC geometry coder, points are first transformed and voxelized into an axis-aligned bounding box before geometry analysis using trisoup or octree scheme. In the trisoup coder, geometry can be represented by a pruned octree plus a surface model. This model approximates the surface in each leaf of the pruned octree using 1 to 10 triangles. In contrast, the octree coder partitions voxelized blocks until sub-cubes of dimension one are reached. First, the coordinates of isolated points are independently encoded to avoid "polluting" the octree coding (Direct Coding Mode -DCM) [START_REF]Inference of a mode using point location direct coding[END_REF]. To encode the occupancy pattern of each octree node, G-PCC introduces many methods to exploit local geometry information and obtain an accurate context for arithmetic coding, such as Neighbour-Dependent Entropy Context [START_REF]Neighbour-dependent entropy coding of occupancy patterns[END_REF], intra prediction [START_REF]Intra mode for geometry coding[END_REF], planar/angular coding mode [START_REF]Planar mode in octree-based geometry coding[END_REF], [START_REF]An improvement of the planar coding mode[END_REF], etc. In this paper, we compare our method against G-PCC lossless geometry coding with octree coding which also targets static point clouds.

In order to deal with the irregular point space, many octreebased lossless PCC methods have been proposed. In [START_REF] Schnabel | Octree-based point-cloud compression[END_REF], the authors proposed an octree-based method which aims at reducing entropy by employing prediction techniques based on local surface approximations to predict occupancy patterns. Recently, more context modeling based approaches are proposed [START_REF] Garcia | Context-based octree coding for point-cloud video[END_REF]- [START_REF] Garcia | Geometry Coding for Dynamic Voxelized Point Clouds Using Octrees and Multiple Contexts[END_REF]. For example, the intra-frame compression method P(PNI) proposed in [START_REF] Garcia | Geometry Coding for Dynamic Voxelized Point Clouds Using Octrees and Multiple Contexts[END_REF] builds a reference octree by propagating the parent octet to all children nodes, thus providing 255 contexts to encode the current octant. Octree coding allows for a progressive representation of point clouds since each level of the octree is a downsampled version of the point cloud. However, a drawback of octree representation is that, at the first levels of the tree, it produces "blocky" scenes, and geometry information of point clouds (i.e., curve, plane) is lost. The authors of [START_REF] Zhu | Lossless point cloud geometry compression via binary tree partition and intra prediction[END_REF] proposed a binary tree based method which analyzes the point cloud geometry using binary tree structure and realizes an intra prediction via the extended Travelling Salesman Problem (TSP) within each leaf node. Instead, in this paper, we employ a hybrid octree/voxel representation to better exploit the geometry information. Besides, the methods in [START_REF] Garcia | Context-based octree coding for point-cloud video[END_REF]- [START_REF] Garcia | Geometry Coding for Dynamic Voxelized Point Clouds Using Octrees and Multiple Contexts[END_REF] produce frequency tables which are collected from the coarser level or the previous frame and must be transmitted to the decoder. Our method predicts voxel distributions in a sequential manner at the decoder side, thus avoiding the extra cost of transmitting large frequency tables.

B. Generative Models and Learning-based Compression

Estimating the data distribution from a training dataset is the main objective of generative models, and is a central problem in unsupervised learning. It has a number of applications, from image generation [START_REF] Theis | Generative image modeling using spatial lstms[END_REF]- [START_REF] Salimans | PixelCNN++: Improving the PixelCNN with Discretized Logistic Mixture Likelihood and Other Modifications[END_REF], to image compression [START_REF] Weinberger | Applications of universal context modeling to lossless compression of gray-scale images[END_REF]- [START_REF] Mentzer | Conditional Probability Models for Deep Image Compression[END_REF] and denoising [START_REF] Chen | Image blind denoising with generative adversarial network based noise modeling[END_REF]. Among the several types of generative models proposed in the literature [START_REF] Goodfellow | Generative adversarial nets[END_REF], auto-regressive models such as PixelCNN [START_REF] Oord | Pixel Recurrent Neural Networks[END_REF], [START_REF] Salimans | PixelCNN++: Improving the PixelCNN with Discretized Logistic Mixture Likelihood and Other Modifications[END_REF] are particularly relevant for our purpose as they allow to compute the exact likelihood of the data and to generate realistic images, although with a high computational cost. Specifically, PixelCNN factorizes the likelihood of a picture by modeling the conditional distribution of a given pixel's color given all previously generated pixels. These conditional distributions only depend on the possible pixel values with respect to the scanned context, which imposes a causality constraint. PixelCNN models the distribution using a neural network and the causality constraint is enforced using masked filters in each convolutional layer. Recently, this approach has also been employed in image compression to yield accurate and learnable entropy models [START_REF] Mentzer | Conditional Probability Models for Deep Image Compression[END_REF]. Our paper explores the potential of this approach for point cloud geometry compression by adopting and extending conditional image modeling and masking filters into the 3D voxel domain.

Inspired by the success in learning-based image compression, deep learning has been widely adopted in point cloud coding both in the octree domain [START_REF] Huang | OctSqueeze: Octree-Structured Entropy Model for LiDAR Compression[END_REF], [START_REF] Biswas | Muscle: Multi sweep compression of lidar using deep entropy models[END_REF], voxel domain [START_REF] Guarda | Point cloud coding: Adopting a deep learning-based approach[END_REF]- [START_REF] Guarda | Point cloud geometry scalable coding with a single end-to-end deep learning model[END_REF] and point domain [START_REF] Yan | Deep autoencoderbased lossy geometry compression for point clouds[END_REF]- [START_REF] Wang | Multiscale point cloud geometry compression[END_REF]. Recently, the authors of [START_REF] Huang | OctSqueeze: Octree-Structured Entropy Model for LiDAR Compression[END_REF] proposed an octree-based entropy model that models the probability distributions of the octree symbols based on the contextual information from octree structure. This method only targets static LiDAR point cloud compression. The extension version for intensity-valued LiDAR streaming data using spatio-temporal relations is proposed in [START_REF] Biswas | Muscle: Multi sweep compression of lidar using deep entropy models[END_REF]. However, these methods target dynamically acquired point clouds, while in this paper we mainly focus on dense static point clouds.

Working in the voxel domain enables to easily extend most 2D tools, such as convolutions, to the 3D space. Many recent 3D convolution based autoencoder approaches for lossy coding [START_REF] Quach | Learning Convolutional Transforms for Lossy Point Cloud Geometry Compression[END_REF]- [START_REF] Guarda | Point cloud geometry scalable coding with a single end-to-end deep learning model[END_REF] compress 3D voxelized blocks into latent representations and cast the reconstruction as a binary classification problem. The authors of [START_REF] Yan | Deep autoencoderbased lossy geometry compression for point clouds[END_REF] proposed a pointnetbased auto-encoder method which directly takes points as input rather than voxelized point cloud. To handle sparse point clouds, recent methods leverage advances in sparse convolution [START_REF] Choy | 4d spatio-temporal convnets: Minkowski convolutional neural networks[END_REF], [START_REF] Graham | Submanifold sparse convolutional networks[END_REF] to allow point-based approaches [START_REF] Huang | 3d point cloud geometry compression on deep learning[END_REF], [START_REF] Wang | Multiscale point cloud geometry compression[END_REF]. For example, the proposed lossy compression method in [START_REF] Wang | Multiscale point cloud geometry compression[END_REF] progressively downscale the point cloud into multiple scales using sparse convolutional transforms. Then, at the bottleneck, the geometry of scaled point cloud is encoded using an octree codec and the attributes are compressed using a learning-based context model. In contrast, in this paper, we focus on dense voxelized point clouds and losslessly encode each voxel using the learned distribution from its 3D context. In addition, we apply this approach in a block-based fashion, which has been successfully employed in traditional image and video coding.

III. PROPOSED METHOD

A. System overview

In this work, we propose a learning-based method for lossless compression of point cloud geometry. We aim at minimizing the encoded rate measured by the number of bits per occupied voxel (bpov) by exploiting the spatial redundancies within point cloud. The general scheme of our method is shown in Figure 1. A point cloud voxelized over a 2 n ×2 n ×2 n grid is known as an n-bit depth PC, which can be represented by an n level octree. In this work, we represent point cloud geometry in a hybrid manner, by combining the octree and voxel domains. We coarsely partition an n-depth point cloud up to level n -6. This allows to coarsely remove most of the empty space in the point cloud. As a result, we obtain a n -6 level octree and a number of non-empty binary blocks v of size 2 6 × 2 6 × 2 6 voxels, which we refer to as resolution d = 64 or simply block 64 (Figure 1 The context is chosen adaptively following a rate optimization algorithm (Section III-C). The high-level octree, partitioning signal, selected context as well as the depth of each block are converted to bytes and signaled to the decoder as side information. We first define a 3D raster scan order that scan voxel by voxel in depth, height and width order. For ease of notation, we index all voxels in block v at resolution d from 1 to d 3 in raster scan order with:

v i = 1, if i th voxel is occupied 0, otherwise. (1)

B. VoxelDNN

Our method losslessly encodes the voxelized point cloud using context-adaptive binary arithmetic coding. Specifically, we focus on estimating accurately a probability model p(v) for the occupancy of a block v composed by d × d × d voxels. We factorize the joint distribution p(v) as a product of conditional distributions p(v i |v i-1 , . . . , v 1 ) over the voxel volume:

p(v) = d 3 Π i=1 p(v i |v i-1 , v i-2 , . . . , v 1 ).
( Each term p(v i |v i-1 , . . . , v 1 ) above is the probability of the voxel v i being occupied given the occupancy of all previous voxels, referred to as a context. Figure 2(a) illustrates such a 3D context. We estimate p(v i |v i-1 , . . . , v 1 ) using a neural network which we dub VoxelDNN.

) (a) 3D voxel context (b) 3D type A mask 2
The conditional distributions in (2) depend on previously decoded voxels. This requires a causality constraint on the VoxelDNN network. To enforce causality, we extend to 3D the idea of masked convolutional filters, initially proposed in PixelCNN [START_REF] Oord | Pixel Recurrent Neural Networks[END_REF]. Specifically, two kinds of masks (A or B) are employed. Type A mask is filled by zeros from the center position to the last position in raster scan order as shown in Figure 2(b). Type B mask differs from type A in that the value in the center location is 1 (colored in red). Type A masks are used in the first convolutional filter to remove the Fig. 3: VoxelDNN architecture, d is the dimension of the input block, masked layers are colored in yellow and blue. A type A mask is applied to the first layer (dashed borders) and type B masks afterwards. 'f64,k7,s1' stands for 64 filters, kernel size 7 and stride 1. Only probabilities of voxels being occupied are kept after the last Softmax layer. connections between all future voxels and the voxel currently being predicted. From the second layer, the value of the current voxel is not used in its spatial position and is replaced by the result of the convolution over previous voxels. As a result, from the second convolutional layer, type B masks are applied which relaxes the restrictions of mask A by allowing the connection from the current spatial location to itself. In order to learn good estimates p(v i |v i-1 , . . . , v 1 ) of the underlying voxel occupancy distribution p(v i |v i-1 , . . . , v 1 ), and thus minimize the coding bitrate, we train VoxelDNN using cross-entropy loss. That is, for a block v of resolution d, we minimize :

H(p, p) = E v∼p(v)   d 3 i=1 -log p(v i )   . (3)
It is well known that cross-entropy represents the bitrate cost to be paid when the approximate distribution p is used instead of the true distribution p. [START_REF] Quach | Learning Convolutional Transforms for Lossy Point Cloud Geometry Compression[END_REF], [START_REF]Improved Deep Point Cloud Geometry Compression[END_REF]. In this lossy context, the motivation behind using focal loss is to cope with the high spatial unbalance between occupied and non-occupied voxels. The reconstructed PC is then obtained by hard thresholding p(v), and the target is thus the final classification accuracy. Conversely, here we aim at estimating accurate soft probabilities to be fed into an arithmetic coder. Figure 3 shows our VoxelDNN network architecture for a block of dimension d. Given the d × d × d input block, VoxelDNN outputs the predicted occupancy probabilities of all input voxels. Our first 3D convolutional layer uses 7×7×7 kernels with a type A mask. Type B masks are used in the subsequent layers. To avoid vanishing gradients and speed up the convergence, we implement two residual blocks [START_REF] He | Deep residual learning for image recognition[END_REF] with 5 × 5 × 5 kernels. Since type A masks are applied at the first layer, identity skip connection of residual block does not violate the causality constraint. Throughout VoxelDNN, the ReLu activation function is applied after each convolutional layer, except in the last layer where we use softmax activation. Using more filters generally increases the performance of VoxelDNN, at the expense of an increase in the number of parameters and computational complexity. After experimenting with various number of filters, we concluded that for input voxel block (d×d×d×1) which only has a single feature, 64 convolutional filters give a good trade-off between complexity and model performance.

C. Multi-resolution encoder and adaptive partitioning

We use an arithmetic coder to encode the voxels sequentially from the first voxel to the last voxel of each block in a generative manner. Specifically, every time a voxel is encoded, it is fed back into VoxelDNN to predict the probability of the next voxel. Note that at this prediction step, all future voxels are filled with zeros. Then, we pass the probability to the arithmetic coder to encode the next symbol.

However, applying this coding process at a fixed resolution d (in particular, on larger blocks) can be inefficient when blocks are sparse, i.e., they contain only a few occupied voxels. This is due to the fact that in this case, there is little or no information available in the receptive fields of the convolutional filters. To overcome this problem, we propose to optimize the block size based on a rate-optimized multiresolution splitting algorithm as follows. We partition a block into 8 sub-blocks recursively and signal the occupancy of sub-blocks as well as the partitioning decision (0: empty, 1: encode as a single block, 2: further partition). The partitioning decision depends on the bit rate after arithmetic coding. If the total bitstream of partitioning flags and occupied sub-blocks is larger than encoding the parent block as a single block, we do not perform partitioning. The details of this process are shown in Algorithm 1. The maximum partitioning level or the maximum number of block sizes is controlled by maxLv and partitioning is performed up to maxLv = 5 corresponding to a smallest block size of 4. Depending on the output bits of each partitioning solution, a block of size 64 can contain a combination of blocks with different sizes. Figure 4

D. Context extension

We have discussed our multi-resolution encoder with multiple block sizes to adapt to the point cloud structure. However, with smaller block sizes, an implicit context model (using the content of the block) will be less efficient because the context may be too small. Therefore, we extend the context of each block to the encoded voxels that are above and on the left of the current voxel (causality constraint). Figure 5 illustrates the context before and after extension. Before extending the context, to encode voxel v c , only voxels from v 1 to v i-1 in Figure 5(a) are considered as contexts. After extending the context to the bigger block, the context is now composed of all voxels in the blue area in Figure 5(b). The white area represent inactive voxels, i.e., not used in Eq. [START_REF] Jang | Video-Based Point-Cloud-Compression Standard in MPEG: From Evidence Collection to Committee Draft [Standards in a Nutshell[END_REF]. Extending the context does not change the partitioning algorithm discussed above, although it might change the optimal selected partitions. Also, the causality is still enforced as long as we use masked filters in our network.

However, extending to a larger context is not always efficient when the extension area is sparse or contains noise, therefore we employ a rate-optimized block extension decision. To limit the computational complexity, we only allow certain combinations of block sizes and extension sizes, as shown in Table I. To encode a block with context extension, in Algorithm 1, we encode a block with all the possible extension sizes and select the best one in terms of bpov. In total, we build 5 models for 5 input sizes which are {128, 64, 32, 16, 8} in the context extension mode.

E. Data augmentation

In order to train more robust probability estimation models and to increase the generalization capabilities of our model, we employ data augmentation techniques specifically suited for PCC. In particular, we observed that methods based on convolutional neural networks are especially sensitive to changes in PC density and acquisition noise. Therefore, in addition to typical rotation and shifting data augmentation used for other PC analysis tasks [START_REF] Wang | Multiscale point cloud geometry compression[END_REF], [START_REF] Alonso | 3d-mininet: Learning a 2d representation from point clouds for fast and efficient 3d lidar semantic segmentation[END_REF], we also consider here alternative techniques, such as downsampling. Note that even though our VoxelDNN operates on voxel domain, to reduce the complexity, all input pipelines process point clouds in the form of x, y, z coordinates before converting into dense block in the final step. Specifically, for each generated block from the training datasets, we rotate them by an angle θ around each x, y, z axis. In addition, to adapt to varying density levels of the test point clouds, we randomly remove points from the original block as well as rotated blocks by controlling the probability of an occupied voxel being kept f s (f s ∈ [0, 1]). Figure 6 shows our data augmentation methods applying on Longdress point cloud from MPEG.

IV. EXPERIMENTAL RESULTS

A. Experimental Setup 1) Training dataset:

We consider point clouds from different and varied datasets, including ModelNet40 [START_REF] Wu | 3D ShapeNets: A deep representation for volumetric shapes[END_REF] which contains 12,311 models from 40 categories and three smaller datasets: MVUB [START_REF] Loop | Microsoft voxelized upper bodies -a voxelized point cloud dataset[END_REF], MPEG CAT1 [START_REF]Common test conditions for PCC[END_REF] and 8i [START_REF] D'eon | 8i Voxelized Full Bodies -A Voxelized Point Cloud Dataset[END_REF], [START_REF]Common test conditions for PCC[END_REF]. We uniformly sample points from the mesh models from ModelNet40 and then scale them to voxelized point clouds with 9 bit precision. To enforce the fairness between the smaller datasets in which we select point clouds for testing, point clouds from MPEG CAT1 are sampled to 10 bit precision as in MVUB and 8i. In addition, we measure the local density ρ of a point cloud, computed as the average portion of occupied voxels in the blocks of size 64, that is:

ρ = 1 N B × Bi∈B 100 × number of points in B i 64 3 (%) (4)
where B is the set of occupied blocks of size 64, and N B is the cardinality of B. The higher the value of ρ is, the denser the point cloud. The selected point clouds, number of frames as well as ρ of the training data are shown in Table II.

To train a VoxelDNN model of size d we divide all selected PCs into occupied blocks of size d × d × d. Table III reports the number of blocks from each dataset for training, with the majority coming from the ModelNet40 dataset. For the models trained with data augmentation, we apply rotation with θ = 45 • on x, y, z axis and then sampling from all blocks with sampling rate f s = [0.7; 0.4]. In total, we augment from each block to 12 variations in terms of density and rotation which significantly increase the volume and diversity of our training set.

2) Test data: We evaluate the performance of our approach on a diverse set of point clouds in terms of spatial density and content type.

B. Performance evaluation and ablation studies

In the following, we evaluate the performance of the proposed approach as well as the impact of its various components. We start with models without data augmentation nor context extension in order to study the optimal maximal partitioning depth for our method and establish a baseline for the evaluation. Next, on top of the best encoder in this experiment (Baseline), we separately add data augmentation (Baseline + DA) and context extension (Baseline + CE). Finally, Baseline + DA + CE incorporates both data augmentation and context extension. We compare our method against the state-of-the-art point cloud compression method G-PCC v12 from MPEG [START_REF] Graziosi | An overview of ongoing point cloud compression standardization activities: video-based (V-PCC) and geometry-based (G-PCC)[END_REF] which has a dedicated lossless geometry mode for static point clouds. We report the number of bits per occupied voxel (bpov) for each test point cloud and the average per dataset.

In all experiments, the high-level octree plus partitioning signal are directly converted to bytes without any compression. For the encoders with context extension, we signal the selected size using two bits (maximum 4 options on block 8). This information is also directly converted to bytes in the bitstream. On average, signaling bits account for 2.44% of the bitstream.

1) Optimal maximum partition depth:

To evaluate the effectiveness of the partitioning scheme, we increase the maximum partitioning level from 1 to 5, corresponding to a minimum block size of 64, 32, 16, 8, and 4. As 3D convolution is not able to efficiently exploit voxel relations on a very small receptive field, we do not train a separate model for block 4. Instead, we use the model trained on blocks of size 8 to predict its probabilities.

Table IV shows the average bpov of our encoder on the 4 test datasets at 4 partitioning levels. The results with 5 partitioning levels are identical to 4 partitioning levels and are not shown in the table. We observe that, as partitioning levels increases, the corresponding gain over single-level also increases. However, adding the 3 rd and 4 th level yields only a slight improvement compared to adding the 2 nd level. This can be explained with Figure 8 showing the percentages of occupied voxels in each partition size. We observe that most voxels are encoded using blocks 64 and 32, while very few voxels are encoded using blocks of smaller size. Adding more partitioning levels enables to better adapt to point cloud geometry, however, this is not often compensated by a bitrate reduction of the sub-blocks, since in the smaller partitions the encoder has access to limited contexts, resulting in less accurate probability estimation. However, there is an increase in the portion of block 32 and 16 on CAT1 and USP compared to MVUB and 8i. This reflects the density characteristics of each dataset: on sparser datasets (CAT1 and USP), the algorithm tends to partition point cloud into smaller blocks to eliminate as much empty space as possible. Based on these observations, we use a maximum of 4 partitioning levels for our baseline codec in later experiments. 2) Comparison with G-PCC: In table V, we report the output bitrate of our methods to compare with MPEG G-PCC. Both our method and G-PCC perform better on dense PCs while having higher rates on sparser PCs. Compared to G-PCC, the Baseline encoder obtains a significant gain -over 29% bitrate reduction on dense point clouds from MVUB and 8i dataset. On CAT1 and USP datasets, our method achieves a comparable rate with G-PCC. In particular, for Arco Valentino and BumbaMeuBoi, the two point clouds having the lowest ρ, our baseline codec yields a rate higher than G-PCC (+7.25% and +5.99%, respectively). For point clouds with high local density levels, our VoxelDNN could efficiently leverage the relations between voxels and predict more accurate probability. In contrast, probability prediction is less accurate on sparser point clouds.

This can be partially solved by adding data augmentation during training. Indeed, by random subsampling the point clouds in the training set, VoxelDNN learns to predict more accurate probabilities when the point cloud is less dense. Baseline + DA yields higher gains over G-PCC on CAT1 and USP compared to Baseline, with average bitrate reductions of about 1.68% and 3.67%, respectively. On the other hand, we observe a small degradation of the performance compared to Baseline for denser datasets, such as MVUB and 8i dataset. This is somehow expected, as data augmentation increases the generalization capability of VoxelDNN, which instead is more adapted to denser PCs in the baseline mode.

The encoder with context extension, Baseline + CE, obtains a better rate on all test point clouds compared to the Baseline encoder, regardless of the density, with an average further bitrate reduction of 4.8% over G-PCC. The cost to be paid for this performance improvement is a higher computational complexity in the encoding process.

The last two columns of Table V show the experiment results for the encoder incorporating both data augmentation and context extension, Baseline + DA + CE. On average, we have a higher gain than Baseline and Baseline + DA because of the Context Extension. As expected, comparing with Baseline + CE, Baseline + DA + CE has increasing gains on CAT1 and USP datasets while obtaining a lower gain on MVUB and 8i datasets. Despite the different performance trends for different densities of the input point clouds, we obtain, on average, a bitrate reduction of 20.66% compared to G-PCC. Note that, in practice, if the characteristics of point cloud to be coded are known in advance, our approach is flexible, in that we could deploy different models targeting a specific application (cultural heritage, tele-immersive conferencing, etc.) and content type to obtain the best compression rate.

3) Effect of PC content and density on coding performance: In order to better understand the performance of our codec for different types of content, we plot in Figure 9 From this figure, we can draw some observations. First, most of blocks are partitioned into 3 levels (smallest dots) and the majority of the remaining blocks are partitioned into 2 or 4 levels. Second, in each point cloud, denser blocks are easier to compress, as mentioned before, due to the better capabilities of convolution to capture spatial relations. On the other hand, our approach becomes inefficient when the blocks are less dense, and the bitrate associated to the very sparse blocks rapidly grows by an order of magnitude compared to the rest. This phenomenon is true for all kinds of contents, although it has a stronger effect when the block density distribution is skewed to the left, such as for Arco Valentino or BumbaMeuBoi, which have the highest bitrates in our experiments.

We can also observe a content-dependence trend in the figure, which appears like a vertical offset for different PCs. Arco Valentino and RomanOilLight overall have higher bpov compared to Phil and Loot with the same number of occupied voxels. This suggests that local density alone is not the only factor affecting the performance of our approach, but that somehow higher-order statistics enter into play. We will speculate more about this behaviour when discussing the bitrate allocation in Figure 11. Further analysis of this trend, as well as how to take better into account the PC characteristics to improve coding performance, are left to future work. 4) Selection of context extension and impact on the partitioning: Figure 10 shows how many times an extended block size is selected in the Baseline + DA + CE experiments. First, it can be seen that in most cases our encoder choose to extend the context to encode the current block, and mostly the immediate larger size is selected. By extending context to exploit geometry information from the neighboring voxels, VoxelDNN can leverage a larger amount of information and predict a better probability. In most cases where the encoder does not extend the context, the blocks are on the border of the volume, corresponding to a mostly empty extending area.

By summing the quantities in each column, we obtain the number of blocks which are encoded using each block size and we observe that large parts of the point cloud are partitioned into block 32 or 16. This is in contrast with the previous observation on baseline experiments where the most frequent partitions are 64 and 32 (Figure 8). This has an intuitive explanation: without context extension, small block sizes of 32 or 16 were insufficient to provide a representative enough context for VoxelDNN in most of the cases, even if they would better adapt to areas with low point density. Conversely, the context extension allows to compensate for the small block dimension and renders these modes competitive. As a result, context extension significantly affects the optimal partitioning and enables VoxelDNN to adapt better to local sparsity while still providing enough contextual information to predict accurate probabilities.

5) Using multiple models for the context: For the multiresolution encoder, instead of using a separate model for each block size, VoxelDNN can use only a single neural network to predict the distribution. Specifically, we place small blocks [START_REF] Garcia | Context-based octree coding for point-cloud video[END_REF][START_REF]Intra mode for geometry coding[END_REF][START_REF]Improved Deep Point Cloud Geometry Compression[END_REF] into a block of size 64 and then use the network for block 64 to predict and extract the corresponding distributions. This method of computing the occupancy distribution is different from Context Extension in that the surrounding voxels are always set to 0. In Table VI, we compare the performance of using a single model with Baseline, which is a multi-models encoder. In this experiment, both encoders have 4 maximum partitioning levels and use the same model 64. On average, by having a separate model for each block size, a multi-model encoder obtains about 1% additional gain over G-PCC compared to the single model encoder. This amount of gain indicates that the bigger VoxelDNN model can predict the conditional distribution on smaller blocks as efficiently as using a separate model for each block size. However, model 64 is trained on blocks of size 64 only, and learns features at that scale. In general, a model trained on small blocks could better capture the context from small input blocks and thus provides a higher gain in some circumstances. 6) Visualization of the bitrate allocation on coded PCs: The bpov heatmaps of 4 point cloud are shown in Figure 11. The blocks in the figures reflect the optimal partitioning obtained by the algorithm. First, we visually confirm what found in Figure 9, i.e., VoxelDNN performs better, i.e., achieves a small bitrate, in the smooth and dense areas of the point cloud. Conversely, in the noisy areas (Phil's hand, Loot's hand), sudden holes (Arco Valentino) or very sparse regions (edges in Arco Valentino, the bottom part of BumbaMeuBoi), which are indicated in red, the performance is worse. We can argue that the density of a point cloud, together with the smoothness and noise characteristics of the content, are among the main factors that influence the performance of VoxelDNN. On the other hand, we can argue that noisy and very sparse areas are intrinsically difficult to code in general, and indeed also the MPEG G-PCC codec requires a large number of bits to encode point clouds such as BumbaMeuBoi and Arco Valentino.

C. Computational complexity analysis

A well-known drawback of auto-regressive generative models such as PixelCNN and VoxelDNN is the sequential generation of the symbol probabilities. This requires to run the network for each voxel, which has a complexity that increases linearly with the number of voxels. Therefore, VoxelDNN has a computational complexity which is 3 orders of magnitures bigger than G-PCC.

Table VII reports the encoding and decoding time of our Baseline and Baseline + CE. Tests are benchmarked on an Intel(R) Xeon(R) Silver 4110 CPU @ 2.10GHz machine with an Nvidia GeForce GTX 2080 GPU and 16 GB of RAM, running Ubuntu 16.04. Our encoding time is highly dependent on the number of blocks and the number of voxels within each block. Besides, the number of modes in the partitioning algorithm and context extension also influence the complexity. The Baseline + CE encoder tries all the extending modes and selects the best one, thus its average encoding time is higher than Baseline -an increase of about 196%. The reason why the encoding time for the Baseline codec is lower than the decoding time is purely implementative: at the encoder it is possible to predict the whole block probabilities in a single batch on a GPU, while in a realistic scenario, at the decoder side the voxels need to be individually decoded. When context extension is enabled, point clouds are partitioned into even smaller blocks, corresponding to a smaller complexity at the decoder, as a smaller number of voxels need to be decoded. On the other hand, the total parameters of each VoxelDNN model corresponds only to about 3.5 MB which is a small-size network in practice. Notice that the bottleneck in our system comes from the adoption of an auto-regressive model, which has the advantage of providing, in principle, an exact likelihood estimation of the data, though at a high computational cost. We are currently investigating the use of alternative generative approaches that avoid sequential probability estimation.

V. CONCLUSIONS AND FUTURE WORK

This paper presents a lossless compression method for point cloud geometry. We extend a well-known auto-regressive generative model initially proposed for 2D images to the 3D voxel space, and we incorporate 3D data augmentation to efficiently exploit the redundancies between points. This approach enables to build accurate probability models for the arithmetic coder. As a result, when using an adaptive partitioning scheme and context extension, our solution outperforms MPEG G-PCC over a diverse set of point clouds.

Our analyses on the performance of the proposed method indicate at least two major avenues for improvement. On one hand, handling low-density point clouds would require to rethink the network architecture to handle sparse input data. On the other hand, a major drawback of VoxelDNN is the high computational cost of sequential probability generation, which we plan to replace in the future by a more efficient generative model.

Fig. 1 :

 1 Fig. 1: Overview of the proposed method. (a): a n bit depth point cloud is partitioned down to the n -6 octree level, yielding occupied blocks of size 64 × 64 × 64. (b): We encode each block of 64 3 voxels as a single block (b1), or divide it into 8 children blocks (b2), depending on the total number of bits of each solution (partitioning level = 2). This procedure is repeated recursively for increasing partitioning levels up to 5. (c): For each occupied block of size d, the context model estimates the distribution of each voxel given the previously encoded voxels.

  (a)). Blocks 64 can be further partitioned at resolution d = {64, 32, 16, 8, 4} corresponding to maximum partitioning level maxLv = {1, 2, 3, 4, 5} as detailed in Section III-C. Figure 1(b) shows the multiresolution encoder with maxLv = 2. A block of size d can be encoded as a single block (b2) or partitioned into 8 sub-cubes (b1). We then encode each non-empty block (blocks in blue in the figure) using the proposed method in the voxel domain (Section III-B) and select the partitioning mode resulting in the smallest bpov. The overview of a single block encoder is shown in Figure 1(c). Our context model predicts the distribution of each voxel given all encoded voxels and pass it to an arithmetic coder to generate the final bitstream.

Fig. 2 :

 2 Fig. 2: (a): Example 3D context in a 5×5×5 block. Previously scanned elements are in blue. (b): 3 × 3 × 3 3D type A mask. Type B mask is obtained by changing center position (marked red) to 1.

Algorithm 1 : 4 if b is empty then 5 child_f 18 total_bit2

 14518 Block partitioning selection Input: block: B, current level: curLv, max level: maxLv Output: partitioning flags: f l, output bitstream: bits 1 Function partitioner(B, curLv, maxLv): 2 f l2 ← 2 ; // encode as 8 child blocks 3 for block b in child blocks of B do = sizeOf (bit2) + len(f l2) × 2; 19 f l1 ← 1; // encode as a single block 20 bit1 ← singleBlockCoder(B); 21 total_bit1 = sizeOf (bit1) + len(f l1) × 2; / * partitioning selection * / 22 if total_bit2 ≥ total_bit1 then 23 return f l1, bit1;

  shows 4 partitioning examples for an encoder with maxLv = 4. Note that VoxelDNN learns to predict the distribution of the current voxel based on previously encoded voxels. As a result, we can use a bigger model size to predict the probabilities for smaller input block size.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Partitioning a block of size 64 into: (a) a single block of size 64, (b): blocks of size 32, (c): 32 and 16, (d): 32, 16 and 8. Non-empty blocks are indicated by blue cubes.

Fig. 6 :

 6 Fig. 6: Example of data augmentation applied on the Longdress point cloud. (a) Original; (b) After removing color attributes; (c),(d),(e) Rotation with θ = 45 • on x, y and z axis; (f),(g) Sampling rate f s = 0.7 and f s = 0.4

Fig. 7 :

 7 Fig. 7: Point clouds in the test set. (a) Phil, (b) Ricardo (c) BumbaMeuBoi (d) RomanOilLight, (e) Redandblack, (f) Loot, (g) Thaidancer (h) Boxer, (i) Frog, (j) Arco Valentino, (k) Shiva.

Fig. 8 :

 8 Fig. 8: Percentage of occupied voxels encoded in each partition size. From top to bottom: block 8, 16, 32, 64. Most of occupied voxel are encoded in block 64 and block 32.

  the average bpov as a function of the percentage of occupied voxels for each block 64 of Phil, Loot, Arco Valentino and BumbaMeuBoi with the Baseline + DA + CE encoder. Notice that each block 64 can be split up to different partition levels, indicated by the size of the dots in the figure. The distribution of the density of blocks 64 is shown in the top panel.

Fig. 9 :

 9 Fig. 9: Performance on block 64 on four test point clouds. Each point corresponds to a block 64 with percentage of occupied voxels (ρ) and bpov (log scale) performance of that block. The size of each point indicates the partitioning level and each partitioning level was sized according to its frequency. Higher points indicate that VoxelDNN is performing worse. The marginal distributions of occupied voxels for each point cloud are on the top of the scatter plot.

Fig. 10 :

 10 Fig. 10: Number of extending block size for each block. (a) Phil, (b) Loot, (c) Arco Valentino, (d) BumbaMeuBoi. Most of the time, the encoder extend the context to neighboring voxels instead of independently encoding a block.

Fig. 11 :

 11 Fig. 11: Output geometry bitrate in bpov per block. (a) Phil, (b) Loot, (c) Arco Valentino, (d) BumbaMeuBoi. The heatmap bar below each subfigure shows the minimum and maximum bpov and the corresponding color.

TABLE I :

 I Extending block size

	Block size	Extending block size
	64	128,64
	32	64,32
	16	64,32,16
	8	64,32,16,8

TABLE II :

 II Training and Testing Point Clouds

	Training Set		Test Set
	Point Cloud	# Fr	ρ	Point Cloud	ρ
	MVUB, 10 bits		MVUB, 10 bits, dynamic upper body
	Andrews	6	1.70	Phil	1.64
	David	5	1.65	Ricardo	1.77
	Sarah	6	1.72
				8i, 10 bits, dynamic full body
	8i, 10 bits		Redandblack	1.49
	Soldier	9	1.51	Loot	1.43
	Longdress	9	1.52	Thaidancer	1.68
				Boxer	1.56
	CAT1, 10 bits
	Facade	1	1.20	CAT1, 10 bits, static cultural heritage
	Egyptian mask	1	0.12	Frog	1.13
	Statue klimt	1	0.89	Arco Valentino	0.45
	Head	1	1.43	Shiva	0.88
	House w/o roof	1	1.21
				USP, 10 bits, static cultural heritage
	ModelNet40, 9 bits		BumbaMeuBoi	0.18
	200 largest PCs	200	1.53	RomanOilLight	0.94

TABLE III :

 III Number of blocks in the training sets of each model.

		MVUB	8i	CAT1 ModelNet40	Total
	Model 128	1516	1101	677	2860	6154
	Model 64	5777	4797	2777	11147	24498
	Model 32	22082	20436	15243	50611	108372
	Model 16	87578	86106	45626	224951	444261
	Model 8	354617 349760 180037	986253	1870667

  All selected point clouds are either used in MPEG Common Test Condition or JPEG Pleno Common Test Condition to evaluate point cloud compression methods. As shown in TableIIthe test PCs can be categorized into four sets: Dense point clouds from 8i Labs. They are also dynamic voxelized point clouds but each sequence contains the full body of a human subject. In the test set, loot USP: an inanimate dataset from the University of São Paulo, Brazil, related to cultural heritage with 10 bits geometry precision[START_REF] Zuffo | University of sao paulo point cloud dataset[END_REF]. BumbaMeuBoi and RomanOil-Light are two selected point clouds from this dataset. PCs from USP dataset have simple shape with smooth surfaces. BumbaMeuBoi is the sparsest PC in our test set with the smallest ρ. Figure7illustrates the test point clouds.3) Training procedure: We train 5 models for 5 input block sizes, i.e., 128, 64, 32, 16, 8. The mini-batch sizes are 1, 8, 64, 128, 128, respectively. Our models are implemented in TensorFlow and trained with Adam[START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] optimizer, a learning rate of 0.001 for 80 epochs on a GeForce RTX 2080 GPU.1

	(1000) and redandblack (1510) are from 8i Voxelized Full
	Bodies (8iVFB v2) [45] while boxer and thaidancer are
	selected and downsampled to 10 bits from 8i Voxelized
	Surface Light Field (8iVSLF) dataset [46].
	• CAT1: static point clouds for cultural heritage and other
	related 3D photography applications [44]. We select
	Arco_Valentino_Dense_vox12, Frog_00067_vox12, and
	Shiva_00035_vox12 from this dataset and downsample
	to 10 bits to validate the performance of our method.
	PCs from this dataset are less dense compared to the
	previous two datasets. Frog_00067 has smoother surfaces
	compared to the other two PCs which contain rough
	surfaces.

• MVUB: Microsoft Voxelized Upper Bodies

[START_REF] Loop | Microsoft voxelized upper bodies -a voxelized point cloud dataset[END_REF]

-a dynamic voxelized point cloud dataset containing five subjects. For testing, we randomly select 2 frames from Phil (frame number 10) and Ricardo (76) sequences which are both very dense (high ρ) with smooth surfaces. • 8i: •

TABLE IV :

 IV Average rate in bpov per dataset at different partitioning levels and the gain over the encoder with 1 partitioning level.

			1 level	2 levels	3 levels	4 levels
	Dataset	Point Cloud	bpov	bpov	Gain	bpov	Gain	bpov	Gain
		Phil	0.8943	0.8295	-7.25%	0.8206	-8.24%	0.8205	-8.25%
	MVUB	Ricardo	0.8109	0.7511	-7.37%	0.7440	-8.25%	0.7440	-8.25%
		Average	0.8256	0.7903	-7.31%	0.7823	-8.25%	0.7823	-8.25%
		Redandblack	0.7920	0.7269	-8,22%	0.7191	-9.20%	0.7190	-9.22%
		Loot	0.7017	0.6347	-9.56%	0.6271	-10.63%	0.6271	-10.63%
	8i	Thaidancer	0.7941	0.7360	-7.32%	0.7298	-8.10%	0.7297	-8.11%
		Boxer	0.6462	0.5960	-7.77%	0.5901	-8.68%	0.5900	-8.70%
		Average	0.7335	0.6734	-8.22%	0.6665	-9.15%	0.6665	-9.16%
		Frog	1.9497	1.8406	-5.60%	1.8216	-6.57%	1.8214	-6.58%
	CAT1	Arco Valentino Shiva	5.4984 3.7964	5.2947 3.6632	-4.52% -3.51%	5.2051 3.6400	-5.33% -4.01%	5.2050 3.6403	-5.34% -4.11%
		Average	3.7482	3.5845	-4.54%	3.5569	-5.31%	3.5556	-5.34%
		BumbaMeuBoi	6.3618	5.8235	-8.46%	5.7305	-9.92%	5.7305	-9.92%
	USP	RomanOilLight	1.8708	1.7157	-5.14%	1.7030	-5.84%	1.7030	-5.84%
		Average	4.0853	3.7696	-6.80%	3.7168	-7.88%	3.7168	-7.88%

TABLE V :

 V Average rate in bpov of proposed method and percentage gains compared with MPEG G-PCC v12 (negative percentages mean bitrate reduction).

			G-PCC	Baseline	Baseline + DA	Baseline + CE	Baseline + DA + CE
	Dataset Point Cloud	bpov	bpov	Gain over	bpov	Gain over	bpov	Gain over	bpov	Gain over
					G-PCC		G-PCC		G-PCC		G-PCC
		Phil	1.1599	0.8205	-29.26%	0.8954	-22.80%	0.7601	-34.47%	0.8252	-28.86%
	MVUB	Ricardo	1.0673	0.7440	-30.29%	0.8235	-22.84%	0.6874	-35.59%	0.7572	-29.05%
		Average	1.1136	0.7823	-29.78%	0.8595	-22.82%	0.7238	-35.01%	0.7912	-28.95%
		Redandblack	1.0893	0.7190	-33.90%	0.7772	-28.65%	0.6645	-39.00%	0.7003	-35.71%
		Loot	0.9524	0.6271	-34.16%	0.6282	-34.04%	0.5766	-39.46%	0.6084	-36.12%
	8i	Thaidancer	0.9990	0.7297	-26.96%	0.7253	-27.40%	0.6769	-32.23%	0.6627	-33.66%
		Boxer	0.9492	0.5900	-37.81%	0.6573	-30.75%	0.5503	-41.94%	0.5906	-37.78%
		Average	0.9975	0.6665	-33.24%	0.6870	-30.21%	0.6171	-38.13%	0.6405	-35.79%
		Frog	1.8990	1.8214	-4.09%	1.7662	-6.99%	1.6971	-10.63%	1.7071	-10.11%
	CAT1	Arco Valentino Shiva	4.8531 3.6716	5.2050 3.6403	+7.25% -0.85%	5.0639 3.5838	+4.34% -2.39%	4.9923 3.4619	+2.87% -5.71%	4,9900 3.5135	+2.82% -4.31%
		Average	3.4746	3.5556	+0.77%	3.7413	-1.68%	3.3838	-2.61%	3.4035	-3.86%
		BumbaMeuBoi	5.4068	5.7305	+5.99%	5.3831	-0.44%	5.3580	-0.90%	5.066	-6.29%
	USP	RomanOiLight	1.8604	1.7030	-8.46%	1.7319	-6.91%	1.6130	-13.30%	1.6231	-12.76%
		Average	3.6336	3.7168	-1.24%	3.5575	-3.67%	3.4855	-7.10%	3.4855	-9.52%

TABLE VI :

 VI Single model and multi-models comparison.

			G-PCC	Single model	Multi-models
		Point Cloud	bpov	bpov	Gain over	bpov	Gain over
					G-PCC		G-PCC
	MVUB	Phil Ricardo Average	1.1599 1.0673 1.1136	0.8312 0.7541 0.7927	-28.33% -29.34% -28.81%	0.8205 0.7440 0.7823	-29.26% -30.29% -29.78%
		Redandblack	1.0893	0.7320	-32.80%	0.7190	-33.99%
		Loot	0.9524	0.6403	-32.77%	0.6271	-34.16%
	8i	Thaidancer	0.9990	0.7305	-26.87%	0.7297	-26.96%
		Boxer	0.9492	0.6008	-36.70%	0.5900	-37.84%
		Average	0.9975	0.6759	-32.24%	0.6665	-33.24%
		Frog	1.8990	1.8433	-2.93%	1.8214	-4.09%
	CAT1	Arco Valentino Shiva	4.8531 3.6716	5.2173 3.6595	+7.50% -0.32%	5.2050 3.6403	+7.25% -0.85%
		Average	3.4746	3.5734	+2.84%	3.5556	+0.77%
	USP	BumbaMeuBoi RomanOiLight Average	5.4068 1.8604 3.6336	5.7501 1.7094 3.7298	+6.34% -8.11% -2.64%	5.7305 1.7030 3.7168	+5.10% -8.46% -1.24%

TABLE VII :

 VII Average runtime (in seconds) of different encoders comparing with G-PCC v12

		G-PCC	Baseline	Baseline + CE
	(Enc)	2.90	298	885
	(Dec)	2.79	672	640

The source code is available at https://github.com/Weafre/VoxelDNN_v2.

He is currently with CNRS/LSS (Laboratoire de Signaux et Systemes, Gif sur Yvette, France), where he developed studies in Signal processing for communications and signal/image processing for multimedia applications, including source/protocol/channel coding/decoding. He is also investigating the connections between communication theory and networking as well as information theory and AI. He has been "directeur de recherches émérite" since March 2019.

He has published more than 100 articles in international journals, more than 300 papers in international conferences, and holds 29 patents. He is a co-author of the book 'Joint Source and Channel Decoding: A cross layer perspective with applications in video broadcasting" which appeared in 2009, Academic Press. He successfully advised or co-advised more than 60 PhD students, and two of them are now fellows of the IEEE.


Farzam Tajdari
email: f.tajdari@tudelft.nl

Mahsa Tajdari

Amin Rezaei

Discrete Time Delay Feedback Control of Stewart Platform with Intelligent Optimizer Weight Tuner

Keywords: Time-Delay, Integral Control, Stewart Platform, Artificial Neural Network, Optimal Solution

In the presence of complicated kinematic and dynamic, we present a generalizable robust control technique for the 6-Degree of Freedom (6DoF) Stewart integrated platform with revolving, time-delayed torque control actuators to achieve faster, and reliable efficiency for parallel control manipulators. The suggested optimal solution involves the construction of a timedelay Linear Quadratic Integral (LQI) controller integrated with an on-line Artificial Neural Network (ANN) as the cost function gain tuner. The controller is formulated to robustly mitigate the nonlinear system's real-time tracking error with large timedelay, which is implemented via ADAMS software. The method is validated through simulation experiments to demonstrate that the developed methodology is practical, optimum, and zero-error convergence.

I. INTRODUCTION

Parallel systems are recommended over serial manipulators for applications such as surgeries, scanning, and 3D printing that require higher speeds and acceleration, greater accuracy, and the lightest weight [START_REF] Taghirad | Parallel robots: mechanics and control[END_REF]. The invention of anatomically adaptable recovery equipment [START_REF] Nomura | Development of six-dof human ankle motion control device using stewart platform structure for fall prevention[END_REF] and the upcoming generation of Ultra Personalized Products and Services (UPPS) [START_REF] Yang | Posture-invariant 3d human hand statistical shape model[END_REF], [START_REF] Tajdari | Image-based modelling for adolescent idiopathic scoliosis: Mechanistic machine learning analysis and prediction[END_REF] are of particular concern due to the intense importance of robot-influenced well-being and human life. Rapid automatic breast scanning [START_REF] Sun | Robot-assisted decompressive laminectomy planning based on 3d medical image[END_REF] is an interesting scanning scenario for the study's potential implementation, considering the fact that it is a very complex procedure due to the flexibility and deformity of breast tissue, particularly in motion scanning. The extremely high precision necessitates a thorough knowledge of parallel robotics, including sensor or actuator time delays and implementation objections [START_REF] Shoham | Bone-mounted miniature robot for surgical procedures: Concept and clinical applications[END_REF]. This examination necessitates a thorough understanding of the vector, complex, and complicated control approaches. Well-known parallel robots with widespread use are derived from the mechanism suggested by Stewart platform [START_REF] Stewart | A platform with six degrees of freedom[END_REF] as the basis of recent parallel robots in the applications mentioned above. The complex approach explains the interaction between internal or external forces and torques, as well as the state of the end-effector controller, in terms of the dynamic equations' main problems. With respect to dynamics, many approaches have been used, including the momentum-based approach [START_REF] Lopes | Dynamic modeling of a stewart platform using the generalized momentum approach[END_REF], the Newton-Euler methodology [START_REF] Dasgupta | A newton-euler formulation for the inverse dynamics of the stewart platform manipulator[END_REF] and the Lagrangian approach [START_REF] Bingul | Dynamic modeling and simulation of Stewart platform[END_REF]. The key aim of the parallel robot control system architecture is to ensure an exact continuation for the target position and orientation of dynamic and static variables in the moving endeffector of the robot [START_REF] Merlet | Parallel robots[END_REF]. Due to their extreme complexity, the robot's dynamics have lower control methods for the rotary Stewart robot while a wide range of control mechanisms are available, for example, optimum feedback robotic control [START_REF] Tourajizadeh | Closed loop optimal control of a stewart platform using an optimal feedback linearization method[END_REF]; backstage adaptive control [START_REF] Huang | Adaptive backstepping tracking control of the stewart platform[END_REF]- [START_REF] Tajdari | Semi-real evaluation, and adaptive control of a 6dof surgical robot[END_REF]; and backstage adaptive control [START_REF] Velasco Pascua | Experimental validation of a sliding mode control for a stewart platform used in aerospace inspection applications[END_REF], [START_REF] Tajdari | Robust control of a 3-dof parallel cable robot using an adaptive neuro-fuzzy inference system[END_REF]. In order to directly exploit the movement of the end-effector, reverse dynamic controls are extremely necessary. In limited studies [START_REF] Yang | Dynamic modeling and decoupled control of a flexible stewart platform for vibration isolation[END_REF], however, rotary torque controls are explored and provide the possibility of further use of recognition technologies, non-linear testing methods and impedance control systems. Moreover, few of them studied time delays (particularly in the operational file of robots) as an undeniable part of robotic controllers in the studies on rotary control Stewart, which highlighted high accuracy. Time-Delay Control (TDC) is generally considered an optimal way to solve this problem. The TDC scheme is also used for large applications such as robot manipulators [START_REF] Jin | Robust control of robot manipulators using inclusive and enhanced time delay control[END_REF], and flexible robots [START_REF] Roy | Adaptive-robust control of a class of el systems with parametric variations using artificially delayed input and position feedback[END_REF]. Time-delay estimation (TDE) errors will eventually be estimated by using time-related controls in functional management approaches. There are also increased TDE errors, which allow the trust of the closed-loop system regulations to be extremely complicated and non-linear (e.g. saturation.) As a general solution, TDC is the main benefit of its model-free functionality.

In the face of a large time lag with the actuators to maximized torques and reduce tracking errors, the objective of this paper is to establish and test a modern, stable torque management strategy for the complex parallel system. In addition, the controller is evaluated by a valid non-linear test bed using the ADAMS model proposed for the intelligent implementation of [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF], [START_REF] Tajdari | Implementation, and intelligent gain tuning feed-back based optimal torque control of a rotary parallel robot[END_REF] and experimentally analyzed using a robot manufactured.

II. EQUATIONS OF MOTION

A. Kinematics equation

The power torques of the motors and their rotational angles are controllable variables on a Stewart base powered by rotat- ing engines. The relationship between the motor variables and the end-effector position is stated in the following equations:

L = f ( X), θ M = g( L), (1)
where vector X = [φ θ ψ x y z] T represents the end-effector's variables which are the rotation angles φ, ψ, and θ of the moving plate about the X, Y and Z axis respectively, and the position of the end-effector center of mass as x, y, and z with respect to the fixed coordinate XY Z, according to Fig. 1(a). Moreover, the vector

L = [L 1 L 2 L 3 L 4 L 5 L 6 ]
T representing the distance between the joints on the fixed base ( P i ) and the corresponding joints on the movable frame ( p i ), and eventually

θ M = [θ M1 θ M2 θ M3 θ M4 θ M5 θ M6 ]
T is the motor angular shown in Fig. 1(a). We've used the findings of the same model in [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF] to find the relationship between Li and θM i, then

θ Mi = -(asind C i A 2 i + B 2 i -acosd B i A 2 i + B 2 i + 180), (2)
where A i , B i , C i ∈ U ( X) and

A i = 2| u i | | L ix |Sin(θ h ) -| L iy |Cos(θ h ) , (3)
B i = 2| u i |(| L ix |Cos(β)Cos(θ h ) + | L iy |Cos(β)Sin(θ h ) -| L iz |Sin(β)), (4)
C i = | l i | 2 -| L i | 2 (5)
in which θ h is the angle between P i and X in global axis, and β is the installation angle between the motor and horizon. The equation explained in [START_REF] Nomura | Development of six-dof human ankle motion control device using stewart platform structure for fall prevention[END_REF] specifies that with a X on the end-effector sensor, the angle value of each motor is defined. What helps to monitor these robots location.

B. Kinetics equations

In this section, the dynamic equations of the Stewart platform powered by six rotary motors are obtained using the Newton-Euler method. The derivations are summarized only in order to illustrate the different dynamic properties of the systems and the effects. The Stewart system with six rotary motors as actuators studied in this paper is shown in Fig. 1 (a). The frame comprises, as can be seen, of a revolving plate as an end-actuator, of a fixed plate as a base and, as an end-actuator manipulator, of six legs connected to six rotary actuators. The circular joints are used to connect the end-effector and the spine to the six thighs. The kinetics equations for end-effector can be written as:

Σ M = Ī α, Σ F = m a ( 6
)
where α is the angular acceleration vector, and a is the linear acceleration vector of the end-effector. In addition, m contains the mass value (m) of the end-effector, and Ī includes the endeffector moment of inertia around x, y, and z axis, then

Ī 0 3×3 0 3×3 m X = M F (7) Ī =   I xx 0 0 0 I yy 0 0 0 I zz   , m =   m 0 0 0 m 0 0 0 m   (8)
and,

M = [M x , M y , M z ] T , F = [F x , F y , F z ] T . (9)
In Fig. 1(b), M and F represents the torques and the exerted forces on the end-effector, respectively. The complex equation should be driven in the following way, when discussing the torque-manipulation of the rotary motors

Ī 0 3×3 0 3×3 m Ẍ = M F = τ 6×6 T (10)
where T = [T 1 T 2 T 3 T 4 T 5 T 6 ] T , and each component of T as T i is the manipulated torque with motor i. The kinetics equations can also be formed using their unit vectors. We can write:

Σ M = Σ r i × F li = Σ p i × e li | F li | = Σ e Mi | F li | (11)
in which, e Mi = p i × e li (i.e., e denotes the unit vector of its corresponding vector). Then, the summation of the exact forces can be obtained by:

Σ F = Σ e li | F li |, (12)
consequently,

M F = e M1 • • • e M6 e l1 • • • e l6    F l1 . . . F l6    (13)
where,

τ 1 =   e M1 • • • e M6 e l1 • • • e l6 .   . (14)
As

F li = e li • e Ni | Ti|
| ui| where e Ni is the unit vector of T i , then Considering,

   F l1 . . . F l6    =     e l 1 • e N 1 | u1| • • • 0 . . . . . . . . . 0 • • • e l 6 • e N 6 | u6|        T 1 . . . T 6    . (15)
τ 2 =     e l 1 • e N 1 | u1| • • • 0 . . . . . . . . . 0 • • • e l 6 • e N 6 | u6|     (16)
one can obtain

Ī 0 3×3 0 3×3 m Ẍ = M F = τ 6×6 T = τ 1 τ 2 T. (17)
Therefore, the ultimate dynamic transfer matrix, from the endeffector to the base is τ = τ 1 τ 2 , which defines the equations of motion based on the end-effector states and the applied torques by the actuators, which enables us to directly implement the torque control method as below:

Ẍ = τ tot T (18)
where,

τ tot = Ī 0 3×3 0 3×3 m -1 τ. (19)

C. Non-linear model

In order to evaluate the controller design, it is vital to have a reasonable non-linear model to characterize the real device. While certain physical restrictions, such as collision of objects, stiffness and elasticity of rigid bodies, and the friction between hard surfaces, cannot not be considered in the MATLAB software, the ADAMS model can be efficiently used [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF] to examine these physical constrains. Then, the 3D machine model that is shown in Fig. 2 can be integrated into ADAMS and the basic constrains can be implemented on the robot model. The given parameters in Fig. 2(a) and Fig. 1(a) are as follows: R up =0.12 mm, R down =0.22 mm, θ 0 =5 • , and θ 1 =15 • . In addition, the dynamic parameters are reported in Table I in which, the system consists of fourteen components of mass and moment of inertia.

III. CONTROLLER DESIGN

Given the fact that the dynamic system in the ADAMS software is developed based on the nonlinear parameters of the real system, it can be seen as a nonlinear Stewart platform; thus, in order to proficiently monitor the system, the proposed controller has to minimize the tracking errors in the presence of nonlinearities [START_REF] Åström | PID controllers: theory, design, and tuning[END_REF]. Besides robustness, input signals optimization of the parallel manipulator is another problem to be tackled in the control design step. In this respect, an optimal control unit for the Linear Quadratic Integral (LQI) is a strong and suitable candidate.

A. State-space equations

The state space and dynamic equations of the Stewart robot are formulated as follows:

ẋ(t) = Āc x(t) + Bc (t)T (t) + dc (t) (20)
where x represents the main states in X, and the speed of the states. Consequently, the error dynamic is defined is the following fashion:

ė(t) = Āc ē(t) + Bc (t)T (t) + dc (t) (21)
where, vector dc depicts the time variant disturbance. ē is the error matrix of the states and is defined as ē = x -x d , where x d denotes the desired values of each states. In addition, matrix Āc includes all the a i,j elements in ( 22) which represents the system dynamics. Finally, matrix Bc , carries b i,j elements in [START_REF] Åström | PID controllers: theory, design, and tuning[END_REF] which produces the interconnections between the applied torques and corresponding states as follows:

a i,j = 1 if i = 2n -1, and j = i + 1 0 otherwise (22) b i,j = τ tot ( i 2 , j) if i = 2n, 0 otherwise (23)
where n = 1, • • • , 6. Since in any real system the presence of disturbances (i.e., uncertainties in dynamic parameters,) is inevitable, an integral controller is used to eradicate constant disturbances [START_REF] Åström | PID controllers: theory, design, and tuning[END_REF]. The key control issue resulting from [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF] can therefore be reformulated by taking into account integral states (z) of S (i.e., as many as the number of links,) described by:

ż(t) = Cc ē(t) (24)
where Cc consists of ci,j denoted in the following way:

ci,j = 1 if i = n, and j = 2i -1 0 otherwise [START_REF] Tajdari | Integrated ramp metering and lane-changing feedback control at motorway bottlenecks[END_REF] hence, by adding the integral states, the resulting dynamic error of the system can be written as bellow:

ė = Ae(t) + B(t)T (t) + d(t) (26)
e = ē z , d = dc 0 S×1 (27)
A = Āc 0 H×S Cc 0 S×S , B = Bc 0 S×F , C = Cc 0 n×S 0 S×2n I S×S . (28)
B. Time-delay control scheme

If we consider a time-delay of η for the system in (26) (η is the summation of all the delays including actuator and sensor), then the time variant form of the system is written as:

ė(t) = Ae(t) + BT (t -η) y(t) = Ce(t). (29)
Accordingly, the solution for the well-known differential equation in ( 29) is presented in [START_REF] Karafyllis | Stabilization of nonlinear delay systems using approximate predictors and high-gain observers[END_REF]:

ė(t) = e At(t-t0) e(t 0 ) + t t0 e A(t-ζ) BT (ζ)dζ (30)
In the discrete time domain and considering only one step delay, t can be considered as

t 0 = k ∆t, t = k ∆t + ∆t. (31)
Generally, it may happen that the resolution of the time step (∆t) cannot cover the delay of η, then η = n∆t + η 0 , and e(k + 1) =e A∆t e(k) + η0 0 e A(∆t-ν) Bdν T (k -1)

+ ∆t η0 e A(∆t-ν) Bdν T (k) (32)
where k∆t -ζ = -ν, and T (k -1) is the delayed output of T (k) for one time step. From now on, the µ number of steps delayed output is shown as T k-µ Which is considered as extra states to be controller as following,

e 1 dis = e T k-1 . (33)
Then, the state space in (32) is formulated as

e 1 dis (k + 1) = A 1 dis e 1 dis (k) + B 1 dis T (k) (34)
where,

A 1 dis = e A∆t η0 0 e A(∆t-ν) Bdν 0 n×(N +S) 0 n×n B 1 dis =   ∆t η0 e A(∆t-ν) Bdν ------- I n×n   , C 1 dis = C 0 S×n (35)
where, n is the number of the system inputs, S is the number of integral states, N denotes the number of elements in matrix ė in [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF], and I is the identity matrix. Now, we explain the discrete system with µ steps delay which can be written as:

e 1 dis (k + 1) = A 1 dis e 1 dis (k) + B 1 dis T k-µ y(k) = C 1 dis e 1 dis (k). (36)
Considering each delay line of the control signal as a single state of the discrete system and the total number of variable in e 1 dis as H = N + S + n, then, the state variable matrix of the discrete system with the new states is

e dis(H+µn)×1 =     
e 1 dis e dis H+1:H+n . . .

e dis H+(µ-1)n+1:H+µn      =      e 1 dis T k-µ . . . T k-1      . (37)
And, the system with µ-step delay is

e dis (k + 1) = A dis e dis (k) + B dis T k y(k) = C dis e dis (k) (38)
where,

A dis =   A 1 dis B 1 dis 0 H×(µ-1)n 0 (µ-1)n×H 0 (µ-1)n×n I (µ-1)n×(µ-1)n 0 n×H 0 n×n 0 n×(µ-1)n   B dis = 0 (H+(µ-1)n)×n I n×n , C dis = C 1 dis 0 S×(H+µn) . (39)
Now, we can define the following quadratic cost function over an infinite time horizon, which is defined for minimization of all states and control inputs:

minJ = ∞ k=0 [e dis (k) T C T dis QC dis e dis (k) + T T (k)RT (k)] (40) where, Q = ω Q I (H+µn)×(H+µn) , R = ω R I n×n . ( 41
)
The Q and R matrices are weight matrices that tend to minimize the error of all states and the control signals. And, ω Q > 0 and ω R > 0 are formulated for weighing matrices.

A Linear Quadratic Regulator (LQR) is used to solve the resultant optimal control problem (40), (41), in which the goal is to stabilize feedback gains via the assumption. Although the stability and detectability parameters are fulfilled according to Hautus-test in [START_REF] Tajdari | Integrated ramp metering and lane-changing feedback control at motorway bottlenecks[END_REF]- [START_REF] Tajdari | Adaptive traffic control at motorway bottlenecks with time-varying fundamental diagram[END_REF] extended to [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF] by the initial resulting method.

C. Controller design

In order to solve the LQI problem, a linear feedback control law is proposed as follows:

T (k) = -Ke dis (k) (42)
where

K = R + B T dis P B dis -1 B T dis P A dis , (43)
P = C T dis QC dis + A T dis P A dis -A T dis P B dis R + B T dis P B dis -1 . (44)
Resulting K in (43) as optimal gain, and the Algebraic Riccati Equation (44) are investigated in [START_REF] Navvabi | New afsmc method for nonlinear system with state-dependent uncertainty: Application to hexapod robot position control[END_REF].

D. Intelligent estimator design

As a function of time variables in the dynamic equation ( 20), B dis (k) is the time variable in (38) according to [START_REF] Åström | PID controllers: theory, design, and tuning[END_REF], which leads to the conclusion that the process behaves differently in various paths. Thus, using intelligent approaches is an achievement in solving non-linearity and maximum efficiency [START_REF] Ghaffari | New fuzzy solution for determining anticipation and evaluation behavior during car-following maneuvers[END_REF]- [START_REF] Li | Clustering discretization methods for generation of material performance databases in machine learning and design optimization[END_REF]. As the cost efficiency can be maximized by well-defining the w Q and w R in (40) it's a path dependence approach to approximate true variables by time.

Where w Q indicates the degree of state error minimization and w R optimizes the torque values. An ANN-based estimator will be defined where the input variables are the density of the functions of the state error and the torque density of the function and the outputs of the estimator are w Q and w R .

The density functions are derived from [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF] which describes the density function for any arbitrary parameter by the following equation:

ρ (k) = Ω ( (k)) = (w in + 1) (k) k i=k-win (i) (45)
where ρ denotes the parameter's density function. When the average of the parameter with window size w in + 1 is significantly less than (k), the ρ value becomes highly sensitive. As a consequence, the ρ that is appropriate is the one that has the least divergence from the average values. This gives the a parameter changes a smooth and rational instigation, according to [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF]. Therefore, the inputs of the intelligent estimator are ρ ē1×6 and ρ T1×6 , and the outputs are w Q1×1 , and w R1×1 as shown in Fig. 3. We used controller output, for w in = 30, and w Q and w R as 0.1, 1, 4 and various quantities of delay (µ) as 0, 5, 10, 20 and 30 intervals including 45 experimental intervals in all. We have repeated the above measurements with frequencies of 0.5 Hz, 2 Hz, 4 Hz, and 6 Hz to consider the frequency effects of the desired values, so that the final datasets are composed of 180, and the sampling time during these tests is 0.1s and the time spent is 10 s. That ends the data set with 18000 samples of the input-output.

Regarding the Fig. 3, the input layer consists of 18 neurons, the hidden layer contains 10 neurons, and the output layer includes 2 neurons. Also, the employed activation functions are linear and sigmoid for the input and the output, and the hidden layer respectively. In addition, the training algorithm of the ANN is a Levenberg-Marquardt back-propagation method based on [START_REF] Kosko | Neural networks and fuzzy systems[END_REF]. A collection of data, including the inputsoutputs, is needed to design an ANN-based estimator. As a result, the estimator's controller (42) is used as a master. Thus the randomly chosen 70% data set for training is considered, while the remaining 30% is used to test the model collected. The closed-loop optimum feedback management diagram of the proposed device displays in Fig. 4 where the ANN estimator changes the controller gains of w Q and w R at each step, depending on error states and torque values.

IV. SIMULATION RESULTS

The results of the simulation are produced by implementing w in = 30, µ = 8, and ∆t = 0.1 s in ADAMS, while the desired variables ( Xd ) are considered to be following:

Xd (t) = φ d (t), θ d (t), ψ d (t), x d (t), y d (t), z d (t) T = π 12 , π 12 , π 12
, 0.02, 0.02, .08

T Sin(t)+ 0, 0, 0, 0, 0, .38

T ( 46
)
where t is time. The corresponding simulation results are obtained by implementing a proposed controller (42). Two different control methodologies (controllers without integral states, integral controllers, intelligent time-delay controller without integral state, and intelligent time-delay integral controller) are compared, according to Fig. 5(a). While the first two controllers are technically stable without any delay in the actuator, as discussed in [START_REF] Tajdari | Intelligent optimal feedback torque control of a 6dof surgical rotary robot[END_REF], they cannot monitor the robot in the presence of a time-delay (cyan and green line) due to the unique dynamics of the robots (resulting from pair instability ( Āc , Bc )). In Section III-C, the other two controllers have a stable output, however, with finite errors. The intelligent integrated time-delay controller ensures flawless monitoring (looking at Fig. 5(a) on the blue lines and refuses errors as seen in Fig. 5 (b). The reduction in error was completed by an integral status zero error function and with well-specified w Q and w R controller time variant parameters recorded in Fig. 6 (a). The smart time-delay controller without any detailed states does, however, cause a substantial error, especially in the crucial points of the required path (i.e. the output of the controller with respect to a x, y and z in Fig. 5(b)). This is a discrepancy even though the controller attempts to generate a periodic motion close to the desired value. To compare the experiments with regard to the monitoring error, a criterion with the accumulated absolute percentage error E ac is added

E ac = 2 H H 2 i=1 t 0 | Xi (t) -Xdi (t)| A mi × 100 dt (47)
where E ac is the sum of all absolute amounts, from the beginning to the end of the simulation. Where i = 1, . . . , H 2 , t end = 10 s, and A mi represents the peak-to-peak magnitude of the desired values ( Xdi ) of the relative vector ( Xi ) (where the desire is constant, A mi is presumed to be 1). The output of the controllers is then defined separately on the basis of the parameters in TABLE II. As for the table, all controllers without the time-delay scheme (the controller with no integral states and the integral controller) failed as expected, at around t = 0.8 s looking at the Figs. 5(a) and 5(b). While timedelay controllers kept the system in the allowed work-space, the time-delay integral controller had a significant advantage over the time-delay controller without integral states, with a 64% improvement. This reflects integral states' ability to reject unknown shocks and uncertainties. In comparison, Fig. 7 displays the torques that have been applied on the device for the two time-delay management scenarios. As the torques in Fig. 7(b) are compared to the torques in Fig. 7(a), the control effort produced by the intelligent time-delay integral controller decreases over time as the errors decrease. In addition, the Integral controller's torques have more trembling (from t = 0 s to t = 5 s) than the torques of the controller without integral states, which is attributable to the integral states' increased sensitivity to the tracking error due to [START_REF] Karafyllis | Stabilization of nonlinear delay systems using approximate predictors and high-gain observers[END_REF]. Due to the significant tracking error, the amplitude of the control signals produced by the intelligent time-delay controller without integral states is increasing.

V. CONCLUSION

In this paper, we investigated an innovative control methodology for a Stewart platform parallel robot with delayed rotary actuator. We derived the equations of motion for the system, and optimally designed the dimensions of the platform. We also implemented a 3D model of the system in ADAMS software. The current dynamic is controlled online via the intelligent time-delay tuning of the LQI controller. The robustness, as well as the adaptive component, i.e., the ANN estimator, enables us for an easy practical implementation, without the necessity of lengthy and costly measurements of the dynamic parameters. Finally, the simulation results showed that the ideal time-delay LQI controller with an intelligent estimator was able to monitor the dynamics of the device and eradicate the error of monitoring given the presence of actual dynamic forces and the enormous delay of the actuators in the ADAMS model. Numerically, we showed that the time-delay intelligent integral controller achieved 64% improvement in terms of accumulative absolute percentage error (E ac ). Further developments include the study of the controller parameters sensitivity, robustness of the method in presence of noise and actuator saturation which will appear in our future scholar.

Fig. 1 :

 1 Fig. 1: Schematic of a Stewart platform. (a) Defined variables and vectors. (b) Dynamic force-torque diagram.

Fig. 2 :

 2 Fig. 2: 3D model of the Stewart platform designed with Solidworks software. (a) Top view. (b) Perspective view.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Structure of the ANN-based intelligent estimator.

Fig. 5 :

 5 Fig. 5: Controlled case performance via the ADAMS model. (a) Main states ( X). (b) Error of Main states ( X -Xd )

Fig. 6 :

 6 Fig. 6: Online estimated of w Q ,and w R . (a) Time-delay controller without integral states. (b) Time-delay integral controller.

Fig. 7 :

 7 Fig. 7: Motor torques. (a) Time-delay controller without integral states. (b) Time-delay integral controller.

TABLE I :

 I Dynamical features of the robot's components.

	Component Dimension	Inertia [Ix,Iy,Iz]	Mass	Quantity
		(m)	(×10 -4 kg.m 2 )	(kg)
	l i	0.404	[16, 16, 0.01]	0.08	6
	u i	0.11	[0.01, 0.58, 0.6]	0.04	6
	End-	Circle	[30, 30, 60]	4	1
	effector	(Rup=0.12)
	Base	Circle	[270, 270, 540]	7	1
		(R down =0.22)

TABLE II :

 II Accumulative absolute percentage error (E ac ) report.

	Controller	Eac
	Controller without integral states	NaN
	Integral controller	NaN
	Time-delay controller without integral states	225
	Time-delay integral controller	80


Olivier Pourret
email: olivier.pourret@unilasalle.fr

Jennifer Middleton

Dasapta Daniel E Ibarra

Erwin Irawan

Ashaki Rouff

Pallavi Anand

Aradhna Tripati

Amy J V Riches

Anthony Dosseto

Dasapta Erwin Irawan

Triple Point section -Elements 1 Diversity among Editorial Boards of Elements and other selected

Original sketch drawn by Dasapta Erwin Irawan

The publication of scientific work is foundational to our disciplines. To ensure equitable publication standards during the global flow of knowledge production, professional societies and publishers must take positive steps to avoid any biases that might hinder the publication of scientific work (e.g., [START_REF] Liévano-Latorre | Pervasive gender bias in editorial boards of biodiversity conservation journals[END_REF]. Biases among editors and reviewers can be unconscious and influenced by different aspects of an author's identity, from country of origin, first language, affiliation, gender, ethnicity, and/or other factors, that could result in challenges to publication rates and visibility in key journal forums for under-Triple Point section -Elements 3 represented groups [START_REF] Lerback | Association Between Author Diversity and Acceptance Rates and Citations in Peer-Reviewed Earth Science Manuscripts[END_REF]. Ensuring that there is diversity in the peer review and publishing process and on editorial boards may help to eliminate bias. Diversity promotes innovation from hypothesis through peer review to final publication (e.g., [START_REF] Hofstra | The Diversity-Innovation Paradox in Science[END_REF], and should be set as a new standard, as shown by

The Royal Society of Chemistry (https://www.rsc.org/new-perspectives/talent/jointcommitment-for-action-inclusion-and-diversity-in-publishing/). Personal identity impacts how we engage with our science; it impacts how we approach a problem, and what we value, study, and write. It influences how we select reviewers, how we review, and ultimately what is successfully published. Therefore, the limited diversity of major editorial boards will act as a barrier to representation of all academic members. The members of editorial boards shape the direction and success of a journal, and influence the authorship of papers and what is published within the journal. And, differences in scientific networks may be a core reason for persistence of implicit bias from editorial boards, particularly with regard to gender [START_REF] Hanson | Age, Gender, and International Author Networks in the Earth and Space Sciences: Implications for Addressing Implicit Bias[END_REF]. Therefore, editorial bias or perceived editorial bias, can exclude certain groups and exacerbate historic inequities regarding under-representation of entire continents within the geoscience literature (e.g., Africa; see [START_REF] North | Out of Africa: The underrepresentation of African authors in high-impact geoscience literature[END_REF]. Hence, for an editorial board to be inclusive and unbiased, it needs to be as diverse as the research community it represents that we know do not exactly represent the general population.

In this piece, we provide an exploration of diversity among editorial boards by presenting data for Elements, in terms of gender and geographic affiliation. We further compare these data with editorial board data from journals that are published by the societies that jointly publish Elements (thereafter termed "journals of participating societies").

Historically (from 2005 to 2021), 19 principal editors (PE) have served with Elements, among them four were women (21%). In the past 10 years, there have been 11 PEs, 3 of which have been women (27%). These numbers, though falling short of gender parity, are representative of the proportion of the mid-to late-career women in the field and are indicative of wider challenges to diversity in our discipline. The executive editor is also considered part of the editorial board: both executive editors have been women. So, Elements editorial board (executive and principal editors) has always had 25% to 50% women at any one time. All serving editors are white and affiliated with institutions in Northern America (n=13, 62%; USA and Canada) or Western Europe (n=8, 38%; UK, France, Germany, Denmark).

Examination of the editorial team among the journals of participating societies of Elements highlights the pervasiveness of editorial lack of diversity within our field (17 journals published and edited by various commercial publishers and learned societies;

Figure 1). Here we identify that as of April 2021, the current editorial boards span a range of gender representation from 7 (50%) men and 7 (50%) women for Geochemical Perspectives Letters to 37 (97%) men and 1 (3%) woman for Journal of Mineralogical and Petrological Sciences. Of the 683 total Editorial Board members, 142 (21%) are women and 539 (79%) are men. These numbers are comparable to the February 2021 Elsevier Benchmark Gender Diversity distribution across "Geochemistry & Planetary Science" (25% women, 74% men and 1% prefer not to disclose) or "Applied Geosciences" portfolio Editors (with 14% women, 86% men). It must be noted that numbers from Elsevier are from an incomplete voluntary survey of the editors yet provide a good indication of gender distribution. Moreover, in some of the journals (i.e., those with larger editorial boards), turnover happens more quickly and thus numbers can fluctuate (up to a few percent in six months).

Figure 1 Gender breakdown of editorial board members from selected journals (data accessed on journal websites on 8 April 2021). Editorial board members were assigned a binary gender using first names and in some cases, based upon the authors' own perceptions and knowledge, n stands for the number of editorial board members. This approach has multiple limitations: (i) it may misgender people and (ii) gender is not binary, non-binary people were not included in this first-stage analysis due to lack of available information.

Triple Point section -Elements 6 As with gender distribution, geographic distribution is also strongly biased (Fig. 2).

Country is based on the affiliation of editors and implicitly may create a bias when an editor originates from a country and moves to another country. Journals like Chemical Geology (36% from Northern America, 41% from Western Europe, 2% from Sub-Saharan Africa, 2% from Northern Africa and Western Asia, 2% from Central and Southern Asia, 5% from Eastern and South-Eastern Asia and 12% from Oceania) and Geochemistry: Exploration, Environment, Analysis (14% from Northern America, 10% from Latin America and the Caribbean, 38% from Western Europe, 10% from Eastern Europe, 5% from Sub-Saharan Africa, 10% from Eastern and South-Eastern Asia and 14% from Oceania) are the more geographically diverse whereas Elements (75% from USA and 25% from Western Europe), Journal of Mineralogical and Petrological Sciences (79% from Japan) and Swiss Journal of Geosciences (79% from Switzerland) are more limited in geographic representation among editors. Such differences, in spite of ongoing efforts by participating societies to progress scientific excellence via improved diversity and inclusion, may be explained by unconscious bias arising from legacy influences and the current management of expectations for editorial roles. Most journals are published in English and editors are expected to be native or fluent in English, therefore, most editorial boards would consist of individuals who originate from countries where English is an official language (e.g., United Kingdom). Moreover, the regional scope of some journals like the Swiss Journal of Geosciences (readership from Switzerland) and Journal of Mineralogical and Petrological Sciences (readership from Japan) can also explain a distribution biased towards a specific region.

Overall, editorial board members are predominantly from Western Europe 39% , Northern America 29%, Eastern and South-Eastern Asia 16% and Oceania 5%. If we look at gender distribution among regional groupings, Latin America and the Caribbean editorial board members are 46% women and 54% men (but represent less than 2% of the total) whereas Northern America and Western Europe have a 23% to 77% and 26% to 74% distribution, respectively. Eastern and South-Eastern Asia distribution is far less balanced Triple Point section -Elements 7 (5% of women and 95% of men). The observed differences in gender distribution likely reflect regional or local progress towards gender equality [START_REF] Hori | Progress and problems of gender equality in Japanese academics and geosciences[END_REF]. The results presented in this article show the same and persisting imbalance of diversity in the editorial boards compared to the research community as a whole as was Triple Point section -Elements 8 identified by [START_REF] Mukasa | Under-representation of women and minority awardees in geoscience societies[END_REF] more than a decade ago. This inequitable representation is exacerbated because our community demographics do not demonstrate gender or racial parity [START_REF] Bernard | No progress on diversity in 40 years[END_REF][START_REF] Pourret | Diversity, Equity, and Inclusion: tackling under-representation and recognition of talents in Geochemistry and Cosmochemistry[END_REF] and women are also underrepresented as first authors relative to their representation in the field of geoscience [START_REF] Pico | First Authorship Gender Gap in the Geosciences[END_REF]. Scientific excellence may suffer as a consequence of this imbalance. The reality is that until there are more women and under-represented groups in the field (which is happening, albeit slowly), or there is an equity lens used for assigning workloads that reflects different service loads, early, mid-career and senior scientists from these under-represented groups are being asked to carry heavy loads in order to fill the diversity quotas for committees and editorial boards. It is not uncommon for invitations to serve on an editorial board to be rejected by white women and scientists of other under-represented groups because they are already too busy with other service activities that may provide more immediate professional benefit.

To improve scientific excellence and diversity, journals could implement the following: (i) Set up a diversity working group that can help identify potential qualified editorial board members and Editors-in-Chief, while also targeting an increase in diversity.

(ii) Editors-in-Chief may invite identified people (see previous point) to join their editorial board when a position is available (no necessity of expansion, but expansion may accelerate the changes). They should emphasise their results and efforts toward diversity/equity/inclusion at the journal's society meetings to educate members about editorship diversity.

(iii) Individual editorial board members may give personal encouragement, and act as mentors, to potential editors who are women or of diverse backgrounds.

(iv) Ask scientists from under-represented groups how they could be supported in order to participate in an editorial board, for example through workloads that are set using principles of equity.

(v) Journals are suggested to present an infographic of diversity of the editorial board and/or the geographical scope of the published articles. This may attract the attention from diverse researchers, as well as raising the awareness of diversity/equity/inclusion in the scientific publishing space.

One-time actions to tackle diversity are not enough. Journals must monitor the impact of new diversity efforts to ensure real change is happening on their boards.

Recent shifts towards more equal gender and geographic representation for Elements, Geochemical Perspectives Letters or Geochimica et Cosmochimica Acta are very encouraging, but this is not the case for all participating societies journals. Achieving representative diversity on editorial boards needs sustained effort. Further, we recommend that the editorial boards of the journals published by participating societies linked to Elements consider prioritizing and establishing a mentoring approach to address negative and unconstructive critique for articles. Additionally, minimizing barriers to publishing is particularly important now, given the unequal impacts of the COVID-19 pandemic on submissions by men and women and geographical location.

Figure 2

 2 Figure 2 Geographic distribution of editorial board members by journal (regional grouping