2434 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 3D Perception With Slanted Stixels on GPU Daniel Hernandez-Juarez , Antonio Espinosa , David Vazquez , Antonio M. Lopez , and Juan C. Moure Abstract—This article presents a GPU-accelerated software design of the recently proposed model of Slanted Stixels, which represents the geometric and semantic information of a scene in a compact and accurate way. We reformulate the measurement depth model to reduce the computational complexity of the algorithm, relying on the confidence of the depth estimation and the identification of invalid values to handle outliers. The proposed massively parallel scheme and data layout for the irregular computation pattern that corresponds to a Dynamic Programming paradigm is described and carefully analyzed in performance terms. Performance is shown to scale gracefully on current generation embedded GPUs. We assess the proposed methods in terms of semantic and geometric accuracy as well as run-time performance on three publicly available benchmark datasets. Our approach achieves real-time performance with high accuracy for 2048  1024 image sizes and 4  4 Stixel resolution on the low-power embedded GPU of an NVIDIA Tegra Xavier. Index Terms—Stereo vision, stixel world, autonomous vehicles, scene understanding, computer vision, embedded systems, GPU acceleration Ç 1 INTRODUCTION as ground-like planes, upright objects or sky, which are the primitive geometric elements found in urban environments. driver assistance systems (ADAS), autono- A DVANCED mous vehicles, robots and other intelligent devices need to understand their environment. Stereo camera sys- This representation transforms millions of disparity pixels to hundreds or thousands of Stixels. At the same time, most scene structures, such as free space and obstacles, which are tems provide geometric (distance) and semantic (classifica- relevant for autonomous driving tasks, are adequately tion) data to estimate both the semantic class and the represented. distance of objects and the free space in a given scene. The The Stixel world can model the scene structures with cer- large amount of low-level per-pixel data is very costly to tain constraints, e.g. sky is above the horizon line and transmit and process and commonly a medium-level repre- objects usually lie or stay on ground. Generally, the geomet- sentation known as the Stixel World [1], [2] is used. It relies ric constraints of a scene are tightened to the vertical direc- on the fact that man-made environments mostly present tion. Hence, the environment can be modeled as a column- horizontal and vertical planar surfaces, like roads, side- wise segmentation of the image with a 3D stick-like shape, walks or soil (horizontal), and buildings, pedestrians or cars i.e. a set of Stixels, c.f. Fig. 1. The segmentation of the image (vertical). This medium-level representation must be com- is estimated by solving a column-wise energy minimization puted in real-time to serve as a building block of higher- problem, taking depth and semantic cues as inputs as well level modules, such as localization and planning. as a priori information that is used to regularize the solution. The Stixel world has been successfully used for repre- The Stixel model has been successfully used for automotive senting traffic scenes, as introduced in [2]. The field of intel- vision applications either to decrease parsing time, increase ligent vehicles has been using this model over the last years accuracy, or both [14], [15], [16], [17], [18], [19], [20]. [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]. The Stixel The Slanted Stixel world model recently proposed in [11], world defines a compact representation of the dense 3D dis- [12] generalizes the original proposal with a more flexible parity data obtained from stereo vision that uses rectangles, plane model that overcomes the previous rather restrictive the so called Stixels, as elements. Stixels are classified either constant depth and constant height assumptions for object and ground Stixels, respectively, and accurately represents  Daniel Hernandez-Juarez is with the SLAMcore Ltd., SE1 1JA, London, U.K. arbitrary kinds of slanted objects and non-flat roads. Basi- E-mail: dhernandez0@gmail.com. cally, the Slanted model defines a plane in the disparity  Antonio Espinosa and Juan C. Moure are with the Universitat Autonoma space defined by two random variables: line slope and inter- de Barcelona, 08193 Bellaterra, Spain. E-mail: {AntonioMiguel.Espinosa, juancarlos.moure}@uab.cat. cept. It has been proved to provide substantially better qual-  David Vazquez is with the Element AI, Montreal, QC H2S 3G9, Canada. ity on scenarios with non-flat roads. E-mail: dvazquez@elementai.com. Stixel segmentation is a highly computationally complex  Antonio M. Lopez is with the Universitat Autonoma de Barcelona, 08193 Bellaterra, Spain, and also with the Computer Vision Center, Universitat optimization problem, and the higher representation quality Autonoma de Barcelona, 08193 Bellaterra, Spain. of the slanted model comes at the price of even higher com- E-mail: antonio@cvc.uab.es. plexity: if the input image contains h rows and w columns, Manuscript received 16 Apr. 2020; revised 13 Mar. 2021; accepted 18 Mar. 2021. the algorithmic complexity is Oðw  h3 Þ. One way to reduce Date of publication 22 Mar. 2021; date of current version 8 Apr. 2021. complexity is to reduce the size of the input images (depth (Corresponding author: Daniel Hernandez-Juarez.) and semantics), both in horizontal and vertical dimensions. Recommended for acceptance by F. Wolf. Digital Object Identifier no. 10.1109/TPDS.2021.3067836 The problem is that increasing the granularity of the Stixels This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2435 Fig. 1. Scene representation of a challenging street environment obtained by our method. Geometric or Depth (left), semantic (center) and 3D (right) representations are shown. In Fig. 1a, color encodes depth from close (red) to far (green). In Fig. 1b, color encodes the semantic class following [21]. (or, equivalently, decreasing the resolution of the Stixel) also used to compute the cost of each Stixel in constant time. The reduces the quality of the segmented representation. definition of Slanted Stixels using two random variables 2048-by-1024-pixel images scaled 4 and 8 times are seg- (slant and intercept) instead of one, disables the usage of mented, respectively, at 0.7 and 6.6 frames per second (fps) the same algorithmic strategies used for the basic Stixels [2]. on a six-core Intel i7-6800K processor [12]. These perfor- For more details, see Section 3.4. mance results do not meet the real-time or energy efficiency We have also developed a completely new massively requirements of autonomous driving applications. Real- parallel design for GPU execution based on some ideas time performance is achieved by further reducing the reso- taken from our previous work [4], which implements the lution of the input images to levels that sacrifice the quality algorithm for the original Stixel model [2]. Compared to [4], of the resulting segmentation. Dedicated hardware designs our proposal makes better use of local (register) and shared (e.g. FPGA or ASIC) may provide faster implementations, memory for storing cost, index and Summed Area tables, but are very inflexible and expensive regarding changes in and avoids some of the synchronization operations. These the algorithms. enhancements are achieved by modifying the distribution Embedded GPU-accelerated systems, like the NVIDIA of work among threads, and by reusing the same memory Jetson and DrivePX platforms, allow low-cost and low- areas at different phases of the execution, see Section 4. energy consumption, real-time Stixel computation. GPUs Table 1 helps comparing our contribution with respect to are very well suited for algorithms exhibiting massive paral- the previous work. The last row in the table describes our lelism, like the dynamic programming techniques used for proposal for a novel cost equation that basically removes Stixel segmentation. However, to achieve competitive per- the usage of a uniform distribution to model the occurrence formance, inefficiencies due to existing data dependencies of disparity measurement outliers by the usage of the confi- that lead to explicit synchronizations, must be reduced by dence of the depth estimation and of the semantic cues. careful work distribution and cooperation, along with a Thanks to this reformulation, our algorithm has lower proper data layout design. computational complexity than the original algorithm for Our goal is to develop a faster algorithm and GPU-accel- Slanted Stixels. Also, the last column in the table describes erated implementation solving the optimization problem the characteristics of the two GPU implementations that we for Slanted Stixel segmentation. We want to process high- have developed and evaluated in this paper. The proposed resolution input images at high speed in order to reach real- cost equation does not only lead to an algorithm that com- time performance with low energy consumption, but with- putes Stixel costs in constant time, but also allows using out sacrificing segmentation accuracy. To this extent, we smaller SATs (Summed Area Tables), which fit into the propose a novel depth measurement model that enables the shared memory of the GPU and contribute to a more effi- application of several algorithmic techniques to reduce the cient execution. computational complexity of Slanted Stixels from Oðw  h3 Þ We have done an in-depth evaluation of the two GPU- to Oðw  h2 Þ. Very briefly, the proposed cost equations can accelerated programs in terms of run-time as well as seman- be formally derived to find an optimal solution of the result- tic and depth accuracy, carried out on several benchmarks. ing weighted least squares problem in closed form, an then Fig. 2 shows some of the main results to illustrate the trade- multiple pre-computed Summed Area Tables (SATs) can be off between accuracy and run-time. Our proposal using the TABLE 1 Comparison of the Different Stixel Models in Terms of Inputs, Complexity, and GPU Implementation Details Model Inputs Complexity Accuracy GPU version Stixel World [2] Disparity w  h2 Bad Implemented in [4], Most SATs in global memory, Cost and Index table in shared memory, Cost computed in constant time Disparity, Disparity confidence, Implemented in this work, Semantic SATs in shared memory, Slanted Stixels [12] w  h3 Good Cost and Index table in registers, Disparity SATs not feasible: Semantic segmentation cost computed in linear time, Very unbalanced work Novel Disparity, Cost Disparity, Disparity confidence, Implemented in this work, Disparity and Semantic SATs in w  h2 Good shared memory, Cost and index table in registers, Cost Measurement [ours] Semantic segmentation computed in constant time 2436 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 stripes allowing to capture the camera data in a polar fash- ion. Also, the Stixel data model is similar to the forward step usually found in occupancy grid maps [7]. However, the Stixel inference method in the image domain presents important differences compared to classical grid-based approaches. Our work builds upon the proposal from [3]: they use semantic cues in addition to depth to extract a Stixel repre- sentation, which is able to provide a rich yet compact repre- sentation of the traffic scene. We also base our method on [12]: the Slanted Stixels model incorporates a novel plane model together with effective priors on the plane parame- ters, and it is able to represent scenes with complex non-flat roads. Fig. 2. This chart illustrates the trade-off between accuracy and run-time. There are some methods [1], [5], [8], [9], that represent We compare Slanted Stixels [12] and our method, both with Stixel reso- lution 4  4 and 8  8. Disparity error corresponds to results for SYN- simplified scene models with a single Stixel per column. THIA-SF [11], and run-time was measured on the NVIDIA Tegra Xavier The advantage of these approaches is that the computa- embedded GPU. Our model with 4  4 Stixel resolution (blue triangle) tional complexity of the underlying algorithms is linear, but shows a good trade-off between accuracy and run-time, being the one they cannot represent some complex scenarios found in the closer to the bottom-left corner. real world, e.g. a pedestrian and a building in the same column. modified cost equation increases the average disparity error A recent work [10] uses edge-based disparity maps to around 3.5 percent for both Stixel resolutions (4x and 8x), compute Stixels. Their method is fast but gives inferior but the resulting program runs on an embedded low-power accuracy compared to the original Stixel model [25]. GPU more than 8 times faster. In this particular H/W sys- Finally, there is some work proposing fast implementa- tem, our GPU implementation of the original Slanted Stixels tions for Stixel computation. The FPGA implementation achieves real-time performance only for Stixel resolution from [17] runs at 25 Hz with an image resolution of 1 mega- 8x8 (red square), while our GPU implementation using the pixels and a Stixel resolution of 5 pixels. novel model achieves real-time performance for Stixel reso- An earlier work [11], [12] proposed a method to rule out lution 4x4 (blue triangle), which improves the disparity the most unlikely Stixel cuts, thereby saving computational error by more than 9 percent. More details are provided in time by applying the Stixel algorithm only to the probable Section 5. Stixel cuts. The methods for generating the over-segmenta- To sum up, we provide two versions for GPU-accelerated tion of the image where restricted to have linear time com- Stixel segmentation that can achieve real-time and low- plexity so that the added run-time is not high. Two methods energy consumption in our embedded system hardware. were analyzed, one based on times series compression and Both versions offer a different trade-off between segmenta- one using a trained neural network. In practice, most of the tion accuracy and running time. In a time-limited and low- time the number of feasible Stixel cuts remaining in the power scenario, the version using our proposed cost mea- over-segmentation is significantly small. But the biggest surement equation ends up providing better segmentation drawback of this approach is that the run-time of the algo- accuracy. rithm is variable and then non-predictable, which is a prob- The remainder of this paper is structured as follows. Sec- lem for building a real-time system. Anyway, this method is tion 2 reviews the state of the art. Section 3 presents the for- orthogonal to our proposal and both could be combined. mulation of Slanted Stixels, while the modification of the A previous GPU-accelerated version of the Stixel’s measurement model and its rationale is detailed in Sec- method [4] only includes the original non-slanted model tion 3.4. Section 4 explains basic concepts for efficient GPU [2]. The GPU version proposed in this work implements a acceleration and describes our proposed GPU-based optimi- richer Stixel method that incorporates more cues, e.g. zations for real-time Stixel computation. Section 5 presents semantic segmentation [3], disparity confidence [25], as well the experiments we carried out and analyzes the accuracy as, a novel depth model for slanted scenes [12]. Compared and execution speed of our proposed method. Finally, we to [4], the usage of local (register) and shared memory is state our conclusions in Section 6. improved by an enhanced arrangement of the Summed Area Tables (SATs) and the cost and index tables, by a better distribution of the work among the parallel threads, and by 2 RELATED WORK a simpler synchronization required during the execution. We will first comment on work proposing different road scene models. Occupancy grid maps are models used to represent the surrounding of the vehicle [17], [22], [23], [24]. 3 THE STIXEL MODEL Typically, a grid in bird’s eye perspective is defined and The Stixel world is a compressed representation of a 3D used to detect occupied grid cells and then, from this infor- scene that preserves its relevant structure. Given that the mation, to extract the obstacles, navigable area, and non- vertical dimension dominates the structure of street envi- observable areas from range data. These grids and the Stixel ronments, the Stixel world segments the image into inde- world both represent the 2D image in terms of column-wise pendent columns composed of stick-like super-pixels with a HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2437 segmentation S given all measurements M from that partic- ular column. Applying the Bayes’ theorem, the posterior probability can be rewritten using the unnormalized likelihood and prior distributions as Z1 P~ ðM j S ÞP~ ðS Þ. In order to avoid numerical problems with small magnitudes of the individ- ual probabilities, the likelihoods are transformed to log-like- lihoods via P ðS j M Þ ¼ eE ðs ;m mÞ , and the MAP estimation problem is then converted to a cost minimization problem, where E ðÞ is the energy (or cost) function. The energy function is the summation of the energies of the whole Stixel segmentation, which can be separated into the likelihood or data term, Edata ðÞ, and the prior term, Eprior ðÞ. X N E ðs ; m Þ ¼ ðEdata ðs i ; m Þ þ Eprior ðsi ÞÞ : (1) i¼1 The likelihood or data term Edata ðÞ rates how well the meas- urements m fit to the overlapping Stixel si . This energy is Fig. 3. Example of the Stixel segmentation and labeling of a column in a further split in a semantic term and a depth term typical scene (on the right). The input disparity measurements (black thin lines) and output Stixels encoded with semantic colors (colored thick Edata ðs i ; m Þ ¼ Edepth ðsi ; d Þ þ wl  Esem ðsi ; l Þ : (2) lines) are shown on the left. Taken from [12]. The parameter wl controls the influence of the semantic data term. The input is provided by a fully convolutional net- 3D planar depth model and semantic labels. There are three workP (FCN) that delivers normalized semantic scores lv ðci Þ structural classes derived exclusively from depth data: with ci lv ðci Þ ¼ 1 for all classes ci at pixels v. The semantic ground (Stixels with a slant similar to the expected ground energy favors semantic classes of the Stixel that fit to the plane), object (almost vertical Stixels, usually lying on the observed pixel-level semantic input [3]. The semantic likeli- ground), and sky (Stixels at infinite distance). Semantic clas- hood term is ses are refinements to those structural classes (e.g. road or t sidewalk are ground classes, whereas building and vehicle X vi are object classes). Prior to the segmentation, the per-pixel Esem ðsi ; l Þ ¼ logðlv ðci ÞÞ: (3) input images are downsized to the desired vertical and hori- v¼vb i zontal Stixel resolution. An example of Stixel segmentation is presented in Fig. 3. The depth term is defined by a probabilistic and genera- The column highlighted in the image on the right is down- tive sensor model Pv ðÞ that considers the accordance of sized, and the disparity measurements (inverse of depth) the depth measurement dv at row v to the depth model of for each Stixel on the column are shown on the left. The Stixel si t resulting Stixel segmentation and labeling are defined by X vi the colored thick lines. Edepth ðs i ; d Þ ¼ log ðPv ðDv ¼ dv j S i ¼ si ÞÞ: (4) The rest of this section defines the mathematical formula- v¼vb i tion of the Stixel model and how to solve the problem of joint optimization through dynamic programming. The last Following Slanted Stixels [12], we use a plane depth subsection presents our proposal for modifying the mathe- model that overcomes the previous rather restrictive con- matical model in order to reduce the computational com- stant depth and constant height assumptions for object and plexity of the problem. ground Stixels, respectively. To this end, we formulate the depth model Dðs i ; vÞ using two random variables defining a plane in the disparity space (slope and intercept) that evalu- 3.1 Mathematical Formulation ates to the disparity in row v via A Stixel column segmentation S consists of an arbitrary number N of Stixels, si , each representing four random vari- Dðs i ; vÞ ¼ bi  v þ ai : (5) ables: the Stixel extent via bottom row Vib and top row Vit , as Note that we assume narrow Stixels and thus can neglect well as its semantic class Ci and depth model Di (slope and one plane parameter, i.e. the roll. intercept). Thereby, the number of Stixels itself is a random The measurement model for disparities is then defined as variable that is optimized jointly during inference. The joint a combination of a Gaussian and a uniform distribution segmentation and labeling problem is carried out indepen- dently for each image column via optimization of the poste-  2 cv ðdv Dðs i ;vÞÞ rior distribution P ðS j M Þ, a Maximum A Posteriori pout 1  pout  s ðs i Þ Pv ðDv j S i Þ ¼ þ e : (6) estimation problem (MAP) defined over a Stixel ZU ZG ðs i Þ 2438 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 The Gaussian distribution models the typical disparity 8 > > Edata ðobk Þ þ Eprior ðobk0 Þ noise and the uniform distribution, weighted by a constant > E ðob0k Þ þ E > 0 0 > prior ðob1 ; ob Þ þ OB k > > data 1 probability for outliers pout , makes the model more robust > > Edata ðob1 Þ þ Eprior ðob1 ; gr Þ þ GR0 k k 0 > < to outliers. The Gaussian sensor noise model is centered at Edata ðobk1 Þ þ Eprior ðobk1 ; sk0 Þ þ SK 0 the expected disparity Dðsi ; vÞ given the depth model of the OBk ¼ min > > ::: Stixel and has confidence cv . ZU and ZG ðs i Þ normalize > > > > E ðob k Þ þ E ðob k k ; ob k1 Þ þ OBk1 the distributions. Similarly to [25], we use the confidence > > data k prior > > E ðob Þ þ Eprior ðobk ; gr Þ þ GRk1 k k k1 of the depth estimates cv to influence the shape of the dis- : data kk Edata ðobk Þ þ Eprior ðobkk ; skk1 Þ þ SK k1 tribution s ðsi Þ. (9) The prior or smoothness term captures the knowledge about the traffic scene, such as, sky Stixels are unlikely below the horizon line, objects tend to be close to the Eq. (8) defines the solution for the base case problem, ground, or there is a small number of objects in the which is the case of one Stixel made by the first single pixel. scene. In order to model the complexity of the segmenta- Eq. (9) indicates how to solve a problem of size k, i.e. how to tion, we include a constant term for each segment to compute the partial solutions OBk , GRk , and SK k , using the favor configurations composed of fewer Stixels. The solutions for smaller problems. We only show the case for Markov property is used to reduce the prior definition to object Stixels, but the other cases are solved similarly. All the pairwise mutual dependencies of each pair of adjacent possible object Stixels ending at position k (and starting at Stixels and the likelihood of the bottom Stixel. Refer to positions from 1 to k) are connected with the last Stixel of [7], [12] for a more comprehensive definition of the the segmentation with minimal cost of the corresponding priors. size, which were previously computed and memorized in We define a prior term for the depth model of Stixels, C. Connections are evaluated for the three Stixel structural Eplane ðs i Þ, that expects the two random variables A, B repre- classes using the prior term. senting the plane parameters of a Stixel to be Gaussian dis- Once the cost table C is completely computed, a back- tributed, i.e. tracking procedure retrieves the resulting Stixel segmenta- tion by starting from the top row of C and computing the !2 !2 a  maci b  mbci successive minimum value Cmin k ¼ minðOBk ; GRk ; SK k Þ. Eplane ðs i Þ ¼ þ log ðZ Þ: (7) s aci s bci 3.3 Reduce the Algorithm’s Complexity Using SATs As shown by Eq. (9), solving a sub-problem of size k This prior favors planes in accordance to the expected 3D requires computing the minimum cost of all the k possible layout corresponding to the particular geometric class ci . E. positions of a cut between Stixels for the 3 possible struc- g. object Stixels are expected to have an approximately con- tural classes. We use this to calculate the time complexity of stant disparity, i.e. mbobject ¼ 0. The expected road slant the algorithm. Since the number of structural classes is con- maground can be set using prior knowledge or by means of an stant and k ranges from 0 to the total number of pixels in a specific method for road surface detection. column, h, then the Stixel segmentation problem for a single column requires Oðh2 Þ steps. The backtracking phase can be 3.2 Algorithm Based on Dynamic Programming done in a linear number of steps, OðhÞ, by creating an index Dynamic Programming (DP) solves a complex problem by table linking each Stixel and the next Stixel with minimum dividing it into simpler sub-problems and storing the partial cost during the DP solving phase. solutions on memory. This way, when the same sub-prob- Each step of the DP process must compute the prior and lem appears, computation time is saved by retrieving the data terms of one single Stixel. The prior term is a function partial solution from memory instead of solving the sub- of the parameters of one or two Stixels, and can be com- problem repeatedly. puted in a constant number of operations. The data term, We apply the DP strategy to compute the column seg- though, depends on the depth, confidence, and semantic mentation with minimum global cost. In order to express class measurements of all the pixels composing the Stixel, the optimization problem as a recursive resolution of and therefore requires a number of operations proportional smaller sub-problems, we use a special notation for the to the Stixel length, which ranges between 1 and h. The chal- three different structural classes: obtb ¼ fvb ; vt ; objectg, lenge is to express the computation of the data term as a grtb ¼ fvb ; vt ; groundg, and sktb ¼ fvb ; vt ; skyg. OBk (respec- constant number of operations. We achieve this goal by pre- tively, GRk and SK k ) refers to the aggregated cost corre- computing partial results derived from the measurement sponding to the optimal Stixel segmentation from data (similarly to [4], [26]) and by a slight modification of position 0 to k of the given column, assuming that the last Eq. (6) that facilitates the parallelization. Stixel is an object (respectively, ground and sky). Given The semantic cost of a Stixel is the summation of the loga- the previous notation, we next show the recursive defini- rithm of the probabilities corresponding to the individual tion of the problem: pixels (c.f. Eq. (3)). We pre-compute the values correspond- ing to each pixel and store them into a Look-Up Table OB0 ¼ Edata ðob00 Þ þ Eprior ðob00 Þ (LUT), one for each semantic class. Then, we pre-compute GR0 ¼ Edata ðgr00 Þ þ Eprior ðgr00 Þ (8) the prefix sum or Summed-Area Table (SAT [27]) of each SK 0 ¼ Edata ðsk00 Þ þ Eprior ðsk00 Þ LUT, i.e. the successive accumulated costs corresponding to HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2439 all the previous pixels. The computation of the semantic The optimal plane parameters of a Stixel si can be substi- data term of a Stixel is then performed in constant time as tuted in Eq. (6) to compute the depth cost of each pixel and follows: then compute the summation of the cost to obtain the over- all cost of the Stixel, Edepth ðs i ; d Þ, as in Eq. (4). This is how Esem ðsi ; l Þ ¼ SATci ðvti Þ  SATci ðvbi  1Þ: (10) the computation is implemented in [11], [12], giving raise to a total algorithm complexity of Oðh3 Þ steps per image col- Notice that the total computation complexity for creating umn, which is very expensive. the SATs corresponding to an image column is OðhÞ, which Eq. (6), in its current form, cannot be formally derived to is lower than the computation complexity of the DP algo- apply the same kind of mathematical and computational rithm: Oðh2 Þ. transformations as the ones done for computing in constant The data term of sky Stixels only depends on the dispar- time the plane parameters. The problem is due to the uni- ity and confidence of the pixel individually, and can be form distribution that was proposed in the original Stixel computed in constant time by using SATs. But the data world, and was critical to model the occurrence of disparity term of ground and object Stixels depends not only on the measurement outliers. Our model, though, includes alterna- measurements but also on the depth model used. tives to soften the effect of those outliers, like the usage of In a previous work [4], [26], they implemented a non- confidence for the disparity measurements (invalid dispar- slanted depth model with a pre-determined constant road ities are modelled as having zero confidence) and the usage slant and a constant depth for objects. For ground Stixels of semantic cues. Our proposal, then, is to remove the uni- they substituted the depth model Dðsi ; vÞ in Eq. (6) and pre- form distribution from the depth model computed the LUTs and the corresponding SATs for each  2 cv ðdv Dðs i ;vÞÞ image column. Object Stixels, though, have a constant model 1  s ðs i Þ that is set as the mean disparity of the Stixel. They solved Pv ðDv j S i Þ ¼ e : (12) ZG ðsi Þ the problem by creating a separate SAT for each possible integer value for Dðsi ; vÞ; i.e. they quantized the mean dis- The logarithm of the previous equation can be computed parities into integer values. This approach proved empiri- in constant time by using multiple pre-computed SATs. An cally to be both accurate and efficient, with time and additional advantage of this computational design versus memory complexities proportional to Oðh  dmax Þ, where the proposal in [4] is that all the required SATs have time dmax is the maximum disparity measurement, which in and memory complexity OðhÞ, instead of Oðh  dmax Þ, and practice is lower than h. that the disparity range is not quantized. If the input image contains w columns, then the time 3.4 Modified Measurement Model for Slanted Stixels complexity for the proposed algorithm is Oðw  h2 Þ. If out- Compared to the model used in [4], the Stixel model consid- lier disparity measurements get a low confidence estima- ered here and described in Section 3.1 is much more elabo- tion, cv , or the semantic data is robust for those outliers, rated. One crucial drawback is that, since the slanted plane then the accuracy provided by the proposed depth model, depth model defined by Eq. (5) depends on two random Eq. (12), will be similar to the accuracy provided by the orig- variables A and B (intercept and slope) and not one, the inal model, Eq. (6). quantization approach is no longer viable, for the time and memory complexity to create the SATs would exceed the 4 MASSIVE PARALLELIZATION work saved. The advantage of the new model is that it incorporates semantic cues and confidence for the disparity This section describes and discusses the massively parallel measurements, that can be used to slightly modify the mea- organization and data layouts designed for the Stixel com- surement model without significantly affecting accuracy. putation pipeline. We first start with a brief explanation of First, we will describe how to compute the optimal the performance-critical elements of a GPU architecture and parameters a; b for a given Stixel in constant time. Next we then follow with a description of the GPU-accelerated will explain the modification in Eq. (6) that allows comput- design and the analysis of the design trade-offs. ing the data term cost in constant time. Similarly to [11], when optimizing for the plane parame- 4.1 GPU Architecture and Performance ters ai ; bi of a certain Stixel si , all other optimization parame- GPUs are massively-parallel devices containing tens of ters are independent of the actual choice of the plane throughput-oriented processing units called streaming multi- parameters, and we can simplify processors (SMs). Compute and memory operations are exe- cuted as vector (SIMD) instructions and are highly argmin E ðs ; m Þ ¼ argmin Estixel ðsi ; m Þ þ Eplane ðs i Þ; pipelined in order to save energy and transistor budged. ai ;bi ai ;bi SMs can execute several vector instructions per cycle, (11) selected from multiple independent execution flows: the higher the available instruction-, vector- and thread-level and minimize the global energy function with respect to the parallelism, the better the pipeline utilization. plane parameters of all Stixels and all geometric classes inde- The CUDA programming model merges vector-level and pendently. By deriving the mathematical expressions we can thread-level parallelism, and allows defining a very large find an optimal solution of the resulting weighted least squares number of potentially concurrent execution instances problem in closed form. The calculation of the solution in con- (called threads) of the same program code. A unique two- stant time relies on the pre-computation of multiple SATs. level identifier (ThrId, CTAid) is used to specialize each 2440 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 thread for a particular data and/or function. A CTA (Cooper- downsizes the inputs, both in the horizontal and vertical ative Thread Array) comprises all the threads with the same dimensions, to produce a more compact representation and CTAid, which run simultaneously and until completion in also to reduce the computational load of the subsequent the same SM, and can share a fast but limited memory stages. Since the downsized 3D output matrix will be later space: the so-called Shared Memory. processed by columns, the output data is transposed on the The CUDA 9.0 specification introduces cooperative groups fly, stored as consecutive columns of memory (column- to dynamically organize groups of threads to perform col- wise) instead of consecutive rows of memory (row-wise). lective operations involving communication and synchroni- Fusing the downsampling and transposition stages saves zation, which enable complex patterns of parallel expensive intermediate reads and writes to global memory. cooperation. The hardware scheduler maps threads in the The left table in Fig. 4 depicts the computational analy- same cooperative group to vector instructions, which are sis of the algorithm. Each input data element must be executed efficiently, specially when the size of the group read once, and must be added to its neighbor elements to matches the hardware warp size (or number of vector provide a mean value written to the output matrix. Since lanes). the amount of input and output data on practical scenar- Each thread has its own private Local Memory space, com- ios is too large to fit into the last-level cache of a GPU, monly mapped to registers by the compiler. A large space of memory operations will be solved on the device memory. Global Memory is public to all execution instances, and is Although the theoretical arithmetic intensity (ratio of mapped into a large-capacity but long-latency device mem- abstract compute operations to memory operations) is ory, which is accelerated using a two-level hierarchy of constant, since device memory accesses are more expen- cache memories. sive than the involved compute operations, the perfor- The parallelization scheme of an algorithm and the data mance of executing this stage on a GPU will be limited by layout determine the available parallelism at the instruction, the performance of the device memory. Since there are vector and thread level and the memory access pattern. A much more memory reads than writes, this analysis large amount of parallelism is required for hiding operation encourages a thread layout aimed at maximizing the read latencies and achieving high resource usage. Additionally, bandwidth from the device memory. efficient memory performance requires that the set of The proposed parallel scheme, depicted in the upper part addresses generated by a group of consecutive threads refer of Fig. 4, distributes the average reduction of the data in to consecutive positions that can be coalesced into a single, image blocks of size s  t (where s and t are the horizontal wider memory transaction. Since the bandwidth of the and vertical Stixel resolution, respectively) to cooperative device memory is often a performance bottleneck, an effi- groups of t threads, with each group operating indepen- cient CUDA code should promote data reuse on the internal dently to calculate a single output value. Each thread first caches, the shared memory, and the registers. accumulates the values corresponding to a column of s pix- els, then the t threads in the group perform a cooperative horizontal reduction, and finally the first thread in the 4.2 Downsampling and Transpose group writes the average result in the transposed position. The input to the Stixel segmentation pipeline is a collection The reduction operation is implemented using shuffle oper- of z dense images of width W and height H (c.f. Fig. 4). The ations when t is a power of two, or else using Shared first image contains the disparity for each pixel, the second Memory. image holds the disparity confidence, and the remaining The CTA size have been set to 256 threads and the SM images contain the probabilities corresponding to each occupancy is 70 percent, limited by the available register semantic class. The first stage in the algorithm pipeline storage. However, the performance bottleneck has been Fig. 4. Downsampling and Transpose: computational analysis and parallel scheme. The width and height of the input images are W and H. There are z channels corresponding to disparity map, disparity confidence and semantic probabilities. The horizontal and vertical Stixel resolution is defined by s and t, respectively. Accordingly, h ¼ Hs and w ¼ Wt . HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2441 empirically measured to be the read bandwidth to the The prefix sum is computed on the Shared Memory and device memory, which approaches between 70 and 90 per- involves a cooperative parallel pattern, requiring communi- cent of the peak bandwidth. The most important perfor- cation and synchronization. We use the parallel scan algo- mance issue is to make consecutive threads (from the same rithm proposed by Harris et al. [28], but modify the original group and from consecutive groups) to read data from con- implementation using register-to-register shuffle instruc- secutive pixels (row-wise) from the device memory, and tions, in order to afford Shared Memory reads and writes. then promote the coalescing of memory read operations. The collective prefix sum operations involve log2 ðhÞ extra computation steps with respect to the serial computation. 4.3 Computation of Summed Area Tables Using less than h threads improves the work-efficiency of As explained in Section 3, our implementation makes exten- the algorithm, and using warpsz = 32 threads is the best sive use of Summed Area Tables (SATs). There are five option, thanks to the fast hardware support for synchroni- SATs per input image column for computing the plane zation and communication at the warp level. However, in parameters a and b determining the depth model of a stixel, practice, since the prefix sum stage involves a very small and three additional SATs for computing the depth term percentage of the total computation load, we do not see any cost for the three structural classes. The semantic term cost performance difference. requires one SAT per semantic class; we use 18 classes in The GPU implementation of the original Stixels, [4], used the experiments presented in Section 5. a very large SAT that did not fit into the Shared Memory The generation of each SAT involves three steps: (1) read and provoked a large amount of accesses to the device values from the input 3D matrix generated in the previous memory that reduced the performance. The proposals pipeline stage; (2) compute the data terms and storing them described in Section 3.3 and Section 3.4 to implement the in a LUT; and (3) calculate the prefix sum of the LUT to gen- Slanted Stixels model reduce the memory requirements for erate the SAT. The arithmetic intensity is constant but rela- the SATs and allows storing them completely into the tively high (c.f. left table on Fig. 5), due to the expensive Shared Memory. mathematical operations that are applied to compute the data term costs. 4.4 Dynamic Programming Stage There are several parallel configurations that are efficient The Dynamic Programming (DP) computation stage, both for this stage. However, we opt to fuse this stage and the fol- on the original model of Slanted Stixels [12] and our pro- lowing stages into the same CUDA kernel in order to save posal, has the higher computational complexity (c.f. left the intermediate reads and writes to global memory, and table on Fig. 5), and for practical cases is the most time-con- reuse data on the Shared Memory. Then, the best thread suming step. Our proposed design exploits the locality of configuration is determined by the computational character- the data accesses to move most memory accesses to the istics of the next stage. Shared Memory and the Local Memory, making the arith- The proposed parallel scheme (upper-left part of Fig. 5) metic intensity proportional to h and, therefore, we can consists of w cooperative groups of h threads. Each coopera- ignore the device memory accesses. However, this stage is tive group generates the 26 SATs corresponding to one the most elusive for massive parallelization. image column and stores them into the Shared Memory. The parallel processing of each input column is simple, Threads read the input (transposed) column data from but not enough to efficiently exploit current GPUs for the Global Memory in a coalesced way, compute the LUT val- image sizes considered in typical applications. The chal- ues in parallel and write them to Shared Memory. lenge is to extract fine-grain parallelism when processing Fig. 5. Stages for generating SATs, Dynamic Programming (DP) computation, and Backtracking, which are fused into the same CUDA kernel. The computational analysis and parallel scheme are shown for the original proposal [12] and our proposal. w, h, z are defined in Fig. 4. vbottom , vtop , dslope , dintercept , and C define the Stixel properties: bottom & top row, depth model, and semantic class. 2442 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 each column, since there are data dependencies and irregu- percent. Moving some data to Global Memory releases lar parallelism that complicate the task. To this end, we space on the Shared and Local Memory and increases the assign a Cooperative Thread Array (CTA) of h threads to a potential thread-level parallelism, but results in a much DP task associated with each column (see Fig. 5). higher instruction count and performance becomes limited The DP recurrence shown in Eq. (9) defines how to calcu- by the device memory latency. late the minimum cost of a problem OBk with k pixels using The computational complexity of the original model of the results computed for smaller problems. The most Slanted Stixels [11], [12] is higher (Oðw  h3 Þ) than that of our straightforward parallel design option (A) is to use k + 1 of novel depth model (Oðw  h2 Þ), described in Section 3.4. Also, the CTA threads to cooperatively compute OBk (and GRk since the original algorithm computes the cost of a Stixel with and SK k ) for each problem size k (0  k < h). An alterna- linear time complexity on the Stixel height, it suffers from a tive option (B) is to assign each CTA thread, i, the task of high load unbalance, which reduces the effective parallelism computing OBi (and GRi and SK i ). Both parallel schemes, and, therefore, the utilization of the GPU resources. A and B: (1) do not balance the computation work evenly; and (2) involve data dependencies that reduce parallelism 4.5 Backtracking and Data Compaction and require additional synchronization. The backtracking step is an inherently sequential process for The first parallel design (A) starts using a single thread each column. As described in Section 3.2, the program navi- and increases the number of running threads progressively. gates back on an index table created during the DP solving Each step requires a cooperative parallel minimum opera- stage and produces a variable-size list of Stixels, c.f. Fig. 5. tion. Option B starts using h threads and decreases the num- The lack of parallelism seems to discourage a GPU imple- ber of active threads on every step of the DP solving mentation, but the time to transfer the index tables to the process. Each step involves a broadcast of the cost values CPU, or even from Shared Memory to Global Memory, is computed by the running thread with minimum identifier. higher than the time to perform the task on the GPU (less This last option is the one selected, and depicted on Fig. 5. than 0.5 percent of the overall execution time). Both parallel options involve multiple reads to consecutive As shown in Fig. 5, we fuse the backtracking stage with positions or to the same position on the 26 SATs. As explained the two previous stages into the same CUDA kernels. The before, the SATs are stored into Shared Memory, which pro- CTA threads copy the index table from local registers to vides efficient accesses. Only option B allows each thread to Shared Memory (reusing the space devoted to the SATs, not hold into the thread-private registers (or Local Memory) its cor- needed in the backtracking stage), and then a single thread responding portion of both the cost table and the index table. processes the index table and generates the final output. A In each iteration, the thread with minimum identifier computes fixed (and conservatively large) amount of Global Memory the final value in the corresponding cost table entry, and then is allocated per column to hold the variable-size lists of Stix- uses the Shared Memory to broadcast that value; a barrier is els. A final and very fast execution kernel is used to compact used to enforce the required synchronization; and finally, that the information into a contiguous region of Global Memory. thread becomes idle. Option A is slower because it requires more synchronization operations and data movements. The performance of this stage is latency-bounded due to 5 EXPERIMENTS the lack of parallelism. There are three causes for the limited This section assesses the accuracy and performance of our parallelism: (1) the relatively high memory requirements on proposal. We first verify that our method maintains the the Shared and Local Memory; (2) the decreasing amount of same accuracy level as the previous Slanted Stixel model independent work as the recurrence loop advances; and (3) [12]. For that purpose, we evaluate on both synthetic and the synchronization barriers between recurrent steps, which real data, c.f. Section 5.1.1, and report quantitative and qual- reduce the effective parallelism. itative results, c.f. Section 5.1.3. We also show and discuss Specifically, each thread holds an average of 26 float the performance advantage of our novel depth model for numbers in the Shared Memory and uses 79 local registers. GPUs and show quantitative results, c.f. Section 5.2. The best thread block configuration contains 256 threads, which requires 19.75 Ki registers and 26 KiB of Shared Memory per block. Both Pascal and Volta CUDA architec- 5.1 Accuracy and Compression Experiments tures provide 64 Ki registers per SM (see Table 4), which 5.1.1 Datasets pose a limit of three 256-thread blocks per SM We use two datasets with real images, KITTI 2015 [29], [30] (3  19:75 ¼ 59:25 Ki registers out of 64). However, the Pas- and Cityscapes [21], and one dataset with synthetic images, cal architecture only provides 64 KiB of Shared Memory, SYNTHIA-SF (SYNTHIA San Francisco) [11]. which allows allocating two blocks of threads, while the The well-known stereo challenge KITTI 2015 contains Volta architecture provides 96 KiB of Shared Memory, and images with sparse disparity ground truth obtained from a allows reaching the limit of 3 blocks of threads. Overall, the laser scanner and semantic segmentation ground truth. City- maximum GPU occupancy is 512 threads out of 2048 (25 scapes is a highly complex dataset with dense annotations of percent) in a Pascal GPU, and 768 out of 2048 threads (37.5 19 classes. SYNTHIA-SF (SYNTHIA San Francisco) is a syn- percent) in a Volta GPU. The effective average GPU occu- thetic dataset that consists of photo-realistic frames rendered pancy is almost half of the peak values, due to the reduction from a virtual city, with precise pixel-level depth and semantic of parallelism in the algorithm (2), and the synchronization annotations. operations (3). Even with this hard limitations, the GPU We evaluate depth accuracy on the training images of computation cores have an utilization between 30 and 50 KITTI (200) and all the images of SYNTHIA-SF (2224). The HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2443 TABLE 2 of the disparity estimates. We generate back the dense dis- Data Compression Measured in Pixels Per Stixel of Our Method parity image from the segmentation obtained from our and [12] method and from [12]. Then, an outlier is a disparity estima- Stixel resolution: 8  8 Stixel resolution: 4  4 tion with an absolute error larger than 3 px or a relative deviation larger than 5 percent compared to the ground Dataset Slanted [12] Ours Slanted [12] Ours truth. KITTI 2015 587 572 254 242 The semantic accuracy is evaluated as the average Inter- Cityscapes 1105 877 475 331 section-over-Union (IoU) over all 19 classes, which is also a SYNTHIA-SF 1439 1379 637 606 standard measure for semantics [33]. We evaluate on three datasets: KITTI 2015 [29], Cityscapes [21] and SYN- Data compression is measured as the average number of THIA-SF [11], c.f. Section 5.1.1. pixels per Stixel, and quantifies the complexity of the obtained representation. semantic accuracy is measured on KITTI, the validation images of Cityscapes (500), and SYNTHIA-SF. 5.1.3 Results We want to remark that the datasets chosen for evalua- The quantitative results of our proposal and baseline as tion are not representative of the world, i.e. lacking some described in Section 5.1.2 are shown in Tables 2 and 3 . very common kind of pedrestians such as women with The first observation, taken from Table 2, is that both var- strollers, children, pets and disabled people. Due to the cost iants provide compact representations of the surrounding, of acquiring more data we evaluate our methods with these with a compression larger than 100 compared to the high datasets but results can not be considered conclusive. resolution input images. Segmentations generated by our proposal are around 5 percent less compact for KITTI and SYNTHIA-SF, and 36 and 45 percent less compact for City- 5.1.2 Experiment Details scapes, which is the more complex scenario. Slanted Stixels [12] serves as baseline for the comparison Second, results from Table 3 indicate that our method with our proposal, c.f. Section 3, because it represents the achieves very similar accuracy results on all datasets, with state-of-the-art results in terms of Stixel accuracy. an increase of less than 3.5 percent of the disparity error As input, we use disparity maps obtained via semi- and a decrease of less than 1 percent of the IoU. We consider global-matching (SGM) [31] and pixel-level semantic labels that this slight degradation of the accuracy results is a small computed by a fully convolutional network (FCN) [32]. The price to pay for the speed improvement that we will show parameters are taken from [3] for fair comparison. For the on the next section. same reason, we use the same FCN architecture from [3]. Finally, a 4 higher Stixel resolution (4  4 versus 8  8) However, FCN weights are not the same, but accuracy decreases the compression of the representation between 2 (IoU) of the input is provided for reference. Similarly, SGM and 3 times but, in return, is more accurate. Note that the implementation details differ between our implementation disparity error of the compact representation is lower than and those of previous works, therefore, input disparity the one of the input generated by the SGM algorithm. The accuracy is also provided in the respective tables. Stixel world model helps removing input noise thanks to Following [12], we use the known camera calibration to the joint inference of semantic and depth data. obtain expected maground and mbground . We assume that objects The observations from the quantitative evaluation are are vertical and set s bobject ! 0; mbobject ¼ 0, because the dis- confirmed also in the qualitative results, c.f. Fig. 6. parity is too noisy for the slanted object model. Sky Stixels Fig. 7 shows the distribution of the disparity errors on the are assumed to be vertical and very far. individual images of the KITTI dataset when using both We use three metrics to evaluate our proposed method in Slanted Stixels and our model. The distributions of the errors terms of depth and semantic accuracy, and also in terms of are very similar, with most of the images containing a small data compression. number of errors and some images containing more errors. The depth accuracy is defined as the same standard met- The last chart shows a histogram of the relative differences ric used to evaluate on KITTI [29], which is the outlier rate of the disparity error for each individual image. The TABLE 3 Accuracy of Our Method Compared to Slanted Stixels [12], Input SGM, and FCN Input Stixel resolution: 8  8 Stixel resolution: 4  4 Metric Dataset SGM FCN Slanted Stixels [12] Ours Slanted Stixels [12] Ours Disp Error (%) KITTI 15 8.51 - 8.53 8.72 7.81 7.93 SYNTHIA-SF 10.26 - 8.55 8.85 7.56 7.83 IoU (%) KITTI 15 - 44.51 43.97 43.63 44.54 44.23 Cityscapes - 68.22 66.87 66.75 67.92 67.78 SYNTHIA-SF - 34.01 33.41 33.39 33.82 33.83 We evaluate on three datasets: KITTI 15 [29], Cityscapes [21] and SYNTHIA-SF [11] using these metrics: Disparity Error (lower is better) and Intersection over Union (higher is better), c.f. Sections 5.1.1, and 5.1.2. Ours is detailed in Section 3.4. Fig. 6. Exemplary outputs on real data (KITTI 2015): RGB Image (left), Slanted Stixels [12] (center) and Our Stixels (right) representations are shown. Color encodes the distance from close (red) to far (green). As we can see, the representation of our proposal is visually similar to the baseline. Fig. 7. Histogram of the average disparity errors per image (KITTI 15 dataset [29] and 4  4 Stixel resolution) when using (a) Slanted Stixels [12] and (b) our model; and (c) histogram of the relative differences of the disparity error for each individual image (Slanted Stixels - Ours, meaning that nega- tive numbers indicate the Slanted Stixels outperforms. differences in the disparity error are relatively small; some- disparity map from stages (stereo matching and FCN) that are times they benefit our proposal and to a somewhat greater expected to be both executed on the GPU (e.g. using SGM extent they harm our proposal. implemented on the GPU [34]). The list of Stixels generated by the computation could be post-processed on the GPU, or is small enough to discard the time for transferring the data to 5.2 Performance Experiments the CPU memory. The main goal of our performance analysis is to evaluate run- time and efficiency on embedded devices such as the NVIDIA Tegra X2 and Tegra Xavier c.f. Table 4. All the metrics are mea- 5.2.1 Results sured using NVIDIA performance tools. We assume the input Figs. 8a and 8b show the performance throughput (frames per data for Stixel segmentation is already into the GPU memory, second, or fps) of the CUDA code implementing both Slanted and we do no add the time for moving these data from the Stixel models (original and ours), on two GPU systems, using CPU memory: Stixel estimation is just a stage in a computer- an image size of 2048  1024 and for both 8  8 and 4  4 vision pipeline that receives the semantic segmentation and Stixel resolution. It is remarkable that real-time rates higher HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2445 TABLE 4 A multi-threaded implementation of the original Slanted Specifications of GPUs Employed in Our Experiments Stixels model was evaluated for the same images and a Stixel resolution of 8 pixels, and reached 6.6 fps on a six-core Intel i7- Tegra X2 Tegra Xavier 6800K CPU [12]. Our implementation reaches 344.3 fps on a architecture Pascal Volta Tegra Xavier, i.e. more that 50 times faster, with a reduced cost clock frequency 1465 MHz 1377 MHz and power envelop (TDP of 30 Watt compared to 140 Watt). number of SMs 4 8 number of cores 256 512 The alternative over-segmentation variant that selects promis- registers per SM 64 Ki 64 Ki ing Stixel cuts by means of a FCN runs at an average of 27.5 fps shared memory per SM 64 KiB 96 KiB on the same six-core CPU. This approach has the drawback device memory size 8 GB 32 GB that the execution time is dependent on the image content (not device memory bandwidth 59.7 GB/s 136.5 GB/s predictable), with a worst-case scenario were all the Stixel cuts L2 cache size 4096 KiB 6144 KiB must be evaluated and runs even slower than the original ver- GFLOPS (single precision) 750 1410 sion. Our GPU-accelerated design runs 12 times faster and TDP 15 W 30 W with predictable times. Fig. 9 presents a breakdown of the elapsed time and the IPC (ratio of machine instructions executed per clock cycle and per SM). For low Stixel resolution (8  8), the time for the common Down-sampling & Transpose stage represents a substantial portion of the total time: 62 percent on the Tegra Xavier, and 47 percent on the Tegra X2. The performance bottleneck of this stage is the GPU memory bandwidth and its execution time is proportional to the size of the original and down-sampled images. Increasing the Stixel resolution makes the time of the Dynamic Programming stage to domi- nate on both GPUs, since the computational complexity of the DP stage with respect to the image height after down- sampling is quadratic for our proposal, while it is cubic for the original Slanted Stixels proposal. Considering only the Dynamic Programming stage, our pro- posal executes 7.2 times (8  8) and 10.9 times (4  4) faster than [12] on the Tegra Xavier. Our implementation achieves higher IPC ratios (1.31 and 1.33 better) on each SM, which means that our approach exhibits more parallelism. But most of the speedup is due to a 5.8 and 8.7 reduction on the total number of machine instructions executed by the GPU (these data can be derived from the results in Fig. 9). This corroborates the better algorithmic scalability of our approach. We now assess the performance differences when using both GPUs. The Tegra Xavier contains 8 Volta SMs (512 cores) running at a slightly lower clock frequency than the 4 Pascal SMs (256 cores) in the older Tegra X2 (c.f. Table 4). This means a potential raw performance advantage of 1:88. The speedup on the execution time of the Dynamic Programming stage is around 6.5 times for both resolutions, which means that our implementation is using the Volta cores more efficiently than Fig. 8. Frame-rate of our method compared to Slanted Stixels [12] for the Pascal cores, partially due to a higher GPU occupancy (see 2048  1024 image resolution on the NVIDIA Tegra Xavier and Tegra X2 Section 4.4), which improves the IPC ratio (from 1.34 to embedded GPUs. 1.48), and partially due to a better low-level codification effi- ciency (less machine instructions to implement the same basic than 90 fps are achieved by the Tegra Xavier GPU for both operations) of around 2.2 times (derived Fig. 9). The speedup Stixel resolutions (see Fig. 8a), and even the older Tegra X2 of the Down-sampling & Transpose stage is around 3.6, closer GPU is able to achieve real-time rates for Stixel resolution of but higher than the 2.3 improvement on the device memory 8  8 pixels (76 fps, c.f. Fig. 8b). On the Tegra Xavier and for a bandwidth (from 59.7 to 136.5 GB/s). Therefore, our proposal low Stixel resolution (8  8), our method is 3:4 faster than achieves very good scalability when ported to the new GPU Slanted Stixels [12] (344.3 fps vs 102.4 fps). A higher resolution Xavier architecture. (4  4) provides more accurate segments (c.f. Table 3), but while our proposal achieves practical frame rates (92.3 fps), the implementation following [12] is 8:5 slower, which 6 CONCLUSION impedes real-time execution. In fact, the bigger the problem We have described and assessed the performance of the first size (either increasing the Stixel resolution or the image size), GPU-accelerated implementation of Slanted Stixels and we the higher the advantage of our proposal. show that our algorithmic proposal is efficient for GPU 2446 IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. 32, NO. 10, OCTOBER 2021 Fig. 9. Breakdown of low-level performance metrics for the two main stages of our method: Down-sampling & Transpose stage (common) and Dynamic Programming stage (ours vs. [12]). IPC is the ratio of machine instructions executed per clock cycle and per SM. Image size is 2048  1024. parallelization. Our proposal achieves real-time perfor- REFERENCES mance for realistic problem sizes, proving that the low- [1] H. Badino, U. Franke, and D. Pfeiffer, “The stixel world - A com- power envelope and remarkable performance of embedded pact medium level representation of the 3D-world,” in Proc. 31st CPU-GPU hybrid systems make them good target platforms DAGM Symp. Pattern Recognit., 2009, pp. 51–60. for most real-time image processing tasks. [2] D. Pfeiffer and U. Franke, “Towards a global optimal multi-layer stixel representation of dense 3D data,” in Proc. British Mach. Vis. The reformulation of the measurement depth model pro- Conf., 2011, pp. 1–12. [Online]. Available: http://dx.doi.org/ posed for Slanted Stixels improves the performance and scal- 10.5244/C.25.51 ability of the original proposal, while slightly reducing [3] L. Schneider et al., “Semantic stixels: Depth is not enough,” in precision. However, in a real environment with run time limi- Proc. IEEE Intell. Veh. Symp., Gotenburg, Sweden, 2016, pp. 110–117. [Online]. Available: http://dx.doi.org/10.1109/IVS.2016.7535373 tations, the shorter execution time with respect to the original [4] D. Hernandez-Juarez , A. Espinosa, J. C. Moure, D. Vazquez, and proposal allows to increase the resolution of the stixels and A. M. L opez, “GPU-accelerated real-time stixel computation,” in then improve the overall accuracy of the segmentation process. Proc. IEEE Winter Conf. Appl. Comput. Vis., Santa Rosa, CA, USA, 2017, pp. 1054–1062. [Online]. Available: https://doi.org/ Compared to the over-segmentation proposal, our approach is 10.1109/WACV.2017.122 more accurate, faster and with predictable run-times. [5] R. Benenson, R. Timofte, and L. J. V. Gool, “Stixels estimation The proposed parallel scheme and data layout for the without depth map computation,” in Proc. IEEE Int. Conf. Comput. irregular computational pattern corresponding to the Vis. Workshops, Barcelona, Spain, 2011, pp. 2010–2017. [Online]. Available: http://dx.doi.org/10.1109/ICCVW.2011.6130495 Dynamic Programming stage follows general optimization [6] M. Cordts, L. Schneider, M. Enzweiler, U. Franke, and S. Roth, rules based on a simple GPU performance model. We have “Object-level priors for stixel generation,” in Pattern Recognit. 36th shown that the parallel implementation scales from a previ- German Conf., M€ unster, Germany, 2014, pp. 172–183. ous generation embedded GPU system to a new generation [7] M. Cordts et al., “The stixel world: A medium-level representation of traffic scenes,” Image Vis. Comput., vol. 68, pp. 40–52, 2017. GPU, and we expect it to scale gracefully on the forthcom- [Online]. Available: http://www.sciencedirect.com/science/ ing GPU architectures. Our parallelization strategy is gen- article/pii/S0262885617300331 eral enough to be applied to similar Dynamic Programming [8] O. Ignat, “Disparity image segmentation for free-space detection,” computational patterns, where parallelism may decrease in Proc. IEEE 12th Int. Conf. Intell. Comput. Commun. Process. (ICCP), Sept. 2016, pp. 217–224. along the processing task. [9] D. Levi, N. Garnett, and E. Fetaya, “StixelNet: A deep convolu- tional network for obstacle detection and road segmentation,” in Proc. British Mach. Vis. Conf., Swansea, U.K., 2015, pp. 109.1–109.12. ACKNOWLEDGMENTS [Online]. Available: http://dx.doi.org/10.5244/C.29.109 The contribution of Daniel Hern andez, Antonio Espinosa, [10] D. A. P. Carrillo and A. Sutherland, “Fast obstacle detection using sparse edge-based disparity maps,” in Fourth Int. Conf. 3D Vis. and Juan C. Moure was focused on the GPU acceleration of (3DV), 2016, pp. 66–72. stixels. The contribution of Daniel Hern andez, David [11] D. Hernandez-Juarez et al., “Slanted stixels: Representing San Vazquez, and Antonio M. L opez was focused on the slanted Francisco’s steepest streets,” in Proc. Mach. Vis. Conf., London, U.K., stixels concept itself from the computer vision viewpoint. 2017, pp. 87.1–87.12. [12] D. Hernandez-Juarez et al., “Slanted stixels: A way to repre- As CVC member, Antonio M. L opez also thank the General- sent steep streets,” Int. J. Comput. Vis., vol. 127, pp. 1643–1658, itat de Catalunya for its CERCA Program and the ACCIO 2019. [Online]. Available: https://doi.org/10.1007/s11263–019- agency. The work of Daniel Hern andez, Antonio Espinosa, 01226-9 [13] T. M. Hehn, J. F. P. Kooij, and D. M. Gavrila, “Instance stixels: Seg- and Juan C. Moure was supported by the Ministerio de menting and grouping stixels into objects,” in Proc. IEEE Intell. Economıa, Industria y Competitividad under Project Veh. Symp (IV), Paris, France, 2019, pp. 2542–2549. [Online]. TIN2017-84553-C2-1-R. The work of Antonio M. L opez was Available: https://doi.org/10.1109/IVS.2019.8814243 supported in part by the Project TIN2017-88709-R [14] R. Benenson, M. Mathias, R. Timofte, and L. J. V. Gool, “Fast stixel computation for fast pedestrian detection,” in Proc. Comput. Vis. - (MINECO/AEI/FEDER, UE), and in part by ICREA under ECCV 2012 Workshops Demonstrations, Florence, Italy, 2012, the ICREA Academia Programme. pp. 11–20. HERNANDEZ-JUAREZ ET AL.: 3D PERCEPTION WITH SLANTED STIXELS ON GPU 2447 [15] M. Enzweiler, M. Hummel, D. Pfeiffer, and U. Franke, “Efficient Daniel Hernandez-Juarez received the BSc stixel-based object recognition,” in Proc. IEEE Intell. Veh. Symp., degree in computer science from Universitat Alcala de Henares, Madrid, Spain, 2012, pp. 1066–1071. [Online]. Autonoma de Barcelona (UAB), in 2014, and the Available: http://dx.doi.org/10.1109/IVS.2012.6232137 MSc degree in computer vision from the UAB, [16] F. Erbs, B. Schwarz, and U. Franke, “Stixmentation - probabilistic UPC, UPF, and UOC, in 2015. He received the stixel based traffic scene labeling,” in Proc. British Mach. Vis. Conf., PhD degree in computer vision, in 2020, from the Surrey, U.K., 2012, pp. 71.1–71.12. [Online]. Available: https:// UAB, he was supervised by Dr. Juan Carlos doi.org/10.5244/C.26.71 Moure and Dr. David Vazquez. His focus is to [17] M. Muffert, N. Schneider, and U. Franke, “Stix-fusion: A probabi- improve 3D perception algorithms and adapt listic stixel integration technique,” in Proc. Can. Conf. Comput. them to GPU devices. He is currently immersed Robot Vis., Montreal, QC, Canada, 2014, pp. 16–23. [Online]. Avail- in the study of the Stixel World. His research able: http://dx.doi.org/10.1109/CRV.2014.11 interests include self-driving cars, deep learning [18] D. Pfeiffer and U. Franke, “Modeling dynamic 3D environments techniques, 3D perception, and embedded by means of the stixel world,” IEEE Intell. Transport. Syst. Mag., systems. vol. 3, no. 3, pp. 24–36, Fall 2011. [Online]. Available: https://doi. org/10.1109/MITS.2011.942207 [19] T. Scharw€ achter, M. Enzweiler, U. Franke, and S. Roth, Antonio Espinosa is currently an associate pro- “Stixmantics: A medium-level model for real-time semantic scene fessor with the Computer Architecture and Oper- understanding,” in Proc. Comput. Vis. - ECCV 2014–13th Eur. Conf., ating Systems Department, Universitat Autonoma Zurich, Switzerland, 2014, pp. 533–548. de Barcelona. During the last ten years, he has [20] X. Li et al., “A new benchmark for vision-based cyclist detection,” participated in several European and national in Proc. IEEE Intell. Veh. Symp. IV, Sweden, 2016, pp. 1028–1033. projects related to computer science, high-perfor- [Online]. Available: https://doi.org/10.1109/IVS.2016.7535515 mance computing, and computational accelerator [21] M. Cordts et al., “The cityscapes dataset for semantic urban scene systems in collaboration with a number of compa- understanding,” in Proc. IEEE Conf. Comput. Vis. Pattern Recognit., nies and research institutions. Las Vegas, NV, USA, 2016, pp. 3213–3223. [Online]. Available: http://dx.doi.org/10.1109/CVPR.2016.350 [22] V. Dhiman, A. Kundu, F. Dellaert, and J. J. Corso, “Modern MAP David Vazquez is currently a fundamental inference methods for accurate and fast occupancy grid mapping research scientist with Element AI, where he on higher order factor graphs,” in Proc. IEEE Int. Conf. Robot. Auto- works on computer vision. He was a postdoctoral mat., 2014, pp. 2037–2044. researcher with Computer Vision Center of Bar- [23] D. Nuss, T. Yuan, G. Krehl, M. Stuebler, S. Reuter, and K. Die- celona (CVC) and Montreal Institute of Learning tmayer, “Fusion of laser and radar sensor data with a sequential Algorithms (MILA) and an assistant professor monte carlo Bayesian occupancy filter,” in Proc. IEEE Intell. Veh. with the Department of Computer Science, Symp. (IV), 2015, pp. 1074–1081. Autonomous University of Barcelona. He is an [24] S. Thrun, “Robotic mapping: A survey,” in Exploring Artificial expert in machine perception for autonomous Intelligence in the New Millenium, G. Lakemeyer and B. Nebel, vehicles and on domain adaptation from simula- Eds. Amsterdam, The Netherlands: Morgan Kaufmann, tion to real-world environments. 2002. [25] D. Pfeiffer, S. Gehrig, and N. Schneider, “Exploiting the power of stereo confidences,” in Proc. IEEE Conf. Comput. Vis. Pattern Recog- Antonio M. Lo  pez is the principal investigator nit., Portland, OR, USA, 2013, pp. 297–304. [Online]. Available: http://dx.doi.org/10.1109/CVPR.2013.45 with Autonomous Driving Lab, Computer Vision Center, University Auto noma de Barcelona [26] D. Pfeiffer, “The stixel world - A compact medium-level represen- tation for efficiently modeling three-dimensional environments,” (UAB). He is also an associate professor with Ph.D. dissertation, Math. Nat. Sci. Faculty II, Hu Univ. Berlin, Ber- Computer Science Department, UAB. His lin, Germany, 2014. research interests include computer vision, com- [27] P. Viola and M. Jones, “Rapid object detection using a boosted cas- puter graphics, machine learning, and autono- cade of simple features,” in Proc. IEEE Conf. Comput. Vis. Pattern mous driving. He has been deeply involved in the Recognit., 2001, pp. I–I. creation of the SYNTHIA dataset and the CARLA [28] M. Harris, S. Sengupta, and J. D. Owens, “Parallel prefix sum open-source simulator, both for democratizing (scan) with CUDA,” GPU Gems, vol. 3, no. 39, pp. 851–876, 2007. autonomous driving research. He is actively [29] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for Autono- working hand-on-hand with industry partners to bring state-of-the-art mous Driving? The KITTI Vision Benchmark Suite,” in Proc. IEEE techniques to the field of autonomous driving. He is granted by the Cata- Conf. Comput. Vis. Pattern Recognit., 2012, pp. 3354–3361. lan ICREA Academia Program. [30] H. A. Alhaija, S. K. Mustikovela, L. Mescheder, A. Geiger, and C. Rother, “Augmented reality meets deep learning for car ins- tance segmentation in urban scenes,” Int. J. Comput. Vis., vol. 126, Juan C. Moure is currently an associate profes- pp. 961–972, 2018. sor with Computer Architecture and Operating [31] H. Hirschm€ uller, “Stereo processing by semiglobal matching Systems Department, Universitat Autonoma of and mutual information,” IEEE Trans. Pattern Anal. Mach. Intell., Barcelona (UAB), where he teaches computer vol. 30, no. 2, pp. 328–341, Feb. 2008. [Online]. Available: http:// architecture, performance engineering, and dx.doi.org/10.1109/TPAMI.2007.1166 parallel programming. He is the author of more [32] J. Long, E. Shelhamer, and T. Darrell, “Fully convolutional net- than 50 papers, and has participated in sev- works for semantic segmentation,” in Proc. IEEE Conf. Comput. Vis. eral European and Spanish projects related to Pattern Recognit., Boston, MA, USA, 2015, pp. 3431–3440. [Online]. high-performance computing. His current Available: http://dx.doi.org/10.1109/CVPR.2015.7298965 research interests include massive parallel [33] M. Everingham, S. M. A. Eslami, L. J. V. Gool, C. K. I. Williams, architectures, programming, and algorithms, J. M. Winn, and A. Zisserman, “The pascal visual object classes mainly focused on computer vision, signal processing, and bioinfor- challenge: A retrospective,” Int. J. Comput. Vis., vol. 111, no. 1, matics applications. pp. 98–136, 2015. [Online]. Available: http://dx.doi.org/10.1007/ s11263–014-0733-5 [34] D. Hernandez-Juarez et al., “Embedded real-time stereo estima- " For more information on this or any other computing topic, tion via semi-global matching on the GPU,” in Proc. Int. Conf. Com- please visit our Digital Library at www.computer.org/csdl. put. Sci., 2016, pp. 143–153. [Online]. Available: http://dx.doi. org/10.1016/j.procs.2016.05.305