It has been seen that a vector field decomposition method called the Helmholtz Hodge Decomposition (HHD) can analyze scalar fields present universally in nature. It aids to reveal complex internal flows including energy flows in interference and diffraction optical fields. A gradient field defined in a region R, can be separated into solenoidal and irrotational components. HHD applied onto Magnetic Resonance Elasticity data can also aid to retain the curl field, while revealing the tissue elasticity in such medical measurements.

Keywords: helmholtz hodge decomposition, phase gradient, magnetic resonance, elastography

HHD: Helmholtz Hodge Decomposition; OAM: Orbital Angular Momentum; MRI: Magnetic Resonance Imaging; MRE: Magnetic Resonance Elastography

A vector field decomposition technique namely Helmholtz Hodge Decomposition (HHD) allows the field to be segregated into a solenoidal (divergence–free part) and an irrotational (curl–free part).^1–9 Many problems in electromagnetism, MRI.⁶ and fluid and smoke simulations.⁷ use this decomposition method to visualize real–time data. HHD aids to represent the homogeneous data explicitly by extracting the critical points like sources, sinks and vortices. It had been applied to polarized vector fields and to reconstruct phase for wavefront distortions. We used HHD on scalar optical fields and studied the Orbital angular momentum (OAM) in diffraction optics that has been reported.^10–15

An identity relates the OAM in an optical field to its phase and amplitude distribution.¹⁶ The vectorial nature of the fields is disregarded in such cases. We had shown the usefulness of the HHD in analyzing all such fields, including the ones obtained in interference optics where a single state of polarization (SOP) is assumed. We constructed a phase gradient field $\nabla \phi$  from a scalar field.^17,18 by using the relation

$\nabla \phi =\frac{\mathrm{Im}[{\psi }^{*}\nabla \psi ]}{I}$   (1)

Where $\psi$ a scalar is field resulting due to interference or diffraction and $I={\psi }^{\ast }\psi$  is the intensity distribution.

$\nabla \varphi$  that directed in the local propagation vector direction normal to the phase contour surfaces carries all the features of the wave. The phase gradient field in a singular beam has a non zero curl.^19–27 Hence, the solenoidal or the curl part that is an explicit component of HHD carries the circulating energy features of the field. The irrotational part reveals solely the spreading of energy, whether diverging or converging.

It is worth noting that different research groups have used different names (Helmholtz, Hodge, Helmholtz–Hodge or Hodge–Helmholtz) for this decomposition. The Helmholtz decomposition theorem suggested the segregation of a vector field, defined on real domains, into the solenoidal and the irrotational components. While the Hodge decomposition talked about a third component that is harmonic and is both solenoidal as well as irrotational. This decomposition theorem was defined for differential forms on Riemannian manifolds. Thus, the Hodge decomposition is the differential form analog of the Helmholtz decomposition in vector analysis. To the best of our knowledge, there is no origin of the Helmholtz–Hodge theorem or some formal merging of the names Helmholtz and Hodge.

The HHD technique that we adopted for decomposition is addressed in the next section. The approach was applied on to random fields in which both positive and negative curvatures were added. Spherical waves with positive and negative divergence, random wave fronts containing vortices with curling phase gradients, a vortex lattice field resulting from the superposition of plane waves were some of the fields that were investigated to test our HHD method. All the results showed clear explicitly segmented data and proved that our method worked really fine with scalar optical fields. We also believe that this method will definitely augment the flow visualization of velocity, pressure and temperature using optical methods. This technique can also be applied to process data fields as in medicine. It was anticipated and later seen that HHD using our least squares approach yields good results when applied to Magnetic Resonance Elastography (MRE) in human brain tissues. There are various methods to measure tissue elasticity that are used worldwide. Majority of them use sophisticated hardware and software and are not really cost effective. But our technique is very simple and does not require any complicated assessments and interpretations.

MRE consists of measurements of mechanical properties of a tissue using the Magnetic Resonance Imaging (MRI). An external actuator vibrates the tissue of interest and MRI data is obtained. This encodes the shear wave propagation into the MRI image. Tissue elasticity is then measured using wave inversion, after isolation of these shears waves. But the problem is ill–posed and extremely sensitive to denoising methods. Some noise solutions to sort out this issue have been used. Retaining the curl field to study tissue elasticity using HHD by our least square method also gives good results. The method has been explained in detail in the next section.

The Helmholtz Hodge Decomposition (HHD) is based on the Helmholtz theorem.^28–30 which states that a vector field which F is on a bounded domain V in R³, and is twice continuously differentiable, and whose divergence $\nabla .F=b(r)$  and curl $\nabla \times F=c(r)$ are known, can be segregated into components $f_1$ and $f_2$  determined by

$F=f_1+f_2\Rightarrow \nabla \varphi +\nabla \times A$   (2)

Where $\varphi \left(r\right)$  and $A\left(r\right)$  are scalar and vector potentials respectively, that can be obtained from the Poisson’s equations

$\varphi (r)=\frac{1}{4\pi }{\displaystyle \underset{V}{\int }\frac{b(r)}{r}}d{v}^{\prime }$    (3)

$A(r)=\frac{1}{4\pi }{\displaystyle \underset{V}{\int }\frac{c(r)}{r}}d{v}^{\prime }$   (4)

These potentials $\varphi \left(r\right)$  and $A\left(r\right)$  allow the field $f$ to be segregated into the curl free and divergence free components.

The boundary conditions imposed in HHD ensure a normal boundary flow on the curl free component and a tangential flow on the divergence free component. Considering $\widehat{n}$  as the outward normal to the boundary Ω this implies that for a unique decomposition,

1. The irrotational component f₁ is normal to the boundary $d\Omega$  of $\Omega$ ,e. ${\overrightarrow{f}}_1\times \widehat{n}=0$  , and
2. The solenoidal component f₂ is parallel to the boundary $d\Omega$  of $\Omega$ , i.e. ${\overrightarrow{f}}_2\cdot \widehat{n}=0$

The method adopted to solve the HHD problem involves minimizing the errors in the terms which are constructed from the initial guesses. $f_1$  and $f_2$  are considered as the initial guesses for the curl free and the divergence free component fields of a vector field $f$ .

The error terms/ residuals $\Vert \nabla \times f_1\Vert$, $\Vert \nabla \cdot f_2\Vert$  and $\Vert f_1+f_2-f\Vert$ are then reduced to a minimum.

In Cartesian coordinate system, the difference operator operating on a scalar function f, is defined as

$d\overrightarrow{f}=\frac{\partial f_x}{\partial x}\widehat{x}+\frac{\partial f_y}{\partial y}\widehat{y}+\frac{\partial f_z}{\partial z}\widehat{z}$   (5)

Using Finite Difference Approximation, the $\partial$ operator can be written as a matrix given by

$\partial =\frac{1}{2}\left[\begin{array}{ccccc}0& -1& 0& \cdots & 0\\ 1& 0& -1& \ddots & ⋮\\ 0& 1& 0& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & -1\\ 0& \cdots & 0& 1& 0\end{array}\right]$    (6)

This is true for 1D but for 3D, it is expanded as

${\partial }_{x,3D}=I_m\otimes I_m\otimes {\partial }_x\approx {\partial }_x$   (7)

where $I_m$ is an $m\times m$ Identity matrix and  is the Kronecker delta product. Similarly,

${\partial }_{y,3D}=I_m\otimes {\partial }_y\otimes I_m\approx {\partial }_y$   (8)

${\partial }_{z,3D}={\partial }_z\otimes I_m\otimes I_m\approx {\partial }_z$   (9)

The curl, in Cartesian coordinate system, is written using finite difference operator matrix as

$\nabla \times f=\left[\begin{array}{ccc}\widehat{x}& \widehat{y}& \widehat{z}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ f_x& f_y& f_z\end{array}\right]=\left[\begin{array}{ccc}0& -{\partial }_z& {\partial }_y\\ {\partial }_z& 0& -{\partial }_x\\ -{\partial }_y& {\partial }_x& 0\end{array}\right]\begin{array}{c}\end{array}\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]$   (10)

Similarly, for divergence,

$\nabla \cdot f=\frac{\partial f_x}{\partial x}+\frac{\partial f_y}{\partial y}+\frac{\partial f_z}{\partial z}=\left[\begin{array}{ccc}{\partial }_x& {\partial }_y& {\partial }_z\end{array}\right]\begin{array}{c}\end{array}\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]$

Thus, the HHD can then be summarized as:

$\left[\begin{array}{cccccc}0& -{\partial }_z& {\partial }_y& 0& 0& 0\\ {\partial }_z& 0& -{\partial }_x& 0& 0& 0\\ -{\partial }_y& {\partial }_x& 0& 0& 0& 0\\ 0& 0& 0& {\partial }_x& {\partial }_y& {\partial }_z\\ I& 0& 0& I& 0& 0\\ 0& I& 0& 0& I& 0\\ 0& 0& I& 0& 0& I\end{array}\right]\begin{array}{c}\end{array}\left[\begin{array}{c}{f}_{1x}\\ f_{1y}\\ f_{1z}\\ f_{2x}\\ f_{2y}\\ f_{2z}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ f_x\\ f_y\\ f_z\end{array}\right]$   (12)

The boundary condition considered is that the fields tend to go to zero at infinity. By applying the above mentioned boundary conditions, the system of equations can be efficiently solved. 12 can be solved as an equation $Px=Q$  where $P$  represents $\left[\begin{array}{cccccc}0& -{\partial }_z& {\partial }_y& 0& 0& 0\\ {\partial }_z& 0& -{\partial }_x& 0& 0& 0\\ -{\partial }_y& {\partial }_x& 0& 0& 0& 0\\ 0& 0& 0& {\partial }_x& {\partial }_y& {\partial }_z\\ I& 0& 0& I& 0& 0\\ 0& I& 0& 0& I& 0\\ 0& 0& I& 0& 0& I\end{array}\right]$ , while $\left[\begin{array}{c}{f}_{1x}\\ f_{1y}\\ f_{1z}\\ f_{2x}\\ f_{2y}\\ f_{2z}\end{array}\right]$  and $\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ f_x\\ f_y\\ f_z\end{array}\right]$  are represented by x and Q respectively.

Since the system of equations is not full rank, appropriate weights are applied to the residuals or error terms in order to get a unique solution. The weight parameters α, β and γ are defined for the curl free, divergence free and the sum residual respectively.

Thus, one minimizes the expression $\Vert W\left(Px-Q\right)\Vert$  where W is a diagonal matrix defined as

$W=diag\left(\begin{array}{cccc}\begin{array}{ccc}\begin{array}{cc}\begin{array}{c}\end{array}[\alpha & \alpha \end{array}& \alpha & \beta \end{array}& \gamma & \gamma & \gamma \end{array}]\begin{array}{c}\end{array}\right)$   (13)

The results obtained with the specific boundary condition are shown in next section. It is to be noted that the fields are assumed to go to zero at infinity. One clearly observes the vanishing of the normal component of the field in the solenoidal part and tangential component in the irrotational part at the boundaries. The boundary conditions imposed ensure a unique and orthogonal decomposition of the original field. The curl free part is the projection of the original field onto the space of solenoidal fields. Similarly, the divergence free part is the projection of the original field onto the space of irrotational fields. This is possible only when proper boundary conditions are satisfied.

HHD has, hitherto, seemed to be a great technique to extract and detect vortices, but we were unable to determine their strength. Vortices of any topological charge, however high it was, appeared similar. Thus, this seemed to be more of a visualization and analyzation technique than a method for measurement of the strength of the singularities.

In this section, the decomposition using HHD method for the some of the scalar fields is demonstrated. The simulation work has been done using Matlab.

The phase of the beam is shown in part (a) in Figure 1 (I) shows a Hodge decomposed beam for a positive spherical beam. The computed phase gradient field is shown in part (b). Part (c) shows the Hodge decomposed divergence free part, while part (d) shows the curl free component. As can be seen, the solenoidal field is zero in the core area far away from the boundary in case of spherical beams. The normal component of the field vanishes near to the boundary and the field lines tend to get parallel here. The irrotational component shows diverging field lines emanating from the center. As expected, the normal component tends to be perpendicular at all points on the boundary. We thus visualize that the Orbital Angular Momentum, which is associated with the circulating phase, is explicitly absent, en masse, in a spherical wave.

Figure 1 The Hodge decomposition applied to (I) a spherical wave with a positive divergence and (II) a vortex lattice field.

Speckle fields have also been decomposed and we envisage that HHD can be used to produce speckle free fields wherever required. An experimental result that was a vortex lattice field was considered next for decomposition and explicitly segregated components obtained, as shown in Figure 1 (II).

One can, on similar lines, decompose any field obtained in interference / diffractive optics, and study its propagation dynamics and other topological features. We envisaged that fields in MRE where one requires the curl waves could also be decomposed using our technique of decomposition. We tried our method on brain acquisition data collected from a lab and obtained results. We saw that the shear waves were more clarified as compared to ones obtained using other methods. It is attributed to the fact that the HHD removes the low–frequency artifact that causes overestimation of wavelengths. This helps to reduce the noise and thus, clarifies shear waves. These shear waves are then inverted to create a mechanical property map that gives an estimation of the brain mechanical stiffness. Figure 2 (I) and (II) show the shear waves obtained using our HHD method and the inverted map to see mechanical property in brain tissues.

Figure 2 The (I) clarified shear waves obtained using HHD and (II) inverted shear waves.

We had established that the Helmholtz Hodge decomposition can be used as a tool to analyze scalar fields and had demonstrated that the HHD can be used to segment the solenoidal and irrotational components in them. This has been solved in the rectangular coordinate system with general boundary conditions. The propagation of optical beams in circular cross–sectional channels is also of interest. The HHD method described above yields important results in the study of propagation of optical beams. The segregated component fields give a lot of insight into the generation and annihilation of optical vortices during propagation. Our HHD method has been applied on to the brain tissues and the clarified shear waves have been studied. These shear waves were more clarified as compared to the ones obtained using other methods. This is due to the fact that the HHD removes the low–frequency artifact that causes an overestimation of wavelengths. This helps to reduce the noise and thus, clarifies shear waves. These shear waves were then inverted to study the brain mechanical stiffness. The state of tissues at any instant of time is thus revealed using our method of wave separation. This can, hence, be used to study the state of tissues in case of treatment of any disease in the human body. An estimation of the state of health of tissues can thus be obtained using the decomposition technique.

Special Thanks to Dr. Eric Barnhill, University of Edinburgh for providing the brain acquisition data file.

None.

1. A Globus, C Levit, T Lasinki. A Tool for Visualizing the topology of Three–Dimensional Vector Fields. In Proceedings of IEEE on Visualization, pp. 1991:33–40.
2. K Polthier, E Preuss Identifying vector field singularities using adiscrete hodge decomposition. In: HC Hege & K Polthier (Eds.), Mathematics III and Visualization, pp. 2002:113–134.
3. Y Tong, S Lombeyda, AN Hirani, et al. Discrete multiscale vector field decomposition. ACM Trans Graphics. 2003;22(3):445–452.
4. F Petronetto, A Paiva, M Lage, et al. Meshless Helmholtz–Hodge decomposition. IEEE Transactions on Visualization and Computer Graphics. 2010;16(2):338–342.
5. FM Denaro. On the application of the Helmholtz–Hodge decomposition in projection methods for incompressible flows with general boundary conditions. Int J Numer Meth Fluids. 2003;43:43–69.
6. I Kaya, AP Santhanam, C Imielinska, et al. Modeling Air–flow in the tracheobronchial tree using Computational Fluid Dynamics. Medical Imaging Computing and Computer Aided Intervention Computational Biomechanics workshop (MICCAI), Greece, pp. 2007:142–151.
7. R Fedkiw, J Stam, HW Jensen. Visual simulation of smoke. Proceedings of the 28th annual conference on Computer graphics and interactive techniques, USA, p. 2001:15–22.
8. AM Stewart Angular momentum of Light. J Mod Opt. 2005;52(8):1145–1154.
9. M Hattori, S Komatsu. An exact formulation of a filter for rotation in phase gradients and its applications to wavefront reconstruction problems. J Mod Opt. 2003;50(11):1705–1723.
10. Monika Bahl, P Senthilkumaran. Helmholtz Hodge decomposition of scalar optical fields. J Opt Soc Am A Opt Img Sci Vis. 2012;29(11):2421–2427.
11. Monika Bahl, P Senthilkumaran. Focal plane internal energy flows of singular beams in astigmatically aberrated low NA systems. Opt Soc Am A Opt Img Sci Vis. 2014;31(9):2046–2054.
12. Monika Bahl, Brijesh Kumar Singh, Rakesh Kumar Singh, et al. Internal energy flows of coma affected singular beams in low NA systems. J Opt Soc Am A Opt Img Sci Vis. 2015;32(4):514–521.
13. P Senthilkumaran, Monika Bahl. Young's experiment with waves near zeros. Opt Exp. 2015;23(9):10968–10973.
14. Monika Bahl, Senthilkumaran. Energy circulations in singular beams diffracted through an isosceles right triangular aperture. Phys Review A. 2015;92:013831.
15. Brijesh Kumar Singh, Monika Bahl, DS Mehta, P Senthilkuamran Study of Internal energy Flows in Dipole Vortex Beams by Knife Edge Test. Optics Communications. 2013;293:15–21
16. G Molina–Terriza, JP Torres, L Torner. Twisted Photons. Nature Physics. 2007;3:305–310.
17. MV Berry. Optical Currents. J Opt A: Pure Appl Opt. 2009;11(9):1464–1475.
18. M Born, E Wolf Principle of optics. (7th edn), Cambridge University Press, India 2002:.
19. L Allen, MJ Padgett. The Poynting vector in Laguerre Gaussian beams and the interpretation of their angular momentum density. Opt Comm. 2000;184:67–71.
20. L Allen, MW Beijersbergen, RJC Spreeuw, et al. Orbital angular momentum of light and the transformation of Laguerre–Gaussian laser modes. Phys Rev A. 1992;45(11):8185–8189.
21. SF Arnold, L Allen, M Padgett. Advances in optical angular momentum. Laser Photon Rev. 2008;2(4):299–313.
22. M Padgett, L Allen. Optical tweezers and spanners. Physics World. 1997;10(9):35–38.
23. MS Soskin, MV Vasnetsov Progress in Optics. (1st edn), Elsevier, Netherlands, 2001:42: 423.
24. AM Yao, MJ Padgett. Orbital angular momentum: origins, behavior and applications. Advances in Optics and Photonics. 2011;3(2):161–204.
25. IV Basistiy, MS Soskin, MV Vasnetsov. Optical wavefront dislocations and their properties. Opt Commun. 1995;119:604–612.
26. IV Basistiy, VY Bazhenov, MS Soskin, et al. Optics of light beams with screw dislocations. Optics Communications. 1993;103(5–6):422–428.
27. P Senthilkumaran, F Wyrowski, H Schimmel. Vortex Stagnation problem in iterative Fourier transform algorithms. Opt Laser Eng. 2005;43:43–56.
28. GB Arfken, HJ Weber, Frank E. Harris Mathematical Methods for Physicists. (7th edn), Elsevier, Netherlands, p. 2005:1220.
29. Finite Difference method.
30. T Mc Graw, T Kawai, I Yassine, et al. Visualizing High–Order Symmetric Tensor Field Structure with Differential Operators. Journal of Applied Mathematics. 2011; 27.