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 [6], and fluid and smoke simulations [7] 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 [16]. 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)

$\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 R3, 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$ and $f_2$  determined by

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

$\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)

1. The irrotational component $f_1$ is normal to the boundary $d\Omega$  of$\Omega$ , i.e. ${\overrightarrow{f}}_1\times \widehat{n}=0$  , and
2. The solenoidal component $f_2$  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$  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)

$\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)

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

${\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)

$\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)

$\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]$                      $$ (11)

$\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)

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.

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.