Gilles Perrot
email: gilles.perrot@univ-fcomte.fr

Stéphane Domas

Raphael Couturier

Raphaël Couturier

How separable median

Keywords: median, filter, separable, GPU Graphical Processing Unit 2 NPP: Nvidia Performance Primitives 3 Compute Unified Device Architecture, Nvidia's programming model. 4 Parallel Ccdf-based Median Filter (Ccdf: Cumulative Complementary Distribution Function). 5 Parallel Register-only Median Filter

Introduction

Median filtering is a mathematically simple technique, first introduced by Tukey in [START_REF] Tukey | Exploratory Data Analysis[END_REF]. Applying a median filter to an input image consists in replacing each pixel value by the median value of its neighbors. Most of the time, the neighbors belong to a n = k × k window around the working pixel. Median filtering has been widely studied since its introduction and several median selection methods have been proposed. These methods can be classified into two main classes: histogram-based and sorting-based algorithms. Histogram-based solutions work by building the histogram of the gray levels of all the pixels inside the filtering window. Then, starting from either the top or the bottom class, the number of pixels in each gray level class are summed up until half the total pixel count of the windows is reached. The current class holds the median value. Sorting-based solutions simply sort the gray level values of the filtering window and select the one located at the middle index in the sorted vector.

Originally, the computational cost, the data-dependent execution time and most of all the computational complexity were seen as drawbacks. Researchers have since addressed these issues and designed many implementations, among which efficient histogram-based median filters featuring predictable runtimes [START_REF] Huang | Two-Dimensional Digital Signal Processing II: Transforms and Median Filters[END_REF][START_REF] Weiss | Fast median and bilateral filtering[END_REF] and even a constant-time algorithm proposed in [START_REF] Perreault | Median filtering in constant time[END_REF], which exploits histogram properties to reduce the computational complexity to a single suppress and add operation for each consecutive window position. More recently, authors have managed to take advantage of the newly opened perspectives offered by modern GPUs to develop CUDA 3 -based filters such as the Branchless Vectorized Median filter (BVM) [START_REF] Chen | High performance median filtering using commodity graphics hardware[END_REF][START_REF] Kachelriess | Branchless vectorized median filtering[END_REF] featuring a significant runtime improvement, along with the histogram-based PCMF 4 median filter [START_REF] Sanchez | Highly parallelable bidimensional median filter for modern parallel programming models[END_REF]. Lately, in a previous work, the sorting-based PRMF 5 implementation [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF] has been proposed and remains the fastest method to date for 3×3 and 5×5 windows. Since then, other solutions have been proposed, that allow higher speed for larger kernels, as it is the case for recent works proposed in [START_REF] Green | Efficient Scalable Median Filtering Using Histogram-Based Operations[END_REF] or [START_REF] Battiato | High Performance Median Filtering Algorithm Based on NVIDIA GPU Computing[END_REF], adapted from the original [START_REF] Perreault | Median filtering in constant time[END_REF]. PRMF uses simplified sorting networks to iteratively eliminate min and max values. Recently in [START_REF] Adams | Fast Median Filters Using Separable Sorting Networks[END_REF], the author worked on the separability of sorting networks and obtained significant speedups, reaching almost the highest throughput of PRMF for 3×3 and 5×5 filters with better scalability for larger sizes.

Median filtering is mainly efficient in reducing noises whose Probability Density Functions (PDFs) show high standard deviation values, which amounts to saying that it is a robust estimator of the average value. It is this very property of robustness that makes median filtering so interesting in signal processing, as it can ignore extreme or aberrant values present in the signal and select the most representative value inside the filtering window. Two-dimensional median filters are derived from a non-linear and non-separable mathematical operator but, as shown by Tukey in [START_REF] Tukey | The ninther, a technique for low-effort robust (resistant) location in large samples[END_REF], the median of medians 6 obtained by two computationally simpler stages of median filtering remains as robust an estimator of the mean value as the exact median and can replace it advantageously, mostly for runtime reasons.

Some authors have already studied the smoothing performances of both full and separable median filters, as in [START_REF] Narendra | A separable median filter for image noise smoothing[END_REF] where coarse approximations are given for the output variances under both uniform and Gaussian noise PDFs assumptions.

In the field of noise reduction in digital images, Tukey's observations can be adapted with small datasets built from the pixels of the sliding window. The first computational stage would then provide one median value for each of the k rows of the window; the second stage would issue the median value of those k median values. However, the median of median has not been widely adopted as a noise reduction technique despite the increase in image sizes and the need for larger filtering windows. The reason is probably its higher output variance value. Instead, researchers have mainly kept trying to compute the full 2D median filter at high speed rates on larger filtering windows.

On GPU architecture, some histogram-based algorithms are quite efficiently implemented and show almost constant speed with respect to the filtering window size, but the memory size required by the histogram forbids to achieve near to peak throughput values. Sorting-based algorithms proved highly efficient for small window sizes, but their complexity does not allow an extended scalability. The fastest implementations of both categories [START_REF] Green | Efficient Scalable Median Filtering Using Histogram-Based Operations[END_REF] and [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF] define a threshold around 13×13, where the histogram-based implementation becomes faster than the sorting-based one.

The contribution of this paper is summarized in this paragraph. First, it can be observed that a separable median filter is most of the time stronger than a full median filter. Second, the output variance difference between both versions of the median filter decreases very abruptly when the window size increases. Third, a study is proposed to show that the noise types which are best reduced by median filters have spiky looking PDFs, such as salt and pepper noise. However, this does not lessen the interest of median filtering for Gaussian noise reduction. Finally, a very efficient implementation on GPU is described and throughputs are given on a Volta V100.

In the following, the median of medians two-stage process is referred to as the "separable" median filter for simplicity and clarity's sake, as opposed to the "full" median filter. In addition, it will be shown that the separable filter is almost always a better choice. For this purpose, section 2 first details the comparison between full and separable median filters. Then, section 3 presents a reminder of the median value selection method used in the PRMF and how is has been efficiently implemented on GPU within the proposed 2D separable median filter. Section 4 then presents the two experimental setups used to evaluate the performances of the proposed implementation and to conduct the statistical comparison of both median filter versions. Finally, the execution time measurements and numerical properties of the compared filters are exposed in section 5 before the conclusion drawn in section 6.

2 Theoretical comparison of separable and full median filter under salt and pepper noise This section proposes a theoretical study that demonstrates the conditions under which the separable filter is stronger than the full version on bi-chromatic images. It exposes a threshold for window size above which a separable filter corrupts less pixels than the separable version, while the number of denoised pixels tends to be similar.

Salt and pepper noise is caused by bad pixels in image sensors or by transmission errors. Its name refers to the visual aspect of the resulting corruption: black and white pixels scattered over the image.

The power of salt and pepper noise is represented by the probability P of each pixel to be corrupted. For bpp 8 gray level images, the PDF can be characterized by the expressions below where p(v | u) is the probability to observe the gray-level value v knowing the ground-truth value u. These expressions are derived from the classical equations of the model, whose description can be found, for example in [START_REF] Ng | A switching median filter with boundary discriminative noise detection for extremely corrupted images[END_REF].

p(v | u) =                          P 2 + (1 -P ) if v = 0 and u = 0 P 2 + (1 -P ) if v = 255 and u = 255 P 2 if v = 0 and u = 0 P 2 if v = 255 and u = 255 (1 -P ) if v = u and u / ∈ {0, 255} 0 elsewhere ∀(v, u) ∈ [0; 255] 2 (1)
To establish which filter provides the strongest denoising, let us assume that the original image is noiseless (none of its pixels is white or black) and that it is corrupted by a salt and pepper noise of probability P , which means that all black or white pixels are noise. Such a case is easily simulated by scaling the gray levels and does not cause too much distortion, especially when dealing with bpp 12 or higher-depth images.

Noisy pixels removed by full 2D median filter

If glmax is the highest gray-level value allowed by the pixel format, the probability for one pixel x of the input image to be corrupted is P (x = 0 or x = glmax) = λ. The full median filtering window is k ×k pixels wide and contains a total amount of n = k 2 pixels. At each position of the filtering window, if one of the extrema values (black or white) occurs at least k 2 /2 + 1 times, the output pixel y will be of that same extremum gray level, and will still remain a noisy pixel. The probability of the output pixel y being left noisy, denoted P f (yn) is then given by the following equation:

P f (yn) = 2 n q= k 2 /2+1 C q n λ 2 q 1 - λ 2 n-q
(2)

Noisy pixels removed by separable median filtering

As for the separable filter, similar considerations still apply, but two 1D filters, 1 × k, then k × 1 are processing each image consecutively.

Here, the threshold is k 2 + 1 and the probability P sh (x h n) of x h being a noisy pixel after the first pass (horizontal by assumption) and before the second pass of 1D median filtering is given by the following expression, analogous to equation 2

P sh (x h n) = 2 k q= k 2 C q k λ 2 q 1 - λ 2 k-q
and the probability of y being noisy is then

P sv (yn) = 2 k q= k 2 C q k P sh (x h n) 2 q 1 - P sh (x h n) 2 k-q (3)
Fig. 1 displays the difference P sv -P f when the noise level ranges from 7% to 50% and filtering window sizes from 3 × 3 to 11 × 11. For larger filtering windows, the difference becomes insignificant compared to the other parameters that will be discussed later on in this paper. For example at a 50% level, the probability of both 11 × 11 median filters to output a noisy pixel at a given position differs by less than 2.6 × 10 -4 point. It becomes obvious that, the number of noisy pixels left by both filters rapidly tends towards extremely close values.

However, median filters do not just alter noisy pixels, but are likely to alter other pixels in the input image, thus generating unwanted corruption (noise).

Actually, computing the number of pixels which would be corrupted through this mechanism, inside a natural image, is far beyond standard computing capabilities. Nevertheless, in order to begin to understand the distortion brought by each of the full and separable median filters, let us make the strong hypothesis of a bi-chromatic image composed of a background color c b and some objects or shapes of color c a . Three illustrative examples of such images are provided at Fig. 2.

Fig. 1 Difference between the probabilities of a pixel being left noisy by full and separable median filters for filtering windows ranging from 3 × 3 to 11 × 11 in presence of salt and pepper noise of power 7% to 50%.

Pixels corrupted by full median filters

One pixel is altered by k × k full median filtering if there are at least k 2 2 pixels of the other color inside the filtering window (and at most n -1).

Consequently, a central pixel of color c a is altered if at least k 2 2 pixels are of color c b inside the k × k filtering window. This color balance is provided by a number s a of possible combinations, given by

s a = n-1 q= k 2 2 C q n-1
The combination count s b that leads to altering a central pixel of color c b follows the same expression and the sum s f ull = s a + s b represents the total number of color combinations where the central pixel will be altered by median filtering. The following equation describes s f ull .

s f ull = 2 n-1 -C r-2 n-1 (4)

Pixels corrupted by separable median filtering

A similar path leads to the count of combinations s sep in the separable case. Assuming that the first 1D filtering is achieved through a 1 × k horizontal window, and that the color of the central pixel P is c a , this pixel will be altered if at least k 2 pixels are c b -colored inside the k×1 vertical window of the second The same applies to c b -colored central pixels and the total combination count s sep is obtained by adding up both subtotals. That leads to the expression of equation 5, which obviously applies whatever the direction of the first filtering step.

s sep = k l= k/2 Q l-1 k C l k (2 k -Q k ) k-l-1 k × R k .l(2 k -Q k ) + Q k (k -l)(2 k-1 -R k )
where

Q k = k q= k/2 C q k R k = k-1 q= k/2 C q k ( 5
)
Fig. 3 displays the ratio s sep /s f ull for filtering windows under 13×13 and shows that both counts are very close and are asymptotically tending to identical values when k exceeds 11.

Judging from the above comparisons, it seems clear that in presence of salt and pepper noise, full median filtering removes a higher number of noisy pixels and also corrupts pixels in fewer situations than separable filtering.

However, a more in-depth analysis reveals that, inside natural images, not all color combinations resulting in the corruption of a given pixel have the same probability to occur. Unfortunately, these probabilities cannot be formally expressed as they depend on each image content. The combinatorial complexity also forbids a numerical evaluation. For example, the smallest filtering window sizes (3×3 and 5×5) would lead to analyse respectively more than 4 × 10 21 and 1.6 × 10 60 combinations.

When the separable median outputs better images than its full version.

For k = 3, the filtering window includes k 2 = 9 pixels, each of them colored with one of the two possible gray levels c a and c b . Thus, it allows a binary representation of the entire filtering window, with the first row holding the k most significant bits and the last row holding the k least significant ones. Among the 2 k 2 = 512 possible binary combinations offered by this representation, separable median will add noise by altering the center pixel (ie. changing its color) for 192 of these combinations and only 186 for the full filter. Moreover, there are 30 combinations where the separable filter alters the central pixel and the full filter does not, while the opposite is observed in only 24 situations. Fig. 4 details all these combinations where only one of the filters alters the central pixel.

Nevertheless, the above figures do not take into account the probability of each above combinations to appear in natural images. At a first glance, one can notice that most of the situations of Fig. 4b seem likely to be found in natural images, as they can represent pieces of simple continuous contours or segments. This is the case for 22 out of the 24 combinations detailed in 4b, only numbers 460 and 481 are not continuous. The situations displayed in Fig. 4a describe more scattered pixels less likely to be found in natural images. Indeed only 6 out of the 30 combinations of Fig. 4a could represent a piece of simple continuous contour or segment (numbers 424, 232, 124, 46, 313 and 43). These observations can lead to the assumption that the full filter alters more than the separable one when processing natural images. This is confirmed by parsing the three images of Fig. 2 and counting the occurrences of each combination. The results are shown in Table 1. Table 1 Total count of pixels in the bi-chromatic image that are altered by only one of the 3 × 3 median filters (full or separable).

The full filter removes a higher number of noisy pixels than the separable one but the gap decreases very rapidly when k increases. Meanwhile, the full filter seems to corrupt images more than the separable version. These two remarks tend to imply that for any given power of noise, there is a threshold value for k above which separable median filtering is better than full median filtering. This is confirmed by the measurements conducted on grayscale images detailed in sections 4 and 5.

Implementing a fast 2D median filter in a separable way

In an earlier work, the 2D full-median filter implementation called PRMF [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF] was proposed. Though extremely fast for small filtering windows, it is outperformed for large window sizes by several other implementations like [START_REF] Chen | High performance median filtering using commodity graphics hardware[END_REF], [START_REF] Sanchez | Highly parallelable bidimensional median filter for modern parallel programming models[END_REF], [START_REF] Perreault | Median filtering in constant time[END_REF] and [START_REF] Green | Efficient Scalable Median Filtering Using Histogram-Based Operations[END_REF]. However, section 2 demonstrates that the separable median filter is as robust as the full version, and could output less noisy images in a wide range of situations, notably for large window sizes. This is because the computational complexity of the full median is higher than that of separable median, so that the expected runtimes of the latter are significantly faster.

The main issue to achieve good performances is the choice and implementation of the selection method used to identify the median value among a collection of gray levels. A second key issue is how to rule redundancy between consecutive positions of the sliding window, since two neighboring pixels share some of the values to be sorted, as shown in Fig. 5. The separable implementation proposed in this paper uses the same optimizations for the selection process as in PRMF [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF]. It consists in a forgetful selection algorithm using thread registers and exploiting the Instruction Level Parallelism (ILP) capability of the GPU. These principles are exposed in sections 3.1 and 3.2. Compared to the full version, the way of ruling the redundancy has been extended and is discussed in section 3.3. It is also the case for the strategies to manage the global memory used to store input, output and intermediate images. They are presented in section 3.4.

Using registers

As register access is at least 20 times faster than all the other memory types available on the GPU, it is natural to try to use thread registers as a means to store temporary data inside our kernels, keeping in mind that from Kepler to Volta architectures [START_REF]Nividia Kepler architecture[END_REF]-[18], each individual thread can use a maximum of 255 registers within the limit of 32K or 64K per thread block, depending on the GPU family. However, a high register usage, even if below the above-mentioned limitations, may result in a loss of performance due to a lower parallelism level inside each block, i.e. less threads actually run in parallel. Consequently, registers need to be used as sparingly as possible, in order to preserve high pixel throughput values. To do so, the forgetful selection algorithm already described in [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF] has been implemented. Its principle is to build an initial list of R k pixels values taken among the k pixels inside the 1D window (horizontal or vertical), then to identify and eliminate (forget) both elements showing the maximum and the minimum values in the list. Finally, one of the values left apart of the original list is included in the current list. This process is repeated until no more values can be included in the list. The remaining element in the list is the global median value. It is important to notice that this algorithm has a fixed and known number of steps, equal to k -k 2 , which implies that all threads have almost the same workload, despite the data dependency of the extrema identification step.

The number R k of elements (and thus registers) in the initial list is chosen as the minimum element count which allows to identify the global median value through the above process. It is obtained by considering the constraint of keeping the global median in the list at each elimination step. This leads to:

R k = k 2 + 1
It is also noticeable that each elimination step uses one register less than the previous one as two elements are eliminated and one element added at each step. Fig. 6a illustrates the forgetful selection process applied to a 7-pixel 1D median filter. The process begins with R 7 = 5 elements and ends after 3 iterations, when there is no more candidate element to add to the current list. This also corresponds to the state where there is only one element left in the list: the median value. The selection of both extrema is similar to the first iteration of a bubble sort. It uses a basic 2-element swapping function to push the minimum (resp. the maximum) value to the beginning (resp. the end) of the list. The implementation of this selection process ensures that the GPU kernel code is free of divergent branches liable to severely impact performances.

Hiding Latencies

Optimizing a GPU kernel also means hiding latencies potentially generated by memory accesses and data dependent instruction calls. Indeed, modern GPUs are able to pipeline instruction processes so as to reduce the average latency of an instruction sequence. This capability is called ILP (Instruction Level Parallelism). As for global memory accesses, when two or more consecutive arithmetic operators manipulate (read or write) independent variables, only the first access generates latency. The massive thread parallelism of CUDAenabled devices helps in hiding some of those latencies transparently. In our parallel code, an analysis of the actual computation performed by each thread, shows that ILP has to be optimized: First, the Instruction Level Parallelism can be increased within the forgetful selection method by re-arranging the instruction sequence as it would be done in a sorting network [START_REF] Batcher | Sorting networks and their applications[END_REF], so as to reduce the data dependency of consecutive instructions and thus to prevent frequent empty pipelines. Fig. 6b shows the scheduling of the first extrema identification step of a 21 × 1 median filter, carried out with R 21 = 12 elements. Each arrow represents one call to the 2-element swapping function: after the call, the starting point symbolizes the lowest value element and the ending arrow points out the highest one. In addition, horizontal dashed lines separate packs of instructions, where the first packs contain instructions with independent variables and the last pack contains the swapping instructions with the largest strides. This last pack may contain some interdependent instructions, depending on the size of the array. After the execution of each pack, the current state of the array is also displayed to help in following the individual moves of each value. In the example of Fig. 6b, 3 packs of independent instructions are defined, and the fourth contains interdependent ones. Considering the whole sequence, the schedule of each instruction is adjusted to maximize the distance (in terms of line of code) between interdependent instructions. It allows to maximize the ILP. In this example, the last pack cannot be optimized, but it is negligible compared to a naive implementation that would not have achieved any ILP at all.

Exploiting the window overlapping

The goal here is to increase computations inside the GPU threads, as the performances of the median filter are limited by the cost of memory accesses. Increasing the computation ratio must be done without impairing the possible parallelism level, i.e. by avoiding to process too many redundant operations at thread level, thus multiplying the registers use. The following design intends to reduce the registers count per thread while increasing the threads computation load. It also allows to reduce the effect of global memory access latencies.

Since each thread processes m input pixels, it has to carry out m selections. Thus, the window overlapping cannot be exploited in the same way as it is achieved within histogram-based solutions. Nevertheless, the window overlapping implies that these m selections are not independent. Thus, an efficient implementation must avoid redundant extrema identifications. This is achieved by choosing the R k elements of the first selection step among those shared by all windows. It only makes sense if the consecutive windows actually share at least R k elements. For example, if m = 2, k must be greater or equal to 5. The first selection steps can then be considered common to all windows. The following selection steps apply on the remaining pixels, which are shared by fewer windows.

This principle was already exploited in the 2D PRMF design, but was limited to m = 2. Indeed, when two 2D consecutive windows overlap, a total of 2 × k pixels are not shared and have to be processed separately, using two Algorithm 1: Code generator for the forgetful selection. The generator actually writes the required number of #define minMax i macros inside the CUDA code, in order to achieve the forgetful selection process. Each minMax i function works on a set of i single variables a m stored in registers. During execution of each minMax i function, a swapping function s(a m , a n ) is called multiple time to rearrange a m and a n in ascending order whenever necessary. input : r, the radius of the (2r + 1) × (2r + 1) pixels filtering window input : (a 0 ...a P ) the set of pixel values stored in registers output: The set of #define minMax i definitions 1 kmin ← r + 2; 2 for i ← kmin to 3 do output: "#define minMax i(a 0 , ..., a (i-1) ) " // Init sorting network loop different groups of registers. For example, if k = 5, there are 20 pixels shared by the two 2D windows. The selection starts using R 5×5 = 14 pixels, and thus registers. After seven selection steps, all shared pixels have been processed and there are only 6 registers that contain relevant values for the next steps. Their values are copied to a new group of 7 registers that is used to process the selection for the second window, while the original group of registers is used for the first window. In total, the process uses 14 + 7 = 21 registers, instead of 14 + 14 = 28 if overlapping is not exploited.

Nevertheless, when processing two 1D windows, using the same k as in 2D, the computation load of a thread is very low. For example, if k = 5, a thread processes only 6 pixels, compared to 30 in 2D. Since R 5 = 4, our selection process uses 7 registers instead of 8 if both windows are processed independently. But even if there is a gain, experiments have demonstrated that it is not sufficient to compensate the weak load of each thread. This is why the method needed to be extended to larger pixel packets for 1D filters. Fig. 7 details the extended method in the case of a 9 × 1 horizontal median processing four pixels per thread (m = 4). A solution that does not take the overlapping into account would use m×R 9 = 4×6 = 24 registers. Our method drastically reduces this number. Each of the m center pixels is displayed in a specific color and is associated to a group of registers printed in the same color. These groups are numbered from 0 (the blue one) to 3 (the yellow one), and registers of group i are named x i * . The numbers in the square pixels are sample gray-level values. The first stage consists in devoting group 0 (blue) for each of the R 9 = 6 pixels shared by all filtering windows. The first minM ax[R k ] (here minMax6) selection stage is then achieved, which frees two registers (the min and the max). The four remaining relevant registers of group 0 are copied to group 1 (red), dedicated to the second center pixel. Both groups are next augmented with one pixel value chosen among the remaining pixels, starting from the innermost ones. Now two R k -1 (here minMax5) selection stages are achieved. This duplication-augmentation-minMax process is then repeated. At each iteration i, group i is created and filled by copy from the relevant registers of group i -1, then for each group, one of the remaining pixels is integrated, before executing the minMax. It goes until four minMax3 selection steps output the four median values. In this example, groups have respectively a size of 6, 5, 4 and 3, leading to a total of 18 registers, instead of 24 for a naive implementation. For m = 4, the number of registers used by each thread is then given by (R k + 12), while it is (R k + 42) for m = 8 (provided k is large enough). It is worth noting that since each thread uses more registers than with m = 1, it may lead to overpass the limit of the register count per thread block. Nevertheless, it is easy to stay under this limit by dividing the block size by m, while preserving the grid size.

The GPU specific access patterns required to ensure coalescence to and from global memory lead to choose packet sizes m as powers of two. The maximum packet size m max is then the greater power of two that remains less than k/2. Above k = 65, all packet sizes can be chosen between 2 and 32.

The case k = 3 cannot be treated in the same way, as no overlapping can be exploited at this size. However, devoting one thread to processing only one pixel does not provide a satisfying level of performance, because of unbalanced computations and communications (memory accesses). Instead, experiments show that the highest throughput is obtained by processing m = 4 pixels per thread, without taking the overlapping into account. It appears that it is also the case for 5 × 5 and 7 × 7, even though sufficient overlapping exists at these particular sizes. That confirms that computation load has a higher impact than the parallelism level, which is not a limiting factor in our design, for small sizes.

Optimizing the sequence of both 1D filtering stages

As already stated, the choice of running the horizontal 1D step first is arbitrary and does not impact the filtering performances. Separable implementations (not only median operators) frequently adopt a symmetrical architecture that allows the same function to execute twice, for both vertical and horizontal stages, separated by an image transposition. The proposed implementation does not follow this scheme and uses one specific GPU kernel for each 1D stage. This is imposed by the memory allocations and alignment constraints that must be fulfilled in order to achieve optimal accesses to and from the GPU global memory.

Three global memory areas are allocated onto the GPU before calling the filtering kernels: d in, d outh and d outv. Assuming an input image of M rows and N columns, stored row after row :

d in contains the input image. It has the same height as M , but is larger.

In order to avoid special processing, the input image must be horizontally zero-padded by adding k 2 columns on both the right and the left sides. To limit non coalescent accesses to global memory, the first column of d in must be 128-Byte aligned. Under these conditions, the width of d in is the closest multiple of 128 above N + k -1. It is filled with zeroes and the input image is stored starting at row 0 and column k 2 . d outh receives the intermediate image output by the horizontal kernel that will then become the input of the vertical kernel. For the same reasons as for d in, the output needs to be vertically zero-padded by adding k 2 rows on both bottom and top sides and the first column of d outh must be 128-Byte aligned. Thus, the height of d outh is M + k -1 and its width is the closest multiple of 128 above N . It is filled with zeroes and the intermediate image is stored starting at row k 2 and column 0. d outv receives the output image. It only requires a 128-Byte alignment.

The first pixel column of the image is written in column 0.

To summarize, the optimization techniques described above allow each thread to use fewer registers while processing more pixels, thus increasing Fig. 7 Reducing register count in a 9×1 register-only median kernel by processing 4 input pixels in parallel (in blue, red, green and orange). The 4 filtering windows are numbered from 1 to 4 and contain gray-level values printed in white. A set of 4 groups of registers x ij is defined by its color and by the number of the associated filtering window as i index. The first forgetful selection step is processed once for the 4 filtering windows, as R 9 = 6 is the exact count of pixels that are shared by all of the 4 windows. During each following step, the leftmost branch of register groups follows the standard forgetful selection process until the median value is output. At the same time, the rightmost branch of register groups is divided after each selection step in order to devote, at the end, one group to each filtering window. In this example, only 18 GPU registers are used by each thread to process four 9 × 1 median values.

the level of parallelism and allowing larger filtering windows without suffering a drastic drop in performances, which occurs when the register file limit is reached. In addition, global memory accesses are coalescent and cache misses are minimized. All of these properties contribute to providing a high level of performance that is demonstrated in the next sections.

Experimental setup

This section describes the software and hardware context to evaluate the applicability of the results of Section 2 to grayscale images, and the performances of our GPU implementation presented in Section 3.

Context of the statistical comparison of full and separable median filters

The approach has been to conduct a statistical comparison over a large number of natural images and odd window filtering sizes. The chosen images are a collection of 10,000 gray-level 8-bit depth Portable Gray Map natural pictures from the BOSSRank database [START_REF]Cover database used to generate the Testing database[END_REF], whose original purpose was to benchmark steganalysers. These images, called cover images are meant to represent original images, i.e. images with no significant and/or visible noise. During the simulations, cover images have been corrupted by salt and pepper noise with different strengths: 7%, 10%, 25% and 50%, which cover a wide range from light to strong power of noise.

Measurements have been conducted to compare the output of both full and separable median filters for filtering sizes ranging from 3×3 to 131×131. For each pair of filtered images, the difference between output images has then been measured through PSNR 7 and MSSIM8 [START_REF] Zhou Wang | Image quality assessment: from error visibility to structural similarity[END_REF] indicators. Absolute values have not been displayed as they have been detailed in the literature during the past decades and would not provide any additional information.

For any given image I in (i, nl) of size H × W pixels, where i ∈ [1 : 10, 000] and the noise level nl ∈ [0; 1], the output images I f ull and I sep are obtained respectively by filtering I in (i, nl) with a k ×k median filter in its full or separable version. For each given set of parameters, the two main indicators ∆ P SN R and ∆ M SSIM are obtained by evaluating the following expressions:

∆ P SN R = P SN R(I f ull , I in ) -P SN R(I sep , I in ) ∆ M SSIM = M SSIM (I f ull , I in ) -M SSIM (I sep , I in ) (6)
In addition to the average ∆ P SN R and ∆ M SSIM values, each pair of output images has been analysed to count how many images are best denoised respectively by the full and the separable median filter. This led to produce a supplementary indicator that helps in deciding which filter version to use, regarding the window size as well as the type and power of the noise. This indicator is simply displayed in the various figures as a label with either of the following texts: "full likely better" or "separable likely better".

Eventually, another batch of measurements has been conducted in order to confirm what the study of the bi-chromatic case has revealed: there are more color combinations where only the separable filter would alter the central pixel, but they are less likely, inside the filtering window, than those where the full would be the only one to alter the central pixel. To avoid mixing noisy and non-noisy pixels, both median filters have been applied on every cover images of the database and the two interesting situations described above have been counted.

To avoid any doubt regarding the simulation results, all measurements have been performed under Matlab R2017b, using standard functions like med-filt2(), psnr() and the original mssim implementation proposed in [START_REF] Zhou Wang | Image quality assessment: from error visibility to structural similarity[END_REF].

Evaluation context

Results have been obtained by averaging runtimes of a variable number of executions. The actual number of runtimes was dynamically set to ensure a total execution time under 0.5 second. This setup avoids any overhead caused by the RCU scheduler on Linux 3.16.0-4-amd64 host platform powered by one Xeon E5620 @2.40 GHz processor and CUDA v9.0. Each kernel has been run by Titan-X9 and V10010 GPUs on 8 bit images of sizes 1,024×1,024 and 8,192×8,192. Cover images of the BOSS database, whose original size is 512×512, have then been resized and corrupted by a medium power salt and pepper noise of 25%.

The pixel throughput value of the GPU kernels was used as a main performance indicator. It does not include transfer times to and from the GPU, as median filter kernels are likely to be part of a complex process fully run on the GPU.

To evaluate the absolute performance of the proposed implementation, the maximum effective pixel throughput that our own GPU/host couple is able to achieve has been measured. Such a peak value helps in deciding on further investigation. We performed this measurement by running dummy kernels that fetch the gray-level of each pixel from GPU's global memory and output it into another area into the GPU's global memory, exclusive of any other instruction. Each thread may process several pixels in order to balance the compute and read/write loads.

An experimental tuning stage shows that processing four pixels by thread achieves the peak pixel throughput of respectively 80 (Titan-X) and 200 (V100) billion pixels per second on an 8192×8192 pixel image. The most significant factor in the speedup brought by V100 is its memory bandwidth of 900 Go/s, compared to 336 Go/s for Titan-X. This speedup of ×2.5 achieved by the Volta V100 against the Titan-X is observed whatever the median kernels executed on both GPU models. That led to display only the results of the most recent V100 GPU model. If needed, readers can easily deduce Titan-X throughputs values by dividing those of V100 by 2.5.

Results

The first part of this section presents the denoising quality obtained for the considered image set, using PSNR and MSSIM metrics. It is followed by a comparison on the pixel corruption generated by both versions. Finally, an analysis of the computational complexity and performance is given.

Comparison of full and separable median filters under salt and pepper noise

Fig. 8 displays the evolution of the average ∆ P SN R and ∆ M SSIM values when k ranges from 3 to 131 on the whole set of images. It should be noted that on average, whatever the noise level, the separable median filter achieves better PSNR and MSSIM values for every window size above or equal to 11×11, and the average differences tend very rapidly towards the pair of values below :

   lim k→min(H,W ) ∆ P SN R = -0.4 dB lim k→min(H,W ) ∆ M SSIM = -0.005
Another interesting fact is detailed by Fig. 9 which displays the percentage of images, among the whole database, that are best denoised by full median filtering. It reveals that, whatever the power of noise, the separable version gives better results above 21×21, and that below 25% of salt and pepper noise, this threshold is downsized to 7×7.

Given the filtering size and the power of noise, Fig. 9 can then provide a rule for selecting the best type of median filter to use, from a statistical point of view: below 50%, on average, the separable filter is probably the best choice.

It is also important to note that, up to a 7% of salt and pepper noise, separable filtering is always better than full filtering.

On the basis of the above statistical study, and with all the precautions that must be taken regarding the data-dependency of median filtering as well as its intrinsically limited scope (though quite large) provided by the 10,000 natural scenes of the database, the conclusion is that in common situations, separable median filtering is likely to achieve higher salt and pepper noise reduction than full median filtering. As for larger filtering windows often used for other reasons than noise reduction, both filters have a very close behavior which is in favor of the much faster separable version. An example is presented in Fig. 10, where image #9850 of the database, randomly chosen, is displayed in two configurations: after corruption by 7% and 25% salt and pepper noise (images a and b). Images c and d show the results after a full 3×3 median filtering stage, while images e and f display the results of the separable 3×3 filtering stage. The visual comparison confirms the statistical results of Fig. 9 and Fig. 8 with very close output images in the 7%-noise case, while the full filter outperforms the separable one in the 25%-noise case.

To illustrate the discussion about the filtering size, Figure 11 shows the results of the same processes applied on a 8192×8192 image corrupted by a 25% salt and pepper noise. For a better visibility, a small square region of interest is displayed, making it easier to compare both output images. In this configuration, a 15×15 filter is well adapted to noise reduction and both filters visually produce very close results.

Comparison of the corruption generated by separable and full filters

The percentages of non-noisy pixels corrupted by only one of the median filters, separable or full, are plotted on Fig. 12. It confirms that in the range of our experiments, the full filter always alters more pixels than the separable version, though the difference tends to decrease when the window size increases. The maximum difference of 6.68 points is observed for the smallest size while the minimum of 2.64 points is reached at the largest size of the experiment, ie. 131×131. It represents between 17,500 and 6,900 pixels for the 512×512 images of the database.

Computational complexity analysis

As presented above, computations are mostly achieved by means of multiple calls to a simple swapping function, that actually swaps its two input parameters if their values are not in ascending order. A good approach to the complexity is given by the number of times each median filter calls this function. It actually represents the worst case, as each call will not necessarily result in an actual swap.

For each processed pixel, the proposed separable filter realizes respectively 3, 10 and 16 calls for the smallest sizes 3×3, 5×5 and 7×7 because no overlapping is exploited. For sizes ranging from 9×9 to 19×19 the packet size is set to m = 4 and the number of calls is given by:

W sep (k) = 18 + 1 2 i=6+ k-9 2 i=6 sw(i) (7)
where

sw(i) = 3 i-1 2 + 1 , i odd 3 i 2 , i even
Above 21×21, the packet size is set to m = 8 and the number of calls is then given by: Fig. 13 Complexity of both PRMF (full) and PRSMF (separable) implementations, for filtering window sizes ranging from 3×3 to 131×131. The complexity is displayed as the total number of calls to the swapping function, according to real executions.

W sep (k) = 64.25 + 1 4

i=12+ k-21 2 i=10 sw(i) (8)
According to Equations 7, 8 and results from [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF], the complexity of processing a single pixel for the proposed separable median filter is O( 14 klog(k) 2 ), which is close to k 3 , as opposed to O( 1 10 k 4 ) for the PRMF full version. It is illustrated in Fig. 13 where the proposed PRSMF implementation is compared to the PRMF regarding the real (i.e. during executions) total number of calls to the swapping function made by both filters. Except for a flat area in the PRSMF results, curves perfectly match the given complexities.

Performances of the proposed implementation

Kernel throughput values of the proposed implementation called PRSMF 11 are presented in Fig. 14 for 1,024×1,024 and 8,192×8,192 pixel images and compared with Nvidia's NPP, Green's new GPU 2D filter [START_REF] Green | Efficient Scalable Median Filtering Using Histogram-Based Operations[END_REF] and our previous PRMF 2D filter [START_REF] Perrot | Fine-tuned High-speed Implementation of a GPU-based Median Filter[END_REF]. Green's new GPU implementation represents one of the most recent and most effective 2D median filter but does not run in a separable way. As we do not have access to the source code, the timings on V100 have been extrapolated from [START_REF] Green | Efficient Scalable Median Filtering Using Histogram-Based Operations[END_REF], where a comparison on Titan-X against PRMF is detailed; we trusted author's results and used the constant ×2.5 speed factor measured between Titan-X and Volta (see 4.2) to produce a relevant comparison. Nvidia's NPP is a closed code library which appears to be currently the only one that allows to run in the separable way with reasonable throughputs.

The maximum throughputs of PRSMF are reached by the 3 × 3 filter, processing m = 4 pixel per thread, with around 125 billion pixels per second on V100. That represents around 62.5% of the peak throughput (for both GPUs). The optimal value for m depends on the window size k but only takes two different values: m = 4 when k < 21 and m = 8 when k ≥ 21. This behavior confirms that the computation load has a higher impact than the parallelism level for relatively small windows, as the smallest size for which m = 8 filters use less registers per pixel than m = 4 filters is k = 33, according to section 3.3.

The performances of PRSMF are between 9 and 24 times faster than the Nvidia NPP library, which seems to implement two different methods depending on the filtering size: below k = 21, the plot looks like a sorting-based implementation while from k = 31 and above, the plot seems to reveal a histogram-based implementation (near constant time). As for Green's implementation, the achieved throughput is almost constant whatever the filtering and image sizes, while PRSMF's throughput increases with the image size. Still, Green's implementation achieves lower throughputs whatever the filtering size, at least up to k = 131 where the lowest speedup brought by PRSMF is around ×5 and reaches around ×200 as its highest value (3×3, 8,192×8,192 pixel image).

It is important to notice that the performances of these kernels on 16-bit images remain the same as when run on 8-bit images.

Additionally, it is worth noting that the proposed implementation also optimizes data transfers by allocating aligned memory on the GPU side and pinned memory on the CPU side. As a reference value on both Titan-X and V100, 9.0 ms are needed to copy a 8,192×8,192, 8 bit depth image from CPU to GPU's global memory, and 5.2 ms are necessary to copy back the same image to CPU's pinned-memory. The same transfer from GPU to CPU, using a classical pageable-memory allocation on CPU side(malloc) will cost 40 ms on both GPU models. Notice that transfer times vary linearly with pixel depth and are thus doubled for 16-bit images with respect to 8-bit images.

Conclusion

In this paper, extensive study and theoretical considerations have proved that statistically, better 2-D median filtering is obtained through a separable version than classical full 2-D filtering. In addition, a very high speed implementation has been designed, which makes it possible to process up to 50 billion pixels per seconds on a very common Titan-X GPU, which represents for example, more than 6000 Ultra High Definition (or 4K) images per second. An unprecedented throughput of 125 billion pixels per second has been achieved on a high-end V100 GPU model.

Let us also note that the throughput values achieved by 3×3 window filters come quite close to the measured peak effective pixel throughput allowed by our development platforms. Hence, the search for further improvement seems to be conditioned by new GPU architecture capabilities.

Indeed, new perspectives have been opened with the release of Volta and Turing family GPUs and not yet investigated in this work because of its very recent availability. For example, independent thread scheduling and collaborative groups could be promising features to be utilized.

The generalization of the comparison of the separable 2D median filter and the full 2D median in presence of additive white Gaussian noise has already been conducted and confirms the average superiority of the separable median (not discussed here due to the limited number of pages). Likewise, it would be relevant to study the same situation with the 3D median filter and other noise types. Finally, we plan to revisit other non-linear filters with the pseudoseparability in mind.

Fig. 2

 2 Fig. 2 Examples of 512 × 512 bi-chromatic images.

Fig. 3

 3 Fig. 3 Ratio between the combination counts of pixels altered by separable and full median filters, for windows ranging from 3 × 3 to 11 × 11

  (a) Separable filter alters the output pixel and full filter does not. (b) Full filter alters output pixel and separable does not

Fig. 4

 4 Fig. 4 Combinations, inside a bi-chromatic image, where the central pixel is altered by one and only one version of the 3×3 median filter. Due to the symmetry, colors ca and c b can be associated to either the white or black pixels. The number under each combination is the decimal value of the binary representation obtained by concatenating the rows from top to bottom (most significant bits first) and considering that a black pixel is associated to the digit 1.

Fig. 5

 5 Fig.5Overlapping windows (framed in red and blue) of two 5 × 1 median filters applied on adjacent pixels colored in red and blue. Filtering windows share 4 pixels.

  Fig. 6 (a) Determination of the Median value by the forgetful selection process.(b) First extrema identification step of the forgetful selection. It begins at the top row with unsorted elements and ends with the minimum value at the first position (left) and the maximum at the last position (right). The last pack of instructions, in light purple, may contains interdependent instructions, while the others contain only independent ones. The state of the array is displayed at the end of each pack.

Fig. 8 Fig. 9 Fig. 10

 8910 Fig.8Relative performance of the full median filter compared to its separable version, on average, over 10,000 natural images for filtering window sizes ranging from 3×3 to 131×131. A positive value indicates that full filtering is better than the separable one.

(a) 10 Fig. 11

 1011 Fig. 11 Example of a 8192×8192 image, corrupted by a 25% salt and pepper noise (a) and then processed by both full and separable 15×15 filters. a) shows the whole image, b) focus on a 681×681 pixel square around the last letters "US" on the right side. Images c) and are filtered using respectively the full and the separable version. For each image, results are given through PSNR(dB)/MSSIM and are computed on the whole images.

Fig. 12

 12 Fig.12Percentage of non-noisy pixels corrupted by only one of the median filter versions and not the other. The percentages represent the total pixel count, over the 10,000 natural images of BOSS database, for filtering window sizes ranging from 3×3 to 131×131.

  (a) 1024×1024 input image. (b) 8192×8192 input image.

Fig. 14

 14 Fig. 14 Pixel throughput value comparison, in million pixels per second, of Nvidia's NPP, Green's and PRMF implementations against the proposed PRSMF, on Volta V100 GPU.

In that chapter, Tukey adresses the issue of statistical analysis on large datasets and discusses a method in which data are sliced into nine-value sets (or 81 values for larger datasets) and where the mean estimation is taken as the median value of the median values of all those sets, i.e. the ninther. The process can be recursive. In our implementation, data slices are mapped to image rows.

Peak Signal to Noise Ratio.

Mean Structural SIMilarity.

Titan-X GPU: Maxwell family, compute capability 5.2, 3,072 cuda cores at 1.24 GHz, 12 GByte of RAM.

V100 SXM2 GPU: Volta family, compute capability 7.0, 5,120 cuda cores at 1.45 GHz, 32 GByte of RAM.

Acknowlegments

This work was partially supported by the EIPHI Graduate School (contract ANR-17-EURE-0002). Computations have been performed on the supercomputer facilities of the "Mésocentre de Franche-Comté".

Source codes of the proposed implementation can be generated online and downloaded for testing purposes at: http://info.iut-bm.univ-fcomte.fr/staff/perrot/convomed/ In addition, full definition images and more measurements can be found at: http://info.iut-bm.univ-fcomte.fr/staff/perrot/separable-median/ Availability of data and material Source codes of the proposed implementation can be generated online and downloaded for testing purposes at: http://info.iut-bm.univ-fcomte.fr/staff/perrot/convomed/ In addition, full definition images and more measurements can be found at: http://info.iut-bm.univ-fcomte.fr/staff/perrot/separable-median/ Code availability Sample source codes are available via the above URL, but the source code of the kernel generator is beyond the scope of this paper and is right protected.

Declarations

Funding

The Graphical Processing Units used in the described experiments have been provided by the Mésocentre de calcul de Franche-Comté, the regional computing facility.

Conflicts of interest/Competing interests

Authors declare that they do not have no conflict of interest of any sort.


Uncertainty analysis and sensitivity of scapulothoracic joint angles to kinematic model parameters

Keywords: shoulder, scapula, clavicle, joint center location, ellipsoid

The purposes of this study were to determine the influence of kinematic model parameter variability on scapulothoracic angle estimates, and to define which parameters of the kinematic model have the largest effect on scapulothoracic angle estimates. Nominal subject-specific kinematic models of nine participants were implemented. Fifteen parameters of the nominal models relative to the clavicle length, ellipsoid, sternoclavicular and acromioclavicular joint centers, and contact point location were altered from -1 cm to 1 cm. Then scapulothoracic angles were computed during four movements using multibody kinematic optimizations for nominal and altered models. The percentage of scapulothoracic angle variance explained by each parameter of the kinematic model was computed using Effective Algorithm for Computing Global Sensitivity Indices. When altering simultaneously the 15 parameters of the kinematic model, scapulothoracic angles varied up to 50°. For all movements and degrees of freedom, the clavicle length significantly explained the largest part of scapulothoracic angle variance (up to 25%, p<0.01). In conclusion kinematic model parameters need to be estimated accurately to avoid any bias in scapulothoracic angle estimates especially in a clinical context. The present sensitivity analysis may also be used as a benchmark for future works focusing on improving shoulder kinematic models.

Introduction

Scapular positioning and motion are of great interest in patients suffering from shoulder problems. Indeed, an alteration in the normal position or motion of the scapula relative to the thorax, named scapular dyskinesis [START_REF] Kibler | Clinical implications of scapular dyskinesis in shoulder injury: the 2013 consensus statement from the 'Scapular Summit[END_REF], may be the cause or consequence of shoulder troubles [START_REF] Kibler | Scapular dyskinesis and its relation to shoulder injury[END_REF]. In this context scapular motions have to be assessed and quantified by clinicians [START_REF] Carbone | Scapular dyskinesis after Latarjet procedure[END_REF].

However, the accuracy of such clinical assessment is limited by soft tissue artefacts, due to the bones sliding underneath the skin [START_REF] Begon | Multibody Kinematics Optimization for the Estimation of Upper and Lower Limb Human Joint Kinematics: A Systematized Methodological Review[END_REF]. Improvement of procedures to estimate accurately scapulothoracic joint motion therefore remains crucial to help clinicians in their diagnosis.

Among procedures depicted to quantify scapulothoracic joint angles, the multibody kinematic optimization (MKO) would be the most promising approach [START_REF] Duprey | Kinematic models of the upper limb joints for multibody kinematics optimisation: An overview[END_REF]. The later consists of optimizing joint angles by minimizing the distance between the positions of the measured markers and the virtual markers attached the kinematic model. For the scapulothoracic joint, the MKO has been improved, on one hand, by enforcing the scapula to stay in contact with the surface of an ellipsoid representing the sliding plane between the scapula and the thorax, and on the other hand, by ensuring a constant clavicle length to avoid non-physiological joint dislocation [START_REF] Seth | A Biomechanical Model of the Scapulothoracic Joint to Accurately Capture Scapular Kinematics during Shoulder Movements[END_REF][START_REF] Prinold | Scaling and kinematics optimisation of the scapula and thorax in upper limb musculoskeletal models[END_REF][START_REF] Michaud | Scapular kinematic reconstruction -segmental optimization, multibody optimization with openloop or closed-loop chains: which one should be preferred?[END_REF]. These kinematic models are thus defined by several parameters, namely, the clavicle length, the location of the sternoclavicular and acromioclavicular joints, the ellipsoid's dimension and location, and the location of the scapulothoracic contact point. To ensure accurate estimation of scapulothoracic joint angles, these parameters demand also subjectspecific calibration. Although several studies proposed methods to calibrate the kinematic model [START_REF] Michaud | Scapular kinematic reconstruction -segmental optimization, multibody optimization with openloop or closed-loop chains: which one should be preferred?[END_REF][START_REF] Michaud | Determining in vivo sternoclavicular, acromioclavicular and glenohumeral joint centre locations from skin markers, CT-scans and intracortical pins: A comparison study[END_REF], subject-or patient-specific configuration of all the model parameters may be timeconsuming [START_REF] Charbonnier | A patient-specific measurement technique to model shoulder joint kinematics[END_REF] or irradiant [START_REF] Dumas | Subject-specific model-derived kinematics of the shoulder based on skin markers during arm abduction up to 180° -Assessment of 4 gleno-humeral joint models[END_REF] when based on imagery and therefore sometimes hard-to-use in a clinical context. In consequence, screening and quantifying the kinematic model parameters presenting the greatest influence on the scapulothoracic angle estimates may be useful to find the best trade-off between accuracy and efficiency in a clinical context.

To the best of our knowledge, few studies [START_REF] Hybois | Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study[END_REF][START_REF] Habachi | Global sensitivity analysis of the kinematics obtained with a multi-body optimisation using a parallel mechanism of the shoulder[END_REF] assessed the sensitivity of the scapulothoracic angle estimates to the kinematic model parameters. When altering the radii and the location of the ellipsoid center by more or less 20% of a nominal model, scapulothoracic joint angles varied about 15° during wheelchair propulsion [START_REF] Hybois | Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study[END_REF]. During arm elevation in the frontal plane, scapulothoracic angles were the most sensitive to the clavicle length followed by the scapulothoracic contact point location and the center of the ellipsoid [START_REF] Habachi | Global sensitivity analysis of the kinematics obtained with a multi-body optimisation using a parallel mechanism of the shoulder[END_REF]. Although these two studies help to determine in what extend scapulothoracic angle estimates are sensitive to kinematic model parameters, some concerns may be raised. First, only one participant and one movement were investigated while the sensitivity analysis may also depend on individual and movement variability. Second, no study has investigated the sensitivity of the scapulothoracic angles to the location of the sterno-and acromioclavicular joints.

The purpose of this study was twofold, 1) to determine the influence of the kinematic model parameters variability on scapulothoracic angle estimates (uncertainty analysis), and 2) to define which parameters of the kinematic model have the largest effect on scapulothoracic angle estimates (sensitivity analysis). It was hypothesized that scapulothoracic angles were affected by alterations of the kinematic model parameters, especially when the clavicle length changed.

Methods

Participants

Nine healthy males (29.2 ± 8.2 yrs, 1.74 ± 0.08 m, 73.7 ± 9.6 kg) volunteered to participate in this study. All participants declared having no upper-limb and shoulder injury during the year prior to the experiment. The study was approved by the ethical committee Ouest V (2018-A02449-46) according to the statements of ethical principles of the Declaration of Helsinki, and the participants signed an informed consent document.

Data collection and Procedure

Participants were fitted with spherical reflective markers located on the thorax (6 markers), right clavicle (2) and right scapula (8) (Fig. 1). For further calibration of the kinematic model, participants were asked to hold five static poses corresponding to 0%, 25%, 50%, 75% and 100% of maximal arm elevation in the scapular plane (defined by a goniometer as 40° from the frontal plane [START_REF] Ludewig | Motion of the shoulder complex during multiplanar humeral elevation[END_REF]). At each pose, a scapula locator was used to identify the location of the acromion posterior edge, trigonum spinae and scapular inferior angle Thereafter, participants performed five movements, namely, three unilateral arm elevation/depression (one in the frontal, one in the scapular and one in sagittal plane), one horizontal unilateral abduction and one ball throw. The later consisted of performing an overarm throw, at submaximal preferred speed (5.63 ± 1.96 m.s -1 ), of a tennis ball into a 1-meter side length square located at 3 meters in front of the participant. Marker trajectories were recorded using a 10-camera optoelectronic system with a sample frequency of 200 Hz (Qualysis, Sweden).

Nominal kinematic model

The nominal kinematic model was calibrated to each participant's anthropometry. The model was composed of three segments (thorax, clavicle and scapula) for which coordinate systems were implemented regarding the International Society of Biomechanics recommendations [START_REF] Wu | ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion--Part II: shoulder, elbow, wrist and hand[END_REF] and. The coordinates of the scapula center recorded with the scapula locator (i.e., barycenter of the acromion posterior edge, trigonum spinae and scapular inferior angle) were used to calibrate the dimension of the ellipsoid (i.e., three radii), while its center was fixed to the thorax origin reference system (i.e., midpoint between the xiphoid process and the eight's vertebra) [START_REF] Blache | In Press) Does scapulothoracic contact point trajectory follow an ellipsoid surface?[END_REF].

Sternoclavicular and acromioclavicular joint locations were defined by palpation in line with Michaud, et al. [START_REF] Michaud | Determining in vivo sternoclavicular, acromioclavicular and glenohumeral joint centre locations from skin markers, CT-scans and intracortical pins: A comparison study[END_REF] recommendations. Based on Seth, et al. [START_REF] Seth | A Biomechanical Model of the Scapulothoracic Joint to Accurately Capture Scapular Kinematics during Shoulder Movements[END_REF], the model could be actuated by six degrees of freedom (DoF) between the ground and the thorax, and five DoF between the scapula and the thorax (i.e., two DoF for the translation of the scapula center along the ellipsoid surface and three DoF for scapular rotations) (See the mathematical description in supplementary file for more details).

Model-derived joint angles were computed for each movement using MKO subject to an equality constraint ensuring a constant clavicle length. Scapulothoracic angles were extracted following the sequence upward/downward rotation, internal/external rotation and anterior/posterior tilt.

Uncertainty and sensitivity analysis

Both uncertainty and sensitivity analysis were performed regarding Pianosi et al. [START_REF] Pianosi | Sensitivity analysis of environmental models: A systematic review with practical workflow[END_REF] recommendations. The first step consisted of applying a variation from -1 cm to 1 cm on fifteen parameters of the nominal kinematic models relative to the ellipsoid (6 parameters), clavicle [START_REF] Prinold | Scaling and kinematics optimisation of the scapula and thorax in upper limb musculoskeletal models[END_REF] and the scapulothoracic contact point (2) (Fig. 2). To that end, a Latin-hypercube sampling procedure was used by varying all the parameters at a time to obtain 10 000 different values per parameter. Then for each participant, each movement and each altered kinematic model, the MKO was run to compute the new scapulothoracic joint angles. If the MKO did not satisfy the equality constraint (Tolerance = 10 -6 m), the simulation was removed from the analysis.

To assess whether the uncertainty and sensitivity analyses were independent of the input sample size, the convergence of the scapulothoracic angle variance was visually analyzed from 100 to 10 000 simulations. The slope of the scapulothoracic angles variance between the 9 000 th and 10 000 th simulation was then computed. The linearity of the relationship between the kinematic model parameters and the scapulothoracic angles was graphically controlled.

Concerning the uncertainty analysis, the overall effect of kinematic model parameters on scapulothoracic angles was represented through two standard deviations to the mean computed from the all simulations of each participant and each movement.

For the first-order sensitivity analysis, the percentage of scapulothoracic angle variance explained by each parameter of the kinematic model was computed using Effective Algorithm for Computing Global Sensitivity Indices [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF]. For each DoF, the proportion of scapulothoracic angle variance explained by each kinematic model parameter was averaged over the time. Then, the influence of the kinematic model parameter variability on the variance explained was assessed using linear mixed models with logit correction. When a main effect was revealed by the linear mixed-model, post-hoc pairwise comparisons were performed using the Tukey's honestly significant difference procedure (R core Team, lme4 package). In addition, SPM random effect analyses [START_REF] Pataky | Generalized n-dimensional biomechanical field analysis using statistical parametric mapping[END_REF] were used to compute alpha parameters, quantifying the slope of the linear relationships between each kinematic model parameter and scapulothoracic angles variation over the time.

Results

Convergence analysis of scapulothoracic angle variance showed that 10 000 simulations were adequate since the slope of the variance between the last 1 000 simulations was inferior to 0.001 (See figure S1 and Table T1 in the supplementary file for more details).

Uncertainty

When altering simultaneously the 15 parameters of the kinematic model from -1 cm to 1 cm, mean scapulothoracic angles variation over the time ranged from 9.8° to 32.6°, 8.9° to 33.1°, 9.7° to 43.6°, 18.9° to 44.6° and 9.5° to 49.5° for arm elevation/depression in the frontal, scapular and sagittal plane, for horizontal abduction and for ball throw respectively (Fig. 345, and table T2 and Fig. S2-S3 in the supplementary file). Furthermore, scapulothoracic angle variations changed regarding the phase of the movement and the DoF, with the highest variabilities observed at the middle or at the end of the movement for arm elevations/depressions and throwing respectively.

Sensitivity

For all movements and DoF, the clavicle length is the parameter that most explained the scapulothoracic angle variance (from 14% to 25% on average, p<0.01), followed by the acromioclavicular joint location in the antero-posterior direction (from 9% to 18% on average, p<0.01) and the sternoclavicular joint location in the medio-lateral direction (from 5% to 15% on average, p<0.01) (Fig. 345, and table T2 and Fig. S2-S3 in the supplementary file).

Discrepancies between movements and DoF were observed concerning the influence of both the sternoclavicular joint location in the antero-posterior direction and the center of the ellipsoid in the medio-lateral and antero-posterior directions on scapulothoracic angle variance (from 3% to 9%, from 3% to 11%, and from 2% to 13% respectively). For all movements and DoF, the location of the scapulothoracic contact point, the ellipsoid vertical radius and the ellipsoid center location in the vertical direction presented the least influence on scapulothoracic angle variance (from 0.1% to 1.6% on average, p<0.01). The scapulothoracic angle variance was not fully explained by the 15 kinematic model parameters. The proportion of variance explained by the interaction effects between the kinematic model parameters ranged on average from 30% to 34%, 12% to 26%, 10% to 27%, 14% to 20% and 7% to 11% for arm elevation/depression in the scapular plane, horizontal abduction, ball throw, and elevation/depression in the frontal and sagittal planes respectively.

Alpha parameters computed with SPM random effect analyses revealed that for arm elevation/depression in the sagittal and scapular planes, and ball throwing, lengthening the clavicle by 1 cm increased scapular downward rotation and posterior tilt on average from 4.1° to 7.1° and from 4.0° to 4.4° respectively. Lengthening the clavicle by 1 cm also increased scapular internal rotation from 3.6° to 4.1° on average during arm elevation/depression in the sagittal plane and horizontal abduction (Fig. 6 to 8 and Fig. S4-S5 in the supplementary file). In addition, for all movements but the horizontal abduction, when moving anteriorly the acromioclavicular joint center or when moving laterally or posteriorly the sternoclavicular joint center by 1 cm, scapular downward rotation increased on average from 3.1° to 5.9°, from 2.9° to 5.6° and from 1.9° to 3.2° respectively, while the scapular posterior tilt increased from 3.1° to 3.7°, from 2.1° to 3.1° and from 1.8° to 2.7° respectively. The linear relationship between the above-mentioned parameters and the scapular internal rotation was only significant for the midpart of the movements. Concerning the horizontal abduction, moving anteriorly the acromioclavicular joint center or moving laterally or posteriorly the sternoclavicular joint center by 1 cm led to increase scapular internal rotation by on average 3.7°, 2.8° and 1.9° respectively.

(Fig. 6 to 8 and Fig. S4-S5 in the supplementary file). Finally, moving laterally the ellipsoid center by 1 cm mainly increased the scapular upward rotation and internal rotation for all movements up to 3.0° and 1.8° respectively. When moving anteriorly the ellipsoid center by 1 cm, upward and internal rotation were significantly increased only for arm elevation/depression in the scapular and frontal planes and during horizontal abduction. (Fig. 6 to 8 and Fig. S4-S5 in the supplementary file).

Discussion

The purposes of this study were first to determine the influence of the kinematic model parameters variability on scapulothoracic angle estimates, and second to define which kinematic model parameters have the largest effect on scapulothoracic angle estimates. The main findings were that altering kinematic model parameters from -1 cm to 1 cm resulted in scapulothoracic angle estimate variations up to 50° depending on the movement studied, and that scapulothoracic angles were the most sensitive to the clavicle length.

Very little attention has been paid so far to the uncertainty analysis of shoulder kinematic model parameters on scapulothoracic angles, since to the best of our knowledge only one study focused on this topic [START_REF] Hybois | Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study[END_REF]. We observed a variability of scapulothoracic angles about twice greater than the one described during wheelchair propulsion. However, we have explored movements with higher range of motion and only ellipsoid parameters were altered in Hybois's study [START_REF] Hybois | Effects of ellipsoid parameters on scapula motion during manual wheelchair propulsion based on multibody kinematics optimization. A preliminary study[END_REF], while other parameters relative to the clavicle and scapulothoracic contact point were also altered in our study. In a clinical context, the results of our uncertainty analysis suggest that clinicians should take caution when a kinematic model is used to assess the evolution of scapular positioning over the time. Indeed, slight differences in the calibration of the kinematic model may lead to drastic changes in scapulothoracic angles, which may be a non-negligible confounding factor when assessing scapular dyskinesis for instance [START_REF] Carbone | Scapular dyskinesis after Latarjet procedure[END_REF].

Therefore, calibrating the kinematic model demands caution when assessment of scapular positioning is performed.

Cadaveric and in vivo studies showed that changing the clavicle length significantly affects scapulothoracic angles [START_REF] Hillen | The effect of experimental shortening of the clavicle on shoulder kinematics[END_REF][START_REF] Su | Evaluation of three-dimensional scapular kinematics and shoulder function in patients with short malunion of clavicle fractures[END_REF]. A previous sensitivity analysis also depicts that clavicle length is the kinematic model parameter that most influences scapulothoracic angles during arm abduction [START_REF] Habachi | Global sensitivity analysis of the kinematics obtained with a multi-body optimisation using a parallel mechanism of the shoulder[END_REF]. Such outcomes are reinforced by our study since we observed that the clavicle length was also the most influential parameter for all investigated movements. In consequence, when scapulothoracic angles are estimated with MKO, the clavicle length of the patient must be estimated accurately. To that end, computed tomography [START_REF] Hoj | Validation of Navigation Ultrasound for Clavicular Length Measurement[END_REF] and conventional radiograph [START_REF] Smekal | Length Determination in Midshaft Clavicle Fractures: Validation of Measurement[END_REF] are the favored option for clinicians but expose the patient to radiations and may be time consuming. Low dose biplane X-rays [START_REF] Dumas | Subject-specific model-derived kinematics of the shoulder based on skin markers during arm abduction up to 180° -Assessment of 4 gleno-humeral joint models[END_REF] can be also considered as an alternative solution. In addition, a procedure based on navigation ultrasound has been suggested as a non-invasive method to assess clavicle length [START_REF] Hoj | Validation of Navigation Ultrasound for Clavicular Length Measurement[END_REF]. Although this approach seems promising, it might be also time consuming in a clinical context. Finally, palpation method has also been proposed to measure clavicle length [START_REF] Stegeman | Measurement of clavicular length and shortening after a midshaft clavicular fracture: Spatial digitization versus planar roentgen photogrammetry[END_REF], nevertheless further study should explore the accuracy of such procedure.

Noninvasive palpation methods used to determine sternoclavicular and acromioclavicular joint center location are subject to errors ranging from 1.5 cm to 3 cm on average [START_REF] Michaud | Determining in vivo sternoclavicular, acromioclavicular and glenohumeral joint centre locations from skin markers, CT-scans and intracortical pins: A comparison study[END_REF]. Besides, we observed that altering the location of these joint centers may lead to variations up to 6°/cm. Consequently, determining sternoclavicular and acromioclavicular joint locations with noninvasive methods, such as palpation, may not be satisfying to accurately estimate scapulothoracic angles. In a clinical context, imagery procedures based on low dose radiation should be preferred [START_REF] Nicholson | Determining 3D scapular orientation with scapula models and biplane 2D images[END_REF]. Although time consuming, the use of ultrasonography coupled with a motion capture system seems to be interesting since it avoids any radiation exposure [START_REF] Hoj | Validation of Navigation Ultrasound for Clavicular Length Measurement[END_REF].

However, further studies demand also to be conducted to assess the accuracy of such a method.

When assessing scapulothoracic angles, most of shoulder kinematic models enforce the scapula to stay in contact with an ellipsoid representing the scapulothoracic sliding plane to avoid any non-physiological scapular winging [START_REF] Duprey | Kinematic models of the upper limb joints for multibody kinematics optimisation: An overview[END_REF]. Several procedures have been suggested to calibrate the ellipsoid (i.e., radii and center location) [START_REF] Prinold | Scaling and kinematics optimisation of the scapula and thorax in upper limb musculoskeletal models[END_REF][START_REF] Michaud | Scapular kinematic reconstruction -segmental optimization, multibody optimization with openloop or closed-loop chains: which one should be preferred?[END_REF][START_REF] Bolsterlee | Modelling clavicular and scapular kinematics: from measurement to simulation[END_REF]. Our results led us to suggest that efforts should focus on the ellipsoid center location especially in the anteroposterior and mediolateral directions. In addition, the vertical position of the ellipsoid center or its vertical radii may not deserve too much attention during the calibration procedure since they barely influence scapulothoracic angle estimates. Although, previous studies depicted promising approaches to calibrate the ellipsoid [START_REF] Michaud | Scapular kinematic reconstruction -segmental optimization, multibody optimization with openloop or closed-loop chains: which one should be preferred?[END_REF][START_REF] Blache | In Press) Does scapulothoracic contact point trajectory follow an ellipsoid surface?[END_REF], further researches should be conducted to calibrate subject-and movement-specific ellipsoid when scapulothoracic angle estimates are of interest.

To that end, MKO including model parameters in the optimization process may be explored [START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF][START_REF] Andersen | A computationally efficient optimisation-based method for parameter identification of kinematically determinate and over-determinate biomechanical systems[END_REF].

A previous study highlighted the relevancy of performing the sensitivity analysis throughout the movement instead of considering time averages since variations may occur over the movement [START_REF] Brizard | A screening method to analyse the sensitivity of a lower limb multibody kinematic model[END_REF]. Although the main findings are discussed above, we also observed that the most influential parameters may vary over the time, meaning that ranking the kinematic model parameters that should be calibrated in priority is not always a straight forward process. In addition, only first order effects have been analyzed, while interactions between parameters may have a non-negligible influence on scapulothoracic angles estimates. The fact that the sensitivity analysis revealed major interaction effects is logical as the kinematic chain represents a close loop between scapula and thorax. Therefore, further investigations focusing on screening approach could be conducted to assess these interaction effects. This study presents some limitations warranting discussion. First, our results can only be applied in the domain of variation investigated (i.e., ± 1cm of alteration of the kinematic model parameters). Nevertheless, in a preliminary study, we attempt to make vary kinematic model parameters from -2 cm to 2 cm, but the optimization procedure did not converge in many simulations. Second, scapulothoracic contact point location barely influenced scapulothoracic angle estimates. However, this result should be taken with caution since different conclusions may have been drawn when more than one contact point is defined or when the contact point is initially located at the medial border of the scapula [START_REF] Prinold | Scaling and kinematics optimisation of the scapula and thorax in upper limb musculoskeletal models[END_REF].

Conclusion

Uncertainty of shoulder kinematic model parameters may drastically influence scapulothoracic angle estimates either for analytical or sport movement. In addition, scapular kinematics are particularly sensitive to clavicle length, and to a lesser extent to sterno-and acromioclavicular joint locations, and ellipsoid parameters. Therefore, it is necessary to keep on developing numerical tools and procedures both accurate and easy-to-use to calibrate shoulder kinematic models when looking for assessing scapulothoracic angles.

Figure 1 -

 1 Figure 1-Location of skin markers on a representative participant.

Figure 2 -

 2 Figure 2 -Parameters of the kinematic model undergoing alterations. 3D coordinates of the ellipsoid center (OE) in the thorax reference system {XT, YT, ZT} (n=3), ellipsoid radii (n=3), 3D coordinates of the sternoclavicular joint (OC) in the thorax reference system (n=3), 3D coordinates of the acromioclavicular joint (OS) in the scapula reference system {XS, YS, ZS} (n=3), the clavicle length being the Euclidian distance between the sternoclavicular and acromioclavicular joints (n=1), 2D scapulo-thoracic contact point (STCS) in the scapular plane {XS, ZS} (n=2).

Figure 3 -

 3 Figure 3 -Nominal (thick black line), mean (dashed line), and ±2 standard deviations from the mean (grey lines) scapulothoracic angles (ST) as a function of the normalized time (from 0% to 100%) computed during arm elevation/depression in the scapular plane for the 10 000 simulations and the nine participants. The color bars represent the average proportion of scapulothoracic angle variance explained by each of the 15 parameters of the kinematic model. EllR for ellipsoid radii, EllC for ellipsoid center, ClavLength for clavicle length, AC for acromioclavicular joint location, SC for sternoclavicular joint location, ScapCon for scapulothoracic contact point location

Figure 4 -

 4 Figure 4 -Nominal (thick black line), mean (dashed line), and ±2 standard deviations from the mean (grey lines) scapulothoracic angles (ST) as a function of the normalized time (from 0% to 100%) computed during horizontal abduction for the 10 000 simulations and the nine participants. The color bars represent the average proportion of scapulothoracic angle variance explained by each of the 15 parameters of the kinematic model. EllR for ellipsoid radii, EllC for ellipsoid center, ClavLength for clavicle length, AC for acromioclavicular joint location, SC for sternoclavicular joint location, ScapCon for scapulothoracic contact point location

Figure 5 -

 5 Figure 5 -Nominal (thick black line), mean (dashed line), and ±2 standard deviations from the mean (grey lines) scapulothoracic angles (ST) as a function of the normalized time (from 0% to 100%) computed during ball throwing for the 10 000 simulations and the nine participants. The color bars represent the average proportion of scapulothoracic angle variance explained by each of the 15 parameters of the kinematic model. EllR for ellipsoid radii, EllC for ellipsoid center, ClavLength for clavicle length, AC for acromioclavicular joint location, SC for sternoclavicular joint location, ScapCon for scapulothoracic contact point location

Figure 6 -

 6 Figure 6 -Mean (line) and standard error of the mean (shadow) of the alpha parameter, representing the slope of the linear relationship between kinematic model parameters and scapulothoracic angle variation, with respect to the normalized time (from 0% to 100%) during arm elevation/depression in the scapular plane. Horizontal lines represent, for each degree of freedom, the part of the curve for which a significant relationship is revealed (p<0.05).

Figure 7 -

 7 Figure 7 -Mean (line) and standard error of the mean (shadow) of the alpha parameter, representing the slope of the linear relationship between kinematic model parameters and scapulothoracic angle variation, with respect to the normalized time (from 0% to 100%) during horizontal abduction. Horizontal lines represent, for each degree of freedom, the part of the curve for which a significant relationship is revealed (p<0.05).

Figure 8 -

 8 Figure 8 -Mean (line) and standard error of the mean (shadow) of the alpha parameter, representing the slope of the linear relationship between kinematic model parameters and scapulothoracic angle variation, with respect to the normalized time (from 0% to 100%) during ball throwing. Horizontal lines represent, for each degree of freedom, the part of the curve for which a significant relationship is revealed (p<0.05).

Compliance with ethical standards

Conflict of interest: The authors declare that they have no conflict of interest.

Informed consent: This research involved Human Participants who signed an informed consent after the study was approved by the local ethical comittee.