# Convolution of convex polygons and a Gaussian

I need to find the closest solutions for convolution of convex polygons/circles with a Gaussian function for computer graphics purposes. I was only able to find solutions for rectangles, like this unit square example:

$$f(x,y) = \theta (1+x) \theta (1-x) \theta (1+y) \theta (1-y) \\ g(x,y) = e^{-x^2-y^2} \\ (f∗g)(x,y) = \frac{1}{4} \pi \left(\text{erf}(1-x)+\text{erf}(1+x)\right) \left(\text{erf}(1-y)+\text{erf}(1+y)\right)$$

The result being exactly what I'm looking for with rectangles.

However I can't find a solution with a sphere defined as $f(x,y) = 1-\theta \left(\sqrt{x^2+y^2}-r\right)$ with $r$ as the radius. When I use $\mathcal{F}^{-1}\left(\mathcal{F}(f)\cdot\mathcal{F}(g)\right)$ in Mathematica it runs for about 2 hours before MathKernel actually crashes. If run separately $\mathcal{F}(f)$ takes about half an hour and gives a pretty long result. I'd like to find an exact answer to the sphere convolution, though I can still approximate it using something like: $$h(x,y) = \frac{1}{2} \left(1-\text{erf}\left(\sqrt{x^2+y^2}-r\right)\right)-\frac{1}{2} \left(1-\text{erf}\left(\sqrt{x^2+y^2}+r\right)\right)$$

So the real problem I need to solve is for triangles and generally for non-rectangular convex polygons. Here's an example of a triangle made from rotated Heaviside step functions mutiplied together:

$$f(x,y) = \theta(x\cos(a_0)-y\sin(a_0)+1) * \theta(x\cos(a_1)-y\sin(a_1)+1) * \\\theta(x\cos(a_2)-y\sin(a_2)+1)$$ with for instance $a_0=1, a_1=3, a_2=5$. Again Mathematica has been of little help here, and I know that approximating the solution using rotated error functions multiplied together gives very wrong results (mostly when the Gaussian is large wrt the whole triangle or with very obtuse or acute angles), so I have no idea what I can do.

Edit: I am now mostly just interested in a solution for triangles, or even just a single corner of a triangle, that would be a good start.

• If this is for computer graphics, why do you actually need the explicit solution? Why not just convolve in physical space numerically? If you just need the numerical value that would be much, much quicker. Nov 18, 2014 at 13:19
• Because it would be very slow (for each pixel you'd have to do hundreds of lookups and operations) and impractical, and in some cases (like features that are too small) it would be very inaccurate. I use this approach to draw perfectly antialiased lines, circles, squares etc, and for each pixel I only need a very few table lookups and a handful of simple integer operations. It's both simple, fast and perfect. Now I'd like to go further than just lines, squares and circles and find the solution (or get near) for other convex polygons as it would allow me to make a pretty complete vector engine. Nov 18, 2014 at 13:54
• Well, good luck with that. For the rectangles you are using heavily that the sides are parallel to the axes. The Fourier transform for a general triangle is a complete mess. You are unlikely to get good closed form formulas that apply universally. What you probably want is to divide into cases where the polygon is small/big compared to the Gaussian function. In the small case carry a table of precomputed values for a large collection of polygons. In the big case build your approximations based on pre-computed values for infinite pie-wedges. Nov 18, 2014 at 14:19
• As a final comment though: if you are doing numerical convolution, you don't do "hundreds of lookups and operators" for "each pixel". You store the graph of the original function $f$ as a giant vector, and convolution is simply achieved by matrix multiplication against a good kernel. Sure the dimensions of your image may be large, but the last I checked, vector addition and multiplication is not that slow. (Your kernel is essentially fixed from the get-go and should be prepopulated...) Nov 18, 2014 at 14:24
• This question and answer may at least give some information on calculating the Fourier transform of a convex polygon. I've just used this result to calculate the closed-form FT of a planar triangle, but I get lost when trying perform the inverse FT of this times a Gaussian. (I can post what I have if desired.) Aug 6, 2017 at 0:53

I'll post a partial answer, since it's all I have and it may help someone to arrive at a full answer.

First off, since different people sometimes use different definitions of the Fourier integral, let me just note that my conventions and notation in this post are such that: $$\hat{f}(\vec{q}) \;=\; FT\left\{f\right\} \;=\; \int d\vec{r}\; f(\vec{r})\; e^{-i \vec{q}\cdot\vec{r}}$$ $$f(\vec{r}) \;=\; {FT}^{-1}(\hat{f}) \;=\; \frac{1}{{(2\pi)}^2}\int d\vec{q}\; \hat{f}(\vec{q})\; e^{+i \vec{q}\cdot\vec{r}}$$ $$(f*g)(\vec{r}) \;=\; \int d\vec{r}'\; f(\vec{r}')\; g(\vec{r}-\vec{r}') \;=\; {FT}^{-1}\left\{\hat{f}(\vec{q})\; \hat{g}(\vec{q})\right\}$$ Here vectors like $\vec{r} = (x,y)$ and $\vec{q} = (q_x,q_y)$ are meant to represent 2D vectors.

Let $f(\vec{r})$ represent a function which is 1 inside of a non-self-intersecting polygon (of any number of sides, and not necessarily convex) and 0 outside. Let $g(\vec{r}) = e^{-{|\vec{r}|}^2}$. It is straightforward to show that $$\hat{g}(\vec{q}) \;=\; \pi\, \exp\left(-\frac{1}{4}{|\vec{q}|}^2\right)\, .$$

We wish to calculate $$(f*g)(\vec{r}) \;=\; \pi\,{FT}^{-1}\left\{\hat{f}(\vec{q})\, \exp\left(-\frac{1}{4}{|\vec{q}|}^2\right)\right\}\, .$$ Here is a derivation of an expression for $\hat{f}(\vec{q})$:

Suppose the locations of the vertices of our $N$-sided polygon are (in either clockwise or counter-clockwise order) $\{\vec{v}_1, ..., \vec{v}_N\}$. Furthermore, let $\hat{n}_{ij}$ represent the outward unit normal to the polygon on the side defined by the sequential vertices $\vec{v}_i$ and $\vec{v}_j$. Finally, let $\Omega$ be the set of points in the plane inside the polygon. Then: \begin{align} \hat{f}(\vec{q}) &= \int_{\Omega}d\vec{r}\; e^{-i \vec{q}\cdot\vec{r}}\\[0.1in] &= -\frac{1}{q^2}\int_{\Omega}d\vec{r}\; \nabla^2 e^{-i \vec{q}\cdot\vec{r}} \qquad\qquad\qquad q = |\vec{q}|\\[0.1in] &= -\frac{1}{q^2}\int_{\Omega}d\vec{r}\; \vec{\nabla}\cdot \left( \vec{\nabla}e^{-i \vec{q}\cdot\vec{r}}\right)\\[0.1in] &= -\frac{1}{q^2}\int_{\partial\Omega}d\ell\; \hat{n}\cdot\vec{\nabla}\; e^{-i \vec{q}\cdot\vec{r}} \qquad\qquad\;\; \text{By the Divergence Theorem}\\[0.1in] &= \frac{i}{q^2}\int_{\partial\Omega}d\ell\; \hat{n}\cdot\vec{q}\; e^{-i \vec{q}\cdot\vec{r}} \end{align} In the last expression above, $\partial\Omega$ represents the perimeter of the polygon, $d\ell$ (a scalar) represents an infinitesimal piece of the length of the perimeter, and $\hat{n}$ is the position-dependent outward normal of the polygon. We now take advantage of the fact that this is a polygon having a discrete boundary made of line segments. The perimeter integral can be divided up into pieces, one for each side of the polygon. On each of these sides, the outward normal vector is constant. Below, the notation $\sum_{cyc\{ij\}}$ will denote a sum over all pairs of adjacent vertices $i$ and $j$. $$\hat{f}(\vec{q}) \;=\; \frac{i}{q^2} \sum_{cyc\{ij\}}\left(\hat{n}_{ij}\cdot\vec{q}\right)\;|\vec{v}_j-\vec{v}_i| \;\int_0^1 ds\; \exp\Bigl\{-i \vec{q}\cdot\bigl[\vec{v}_i + s(\vec{v}_j-\vec{v}_i)\bigr]\Bigr\}$$ Here $|\vec{v}_j-\vec{v}_i|$ is the length of the side between vertices $i$ and $j$, and we have made a change of variables so that $d\ell = |\vec{v}_j-\vec{v}_i|\, ds$. The integral along a given side proceeds along the line from $\vec{v}_i$ to $\vec{v}_j$ as $s$ varies from 0 to 1.

If we now perform the simple integral in $s$, and rearrange things a bit, we arrive at the final expression: $$\hat{f}(\vec{q}) \;=\; -\frac{1}{q^2} \sum_{cyc\{ij\}} |\vec{v}_j-\vec{v}_i|\; \frac{\hat{n}_{ij}\cdot\vec{q}}{(\vec{v}_j-\vec{v}_i)\cdot\vec{q}}\; \left(e^{-i \vec{q}\cdot\vec{v}_j} \;-\; e^{-i \vec{q}\cdot\vec{v}_i}\right)$$

A few things to note about this expression: 1) It works for any non-self-intersecting planar polygon, not just triangles. 2) It doesn't matter whether the vertices are numbered going clockwise or counter-clockwise, as long as they are in order. 3) In the limit as $q\rightarrow 0$, we have (from the first line of the derivation above): $\hat{f}(\vec{0}) = \int_{\Omega}d\vec{r} = |\Omega|$, i.e. the area of the polygon.

Putting all of this together with the original question, we have: \begin{align} (f*g)(\vec{r}) &\;=\; \pi\,{FT}^{-1}\left\{\hat{f}(\vec{q})\, \exp\left(-\frac{1}{4}{|\vec{q}|}^2\right)\right\}\\[0.1in] &\;=\; -\frac{\pi}{{(2\pi)}^2} \sum_{cyc\{ij\}} |\vec{v}_j-\vec{v}_i|\; \int d\vec{q}\; \frac{1}{q^2}\; \frac{\hat{n}_{ij}\cdot\vec{q}}{(\vec{v}_j-\vec{v}_i)\cdot\vec{q}}\; \left(e^{-i \vec{q}\cdot\vec{v}_j} \;-\; e^{-i \vec{q}\cdot\vec{v}_i}\right)\; e^{-q^2/4}\; e^{+i \vec{q}\cdot\vec{r}} \end{align} Unfortunately, this is where I get stuck. I don't know how to do the final integral.