# -*- coding: utf-8 -*-
"""
Created on Mon Feb 15 16:07:59 2021
@author: Anna Rauch
This script contains a set of general functions/methods usually related to
linear algebra and connectivty issues. These functions are independent from
the TIE-method, yet are used by it.
"""
import math
import numpy as np
import platform
import tkinter as tk
from tkinter import messagebox
[docs]
def angle2normal(azim, dip):
"""
Calculate the normal vector of a plane (or set of planes) defined with orientational angles in degrees.
Parameters
----------
azim : float or array_like
Angles of azimuth (dip azimuth) in degrees.
dip : float or array_like
Dip angles of the plane(s) in degrees.
Returns
-------
numpy.ndarray
Normal vector of the plane(s) in the form [normalx, normaly, normalz].
"""
dip = np.array(dip) * np.pi / 180
azim = np.array(azim) * np.pi / 180
if np.size(dip) == 1:
dip = [dip]
azim = [azim]
normalz = [np.cos(dp) for dp in dip]
normalx = [np.sin(dip[ind]) * np.sin(azim[ind]) for ind in range(np.size(dip))]
normaly = [np.sin(dip[ind]) * np.cos(azim[ind]) for ind in range(np.size(dip))]
return np.array([normalx, normaly, normalz]).T
[docs]
def angle2vect(trend, plunge):
"""
Calculate the directional vector (length = 1) from a line defined by angles - trend and plunge.
Parameters
----------
trend : float or array_like
Angles of trend (plunge azimuth) in degrees.
plunge : float or array_like
Plunge angles of the directional vector in degrees.
Returns
-------
tuple of numpy.ndarray
Coordinates of the oriented vector in the form (vx, vy, vz).
"""
# Convert angles to radians
trend = np.array(trend) * np.pi / 180
plunge = np.array(plunge) * np.pi / 180
# Handle single values or arrays
if np.size(trend) == 1:
trend = [trend]
plunge = [plunge]
# Calculate coordinates of the oriented vector
vz = [-np.sin(pl) for pl in plunge]
vx = [np.sin(trend[k]) * np.cos(plunge[k]) for k in range(np.size(trend))]
vy = [np.cos(trend[k]) * np.cos(plunge[k]) for k in range(np.size(trend))]
# Return the coordinates as a tuple of numpy arrays
return np.array(vx).T, np.array(vy).T, np.array(vz).T
[docs]
def angleBtwVec(v1, v2):
"""
Calculate the small angle between two directional oriented vectors.
Parameters
----------
v1 : array_like
First vector in the form (x, y, z).
v2 : array_like
Second vector in the form (x, y, z).
Returns
-------
float
Angle between the two vectors in degrees.
"""
if len(np.shape(v1)) > 1:
v1n = np.array([j / np.linalg.norm(j) for j in v1])
v2n = np.array([k / np.linalg.norm(k) for k in v2])
dotv = np.array(
np.round([np.dot(v1n[l], v2n[l].T) for l in range(len(np.shape(v1)))], 5)
)
angle = np.array([math.acos(dt) for dt in dotv])
else:
v1n = v1 / np.linalg.norm(v1)
v2n = v2 / np.linalg.norm(v2)
angle = math.acos(np.round(np.dot(v1n, v2n.T), 5))
return angle / np.pi * 180
[docs]
def azimRot(dip, trend, plunge):
"""Azimuth according to rotation axis and dip.
Calculates the azimuth of a plane with a certain dip
that goes through a rotation axis defined with trend and plunge angles.
Parameters
----------
dip : float
Angles of plane dip (angles in degrees).
trend : float
Angles of trend (plunge azimuth) of the rotational axis.
plunge : float
Plunge of the rotational axis (angles in degrees).
Returns
-------
float
Angles of plane dip azimuths (angles in degrees) that pass through the axis.
"""
trend = trend / 180 * np.pi
dip = dip / 180 * np.pi
plunge = plunge / 180 * np.pi
azim = trend - np.pi / 2 + math.asin(math.tan(plunge) / math.tan(dip))
return azim / np.pi * 180
[docs]
def dipRot(azim, trend, plunge):
"""Dip According To Rotation Axis And Azimuth.
Calculates the dip of a plane with a certain dip azimuth
that passes through a rotation axis defined with trend and plunge angles.
Parameters
----------
azim : float
Angles of plane dip azimuth (angles in degrees) (no strike).
trend : float
Angles of trend (plunge azimuth) of the rotational axis.
plunge : float
Plunge of the rotational axis (angles in degrees).
Returns
-------
float
Angles of plane dips (angles in degrees) that pass through the axis.
"""
trend = trend / 180 * np.pi
azim = azim / 180 * np.pi
plunge = plunge / 180 * np.pi
dip = math.atan(math.tan(plunge) / math.cos(azim - trend))
return dip / np.pi * 180
[docs]
def distance(P1, P2):
"""Calculate the shortest distance between two points in space.
Parameters
----------
P1 : list
Coordinates of the first point [x, y, z].
P2 : list
Coordinates of the second point [x, y, z].
Returns
-------
float
Distance between two points.
"""
return ((P1[0] - P2[0]) ** 2 + (P1[1] - P2[1]) ** 2) ** 0.5
[docs]
def greatCircle(azim, dip):
"""Extract coordinates of the great circle of a given plane orientation.
Parameters
----------
azim : float
Angle of azimuth (dip azimuth) of the plane.
dip : float
Dip angle of the plane.
Returns
-------
list
Coordinates of great circle points [x_stereo, y_stereo].
"""
dip = np.array(dip) * np.pi / 180
azim = -np.array(azim) * np.pi / 180 + np.pi / 2
N = 51
psi = np.linspace(0, np.pi, N)
radint = math.tan(dip) * np.array([math.sin(ps) for ps in psi])
radip = np.array([math.atan(radi) for radi in radint])
rproj = np.array([math.tan((np.pi / 2 - radp) / 2) for radp in radip])
x1 = rproj * np.array([math.sin(ps) for ps in psi])
y1 = rproj * np.array([math.cos(ps) for ps in psi])
x_stereo = np.array(
[x1[k] * math.cos(azim) + y1[k] * math.sin(azim) for k in range(np.size(x1))]
)
y_stereo = np.array(
[x1[k] * math.sin(azim) - y1[k] * math.cos(azim) for k in range(np.size(x1))]
)
return [x_stereo, y_stereo]
[docs]
def neighborPoints(j, rows, columns, connectivity):
"""Define array containing all neighbor indexes of a certain point in a matrix.
Parameters
----------
j : int
Index in the matrix of the point being analyzed (flattened matrix).
rows : int
Number of rows in the matrix.
columns : int
Number of columns in the matrix.
connectivity : str
Type of neighbor connectivity (8-connectivity or 4-connectivity).
Returns
-------
numpy.ndarray
Array with indexes of neighbors.
"""
lowerleft = (columns * rows) - columns
lowerright = columns * rows - 1
if connectivity == 4:
neigh = np.array([j + 1, j - 1, j - columns, j + columns])
if np.remainder(j + 1, columns) == 0:
neigh = np.array([j - 1, j - columns, j + columns])
if np.remainder(j, columns) == 0:
neigh = np.array([j + 1, j - columns, j + columns])
if j < columns:
neigh = np.array([j - 1, j + 1, j + columns])
if j > lowerleft:
neigh = np.array([j - 1, j + 1, j - columns])
if j == 0:
neigh = np.array([j + 1, j + columns])
if j == columns - 1:
neigh = np.array([j - 1, j + columns])
if j == lowerleft:
neigh = np.array([j + 1, j - columns])
if j == lowerright:
neigh = np.array([j - 1, j - columns])
if connectivity == 8:
neigh = np.array(
[
j + 1,
j - 1,
j - columns,
j + columns,
j - columns + 1,
j - columns - 1,
j + columns + 1,
j + columns - 1,
]
)
if np.remainder(j + 1, columns) == 0:
neigh = np.array(
[j - 1, j - columns, j + columns, j - columns - 1, j + columns - 1]
)
if np.remainder(j, columns) == 0:
neigh = np.array(
[j + 1, j - columns, j + columns, j - columns + 1, j + columns + 1]
)
if j < columns:
neigh = np.array(
[j - 1, j + 1, j + columns, j + columns + 1, j + columns - 1]
)
if j > lowerleft:
neigh = np.array(
[j - 1, j + 1, j - columns, j - columns + 1, j - columns - 1]
)
if j == 0:
neigh = np.array([j + 1, j + columns, j + columns + 1])
if j == columns - 1:
neigh = np.array([j - 1, j + columns, j + columns - 1])
if j == lowerleft:
neigh = np.array([j + 1, j - columns, j - columns + 1])
if j == lowerright:
neigh = np.array([j - 1, j - columns, j - columns - 1])
return neigh
[docs]
def neighborPointsD(j, rows, columns, dp):
"""Define array containing the indexes surrounding the index j at a given
distance dp in a matrix [rows,columns] (8-connectivity).
Parameters
----------
rows : int
Number of rows in the matrix.
columns : int
Number of columns in the matrix.
j : int
Index in matrix of the point being analyzed (flattened matrix).
dp : int
Pixel number at the wanted distance (1 being the closest neighbor possible).
Returns
-------
numpy.ndarray
Array with indexes of neighbors.
"""
neigh = np.array([j + dp, j - dp, j - dp * columns, j + dp * columns])
for i in range(1, dp + 1):
neigh2 = np.array(
[
j - dp * columns + i,
j - dp * columns - i,
j + dp * columns + i,
j + dp * columns - i,
]
)
neigh = np.concatenate((neigh, neigh2))
for i in range(1, dp + 1):
neigh2 = np.array(
[
j - (dp - i) * columns + dp,
j - (dp - i) * columns - dp,
j + (dp - i) * columns + dp,
j + (dp - i) * columns - dp,
]
)
neigh = np.concatenate((neigh, neigh2))
remj = np.remainder(j, columns)
remn = [np.remainder(n, columns) for n in neigh]
neigh = neigh.compress((remn <= remj + dp).flat)
remn = [np.remainder(n, columns) for n in neigh]
neigh = neigh.compress((remn >= remj - dp).flat)
neigh = neigh.compress((neigh < rows * columns).flat)
neigh = neigh.compress((neigh >= 0).flat)
neigh = np.unique(neigh)
return neigh
[docs]
def nmbNeighbors(binaryimage):
"""Extract a matrix (of size of the binary image) that contains the
number of positive neighbors for each pixel of the binary image.
Parameters
----------
binaryimage : numpy.ndarray
Matrix of zeros and ones (or TRUE's and FALSE's).
Returns
-------
numpy.ndarray
Matrix with the number of positive neighbors.
"""
bi = binaryimage.flatten()
i_true = (bi == 1).nonzero()[0]
nmbN = np.zeros(np.shape(bi))
for i in i_true:
neigh = neighborPoints(i, np.shape(binaryimage)[0], np.shape(binaryimage)[1], 8)
nmbN[i] = np.size((bi[neigh] == 1).nonzero())
nmbN_mat = np.reshape(nmbN, np.shape(binaryimage))
return nmbN_mat
[docs]
def normal2angle(normal):
"""Calculates the orientation (in dip azimuth and dip) of the plane
according to its normal (in x, y, z).
Parameters
----------
normal : numpy.ndarray
Normal vector: n = [x, y, z].
Returns
-------
Tuple[float, float]
Azimuth (dip azimuth) and dip of the plane expressed in angles.
"""
if np.ndim(normal) > 1:
normal = normal[0]
if normal[2] < 0:
normal = normal * np.array([-1, -1, -1])
if normal[1] == 0:
normal[1] = 0.000001
zvector = np.array([0, 0, 1])
dip = math.atan2(np.linalg.norm(np.cross(normal, zvector)), np.dot(normal, zvector))
dip = dip / np.pi * 180
azim = math.atan(abs(normal[0]) / abs(normal[1]))
azim = azim / np.pi * 180
if normal[0] > 0 and normal[1] <= 0:
azim = 180 - azim
if normal[0] <= 0 and normal[1] < 0:
azim = 180 + azim
if normal[0] <= 0 and normal[1] > 0:
azim = 360 - azim
return [azim, dip]
[docs]
def plungeRot(azim, dip, trend):
"""Calculates the plunge of an axis with a given trend and through which a
plane with a given orientation goes through.
Parameters
----------
azim : float
Angles of plane orientation (angles degree °).
dip : float
Angles of plane dip (angles degree °).
trend : float
Angles of trend (plunge azimuth) (angles degree °).
Returns
-------
float
Angles of axis plunge (angles degree °) that hosts the given plane.
"""
trend = trend / 180 * np.pi
azim = azim / 180 * np.pi
dip = dip / 180 * np.pi
plunge = math.atan(math.tan(dip) * math.cos(azim - trend))
return plunge / np.pi * 180
[docs]
def sortLine(ind, matSize):
"""Sorts a vector of points to connect them into a line and attribute
an order (2D only).
Parameters
----------
ind : array_like
Indexes of points in a matrix.
matSize : tuple
Size/shape of matrix.
Returns
-------
newi : ndarray
Newly sorted/ordered indexes.
"""
matx = np.arange(0, matSize[1])
maty = np.arange(0, matSize[0])
[x, y] = np.meshgrid(matx, maty)
x = x.flatten()
y = y.flatten()
mat = np.zeros(np.shape(x))
mat[ind] = 1
mat = np.reshape(mat, (matSize[0], matSize[1]))
newi = np.zeros(np.shape(ind))
start = []
for i in ind:
neigh = neighborPoints(i, matSize[0], matSize[1], 8)
npos = np.intersect1d(neigh, ind)
if np.size(npos) == 1:
start = i
break
if np.size(start) == 0:
start = ind[0]
neigh = neighborPoints(start, matSize[0], matSize[1], 8)
npos = np.intersect1d(neigh, ind)
newi[-1] = npos[0]
ind = np.array(ind).compress((np.array(ind) != npos[0]).flat)
count = 0
newi[count] = start
count = count + 1
ind = np.array(ind).compress((np.array(ind) != start).flat)
while np.size(ind) > 0:
d = [distance([x[start], y[start]], [x[ii], y[ii]]) for ii in ind]
dimin = [di for di in range(np.size(d)) if d[di] == np.min(d)]
if len(dimin) > 1:
for di in dimin:
neigh = neighborPointsD(ind[di], matSize[0], matSize[1], 8)
npos = np.intersect1d(neigh, ind)
di_neigh = np.array(dimin).compress((np.array(dimin) != di).flat)
npos = np.intersect1d(neigh, ind[di_neigh])
if len(npos) == 0:
dimin = di
break
newi[count] = ind[dimin]
count = count + 1
start = ind[dimin]
ind = np.delete(ind, dimin)
return newi
[docs]
def stereoLine(trend, plunge):
"""Extracts coordinates of a projected line on a stereonet.
Parameters
----------
trend : float
Angle of trend (plunge azimuth) of the line.
plunge : float
Plunge of the line.
Returns
-------
[x_stereo, y_stereo] : list
Coordinates of point(s) in a stereonet.
"""
trend = np.array(trend / 180 * np.pi)
plunge = np.array(plunge / 180 * np.pi)
radius = np.tan(np.pi * (np.pi / 2 - plunge) / (2 * np.pi))
x_stereo = np.sin(trend) * radius
y_stereo = np.cos(trend) * radius
if trend.compress((trend > np.pi / 2).flat) and trend.compress(
(trend <= np.pi).flat
):
trend = np.pi - trend
x_stereo = np.sin(trend) * radius
y_stereo = -np.cos(trend) * radius
if trend.compress((trend > np.pi).flat) and trend.compress(
(trend <= np.pi / 2 * 3).flat
):
trend = trend - np.pi
x_stereo = -np.sin(trend) * radius
y_stereo = -np.cos(trend) * radius
if trend.compress((trend > np.pi / 2 * 3).flat):
trend = 2 * np.pi - trend
x_stereo = -np.sin(trend) * radius
y_stereo = np.cos(trend) * radius
return [x_stereo, y_stereo]
[docs]
def vect2angle(v):
"""Extracts from vectors (vx, vy, vz) its axis directions in angles (geology style).
Parameters
----------
v : list
Vector: v = [vx, vy, vz].
Returns
-------
trend, plunge : tuple
Angles of trend (plunge azimuth) and plunge of the vector.
"""
if np.ndim(v) > 1:
v = v[0]
if v[2] > 0:
v = v * np.array([-1, -1, -1])
if v[1] == 0:
trend = 90
else:
trend = abs(math.atan(v[0] / v[1]))
trend = trend / np.pi * 180
plunge = math.atan(abs(v[2]) / np.sqrt(v[0] ** 2 + v[1] ** 2))
plunge = plunge / np.pi * 180
if v[0] > 0 and v[1] <= 0:
trend = 180 - trend
if v[0] <= 0 and v[1] <= 0:
trend = 180 + trend
if v[0] <= 0 and v[1] > 0:
trend = 360 - trend
return trend, plunge
[docs]
def Mbox(title, text, style):
"""Creates a message box.
Parameters
----------
title : str
Title/name of the box.
text : str
Message text.
style : int
Style of text (bold, fontsize, ...).
Returns
-------
None
Message box appears.
"""
if platform.system() == "Windows":
import ctypes
ctypes.windll.user32.MessageBoxW(0, text, title, style)
else: # Assuming other platforms support tkinter
root = tk.Tk()
root.withdraw() # Hide the main window
messagebox.showinfo(title, text)
root.destroy()