2

Say I have a vector field with non-zero curl, therefore the potential function depends on the path I choose to integrate. In this paper the authors proposed to project the vector field into a gradient field of Fourier basis functions in order to make the potential function independent of the path of integration.

I understand the idea but cannot understand the math of what that projection means. What I think the authors are doing is to approximate the given data using Fourier basis functions (that's why I was asking about Fourier series in another thread) in order to have a smooth function that best fits the given points and complies with the zero-curl of a gradient field. What I don't see is how many basis functions they use in order to approximate the data points, it could be $y_i = C_0 + C_1e^{jwx_i}$ or $y_i = C_0 + C_1e^{jwx_i} + C_2e^{2jwx_i}$ or any other. Moreover, it seems to me that they compute the FFT of the vector field and use all the coefficients of it, so, the function will fit all the points but still will be non-integrable.

I'm sure I am missing something and mixing maybe some concepts, can someone help me unwrapping it all to understand what the idea of this projection is?

Community
  • 1
BRabbit27
  • 1,029
  • 2
  • 11
  • 20

2 Answers2

2

In the matrix case, the least square solution of an overdetermined inconsistent system $Ax=b$ requires the solution of $(A^TA)x=A^Tb$. The projection $A(A^TA)^{-1}Ab$ reduces $b$ to the closest vector that is a consistent right side.

As in the matrix case, the solution of the overdetermined inconsistent system $∇z=V$ is reduced to a solvable problem by applying the adjoint operator, here $$∇⋅ V=∇⋅ ∇z=Δz.$$

If you represent the vector field $V$ by $$V(x)=\sum_{i\in I}C(ω_i)e^{jω_ix}$$ then the divergence of this vector field is $$∇⋅ V(x)=\sum_{i\in I}j(ω_i⋅ C(ω_i))e^{jω_ix}$$ and if a solution is tried that has the form $$z(x)=\sum_{i\in I}A(ω_i)e^{jω_ix},$$ one finds for the Laplacian $$Δz(x)=-|ω_i|^2\sum_{i\in I}A(ω_i)e^{jω_ix}.$$

Putting this together by comparing coefficients of the same basis functions, the solution has the form $$ A(ω_i)=-j\frac{ω_i⋅ C(ω_i)}{ω_i⋅ω_i} $$ which is formula 21 of the paper, but without (explicitly) using the orthogonality.

Note that $x$, $ω_i$ and $C(ω_i)$ are two- or higher dimensional.

You can prune the coefficients of the Fourier series by removing all that fall below some threshold. Orthogonality tells you that you can control the error of that procedure by sorting the magnitudes of the coefficients or coefficient vectors, start the pruning with the smallest and add the squares of the coefficients.

Lutz Lehmann
  • 6,044
  • 1
  • 17
  • 25
1

To give some background to the answer given by LutzL, here is how you'd set up your linear system to start with. You want to find an image $L$ whose gradient $\nabla L$ is as close as possible to your target vector field $T=(T_x, T_y)$. In other words:

$$ \arg\min_{L} \left\|\nabla L - T \right\| $$

You'll actually end up solving the following problem:

$$ \Delta L = \frac{\partial T_x}{\partial x} + \frac{\partial T_y}{\partial y} $$

Where $\Delta$ is the Laplacian operator. Here's why. Let's start from the beginning using matrix algebra (with a slightly imperfect notation):

$$ \arg\min_L \left\| \begin{bmatrix} D_x \\ D_y \end{bmatrix} L - \begin{bmatrix} T_x \\ T_y \end{bmatrix} \right\|. $$

Where $D_x$ and $D_y$ are $x$ and $y$ derivative operator matrices, $L$ is your image represented as a column vector, and $T_x$ and $T_y$ are also in vector form.

You then proceed minimising the above as you would with any least-squares problem, which gives you the linear system:

$$ (D_x D_x + D_y D_y)L = D_x T_x + D_y T_y $$

where indeed $D_x D_x + D_y D_y$ is the matrix form of the Laplacian operator $\Delta$ used above.

Turns out this is a circulant matrix, which has the property that the linear system can be solved using the FFT as described by LutzL.

Roberto
  • 11
  • 1