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.
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: \begin{equation} \hat{f}(\vec{q}) \;=\; FT\left\{f\right\} \;=\; \int d\vec{r}\; f(\vec{r})\; e^{-i \vec{q}\cdot\vec{r}} \end{equation} \begin{equation} 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}} \end{equation} \begin{equation} (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\} \end{equation} 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 \begin{equation} \hat{g}(\vec{q}) \;=\; \pi\, \exp\left(-\frac{1}{4}{|\vec{q}|}^2\right)\, . \end{equation}
We wish to calculate \begin{equation} (f*g)(\vec{r}) \;=\; \pi\,{FT}^{-1}\left\{\hat{f}(\vec{q})\, \exp\left(-\frac{1}{4}{|\vec{q}|}^2\right)\right\}\, . \end{equation} 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$. \begin{equation} \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\} \end{equation} 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: \begin{equation} \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) \end{equation}
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.