Algorithms

In this chapter, we describe the algorithms used in IGMAS+ for the calculation of potential fields interpolation, triangulation, and isosurface extraction.

Anomaly Calculation

Gravitational Constant

The default of the gravitational constant used for any gravitational calculation is $G$ = 6.67384 $\cdot$ 10$^{-11}$ m$^3$ kg$^{-1}$ s$^{-2}$.

Due to compatibility reasons, mantissa of this constant can be changed by user via the property Triangle Kernel of the Model entry in the Object Tree, in the Property Editor Tab:

Warning

Modification of the gravitational constant is valid for the current session only. Be careful, you might turn the gravitational anomalies into nonsense!

Calculation of Anomalies of Triangulated Polyhedrons

This section lists all available algorithms for calculation of the effect for triangulated polyhedrons and some background on how to improve performance for interactivity.

For interactive work it is essential that after model changes the recalculation of the model is done immediately. It is important that this performance is achieved not only when a single point in the model is altered, but also when bigger parts are changed at once (e.g., several triangles). The basis for a fast recalculation is a changed-only recalculation (fortunately possible in gravity and magnetics):

1. Identify parts of the model which have been changed
2. Subtract these effects from the actual field
3. Add the newly calculated effects.

One important aspect is the speed of the calculation of the gravity effect for triangles (single- and multi-z surfaces). This can be achieved, for example, by parallelizing calculations on all available CPU-cores and/or GPU via OpenCL (Alvers et al. 2014¹). Another way of speeding up the calculations is to approximate them. Gaussian Quadrature can approximate the exact calculation of the surface integrals over the triangles, if certain conditions hold. Approximations have to be applied carefully in order to prevent errors. For the recalculation of the model (after model changes) this is often less critical. The introduced error of subtracted and newly added calculation are quite often very similar and cancel out each other at least partly. Parallelization and approximation can obviously be combined. Depending on the available hardware one should decide which algorithms should be applied.

IGMAS+ offers 6 different algorithms for calculation of the effect for triangulated polyhedrons:

1. Triangle Kernel (Multicore): Multicore implementation of the algorithm of Götze and Lahmeyer (1988)².

   Note

   Triangle Kernel (OpenCL) algorithm is the default option.

2. Triangle Kernel (OpenCL): OpenCL implementation of the algorithm of Götze and Lahmeyer (1988)².

   Warning

   Triangle Kernel (OpenCL) algorithm requires double precision, check if your graphics card supports it.

   Warning

   Calculation of the following potential fields is not supported with Triangle Kernel (OpenCL):

   • Magnetic: $MAG_{totr}$
   • Magnetic gradients: $M_{xx}$, $M_{yy}$, $M_{zz}$, $M_{xy}$, $M_{zx}$, $M_{zy}$
   • Geoid

3. Gaussian (Approximation) Quadrature (Multicore): Multicore implementation of triangle approximation using 3-points Gaussian Quadrature.

   Warning

   Calculation of the following potential fields is not supported with Gaussian (Approximation) Quadrature (Multicore):

   • Magnetic gradients: $M_{xx}$, $M_{yy}$, $M_{zz}$, $M_{xy}$, $M_{zx}$, $M_{zy}$
   • Geoid

4. Gaussian (Approximation) Quadrature (OpenCL): OpenCL implementation of triangle approximation using 3-points Gaussian Quadrature.

   Warning

   Gaussian (Approximation) Quadrature (OpenCL) algorithm requires double precision, check if your graphics card supports it.

   Warning

   Calculation of the following potential fields is not supported with Gaussian (Approximation) Quadrature (OpenCL):

   • Magnetic: $MAG_{totr}$
   • Magnetic gradients: $M_{xx}$, $M_{yy}$, $M_{zz}$, $M_{xy}$, $M_{zx}$, $M_{zy}$
   • Geoid

5. Mix exact and approximated (Multicore): This algorithm mixes exact and approximate calculations of Gaussian Quadrature using Multicore.
   Two options are taken into account to select the exact calculation of the effect triangles:

   • a) a threshold for maximum triangle areas (in model units)
   • b) an error threshold for automated error estimation (in mGal), minDepth (in model units).
     That means for deeper parts of the model the calculation can be approximated by Gaussian Quadrature.

   Warning

   Calculation of the following potential fields is not supported with Mix exact and approximated (Multicore):

   • Magnetic gradients: $M_{xx}$, $M_{yy}$, $M_{zz}$, $M_{xy}$, $M_{zx}$, $M_{zy}$
   • Geoid

6. Spherical Triangle Kernel (Multicore): Multicore implementation of the algorithm of Götze and Lahmeyer (1988)² on a sphere.
   IGMAS+ will subdivide large triangles into smaller triangles so that the geometry can be projected onto the curvature of the Earth without accuracy problems.
   The user has to define the maximum length of the triangles (in model units, the default value is 10).

   Note

   Numerical experiments indicate that a maximum length of 10 km gives satisfying results.

   Note

   The coordinate projection of the project is used for spherical calculations:

   • if the projection is known (e.g. UTM, GK, EPSG), it will be used to exactly calculate the 3D position of the vertices and stations
   • if the projection is unknown, a simple sphere will be assumed having a radius of 6371 km.

   Warning

   Interactive modelling is not supported with the Spherical Triangle Kernel (Multicore) algorithm.

Change Triangle Kernel Algorithm

To change the Triangle Kernel algorithm in IGMAS+ use Property Editor Tab of the Model entry from the Object Tree, click ... button near Algorithm property in the Triangle Kernel category:

Alternatively, open the wizard using Research → Wizard [Triangle algorithm].

Select between the available options but take into account the warnings and instructions for each of the algorithm:

Warning

Take into account the warnings and instructions listed in the window for each algorithm.

Change OpenCL Configuration

OpenCL algorithms require double precision calculation, however not all graphics card support it.
To change it, use Research → Plugin → OpenCL Configuration:

Notice the difference between different graphics cards: Intel Iris Xe doesn't support double precision operations, while NVIDIA RTX A2000 supports it:

Intel Iris XeNVIDIA RTX A2000

Calculation of Invariants, Horizontal Gradient and Horizontal Directive Tendency

The invariants $Inv_0$, $Inv_1$, $Inv_2$ are combinations of gravity gradients components, which are the second derivatives of the potential.
Interpretation of invariants can give more information about the high-frequency part of the anomaly field.
Calculations of invariants and gradients are based on Pedersen & Rasmussen (1990)³.

$$\begin{split} Inv_0 &= G_{xx} + G_{yy} + G_{zz}\\ Inv_1 &= G_{xx} G_{yy} + G_{yy} G_{zz} + G_{xx} G_{zz} - G_{xy}^{2} - G_{zy}^{2} - G_{zx}^{2}\\ Inv_2 &= G_{xx} (G_{yy} G_{zz} - G_{yz}^{2}) + G_{xy} (G_{yz} G_{xz} - G_{xy} G_{zz})\\ & + G_{xz} (G_{xy}G_{yz} - G_{xz} G_{yy})\\ \end{split}$$

The horizontal gradient and the horizontal directive tendency are given by:

$$\begin{split} HG_z &= \sqrt{G_{zx}^{2} + G_{zy}^{2}}\\ HDT &= \sqrt{(G_{xx} - G_{yy})^{2} + (2G_{xy})^{2}}\\ \end{split}$$

Calculation of Anomalies of Voxels

Anomalies of Mass Points (Homogeneous Spheres)

Each voxel of the voxel cube is approximated by a mass point - a sphere with its volume being identical to the volume of the voxel.

Let us define the following variables:

• $R$ – radius of the sphere
• $\rho$ – mass density of the sphere
• $r = \sqrt{x^2 + y^2}$ – horizontal distance between the centre of the sphere and the station
• $x$, $y$, $z$ – distance components between the centre of the sphere and the station
• $G$ – gravitational constant

The following equations are used for the calculation of the voxel effects:

Components of the gravity (all to be multiplied by the gravitational constant):

$$g_x = \frac{4}{3} \pi \rho R^3 \frac{x}{(r^2 + z^2)^{3/2}}$$

$$g_y = \frac{4}{3} \pi \rho R^3 \frac{y}{(r^2 + z^2)^{3/2}}$$

$$g_z = \frac{4}{3} \pi \rho R^3 \frac{-z}{(r^2 + z^2)^{3/2}}$$

Gradients of the gravity (all to be multiplied by the gravitational constant):

$$V_{xx} = \frac{4}{3} \pi \rho R^3 \frac{r^2 + z^2 - 3x^2}{(r^2 + z^2)^{5/2}}$$

$$V_{yy} = \frac{4}{3} \pi \rho R^3 \frac{r^2 + z^2 - 3y^2}{(r^2 + z^2)^{5/2}}$$

$$V_{zz} = \frac{4}{3} \pi \rho R^3 \frac{r^2 - 2z^2}{(r^2 + z^2)^{5/2}}$$

$$V_{xy} = \frac{4}{3} \pi \rho R^3 \frac{-3xy}{(r^2 + z^2)^{5/2}}$$

$$V_{zx} = \frac{4}{3} \pi \rho R^3 \frac{3xz}{(r^2 + z^2)^{5/2}}$$

$$V_{zy} = \frac{4}{3} \pi \rho R^3 \frac{3yz}{(r^2 + z^2)^{5/2}}$$

Induced magnetic field:

$$M_x = \frac{4}{3} \pi \sigma |H| R^3 \frac{3x(H_x x + H_y y + H_z z) - H_x (r^2 + z^2)}{(r^2 + z^2)^{5/2}}$$

$$M_y = \frac{4}{3} \pi \sigma |H| R^3 \frac{3y(H_x x + H_y y + H_z z) - H_y (r^2 + z^2)}{(r^2 + z^2)^{5/2}}$$

$$M_z = \frac{4}{3} \pi \sigma |H| R^3 \frac{3z(H_x x + H_y y + H_z z) - H_z (r^2 + z^2)}{(r^2 + z^2)^{5/2}}$$

$$M_{total} = \frac{4}{3} \pi \sigma |H| R^3 \frac{3(H_x x + H_y y + H_z z)^2 - (r^2 + z^2)}{(r^2 + z^2)^{5/2}}$$

where:

• $H_x$, $H_y$, $H_z$ – direction of the external field (components of the unit vector)
• $|H|$ – magnitude of the external field
• $\sigma$ – magnetic susceptibility

Interpolation

Block Average Filter

The Block Average Filter is applied in when importing horizons with regularly spaced points.
If the point data are highly oversampled (i.e. have too high number of points), they have to be down-sampled during import.

For that, each imported point is assigned to the nearest regular vertex. Finally, for each vertex the average of all 3 coordinates ($x$, $y$ and $z$) is calculated, and only these average values are stored and used.

Note

The averaging within the vertex surrounding area is calculated successively during import procedure, so that the required memory does not depend on the number of the points to be imported.

Mundry interpolation

The interpolation algorithm based on the work of Mundry (1970)⁴ is used when importing horizons with irregular distribution of points.

Triangulation

Triangulation is a subdivision of a planar object into triangles (or simplices in a higher-dimension geometry).

Triangle Orientation

A triangle has right and a left hand side, depending on the order of the vertex definition:

Each triangulated surface has one body on its right and another body on its left hand side.

Remember the Right Hand Rule:

If the fingers follow the order of the vertices (1 → 2 → 3), the thumb shows the direction of the positive normal (the right hand side).

Right Hand Rule. Image is taken from the GEO1004 course of TU Delft.

Note

Usually the user should not be concerned about the triangle orientation - it will be chosen correctly by IGMAS+.
However, if a triangulation is imported from another program (e.g. GOCAD), it might be necessary to think about the orientation.

The Object Tree shows the orientation of each interface: it shows two body names for each single interface, which are separated by the body separator < >. The first name (left) specifies the name of the body on the left hand side, the second name (right) the body on the right hand side:

Extraction of Isosurfaces from Voxel Cubes

A voxel cube may be transformed into triangulated isosurfaces using the Voxel Cube Isofurface Extraction Wizard. Use icon to open it:

You can set the limits for the voxel cell values to be included (Lower Value for lower limit and Upper Value for upper limit).
It is also possible to choose, whether the resulting isosurface bodies are closed at the model boundaries or not (Close the body on the sides checkbox).

Using Mesh Simplification checkbox helps reducing the number of triangles after geometry extraction.

Note

The function uses the Marching Cubes algorithm.

References

---

1. M. R. Alvers, H.-J. Götze, L. Barrio-Alvers, S. Schmidt, B. Lahmeyer, and C. Plonka. A novel warped-space concept for interactive 3D-geometry-inversion to improve seismic imaging. First Break, 2014. doi:10.3997/1365-2397.32.4.74375. ↩

2. H.‐J. Götze and B. Lahmeyer. Application of three‐dimensional interactive modeling in gravity and magnetics. Geophysics, 53(8):1096–1108, 1988. doi:10.1190/1.1442546. ↩↩↩

3. L. B. Pedersen and T. M. Rasmussen. The gradient tensor of potential field anomalies: some implications on data collection and data processing of maps. GEOPHYSICS, 55(12):1558–1566, December 1990. doi:10.1190/1.1442807. ↩

4. E Mundry. Zur automatischen Herstellung von Isolinienplänen. BEIH. GEOL. JB., pages 77–93, 1970. ↩