# Fast Gaussian Filter Approximations Comparison on SIMD Computing Platforms
## 4 Approximations
### FIR Approximations
#### 4.1.4 Piecewise Parabolic Approximation  
This approach is based on the piecewise polynomial kernel approximation. The convolution can be decomposed into a sum of convolutions, each of which is different from zero only in a part of the window. Such a convolution is equivalent to a convolution with  
Figure 3: Approximation of the Gaussian with \(\sigma=10\) by running sums method for \(k=3,4,5\).  
a reduced shifted window whose center no longer coincides with the target pixel. This method will be abbreviated below as PcParab.  
Each polynomial can be represented as the sum of terms; i.e., the problem is reduced to calculating convolutions with a constant, \(x\), \(x^{2}\), etc. The important thing is that the product of the number of sites in the piecewise approximation and the number of terms is not large. For approximating the Gaussian kernel, the idea is to divide the segment \([-3\sigma;\,3\sigma]\) into three non-intersecting sub-segments \(p_{k}\), \(k=1,2,3\) such that \(\cup p_{k}=[-3\sigma;\,3\sigma]\) and approximate the parts of Gaussian kernel on each sub-segment by a parabola. The parabola can be described by a quadratic equation of the form \(y=ax^{2}+bx+c\) with real coefficients \(a,b\), and \(c\). By choosing three different points on each sub-segment \(p_{k}\), we can solve a system of three equations  
\[a_{k}\,x_{i}^{2}+b_{k}\,x_{i}+c_{k}=g_{\sigma}[x_{i}],\ \ i=1,2,3,\ \ x_{i}\in p _{k}, \tag{22}\]  
to determine the coefficients \(a_{k},b_{k},c_{k}\) of the corresponding parabola. An example of the obtained approximation is provided in Figure 4.  
In our experiments, we divided the segment at the Gaussian inflection points \(\pm\,\sigma\) and took the ends and the middle of the sub-segment to construct the parabola. Alternatively, partition points can be found using the least squares method.  
To compute convolutions with terms fast, we use the following technique. Convolution (4) with constant function \(h[n]\equiv 1\) is defined as  
\[y^{C}[i]=\sum_{j=i-M}^{i+M}x[j],\quad i=M,\ldots,N-M-1. \tag{23}\]  
Hence, we have  
\[y^{C}[i+1]=y^{C}[i]-x[i-M]+x[i+M+1]. \tag{24}\]  
First, the value for the left-most pixel \(i=M\) is computed according to (23), and then, given (24), the values for each \(i=M+1,\ldots,N-M-1\) are obtained with \(1\) additive and \(1\) multiplicative operation.  
For the convolution with \(h[n]=n\), we can rewrite (4) as follows:  
\[y^{F}[i+1]=y^{F}[i]-M\cdot\big{(}x[i-M]+x[i+M+1]\big{)}+\underbrace{\sum_{j=i -M+1}^{i+M}x[j]}_{y^{C}[i]-x[i-M]}, \tag{25}\]  
Figure 4: Approximation of Gaussian with \(\sigma=10\) by three parabolas.  
so the result is  
\[y^{F}[i+1]=y^{F}[i]+y^{C}[i]-\left(M+1\right)\cdot x[i-M]-M\cdot x[i+M+1] \tag{26}\]  
and can be easily obtained from the already computed values for \(y^{C}\). All convolutions with \(n\) are computed in this way. Furthermore, to convolve with \(h[n]=n^{2}\), we use the following representation:  
\[y^{F2}[i+1]=y^{F2}[i]-M^{2}\cdot\left(x[i-M]+x[i+M+1]\right)+2\underbrace{ \sum_{j=i-M+1}^{i+M}(i-j)\,x[j]}_{y^{F}[i]-M\cdot x[i-M]}+\underbrace{\sum_{j =i-M+1}^{i+M}x[j]}_{y^{C}[i]-x[i-M]}, \tag{27}\]  
so the final expression comes as follows:  
\[y^{F2}[i+1]=y^{F2}[i]+2\cdot y^{F}[i]+y^{C}[i]-(M+1)^{2}\cdot x[i-M]-M^{2} \cdot x[i+M+1]. \tag{28}\]  
The full result of convolution with a parabola \(h[n]=a_{k}\,n^{2}+b_{k}\,n+c_{k}\) therefore becomes  
\[y[i] =\sum_{j=i-M}^{i+M}\left(a_{k}(i-j)^{2}+b_{k}(i-j)+c_{k}\right)x[ j]= \tag{29}\] \[=a_{k}\sum_{j=i-M}^{i+M}(i-j)^{2}\,x[j]\ +\ b_{k}\sum_{j=i-M}^{i+M}(i-j)\,x[ j]\ +\ c_{k}\sum_{j=i-M}^{i+M}x[j]=\] \[=a_{k}\,y^{F2}[i]+b_{k}\,y^{F1}[i]+c_{k}\,y^{C}[i],\]  
which can be computed in 11 additions and 8 multiplications according to (24), (26), and (28). On each of three sub-segments \(p_{k}\), \(k=1,2,3\), we approximate the Gaussian kernel by a fragment of a parabola, so each sub-segment corresponds to a limited range \(j\) in (29). Initialization requires \(3K\) additions and \(2K\) multiplications in total, where \(K=2M+1\) is the window size.