A New Direct Connected Component Labeling And Analysis . - Nvidia

A new Direct Connected Component Labelingand Analysis Algorithm for GPUsArthur Hennequin1,2 , Lionel Lacassagne1LIP6, Sorbonne University, CNRS, FranceLHCb experiment, CERN, Switzerland 21GTC 2019 March 21st1 / 29

What are Connected Component Labeling and Analysis ?Connected Components Labeling (CCL) consists in assigning a unique number(label) to each connected component of a binary imageConnected Components Analysis (CCA) consists in computing some featuresassociated to each connected component like the bounding box [xmin ,xmax ] x[ymin ,ymax ], the sum of pixels S, the sums of x and y coordinates Sx, Sy1gray level imagebinary level image(segmentation by(motion detection)2connected componentlabelingconnected componentanalysis seems easy for a human being that has a global view of the image but, ill-posed problem: the computer has only a local view around a pixel(neighborhood) important in computer vision for pattern recognition, motion detection .2 / 29

Two classes of CCL algorithms multi-pass iterative algorithmsIIIcompute the local positive min over a 3 3 neighborhooduntil stabilization : the number of iterations depends on the datanot predictable, nor suited for embedded systems two-pass direct algorithmsIIIIfirst pass temporary label creation and equivalence buildingneed an equivalence table to memorize the connectivity between labelsthen transitive closure of the tree associated to the equivalence tablesecond pass label relabeling on CPU, scalar algorithms are all direct and can be parallelized on SIMD CPU, until 2019, all SIMD algorithms are iterative, except 1 on GPU, until 2018, all algorithms are iterative, except 3Why so few direct algorithms on GPU and SIMD ? because extremely complex to design (not suited for SIMD nor GPU)3 / 29

Direct algorithms are based on Union-Find structureAlgorithm 1: Rosenfeld labeling algorithmAlgorithm 2: Find(e,T )for i 0 : h 1 dofor j 0 : w 1 doif I [i][j] 6 0 thene1 E [i 1][j]e2 E [i][j 1]if (e1 e2 0) thenne ne 1ex neelser1 Find(e1 , T )r2 Find(e2 , T )ex min (r1 , r2 )if (r1 6 0 and r1 6 ex ) then T [r1 ] exif (r2 6 0 and r2 6 ex ) then T [r2 ] exwhile T [e] 6 e doe T [e]elseex 0E [i][j] exreturn e // the root of the treeAlgorithm 3: Union(e1 ,e2 ,T )r1 Find(e1 , T )r2 Find(e2 , T )if (r1 r2 ) thenT [r2 ] r1elseT [r1 ] r2Algorithm 4: Transitive Closurefor i 0 : ne doT [e] T [T [e]]Parallel algorithms do: sparse addressing scatter/gather SIMD instructions (AVX512/SVE) concurrent min computation recursive atomic min instruction (CUDA)4 / 29

Classic direct algorithm: Rosenfeld (1966)Rosenfeld algorithm is the first 2-pass algorithm with an equivalence table when two labels belong to the same component, an equivalence is created andstored into the equivalence table T for example, there is an equivalence between 2 and 3 (stair pattern) and between 4and 2 (concavity pattern) stair and concavity are the only two patterns generator of equivalence here, background in gray and foreground in white11111000101100010111011111111111binary image of pixels1111100010440004032203222222image of labelse 0 1 2 3 4T[e] 0 1 2 2 2equivalence table2222predecessorpixelsp1predecessorlabelse1p2 pxe2 excurrent pixelcurrent labelimage of pixelsimage of labels111110001022000202220222222212 ex121 1 exstairconcavitypatterns generatorof equivalence2222image of labels after relabeling3124equivalence trees5 / 29

Parallel State-of-the-art Parallel Light Speed Labeling[1](L. Cabaret, L. Lacassagne, D. Etiemble) (2018)IIIparallel algorithm for CPUbased on RLE (Run Length Encoding) to speed up processing and savesmemory accessescurrent fastest CCA algorithm on CPU Distanceless Label Propagation[2](L. Cabaret, L. Lacassagne, D. Etiemble) (2018)Idirect CCL algorithm for GPU Playne-Equivalence[3](D. P. Playne, K.A. Hawick) (2018)IIdirect CCL algorithm for GPU (2D and 3D versions)based on the analysis of local pixels configuration to avoid unnecessary andcostly atomic operations to save memory accesses.6 / 29

Equivalence merge function & concurrency issueThe direct CCL algorithms rely on Union-Find to manage equivalences.A parallel merge operation can lead to concurrency issues:314414314411443411444142444212341244 1st example (top-left): no concurrency, T[3] 1, T[4] 1 2nd example (top-right): no concurrency, T[3] 1, T[4] 2 3rd example (bottom-left): non-problematic concurrency, T[4] 1, T[4] 1 4th example (bottom-right): concurrency issue, T[4] 1, T[4] 2II4 can’t be equal to 1 and 2 4 has to point to 1 and 2 has to point to 1 too.7 / 29

Equivalence merge function (aka recursive Union)The merge function, introduced by Playne and Hawick, solves the concurrencyissues by iteratively merging labels using atomic operationsAlgorithm 5: merge(L, e1 , e2 )while e1 6 e2 and e1 6 L[e1 ] doe1 L[e1 ]// root of e1while e1 6 e2 and e2 6 L[e2 ] doe2 L[e2 ]// root of e2while e1 6 e2 doif e1 e2 then swap(e1 , e2 )e3 atomicMin(L[e1 ], e2 )if e3 e1 then e1 e2else e1 e3// recursive minBy definition, e3 L[e1 ], so: if e3 e1 : no concurrent write, update of L is successful, terminates the loop if e3 e1 : concurrent write, L was updated by another thread, need tomerge e3 and e28 / 29

Hardware Accelerated algorithm : HA4Analysis of state-of-the-art weaknesses: vertical borders (non-coalescent memory accesses) expensive atomic operationsAnalysis of state-of-the-art strengths: equivalence table embedded in the image (Cabaret, Playne) merge function (Komura [4] Playne) segments labeling (Light Speed Labeling) necessary condition to merge two equivalence trees (Playne)Figure 1: All possible 4 pixels configurations. Only (f) need to merge labels. (Playne)9 / 29

Hardware Accelerated: HA4The algorithm is divided into 3 kernels: strip labeling: the image is split intohorizontal strips of 4 rows. Each strip isprocessed by a block of 32 4 threads(one warp per row). Only the head ofsegment is labeled border merging: to merge the labels onthe horizontal borders between strips relabeling / features computation: topropagate the label of each segment tothe pixels or to compute the featuresassociated to the connected components10 / 29

Example – Strip labeling initialization (Step #0)The 8 8 image is divided into 2 strips of 8 4 pixels, warp size 8Initial strip labeling: only the head of each segment (start node)is labeled with an unique label equal to its linear address: L[k] k with k y width x warning: label numbering starts at 0, not 100 01 2 3 41 8183 24260 32342 483 567122 161 405 662043475462(a) Initialization11 / 29

Example – Strip labeling (Step #1)After initialization: detection of merging nodes using necessary conditions in each thread update of start nodes onlyStrips’ segments are now labeled00 01 2 3 41 85 66701 2 3 4121 020182 812263 16180 32340 32323 5643471 32542 40623 48(b) Strip labeling063 242 48762 161 405 60 01234343286403416124843201856475447542662(c) Strip labeledHere, a CC spanning over several strips is represented by 3 disjoint trees of labels12 / 29

Example – Border merging (Step #2)Same merging operations on border nodes only. All the segments are correctlylabeled. A CC spanning to several strips is represented by 1 tree.01 2 3 40 05 600 0761 02 8123 16180 32321 3261 0122 8123 16180 03234472 403 48543 48(d) Border merging01234344754086403416124843267(e) Border merged32185 6061 32342 40201 2 3 45647546232864034161248432018265647546213 / 29

Example – Re-Labeling / Analysis (Step #3)In the final step only, each start node (blue) flattens its equivalence tree to Label the image: broadcast the label to the whole segment to Analyse the image: accumulate features into global memory using atomicsexample of features associated to segment [x0 , x1 [ at line y :I00 0S x1 x0 ,1 2 3 40 01 0Sy S y0 ,0Sx 1 2 3 40 05 60700 000001 000002 0012[x1 (x1 1) (x0 (x0 1)]5 6070000000002 003 00000003 00000000 00000000 00000000001 0001 002 00000002 00000003 00000003 0000000FindRoot5662685412160484718204340 32 3426Relabeling14 / 29

Implementation details: Grid-stride loop first weakness of previous GPU algorithms is the vertical border merging: thenon-coalescent memory accesses are slower we used the grid-stride loop [5] design pattern to divide the image in stripsinstead of tileskernel Classic(width)x blockDim.x blockIdx.x threadIdx.xif x width then// do stuff.kernel Grid stride loop(width)for x threadIdx.x to width by blockDim.x do// do stuff.Benefits: thread reuse: less thread creation. Helps to amortize the cost of threadcreation/destruction thread context is preserved: the loop ensures that pixels are processed in aspecific order and allows to reuse previously computed values15 / 29

Implementation details: horizontal data exchangeAll threads working on the same row are from the same warp, CUDA Warp-LevelPrimitives [6] can be used to directly exchange data from threads registers ballot sync primitive returns a 32-bit bitmask based on the value of aboolean within each thread (1 bit per thread) shfl sync primitive exchanges a 32-bit value between any pair of threads ina warp. Each thread specifies a thread ID to read and a value to share16 / 29

Implementation details: segments each thread needs to find its distance to the segment’s start node distance to the end is also needed for features computation bitwise operations can accelerate the computation of these distances (tx thread number)01 2 3 45 67pixels 01 1 1 00 10start distance0 1 20end distance3 2 111operator start distance(pixels, tx)return clz( (pixels (32 tx))) // clz Count Leading Zerosoperator end distance(pixels, tx)return ffs( (pixels (tx 1))) // ffs Find First Set17 / 29

Implementation details: vertical data exchange classic way of optimizing memory accesses: copying data from global toshared memory shared memory is divided in 32 banks: same bank memory accesses atdifferent addresses get serialized [7]txpixelsy01234567sharedmemorypixelsy-118 / 29

Implementation details: vertical data exchange for each row, we store the bitmasks of the 32 neighbor pixels in differentbanks store: no serialization, load: lsy-119 / 29

One final optimization. two pixels directly next to each other either belong to the same segment orhave a different color we can assign a thread two pixels instead of one. 32-bit 64-bit bitmask: modified distance operators. new version: HA464tx010121 0 1 031 041 0operator start distance64(pixels, tx)b get bit tx of pixelstxb tx breturn clzll( (pixels (64 txb)))operator end distance64(pixels, tx)b get bit tx of pixelstxb tx breturn ffsll( (pixels (txb 1)))20 / 29

Benchmark of CCL and CCA algorithms random 2048x2048 (2k) images of varying density (0% - 100%), granularity(1 - 16, granularity 4 close to natural image complexity) percolation threshold: transition from many smalls CCs to few larges CCsII8C: density 45%4C: density 64%21 / 29

Comparison of CCL algorithms on Jetson TX2 comparison with 2state-of-the-art algorithms[Playne, Cabaret] Cabaret and Playne losetime updating all thetemporary labels thanks to the use ofsegments, HA4’s processingtime decreases after thepercolation thresholdd 64%(a) Playne(b) Cabaret(c) HA432 (ccl)(d) HA464 (ccl) HA464 is 2 faster inaverage than Playne andCabaret CCL throughput: 1.2 Gpx/s(HA464 , 2k, g 4)22 / 29

Comparison of CCA algorithms on Jetson TX2 HA464 CCA: labeling kernel is replaced by on-the-fly analysis kernel other algorithms: features computation kernel after relabeling kernel 7 features: S, Sx, Sy, xmin , ymin , xmax , ymax 1.1 Gpx/s (HA464 , 2k, g 4)(a) Playne(b) Cabaret(a) HA432(b) HA46423 / 29

Performance of CCL on Jetson AGX & V100Latest results on Volta architecture: AGX: 4.6 Gpx/s (HA464 , 2k, g 4) V100: 27.0 Gpx/s (HA464 , 2k, g 4)(a) HA432 Jetson AGX(b) HA464 Jetson AGX(c) HA432 V100(d) HA464 V10024 / 29

Performance of CCA on Jetson AGX & V100Latest results on Volta architecture: AGX: 3.4 Gpx/s (HA464 , 2k, (S, Sx, Sy, xmin , ymin , xmax , ymax ), g 4) V100: 14.9 Gpx/s (HA464 , 2k, (S, Sx, Sy, Sx2 , Sy2 ), g 4)(a) HA432 Jetson AGX(b) HA464 Jetson AGX(c) HA432 V100(d) HA464 V10025 / 29