of Infinite Input Extensions Diego Nehab Andr Maximo IMPA GE - - PowerPoint PPT Presentation

Parallel Recursive Filtering

• f Infinite Input Extensions

Diego Nehab André Maximo IMPA GE Global Research

GTC 2017

Linear time-invariant filters

filter

Linear

• Invariant to scale

filter scale filter scale

Linear

• Invariant to addition

add add filter filter filter

Time invariant

• Invariant to shift

filter t t shift filter t t t shift t

Convolution

• The convolution of sequences 𝐜 and 𝐢 is a sequence 𝐛

𝐛 = 𝐜 ∗ 𝐢 𝑏𝑗 = ෍

𝑘=−∞ ∞

𝑐

𝑘ℎ𝑗−𝑘

𝐜, 𝐢, 𝐛: ℤ → ℝ 𝐜 𝐢 𝐢 𝐢 𝐛 = 𝐜 ∗ 𝐢 𝐢

Linear shift-invariant filters are convolutions

filter’s impulse response unit impulse filter convolve filter

∗

Examples

∗ ∗

Outline of talk

• Introduction
• Recursive filters are very useful
• Initialization at the boundaries is an important problem
• Exact recursive filtering of infinite input extensions
• Closed-form formulas available for the first time
• Enable simple and effective algorithms
• Parallelization
• Fastest recursive filtering algorithms to date
• First to filter infinite extensions exactly

General model for convolutions

• Linear difference equations

𝐁𝐴 = 𝐂𝐱 ⋱ ⋱ ⋱ ⋱ ⋱ 𝑏0 𝑏−1 ⋮ 𝑏−𝑠 ⋱ 𝑏1 𝑏0 𝑏−1 ⋮ 𝑏−𝑠 ⋱ ⋮ 𝑏1 𝑏0 𝑏−1 ⋮ ⋱ 𝑏𝑠 ⋮ 𝑏1 𝑏0 𝑏−1 ⋱ 𝑏𝑠 ⋮ 𝑏1 𝑏0 ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 𝑨𝑗−2 𝑨𝑗−1 𝑨𝑗 𝑨𝑗+1 𝑨𝑗+2 ⋮ = ⋱ ⋱ ⋱ ⋱ ⋱ 𝑐0 𝑐−1 ⋮ 𝑐−𝑒 ⋱ 𝑐1 𝑐0 𝑐−1 ⋮ 𝑐−𝑒 ⋱ ⋮ 𝑐1 𝑐0 𝑐−1 ⋮ ⋱ 𝑐𝑒 ⋮ 𝑐1 𝑐0 𝑐−1 ⋱ 𝑐𝑒 ⋮ 𝑐1 𝑐0 ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 𝑥𝑗−2 𝑥𝑗−1 𝑥𝑗 𝑥𝑗+1 𝑥𝑗+2 ⋮ constant diagonals 𝑏𝑠𝑨𝑗−𝑠 + ⋯ + 𝑏1𝑨𝑗−1+𝑏0𝑨𝑗+𝑏−1𝑨𝑗+1+ ⋯ + 𝑏−𝑠𝑨𝑗+𝑠 = 𝑐𝑒𝑥𝑗−𝑒 + ⋯ + 𝑐1𝑥𝑗−1+𝑐0𝑥𝑗+𝑐−1𝑥𝑗+1+ ⋯ + 𝑐−𝑒𝑥𝑗+𝑒 filter 𝐴 𝐱 𝐴 𝐴

finite impulse response support (FIR) ⋮ 𝑦𝑗−2 𝑦𝑗−1 𝑦𝑗 𝑦𝑗+1 𝑦𝑗+2 ⋮ = ⋱ ⋱ ⋱ ⋱ ⋱ 𝑐0 𝑐−1 ⋮ 𝑐−𝑡 ⋱ 𝑐1 𝑐0 𝑐−1 ⋮ 𝑐−𝑡 ⋱ ⋮ 𝑐1 𝑐0 𝑐−1 ⋮ ⋱ 𝑐𝑡 ⋮ 𝑐1 𝑐0 𝑐−1 ⋱ 𝑐𝑡 ⋮ 𝑐1 𝑐0 ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 𝑥𝑗−2 𝑥𝑗−1 𝑥𝑗 𝑥𝑗+1 𝑥𝑗+2 ⋮

• Decompose into direct and recursive parts

𝐲 𝐴 recursive 𝑦𝑗 = 𝑐𝑡𝑥𝑗−𝑡 + ⋯ + 𝑐1𝑥𝑗−1+𝑐0𝑥𝑗+𝑐−1𝑥𝑗+1+ ⋯ + 𝑐−𝑡𝑥𝑗+𝑡 direct part is what we think of as convolution it is like a matrix multiplication 𝐲 = 𝐂𝐱 𝑃(𝑡 𝑜)

General model for convolutions

𝑏𝑠𝑨𝑗−𝑠 + ⋯ + 𝑏1𝑨𝑗−1+𝑏0𝑨𝑗+𝑏−1𝑨𝑗+1+ ⋯ + 𝑏−𝑠𝑨𝑗+𝑠 = 𝑦𝑗 recursive part is the inverse of a convolution it is like a linear system 𝐁𝐴 = 𝐲 infinite impulse response support (IIR) ⋱ ⋱ ⋱ ⋱ ⋱ 𝑏0 𝑏−1 ⋮ 𝑏−𝑠 ⋱ 𝑏1 𝑏0 𝑏−1 ⋮ 𝑏−𝑠 ⋱ ⋮ 𝑏1 𝑏0 𝑏−1 ⋮ ⋱ 𝑏𝑠 ⋮ 𝑏1 𝑏0 𝑏−1 ⋱ 𝑏𝑠 ⋮ 𝑏1 𝑏0 ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 𝑨𝑗−2 𝑨𝑗−1 𝑨𝑗 𝑨𝑗+1 𝑨𝑗+2 ⋮ = ⋮ 𝑦𝑗−2 𝑦𝑗−1 𝑦𝑗 𝑦𝑗+1 𝑦𝑗+2 ⋮ direct 𝐱 filter 𝐴 𝐱 𝐴 𝐴 direct 𝐱 𝐴 recursive 𝐁𝐴 = 𝐂𝐱 𝐲 = 𝐂𝐱 𝐁𝐴 = 𝐲

𝐅𝐴 = √𝑏0 ⋱ ⋱ ⋱ ⋱ 1 𝑓1 ⋮ 𝑓𝑠 1 𝑓1 ⋮ 𝑓𝑠 1 𝑓1 ⋮ ⋱ 1 𝑓1 ⋱ 1 ⋱ ⋱ 𝐴 = 𝐳 anticausal 𝐴

General model for convolutions

• Decompose recursive part into causal and anticausal passes

direct 𝐱 𝐲 causal 𝐳 anticausal 𝐴 𝑧𝑗 = 1 √𝑏0 𝑦𝑗 − 𝑒1𝑧𝑗−1 − ⋯ − 𝑒𝑠𝑧𝑗−𝑠 causal part is forward-substitution 𝐄𝐳 = 𝐲 𝑃(𝑠 𝑜) 𝑨𝑗 = 1 √𝑏0 𝑧𝑗 − 𝑓1𝑨𝑗+1 − ⋯ − 𝑓𝑠𝑨𝑗+𝑠 anticausal is back-substitution 𝐅𝐴 = 𝐳 𝑃(𝑠 𝑜) 𝐲 causal 𝐄𝐳 = √𝑏0 ⋱ ⋱ 1 ⋱ 𝑒1 1 ⋱ ⋮ 𝑒1 1 𝑒𝑠 ⋮ 𝑒1 1 𝑒𝑠 ⋮ 𝑒1 1 ⋱ ⋱ ⋱ ⋱ 𝐳 = 𝐲 𝐴 recursive 𝐲 𝐁𝐴 = 𝐲 𝐄𝐳 = 𝐲 𝐅𝐴 = 𝐳 𝐲 = 𝐂𝐱

recursive filter

Downsampling with cardinal cubic B-splines

prefiltered post-processed output input Catmull-Rom t cardinal cubic B-spline t

infinite support

cubic B-spline t prefiltered [Nehab & Hoppe 2014]

input direct convolution 𝑡 = 2 𝜏

2 𝜏 𝑜 operations

Fast image blur

Blur with FIR filter (given 𝜏)

causal pass 𝑠 = 3 anticausal pass 𝑠 = 3

6𝑜 operations Blur with IIR filter (any 𝜏) Blur with FFT 𝑜 log 𝑜 operations

[van Vliet et al. 1998]

? ?

What do near input boundaries?

Infinite input extensions

repeat periodically filter reflect periodically filter constant padding filter

periodic repetition clamp to border

Tileable textures

• Textures designed to be tiled in a certain way
• Filtering must respect the periodicity
• riginal texture

filtered, tiled, and shifted

Dealing with boundaries in practice

• In the frequency domain
• DFT/DCT imply infinite extensions
• Computations are exact even for IIR filters
• In the time domain
• Direct convolution can decide out of bounds input arbitrarily
• Recursive filters must define out of bounds outputs!

even more wasted computation filtered even more wasted memory padded filtered more wasted computation padded more wasted memory wasted memory padded wasted computation filtered

Approximation by input padding

input approximation

• utput

amount of padding depends

• n impulse response ”support”

Outline of talk

• Introduction
• Recursive filters are very useful
• Initialization at the boundaries is an important problem
• Exact recursive filtering of infinite input extensions
• Closed-form formulas available for the first time
• Enable simple and effective algorithms
• Parallelization
• Fastest recursive filtering algorithms to date
• First to filter infinite extensions exactly

finite output

Exact recursive filtering of finite input

finite input 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 input extension 𝐲–3 𝐲–1 𝐲–0 𝐲–2 𝐲–4 𝐲–5

…

𝐲8 𝐲10 𝐲7 𝐲9 𝐲11 𝐲12 … input extension 𝐴8 𝐴10 𝐴9 𝐴11 𝐴12 … 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 𝐴7 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 𝐳8 𝐳10 𝐳7 𝐳9 𝐳11 𝐳12 … 𝐳0

…

𝐳–2 𝐳–1 𝐳–3 𝐳–4 𝐳–5

finite output

Exact recursive filtering of finite input

• Instead, first obtain the initial feedbacks
• To initialize causal pass
• To initialize subsequent anticausal pass
• How to obtain these feedbacks in closed form?
• Depends on the choice of infinite input extension

initial causal feedback initial anticausal feedback finite input input extension 𝐲–3 𝐲–1 𝐲–0 𝐲–2 𝐲–4 𝐲–5

…

𝐲8 𝐲10 𝐲7 𝐲9 𝐲11 𝐲12 … input extension 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 𝐴7 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 𝐳0

Constant padding

𝐳8 𝐳10 𝐳7 𝐳9 𝐳11 𝐳12 … series finite input 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 input extension 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 anticausal 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 causal 𝐲0 𝐲0 𝐲0 𝐲0 𝐲0 𝐲0

…

𝐳–2 𝐳–1 𝐳–3 𝐳–4 𝐳–5

…

infinite series 𝐳0 𝐳0 = 𝑵1 𝐲0 initial causal feedback 𝐴8 𝐴10 𝐴9 𝐴11 𝐴12 … infinite series 𝐴7 = 𝑵2 𝐴7 + 𝑵3 𝐳6 𝐲7 initial anticausal feedback 𝐲7 𝐲7 𝐲7 𝐲7 𝐲7 𝐲7 … input extension 𝑵1 = 𝑻𝐺ഥ 𝑩𝐺 𝑻𝐺 = 𝑱 − 𝑩𝐺

𝑠 −1

𝑵2 = 𝑻𝑆𝐺𝑩𝐺

𝑠

𝑻𝑆𝐺 − 𝑩𝑆

𝑠 𝑻𝑆𝐺𝑩𝐺 𝑠 = ഥ

𝑩𝑆 𝑵3 = 𝑻𝑆ഥ 𝑩𝑆 − 𝑻𝑆𝐺𝑩𝐺

𝑠 𝑵1

𝑻𝑆 = 𝑱 − 𝑩𝑆

𝑠 −1

finite output precomputed 𝑠 × 𝑠 matrix precomputed 𝑠 × 𝑠 matrices

ሶ 𝐳2 ሶ 𝐳4 ሶ 𝐳5 ሶ 𝐳6 ሶ 𝐳1 ሶ 𝐳3 1st causal 2nd causal 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 2nd anticausal 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 ሶ 𝐴5 ሶ 𝐴3 ሶ 𝐴2 ሶ 𝐴1 ሶ 𝐴6 ሶ 𝐴4 1st anticausal

Periodic repetition

finite input 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 = 𝑵5 𝐴1 ሶ 𝐴1 𝐴1 initial anticausal feedback 𝐳6 = 𝑵4 ሶ 𝐳6 𝐳6 initial causal feedback 𝑵4 = 𝑱 − 𝑩𝐺

𝑜 −1

𝑵5 = 𝑱 − 𝑩𝑆

𝑜 −1

finite output input extension 𝐲3 𝐲5 𝐲6 𝐲4 𝐲2 𝐲1

…

𝐲2 𝐲4 𝐲1 𝐲3 𝐲5 𝐲6 … input extension

ሷ 𝐴5 ሷ 𝐴3 ሷ 𝐴2 ሷ 𝐴1 ሷ 𝐴6 ሷ 𝐴4 1st anticausal ሶ 𝐳2 ሶ 𝐳4 ሶ 𝐳5 ሶ 𝐳6 ሶ 𝐳1 ሶ 𝐳3 1st causal 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 2nd causal

Periodic reflection

finite input 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 ሶ 𝐳6 𝐳12 = 𝑵8 + 𝑵9 ሷ 𝐴1 ሶ 𝐳6 𝐳12 = 𝑵6 + 𝑵7 𝑵7 = 𝑱 − 𝑩𝐺

2ℎ −1𝑳 ഥ

𝑩𝑆

−1 𝑱 − 𝑩𝐺 𝑠𝑩𝑆 𝑠

𝑵6 = 𝑱 − 𝑩𝐺

2ℎ −1 𝑩𝐺 ℎ − 𝑳𝑩𝑆 ℎ ഥ

𝑩𝐺

−1𝑩𝐺 𝑠 ഥ

𝑩𝑆 𝑵8 = 𝑳 − 𝑩𝑆

𝑠 −𝟐ഥ

𝑩𝑆 𝑵9 = 𝑳 − 𝑩𝑆

𝑠 −𝟐ഥ

𝑩𝑆𝑩𝐺

ℎ

𝐳12 initial causal feedback initial anticausal feedback finite output 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 2nd anticausal input extension

… …

input extension

Is this correct? Is this fast?

• All required 𝑠 × 𝑠 matrices 𝑵1to 𝑵9 exist and can be precomputed
• Proofs in the paper assume only filter stability
• Exact filtering for all infinite extensions takes 𝑃(𝑠 𝑜)
• May require twice as much computation
• Real advantage comes with parallelization
• No additional cost

Outline of talk

• Introduction
• Recursive filters are very useful
• Initialization at the boundaries is an important problem
• Exact recursive filtering of infinite input extensions
• Closed-form formulas available for the first time
• Enable simple and effective algorithms
• Parallelization
• Fastest recursive filtering algorithms to date
• First to filter infinite extensions exactly

GPU

• • •

Threads Threads Shared Memor y Shared Memor y

Global Memory

Modern GPU

• Many multiprocessors each supporting many hardware threads
• 10k threads is typical
• On-chip shared memory within each multiprocessor
• 48k to be shared by threads local to multiprocessor
• Global memory with high throughput but high latency
• High throughput but high latency
• Challenge is hide latency and keeping all cores busy

Multiprocessor Multiprocessor

𝐲13 𝐲15 𝐲16 𝐲12 𝐲14 𝐲1 𝐲3 𝐲4 𝐲5 𝐲2

Independent rows then columns

input 𝐲7 𝐲9 𝐲10 𝐲11 𝐲6 𝐲8

• utput

𝐳1 𝐳3 𝐳4 𝐳5 𝐳2 𝐳13 𝐳15 𝐳16 𝐳12 𝐳14 𝐳7 𝐳9 𝐳10 𝐳11 𝐳6 𝐳8 causal 𝐴1 𝐴2 𝐴3 𝐴4 𝐴13 𝐴15 𝐴14 𝐴16 𝐴8 𝐴7 𝐴6 𝐴5 𝐴10 𝐴11 𝐴9 𝐴12 anticausal [Ruijters and Thevenaz 2010]

Summary of parallel algorithms

𝑠 is the filter order 𝑞 is the number of processors ℎ is the image height 𝑥 is the image width

Algorithm Step complexity Threads Bandwidth independent rows then columns 4𝑠 ℎ𝑥/𝑞 ℎ, 𝑥 8ℎ𝑥

7 1 3 5 Throughput (GiP/s) 1st order Input size (pixels) 642 1282 2562 5122 10242 20482 40962

[Nehab, Maximo, Lima & Hoppe 2011]

independent rows then columns

Performance on Fermi (GF100)

correct feedbacks 𝐴9 𝐴1 𝐴5 𝐴13 ሶ 𝐴9 ሶ 𝐴1 ሶ 𝐴5 ሶ 𝐴13 𝐴10 𝐴11 𝐴9 𝐴12 2nd anticausal 𝐳13 𝐳15 𝐳16 𝐳14 2nd causal ሶ 𝐳13 ሶ 𝐳15 ሶ 𝐳16 ሶ 𝐳14 1st causal ሶ 𝐳4 ሶ 𝐳8 ሶ 𝐳12 ሶ 𝐳16 ሶ 𝐳10 ሶ 𝐳12 ሶ 𝐳9 ሶ 𝐳11 1st causal 𝐳12 𝐳9 𝐳10 𝐳11 2nd causal ሶ 𝐳5 ሶ 𝐳7 ሶ 𝐳6 ሶ 𝐳8 1st causal 𝐳4 𝐳16 𝐳12 𝐳8 correct feedbacks 𝐳5 𝐳7 𝐳6 𝐳8 2nd causal 1st causal ሶ 𝐳1 ሶ 𝐳3 ሶ 𝐳4 ሶ 𝐳2 𝐳3 𝐳1 𝐳4 𝐳2 2nd causal 𝐲13 𝐲15 𝐲16 𝐲14 𝐲1 𝐲3 𝐲4 𝐲2 𝐲12 𝐲9 𝐲10 𝐲11 𝐲5 𝐲7 𝐲6 𝐲8 ሶ 𝐴1 ሶ 𝐴2 ሶ 𝐴3 ሶ 𝐴4 1st anticausal ሶ 𝐴13 ሶ 𝐴15 ሶ 𝐴14 ሶ 𝐴16 1st anticausal ሶ 𝐴8 ሶ 𝐴7 ሶ 𝐴6 ሶ 𝐴5 1st anticausal ሶ 𝐴10 ሶ 𝐴11 ሶ 𝐴9 ሶ 𝐴12 1st anticausal 𝐴1 𝐴2 𝐴3 𝐴4 2nd anticausal 𝐴13 𝐴15 𝐴14 𝐴16 2nd anticausal 𝐴8 𝐴7 𝐴6 𝐴5 2nd anticausal ሶ 𝐴𝑐𝑗 𝐴𝑐𝑗 𝐴𝑐(𝑗+1) = 𝑩𝑆

𝑐

+

Split rows and columns into blocks

• utput

input 𝐲13 𝐲15 𝐲16 𝐲14 𝐲1 𝐲3 𝐲4 𝐲2 𝐲12 𝐲9 𝐲10 𝐲11 𝐲5 𝐲7 𝐲6 𝐲8 ሶ 𝐳𝑐𝑗 𝐳𝑐𝑗 𝐳𝑐(𝑗−1) = 𝑩𝐺

𝑐

+ [Sung & Mitra 1986; Nehab, Maximo, Lima & Hoppe 2011]

Summary of parallel algorithms

𝑠 is the filter order 𝑞 is the number of processors ℎ is the image height 𝑥 is the image width

Algorithm Step complexity Threads Bandwidth independent rows then columns 4𝑠 ℎ𝑥/𝑞 ℎ, 𝑥 8ℎ𝑥 + split rows and columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 9ℎ𝑥

7 1 3 5 Throughput (GiP/s) 1st order Input size (pixels) 642 1282 2562 5122 10242 20482 40962

[Nehab, Maximo, Lima & Hoppe 2011]

independent rows then columns + split rows and columns

Performance on Fermi (GF100)

correct feedbacks 𝐴9 𝐴1 𝐴5 𝐴13 ሷ 𝐴10 ሷ 𝐴11 ሷ 𝐴9 ሷ 𝐴12 1st anticausal 𝐴10 𝐴11 𝐴9 𝐴12 2nd anticausal ሷ 𝐴9 ሷ 𝐴1 ሷ 𝐴5 ሷ 𝐴13 𝐴13 𝐴15 𝐴14 𝐴16 2nd anticausal ሷ 𝐴13 ሷ 𝐴15 ሷ 𝐴14 ሷ 𝐴16 1st anticausal 𝐳15 𝐳13 𝐳16 𝐳14 2nd causal 𝐳12 𝐳9 𝐳10 𝐳11 2nd causal 𝐳5 𝐳7 𝐳6 𝐳8 2nd causal 𝐳3 𝐳1 𝐳4 𝐳2 2nd causal ሶ 𝐳12 ሶ 𝐳9 ሶ 𝐳10 ሶ 𝐳11 1st causal 1st causal ሶ 𝐳1 ሶ 𝐳3 ሶ 𝐳4 ሶ 𝐳2 ሶ 𝐳13 ሶ 𝐳15 ሶ 𝐳16 ሶ 𝐳14 1st causal ሶ 𝐳5 ሶ 𝐳7 ሶ 𝐳6 ሶ 𝐳8 1st causal 𝐳4 𝐳16 𝐳12 𝐳8 correct feedbacks ሷ 𝐴8 ሷ 𝐴7 ሷ 𝐴6 ሷ 𝐴5 1st anticausal ሷ 𝐴1 ሷ 𝐴2 ሷ 𝐴3 ሷ 𝐴4 1st anticausal 𝐴8 𝐴7 𝐴6 𝐴5 2nd anticausal 𝐴1 𝐴2 𝐴3 𝐴4 2nd anticausal ሶ 𝐳4 ሶ 𝐳8 ሶ 𝐳12 ሶ 𝐳16

Overlap causal-anticausal processing

• utput

input ሶ 𝐳𝑐𝑗 𝐳𝑐𝑗 𝐳𝑐(𝑗−1) = 𝑩𝐺

𝑐

+ [Nehab, Maximo, Lima & Hoppe 2011] ሷ 𝐴𝑐𝑗 𝐴𝑐𝑗 𝐴𝑐(𝑗+1) = 𝑩𝑆

𝑐

+ 𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝐳𝑐(𝑗−1) + 𝐲13 𝐲15 𝐲16 𝐲14 𝐲1 𝐲3 𝐲4 𝐲2 𝐲12 𝐲9 𝐲10 𝐲11 𝐲5 𝐲7 𝐲6 𝐲8

Summary of parallel algorithms

𝑠 is the filter order 𝑞 is the number of processors ℎ is the image height 𝑥 is the image width

Algorithm Step complexity Threads Bandwidth independent rows then columns 4𝑠 ℎ𝑥/𝑞 ℎ, 𝑥 8ℎ𝑥 + split rows and columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 9ℎ𝑥 + overlap causal with anticausal ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 5ℎ𝑥

7 1 3 5 Throughput (GiP/s) 1st order Input size (pixels) 642 1282 2562 5122 10242 20482 40962

[Nehab, Maximo, Lima & Hoppe 2011]

independent rows then columns + split rows and columns + overlap causal and anticausal

Performance on Fermi (GF100)

Overlap causal, anticausal, rows, columns

[Nehab, Maximo, Lima & Hoppe 2011]

Stage 1

ഺ 𝐯𝑛,𝑜 𝐯𝑛,𝑜 𝐯𝑛,𝑜−1

= 𝑩𝐺

𝑐 𝑢 +

𝑩𝑆𝐹 + 𝑩𝑆𝐶𝑩𝐺𝐶 𝑈 𝑩𝐺𝐶 𝑢 +

𝐴𝑛+1,𝑜 𝐳𝑛−1,𝑜

𝑩𝑆

𝑐 𝑢 +

𝐴𝑛+1,𝑜 𝐰𝑛,𝑜 𝐰𝑛,𝑜+1

= + 𝑩𝑆𝐶𝑩𝐺𝑄 𝑩𝑆𝐹

𝐯𝑛,𝑜−1 𝐳𝑛−1,𝑜

𝐼 𝑩𝑆𝐶 𝑩𝐺𝐶 𝑢 +

ሷሷ 𝐰𝑛,𝑜

𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝑢 +

Stage 2

ሶ 𝐳𝑛,𝑜 𝐳𝑛,𝑜 𝐳𝑛−1,𝑜

= 𝑩𝐺

𝑐

+

ሷ 𝐴𝑛,𝑜 𝐴𝑛,𝑜 𝐴𝑛+1,𝑜

= 𝑩𝑆

𝑐

+ 𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝐳𝑛−1,𝑜 +

Stage 3

Stage 3

• utput

Summary of parallel algorithms

𝑠 is the filter order 𝑞 is the number of processors ℎ is the image height 𝑥 is the image width

Algorithm Step complexity Threads Bandwidth independent rows then columns 4𝑠 ℎ𝑥/𝑞 ℎ, 𝑥 8ℎ𝑥 + split rows and columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 9ℎ𝑥 + overlap causal with anticausal ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 5ℎ𝑥 + overlap rows with columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 3ℎ𝑥

7 1 3 5 Throughput (GiP/s) 1st order Input size (pixels) 642 1282 2562 5122 10242 20482 40962

[Nehab, Maximo, Lima & Hoppe 2011]

+ overlap rows and columns independent rows then columns + split rows and columns + overlap causal and anticausal

Performance on Fermi (GF100)

+ overlap rows and columns

[Nehab, Maximo, Lima & Hoppe 2011]

Input size (pixels) 642 1282 2562 5122 10242 20482 40962 2nd order 5 1 2 3 Throughput (GiP/s) 4

• verlap causal and anticausal

Performance on Fermi (GF100)

ഺ 𝐯𝑛,𝑜 𝐯𝑛,𝑜 𝐯𝑛,𝑜−1

= 𝑩𝐺

𝑐 𝑢 +

𝑩𝑆𝐹 + 𝑩𝑆𝐶𝑩𝐺𝐶 𝑈 𝑩𝐺𝐶 𝑢 +

𝐴𝑛+1,𝑜 𝐳𝑛−1,𝑜 𝐳𝑛−1,𝑜

𝑩𝑆

𝑐 𝑢 +

𝐴𝑛+1,𝑜 𝐰𝑛,𝑜 𝐰𝑛,𝑜+1

= + 𝑩𝑆𝐶𝑩𝐺𝑄 𝑩𝑆𝐹

𝐯𝑛,𝑜−1

𝐼 𝑩𝑆𝐶 𝑩𝐺𝐶 𝑢 +

ሷሷ 𝐰𝑛,𝑜

𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝑢 + 𝑩𝑆

𝑐 𝑢 +

𝐰𝑛,𝑜 𝐰𝑛,𝑜+1

=

𝐯𝑛,𝑜−1

𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝑢 + ത

𝐰𝑛,𝑜 𝐯𝑛,𝑜 𝐯𝑛,𝑜−1

= 𝑩𝐺

𝑐 𝑢 + ഥ

𝐯𝑛,𝑜

New trick for better performance

ሶ 𝐳𝑛,𝑜 𝐳𝑛,𝑜 𝐳𝑛−1,𝑜

= 𝑩𝐺

𝑐

+

ሷ 𝐴𝑛,𝑜 𝐴𝑛,𝑜 𝐴𝑛+1,𝑜

= 𝑩𝑆

𝑐

+ 𝐼 𝑩𝑆𝐶 𝑩𝐺𝑄 𝐳𝑛−1,𝑜 + sequentially find per-block column feedbacks sequentially find per-block row feedbacks

ഥ 𝐯𝑛,𝑜 ഺ 𝐯𝑛,𝑜

= 𝑩𝑆𝐹 + 𝑩𝑆𝐶𝑩𝐺𝐶 𝑈 𝑩𝐺𝐶 𝑢 +

𝐴𝑛+1,𝑜 𝐳𝑛−1,𝑜 ത 𝐰𝑛,𝑜 𝐴𝑛+1,𝑜

= + 𝑩𝑆𝐶𝑩𝐺𝑄 𝑩𝑆𝐹

𝐳𝑛−1,𝑜

𝐼 𝑩𝑆𝐶 𝑩𝐺𝐶 𝑢 +

ሷሷ 𝐰𝑛,𝑜

new fully parallel intermediate stage (exactly the same as column processing)

Summary of parallel algorithms

Algorithm Step complexity Threads Bandwidth independent rows then columns 4𝑠 ℎ𝑥/𝑞 ℎ, 𝑥 8ℎ𝑥 + split rows and columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 9ℎ𝑥 + overlap causal with anticausal ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 5ℎ𝑥 + overlap rows with columns ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 3ℎ𝑥 + new trick ≈ 8𝑠 ℎ𝑥/𝑞 ℎ𝑥/𝑐 ≈ 3ℎ𝑥

𝑠 is the filter order 𝑞 is the number of processors ℎ is the image height 𝑥 is the image width

13 1 3 5 7 9 11 642 81922 1282 2562 5122 10242 20482 40962 Input size (pixels) Throughput (GiP/s)

• verlap causal and anticausal

+ overlap rows and columns + new trick

Performance on Kepler (GK110)

3rd order 4th order

642 81922 1282 2562 5122 10242 20482 40962 17 1 3 15 13 11 9 7 5 Input size (pixels) Throughput (GiP/s)

Performance on Pascal (GP102)

• verlap causal and anticausal

+ overlap rows and columns + new trick 5th order

ሶ 𝐳2 ሶ 𝐳4 ሶ 𝐳5 ሶ 𝐳6 ሶ 𝐳1 ሶ 𝐳3 1st causal ሷ 𝐴5 ሷ 𝐴3 ሷ 𝐴2 ሷ 𝐴1 ሷ 𝐴6 ሷ 𝐴4 1st anticausal 𝐴5 𝐴3 𝐴2 𝐴1 𝐴6 𝐴4 2nd anticausal 𝐳2 𝐳4 𝐳5 𝐳6 𝐳1 𝐳3 2nd causal

Back to infinite extensions

entire finite input 𝐲2 𝐲4 𝐲5 𝐲6 𝐲1 𝐲3 ሶ 𝐳6 𝐳12 = 𝑵8 + 𝑵9 ሷ 𝐴1 ሶ 𝐳6 𝐳12 = 𝑵6 + 𝑵7 𝑵7 = 𝑱 − 𝑩𝐺

2ℎ −1𝑳 ഥ

𝑩𝑆

−1 𝑱 − 𝑩𝐺 𝑠𝑩𝑆 𝑠

𝑵6 = 𝑱 − 𝑩𝐺

2ℎ −1 𝑩𝐺 ℎ − 𝑳𝑩𝑆 ℎ ഥ

𝑩𝐺

−1𝑩𝐺 𝑠 ഥ

𝑩𝑆 𝑵8 = 𝑳 − 𝑩𝑆

𝑠 −𝟐ഥ

𝑩𝑆 𝑵9 = 𝑳 − 𝑩𝑆

𝑠 −𝟐ഥ

𝑩𝑆𝑩𝐺

ℎ

𝐳12 initial causal feedback initial anticausal feedback entire finite output input extension

… …

input extension ሶ 𝐳6 ሷ 𝐴1

Back to infinite extensions

• Side effect of 2nd stage of parallel algorithm
• Just what we need to compute exact initial feedbacks!

13 1 3 5 7 9 11 642 81922 1282 2562 5122 10242 20482 40962 Input size (pixels) Throughput (GiP/s) 1st order (cubic B-spline interpolation) [Chaurasia et al. 2015] [Nehab et al. 2011] padding exact reflect exact repeat exact clamp to border

Performance on Kepler (GK110)

conv cuFFT [Chaurasia et al. 2015] [Nehab et al. 2011] 642 81922 1282 2562 5122 10242 20482 40962 7 1 3 5 Input size (pixels) Throughput (GiP/s) 3rd order (2D Gaussian blur with 𝜏 = 𝑜/6)

Performance on Kepler (GK110)

padding exact reflect exact repeat exact clamp to border

Summary

• Introduction
• Recursive filters are very useful
• Initialization at the boundaries is an important problem
• Exact recursive filtering of infinite input extensions
• Closed-form formulas available for the first time
• Enable simple and effective algorithms
• Parallelization
• Fastest recursive filtering algorithms to date
• First to filter infinite extensions exactly

Thank you!

• Please download and use our code

https://github.com/andmax/gpufilter

• Questions?