Quick and simple discrete 2D Helmholtz-Hodge Decomposition using FFTs?

Question

For a silly screen saver I'm trying to develop, I'd like to randomly generate a divergence-free 2D array of 2D vectors, and then use it to generate a line integral convolution plot. I've heard$^1$ that one way to do this is to generate random noise, and then project out the solenoidal component of its Helmholtz-Hodge decomposition. To do that, I tried using the following reasoning:

A function $\mathbf{f}:\mathbb{R}^2\rightarrow\mathbb{R}^2$ has Helmholtz-Hodge decomposition$^2$ $$\mathbf{f}=\mathbf{h}+\nabla\phi+J\nabla\psi$$ where $$J=\pmatrix{0&-1\\1&0}$$ and where $\phi,\psi$ are scalar functions. For now, assume the harmonic component $\mathbf{h}$ vanishes.

In Fourier space, this becomes $$\mathcal{F}\mathbf{f}=-i\mathbf{k}\hat{\phi}-iJ\mathbf{k}\hat{\psi}$$ and we can define a solenoidal projection operator on Fourier space as $$P=I-\frac{\mathbf{k}\otimes\mathbf{k}}{\mathbf{k}\cdot\mathbf{k}}$$ which projects a function on its solenoidal component, via $$\mathcal{F}^{-1}P\mathcal{F}\mathbf{f}=J\nabla\psi.$$

I then tried to implement this in Mathematica, applying it to a $21\times 21\times 2$ random array. First I generate the random array, and take apply the FFT to each of its two components:

arr = RandomReal[{-1, 1}, {2, 21, 21}];
fArr = Fourier /@ arr;

I then define $\mathbf{k}$ as a function of array index:

k[k1_, k2_] := Mod[{k1 - 1, k2 - 1}, 21, -10]/21;

Then I perform the projection on the Fourier components (the singularity at $\mathbf{k}=0$ is left alone using an If statement):

dat = Transpose[
   Table[If[k1 == 1 && k2 == 1, fArr[[;; , k1, k2]], 
     fArr[[;; , k1, k2]] - 
      k[k1, k2] (k[k1, k2].fArr[[;; , k1, k2]])/(k[k1, k2].k[k1, 
           k2])], {k1, 21}, {k2, 21}], {2, 3, 1}];

Then I iFFT the two components:

projArr = InverseFourier /@ dat;

This gives a purely real array, and I would naively expect the result to be an approximation of $J\nabla\psi$. My question is:

• In what sense does the result approximate $J\nabla\psi$?

Supposedly the Helmholtz-Hodge decomposition of 2D data is a nontrivial task, as Chris Beaumont's HH_DECOMP routine is supposed to use FFTs to perform Helmholtz-Hodge decomposition, but he also says (in the comments at the top of the code) that the method seems inaccurate. Likewise, there are more complicated variational methods to perform Helmholtz-Hodge decompositions$^3$ of 2D data, which would seem to suggest that the simpler FFT method is somehow inadequate. Why? What does the FFT method get wrong? And is it wrong to assume the harmonic component vanishes for my random noise?

(1): Stable Fluids, Jos Stam.

(2): Feature Detection in Vector Fields Using the Helmholtz-Hodge Decomposition, Alexander Wiebel, page 12.

(3): Discrete Multiscale Vector Field Decomposition, Yiying Tong.

Answer 1

The easiest way to generate a divergence-free 2D function is to use a streamfunction, $\psi$. Where

$$ u=\frac{\partial\psi}{\partial y}\quad \rm{and}\quad v=-\frac{\partial \psi}{\partial x}$$

Once you generate a random $\psi$, you can then use your Fourier representation (or finite differences) to compute the two derivatives. I haven't tried this kind of thing with truly random inputs, though. You might benefit by starting with a smoother but randomized input (like a sum of sine functions with random frequencies and coefficients).

Answer 2

You asked several questions in your post, so while Bill has provided a solution to the underlying problem, I feel someone should say some things about your stated questions too.

In what sense does the result approximate $J\nabla\psi$?

If you interpret the input array as the alias-free sampling of a periodic, bandlimited vector field $\mathbf f=\mathbf h+\nabla\phi+J\nabla\psi$, then the result is precisely the sampling of $\mathbf h+J\nabla\psi$. After all, the FFT of the sampled array is essentially the Fourier series of the underlying periodic function (assuming sufficient sampling).

Note that since $\mathbf f$ is periodic, $\mathbf h$ is merely a constant vector field.

Chris Beaumont's HH_DECOMP routine is supposed to use FFTs to perform Helmholtz-Hodge decomposition, but he also says (in the comments at the top of the code) that the method seems inaccurate.

I'm not sure what's going on there because I can't read IDL, but it looks like the test vector fields aren't periodic?

Likewise, there are more complicated variational methods to perform Helmholtz-Hodge decompositions of 2D data, which would seem to suggest that the simpler FFT method is somehow inadequate. Why?

I think you're worrying too much. On the second page of the paper you cite, they state, "Discrete Helmholtz-Hodge decomposition on regular grids has already been used in graphics (see [25, 10] for instance) and is relatively straightforward to implement with a finite-difference approach. It is, however, much harder to design a practical and accurate method for arbitrary grids."

What does the FFT method get wrong?

If you have periodic boundary conditions, nothing. In fact, under those conditions, spectral methods like this one give you exponential convergence as compared to the polynomial convergence of finite difference methods. They're harder to apply to arbitrary geometries, though.

And is it wrong to assume the harmonic component vanishes for my random noise?

Technically, yes, but in the periodic case the only possible harmonic vector field is a constant. So you're doing fine by leaving the DC component alone.