1
$\begingroup$

The answer to Convolute a gaussian kernel with a large array of off-grid centroids without looping? (how to make "A Thousand (Gaussian) Points of Light" ) involves summing a 3D array over its stacked direction. This question however is about convolution using continuous functions rather than summation.

The answer there mentions

a Fredholm integral equation of the first kind

and that may work here as well.


I'd like to know how to convolute $f(x, y)$ with a parametric shape; a 1D distribution along a parametric path as defined by $g(t), \ x(t), \ y(t)$ in Python, resulting in a 2D array of $f * g$.

A simple, illustrative example of this would be:

$$g(t) = \frac{\Gamma^2}{(t-t_0)^2 + \Gamma^2} \ \ \ \text{ (the 1D function)}$$

$$x(t), \ y(t) = a\cos(t), \ a\sin(t) \ \ \ \text{ (over the 2D path)}$$

$$f(x, y) = \exp(-(x^2+y^2)/2 \sigma^2) \ \ \ \text{ (the 2D kernel)}$$

If $a=1, \ t_0 = \pi/3, \ \Gamma = \pi/30 \text{ and } \sigma=0.05$ then I can approximate it brute force as follows; I've made the parametric variable $t$ step size large intentionally to illustrate the sum. Making it sufficiently densely spaced to reduce the error in the result and adding the appropriate normalization (omitted here) is the brute-force solution, I'm wondering if I'm missing a simpler solution that might be better in cases where the parametric function is longer and more complicated.

note: "Better" in this case means faster, as in running mostly in compiled code of some Python method or library; this calculation will be done perhaps 100 times each time a user touches a slider in an interactive simulation.

result of convolution of a Gaussian with a 1D Lorentzian along a parametric path

import numpy as np
import matplotlib.pyplot as plt

t_0 = np.pi/3
gamma = np.pi/20
a = 1.
sigma = 0.05

all_t = np.linspace(0, 2*np.pi, 50) # big steps to make it bumpy
all_g = gamma**2 / ((all_t - t_0)**2 + gamma**2)

y, x = np.mgrid[-1:2:501j, -1:2:501j]

result = np.zeros_like(y)
for (yc, xc), g in zip(zip(*[a * f(all_t) for f in (np.sin, np.cos)]), all_g):
    result += g * np.exp(-((x-xc)**2 + (y-yc)**2)/(2*sigma**2))

extent = [x.min(), x.max(), y.min(), y.max()]
plt.imshow(result, origin='lower', extent=extent)
plt.show()
    
$\endgroup$
6
  • 1
    $\begingroup$ When you say you will be running this 100s of times interactively, what features are being varied with a slider? If the functions will remain mostly the same and it's just the path changing, you may be able to evaluate the convolution symbolically and then just evaluate over the new path. It would have a high startup cost to symbolically determine the convolution, but interactively evaluating it should then be fairly cheap. $\endgroup$
    – Tyberius
    Jul 18 at 19:36
  • $\begingroup$ Yes, for my application 99% of the time it will indeed be the path, that will be my primary concern. And certainly at first the paths will be circular arcs as shown and $f$ and $g$ will always be simple analytical functions as shown if that helps. $g$ could even be a "top-hat" shape $H(t_0-t+\Gamma) - H(t_0-t-\Gamma)$ to start with (where $H$ is the step function). $\endgroup$
    – uhoh
    Jul 18 at 22:54
  • $\begingroup$ @Tyberius I'm interested if there's a general solution, but wow if there is a specific solution with some constraints that might get me started in my current situation that would be very helpful. I'll try to look into this today. $\endgroup$
    – uhoh
    Jul 18 at 22:55
  • $\begingroup$ I think I misinterpreted the conditions of the problem. If $g$ is only known relative to the path, rather than a general function defined over all $(x,y)$, then I don't think there is a way to precompute the convolution and then just evaluate it over the path. I suppose this is trying to solve a different problem, but you may be able to glean some ideas about improving performance from the technique of Line integral convolution. $\endgroup$
    – Tyberius
    Jul 19 at 16:07
2
$\begingroup$

Not stating that this approach is faster, but maybe it inspires you or someone else. The conventional approach in here was certainly faster compared to my implementation of the approach outlined below. My goal was to somehow combine the parameterized path $(\varphi_1(t),\varphi_2(t))$ and $g(t)$ into a single entity.

Taking the parameterized path as $\varphi_1(t)$ and $\varphi_2(t)$, the expression for the convolution is: \begin{align} q(x,y) &= \int_{0}^{2\pi} e^{-b((x-\varphi_1(t))^2 + (y-\varphi_2(t))^2)}g(t) \text{d}t\\ &= e^{-b(x^2+y^2)}\int_{0}^{2\pi}e^{2bx\varphi_1(t)}e^{2by\varphi_2(t)}g(t)e^{-b(\varphi_1^2(t)+\varphi_2^2(t))} \text{d}t \end{align} Then defining $$\hat{g}(t) = g(t)e^{-b(\varphi_1^2(t)+\varphi_2^2(t))}$$ and using Taylor expansions of the exponential functions \begin{align} q(x,y) &= e^{-b(x^2+y^2)}\sum_{h=0}^{\infty}\sum_{k=0}^{\infty} \frac{(2bx)^k}{k!}\frac{(2by)^h}{h!}\int_{0}^{2\pi}\varphi_1^k(t)\varphi_2^h(t)\hat{g}(t)\text{d}t \end{align} The integral can be rearranged into a matrix $G$ with components $G_{kh}$ \begin{align} G_{kh}=\int_{0}^{2\pi}\varphi_1^k(t)\varphi_2^h(t)\hat{g}(t)\text{d}t \end{align} Thus, the path $(\varphi_1(t),\varphi_2(t))$ and $g(t)$ are encoded in $G$.

The double sum can be rewritten as: \begin{align} \hat{x}^TG\hat{y} = \sum_{h=0}^{H}\sum_{k=0}^{K} \frac{(2bx)^k}{k!}\frac{(2by)^h}{h!}\int_{0}^{2\pi}\varphi_1^k(t)\varphi_2^h(t)\hat{g}(t)\text{d}t \end{align} With $\hat{x}$ and $\hat{y}$ being a vectors with components corresponding to different $k$ and $h$, i.e. $$\hat{x}_k = \frac{(2bx)^k}{k!}\,\,\,\,\,\,\,\,\,\,k=0,1,2,3,4...$$ $$\hat{y}_h = \frac{(2by)^h}{h!}\,\,\,\,\,\,\,\,\,\,h=0,1,2,3,4...$$ For proper computation of the components, Stirling's formula may be used.

Then the convolution $q(x,y)$ can be written as $$q(x,y)=\hat{x}^TG\hat{y} e^{-b(x^2+y^2)}$$ I tried it out in Matlab and it works, though it's rather slow. Maybe it can be done faster.

$\endgroup$
1
  • $\begingroup$ Oh this looks like it's going to be fun! I'll dig in and take it for a spin. Thanks. $\endgroup$
    – uhoh
    Aug 15 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.