Skip to content
Snippets Groups Projects
Commit e59c68a3 authored by hanste's avatar hanste
Browse files

added initial commit of relational shape measure

parents
No related branches found
No related tags found
No related merge requests found
Showing
with 1285 additions and 0 deletions
### Brief Description
The Relational Shape Measure (https://arxiv.org/abs/2008.03927) is a measuring method to measure the shape of an object or a group of objects as a function to the distance to some reference object. The method uses the Hausdorff measures of all possible intersections of the object and a parameterized dilation of the reference object. That is volume, outer surface area, the surface area of intersection, length of intersection contour.
### Use Cases
The method can answer questions like:
- Are the objects distributed more or less closely near the reference object.
- Are the objects larger or smaller near a reference object.
- Does the object have large facing surfaces towards the reference object.
### How to use it
To use the method it requires you to have vertices and faces for the two exterior measures and vertices and tetrahedra for the interior measures. It also requires the distance of each vertex to the reference object. These can be generated by something like mesh-to-point algorithms or interpolation from a distance field from a distance transform method. See the examples.py for a few examples on synthetic data.
\ No newline at end of file
import numpy as np
from .Rotations import axisAngleRot
def RandomRotationMatrix1():
rv = np.random.rand(3)
pi = np.pi
R = np.zeros((3,3))
theta = rv[0] * pi * 2
phi = rv[1] * pi * 2
z = rv[2] * 2.0
r = np.sqrt(z)
Vx = np.sin(phi) * r
Vy = np.cos(phi) * r
Vz = np.sqrt(2.0 - z)
st = np.sin(theta)
ct = np.cos(theta)
Sx = Vx * ct - Vy * st
Sy = Vx * st + Vy * ct
R[0,0] = Vx * Sx - ct
R[0,1] = Vx * Sy - st
R[0,2] = Vx * Vz
R[1,0] = Vy * Sx + st
R[1,1] = Vy * Sy - ct
R[1,2] = Vy * Vz
R[2,0] = Vz * Sx
R[2,1] = Vz * Sy
R[2,2] = 1.0 - z
return R
def RandomRotationMatrix2():
rv = np.random.rand(3)
pi = np.pi
phi = rv[0]*pi*2
theta = np.arccos(2*rv[1] - 1)
thetaAroundVector = rv[2]*pi*2
x = np.array([1,0,0])
y = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
angleBetween = np.arccos(np.dot(x,y)/(np.linalg.norm(x)*np.linalg.norm(y)))
#print angleBetween
cross = np.cross(x,y)
#print cross
cross = cross/np.linalg.norm(cross)
R1 = axisAngleRot(cross,angleBetween)
R2 = axisAngleRot(y,thetaAroundVector)
#print np.linalg.det(R1),np.linalg.det(R2)
return R2.dot(R1)
def RandomRotationMatrix():
return RandomRotationMatrix1()
def RandomRotatedUnitVector3d():
x = np.array([1.0, 0.0, 0.0])
R = RandomRotationMatrix()
return R.dot(x)
def randomDirection():
x = np.random.rand(3,1)*2 - 1
l = np.linalg.norm(x)
while abs(l-1)>0.05:
x = np.random.rand(3,1)*2 - 1
l = np.linalg.norm(x)
x = x / l
return x
def Euler_zxz(alpha, beta, gamma):
# Proper euler extrinsic coordinates (z-x-z)
R = np.dot(np.dot(rotk(gamma), roti(beta)), rotk(alpha))
return R
def InvEuler_zxz(R):
alpha = np.arccos( R[2,1] / np.sqrt(1-R[2,2]**2))
beta = np.arccos( R[2,2])
gamma = np.arccos(-R[1,2] / np.sqrt(1-R[2,2]**2))
return alpha, beta, gamma
def roti(a):
a = a[0]
cos = np.cos
sin = np.sin
Ri = [[1, 0, 0],
[0, cos(a), -sin(a)],
[0, sin(a), cos(a)]]
return np.array(Ri)
def rotj(b):
b = b[0]
cos = np.cos
sin = np.sin
Rj = [[ cos(b), 0, sin(b)],
[ 0, 1, 0],
[-sin(b), 0, cos(b)]]
return np.array(Rj)
def rotk(c):
c = c[0]
cos = np.cos
sin = np.sin
Rk = [[cos(c), -sin(c), 0],
[sin(c), cos(c), 0],
[ 0, 0, 1]]
return np.array(Rk)
def randomCoordinateSystem():
x1 = randomDirection()
x2 = randomDirection()
while abs(np.dot(np.transpose(x1),x2)) > 0.1:
x2 = randomDirection()
x2 = x2 - (np.dot(np.transpose(x1),x2)) * x1
x2 = x2/np.linalg.norm(x2)
x3 = np.cross(x1.T, x2.T).T
alpha, beta, gamma = InvEuler_zxz(np.stack([x1,x2,x3], axis=1))
return alpha, beta, gamma, x1, x2, x3
def jonRandomRotation():
alpha, beta, gamma, x1, x2, x3 = randomCoordinateSystem()
R = Euler_zxz(alpha, beta, gamma)
return R
if __name__ == '__main__':
import time
t1 = time.time()
for i in range(10000):
RandomRotationMatrix1()
r1_time = time.time() - t1
t2 = time.time()
for i in range(10000):
RandomRotationMatrix2()
r2_time = time.time() - t2
print(r1_time, r2_time)
import numpy as np
def axisAngleRot(x,t):
ux = x[0]
uy = x[1]
uz = x[2]
R = np.zeros((3,3))
ct = np.cos(t)
st = np.sin(t)
omct = (1 - ct)
R[0,:] = np.array([ct + ux**2*omct, ux*uy*omct - uz*st, ux*uz*omct + uy*st])
R[1,:] = np.array([uy*ux*omct + uz*st, ct + uy**2*omct, uy*uz*omct - ux*st])
R[2,:] = np.array([uz*ux*omct - uy*st, uz*uy*omct + ux*st, ct + uz**2*omct])
return R
\ No newline at end of file
from .RandomRotations import RandomRotationMatrix
\ No newline at end of file
File added
File added
File added
import numpy as np
import matplotlib.pyplot as plt
from RandomRotations import RandomRotationMatrix
from RPSM import measure_interior, measure_exterior
def Test_Cube_Integration_Interior(do_plots):
r = np.linspace(0, 2, 100, endpoint=True)
T1 = np.array([[0,2,3,7],
[0,2,6,7],
[0,4,6,7],
[0,4,5,7],
[0,1,5,7],
[0,1,3,7]], dtype=np.int32)
T2 = np.array([[2,4,6,7],
[0,1,2,4],
[1,2,4,7],
[1,4,5,7],
[1,2,3,7]], dtype=np.int32)
if do_plots:
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
for i in range(100):
V = np.array([[0,0,0],
[0,0,1],
[0,1,0],
[0,1,1],
[1,0,0],
[1,0,1],
[1,1,0],
[1,1,1]], dtype=float)
D = V[:,0] + 0.5
R = RandomRotationMatrix()
V = np.dot(V, R.T)
ra, volume1, area1 = measure_interior(V, T1, D, r, adaptive=True)
ra, volume2, area2 = measure_interior(V, T2, D, r, adaptive=True)
if do_plots:
ax1.set_title('Volume')
ax1.plot(ra, volume1, label='type 1')
ax1.plot(ra, volume2, ls='--', label='type 2')
ax2.plot(ra, area1, label='type 1')
ax2.plot(ra, area2, ls='--', label='type 2')
if do_plots:
plt.show()
def Test_Cube_Integration_Exterior(do_plots):
r = np.linspace(0, 2, 100, endpoint=True)
F1 = np.array([[0,1,2],
[1,2,3],
[0,1,4],
[1,4,5],
[0,2,4],
[2,4,6],
[1,3,5],
[3,5,7],
[2,3,6],
[3,6,7],
[4,5,6],
[5,6,7]], dtype=np.int32)
F2 = np.array([[0,2,3],
[0,1,3],
[0,1,5],
[0,4,5],
[0,2,6],
[0,4,6],
[1,3,7],
[1,5,7],
[2,3,7],
[2,6,7],
[4,5,7],
[4,6,7]], dtype=np.int32)
if do_plots:
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
for i in range(100):
V = np.array([[0,0,0],
[0,0,1],
[0,1,0],
[0,1,1],
[1,0,0],
[1,0,1],
[1,1,0],
[1,1,1]], dtype=float)
D = V[:,0] + 0.5
R = RandomRotationMatrix()
V = np.dot(V, R.T)
ra, area1, circ1 = measure_exterior(V.copy(), F1, D.copy(), r, adaptive=True)
ra, area2, circ2 = measure_exterior(V.copy(), F2, D.copy(), r, adaptive=True)
if do_plots:
ax1.plot(ra, area1)
ax1.plot(ra, area2)
ax2.plot(ra, circ1)
ax2.plot(ra, circ2)
if do_plots:
plt.show()
def Test_Random_Integrations(do_plots):
np.random.seed(42)
N = 25 # num vertices
M = N//3 # num faces / tetrahedra
R = 15 # max r value (distance) to measure
K = 100 # num reading points
## EXTERIOR TEST
V = np.random.uniform(0, 10, (N, 3))
F = np.random.randint(0, N, (M, 3))
D = np.random.uniform(0, 10, N)
r = np.linspace(0, R, K)
_, mu10, mu11 = measure_exterior(V, F, D, r)
r_adapt, mu10_adapt, mu11_adapt = measure_exterior(V, F, D, r, adaptive=True)
if do_plots:
plt.subplot(1,2,1)
plt.title('volume')
plt.plot(r, mu10, label='fixed r')
plt.plot(r_adapt, mu10_adapt, label='adaptive r')
plt.legend()
plt.subplot(1,2,2)
plt.title('inner surface cut')
plt.plot(r, mu11, label='fixed r')
plt.plot(r_adapt, mu11_adapt, label='adaptive r')
plt.legend()
plt.show()
## INTERIOR TEST
V = np.random.uniform(0, 10, (N, 3))
T = np.random.randint(0, N, (M, 4))
D = np.random.uniform(0, 10, N)
r = np.linspace(0, R, K)
_, mu00, mu01 = measure_interior(V, T, D, r)
r_adapt, mu00_adapt, mu01_adapt = measure_interior(V, T, D, r, adaptive=True)
if do_plots:
plt.subplot(1,2,1)
plt.title('volume')
plt.plot(r, mu00, label='fixed r')
plt.plot(r_adapt, mu00_adapt, label='adaptive r')
plt.legend()
plt.subplot(1,2,2)
plt.title('inner surface cut')
plt.plot(r, mu01, label='fixed r')
plt.plot(r_adapt, mu01_adapt, label='adaptive r')
plt.legend()
plt.show()
if __name__ == '__main__':
do_plots = True
import time
t0 = time.time()
Test_Cube_Integration_Interior(do_plots)
Test_Cube_Integration_Exterior(do_plots)
Test_Random_Integrations(do_plots)
t1 = time.time()
print(t1 - t0, 's')
\ No newline at end of file
###############################################################################
# #
# Written by: Hans Jacob Teglbjaerg Stephensen #
# Date: Nov 10th 2020 #
# #
# Please reference on use: https://arxiv.org/abs/2008.03927 #
# MIT license #
# #
###############################################################################
import numpy as np
from .TriMeshIntegrator import TriMeshIntegrator
from .TetMeshIntegrator import TetMeshIntegrator
def measure_exterior(V, F, D, r, adaptive=False):
"""
Computes the interior volume and interface area of a mesh in a given
distance map
Parameters
----------
V : (N, 3) array_like
Input vertex array.
F : (M, 3) array_like
Input triangular face array.
D : (N) array_like
Distance for each vertex
r : (K) array_like
Distance steps for which to measure the input elements
Returns
__________
mu10 : (K) array_like
The measure of surface area for each of the given distances
mu11 : (K) array_like
The measure of circumference of interface area for each of the given
distances
Examples
--------
>>> import numpy as np
>>> from ReHaMeasure import measure_interior
>>> np.random.seed(42)
>>> V = np.random.uniform(0, 10, (100,3))
>>> F = np.random.randint(0, 100, (100,3))
>>> D = np.random.uniform(0,10,100)
>>> r = np.linspace(0,15,100)
>>> mu10, mu11 = measure_interior(V, F, D, r)
"""
V = V.astype(float)
D = D.astype(float)
return TriMeshIntegrator(V, F, D, r, adaptive)
def measure_interior(V, T, D, r, adaptive=False):
"""
Computes the interior volume and interface area of a mesh in a given
distance map
Parameters
----------
V : (N, 3) array_like
Input vertex array.
T : (M, 4) array_like
Input tetrahedral element array.
D : (N) array_like
Distance for each vertex
r : (K) array_like
Distance steps for which to measure the input elements
Returns
__________
mu00 : (K) array_like
The measure of volume for each of the given distances
mu01 : (K) array_like
The measure of interior interface area for each of the given distances
Examples
--------
>>> import numpy as np
>>> from ReHaMeasure import measure_interior
>>> np.random.seed(42)
>>> V = np.random.uniform(0, 10, (100,3))
>>> T = np.random.randint(0, 100, (100,4))
>>> D = np.random.uniform(0,10,100)
>>> r = np.linspace(0,15,100)
>>> mu00, mu01 = measure_interior(V, T, D, r)
"""
V = V.astype(float)
D = D.astype(float)
return TetMeshIntegrator(V, T, D, r, adaptive)
###############################################################################
# #
# Written by: Hans Jacob Teglbjaerg Stephensen #
# Date: Nov 10th 2020 #
# #
# Please reference on use: https://arxiv.org/abs/2008.03927 #
# MIT license #
# #
###############################################################################
import numpy as np
numba_available = True
try:
from numba import njit
except ImportError:
numba_available = False
def if_njit(function, **kwargs):
if not numba_available:
return function
else:
return njit()(function)
TOL = 1e-10
@if_njit
def copy_verts(verts, v_indices, ordering):
t_verts = np.empty(12, dtype=verts.dtype)
for j in range(4):
vn = v_indices[ordering[j]]
for i in range(3):
t_verts[3*j + i] = verts[3*vn + i]
return t_verts
@if_njit
def get_sorted_indices(dists):
# make a python version at check implementation
ordering = np.array([0,1,2,3], dtype=np.int32)
for i in range(3):
for j in range(3-i):
if (dists[j] > dists[j+1]):
dists[j], dists[j+1] = dists[j+1], dists[j]
ordering[j], ordering[j+1] = ordering[j+1], ordering[j]
return dists, ordering
@if_njit
def rotx(vy, vz, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vz = vy * st + vz * ct
vy = vy * ct - vz * st
vz = tmp_vz
return vy, vz
@if_njit
def roty(vx, vz, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vz = ct * vz - st * vx
vx = vx * ct + st * vz
vz = tmp_vz
return vx, vz
@if_njit
def rotz(vx, vy, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vy = vx * st + vy * ct
vx = vx * ct - vy * st
vy = tmp_vy
return vx, vy
def get_linear_params(verts, dists):
A = np.empty((4,4))
v = np.empty(4)
for i in range(4):
v[i] = dists[i]
for j in range(4):
if j < 3:
A[i, j] = verts[i*3 + j]
else:
A[i, j] = 1.0
if np.abs(np.linalg.det(A)) < TOL:
# this will be caught later
return 0., 0., 0.
x = np.linalg.solve(A, v)
a, b, c = x[:3]
return a, b, c
@if_njit
def full_tri_area(verts):
# cross product and multiply by last vector
a = np.array([verts[3] - verts[0], verts[4] - verts[1], verts[5] - verts[2]])
b = np.array([verts[6] - verts[0], verts[7] - verts[1], verts[8] - verts[2]])
cx = a[1] * b[2] - b[1] * a[2]
cy = b[0] * a[2] - a[0] * b[2]
cz = a[0] * b[1] - b[0] * a[1]
return np.sqrt(cx**2 + cy**2 + cz**2) / 2.0
@if_njit
def full_tet_volume(verts):
# cross product and multiply by last vector
a = np.array([verts[3] - verts[0], verts[4] - verts[1], verts[5] - verts[2]])
b = np.array([verts[6] - verts[0], verts[7] - verts[1], verts[8] - verts[2]])
c = np.array([verts[9] - verts[0], verts[10] - verts[1], verts[11] - verts[2]])
cm_x = (b[1] * c[2] - c[1] * b[2]) * a[0]
cm_y = (c[0] * b[2] - b[0] * c[2]) * a[1]
cm_z = (b[0] * c[1] - c[0] * b[1]) * a[2]
return np.abs(cm_x + cm_y + cm_z) / 6.0
def transform_tet(verts, dists):
# Normalize so first vertice is at the origin
for j in range(1,4):#(int j = 1; j < 4; ++j)
for i in range(3): #(int i = 0; i < 3; ++i)
verts[j*3 + i] -= verts[i]
for i in range(3): # (int i = 0; i < 3; ++i)
verts[i] = 0
a, b, c = get_linear_params(verts, dists)
#print('Python', a, b, c)
gradient_norm = np.sqrt(a**2 + b**2 + c**2)
scale = gradient_norm
t = -np.arctan2(c, b)
for i in range(1,4):
# rotx (x is unchanged, rot y and z)
verts[i*3+1], verts[i*3+2] = rotx(verts[i*3+1], verts[i*3+2], t)
b, c = rotx(b, c, t)
t = -np.arctan2(b, a)
for i in range(1,4):
# rotz (z is unchanged, rot x and y)
verts[i*3], verts[i*3+1] = rotz(verts[i*3], verts[i*3+1], t)
# Scale to h scale
for i in range(4):
verts[i*3] = verts[i*3] * gradient_norm + dists[0]
return verts, scale
@if_njit
def h_midpoint(p0, p1, h):
p01 = np.empty(3)
K = np.array([p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]])
h0 = p0[0]
h1 = p1[0]
a = h0
b = 1.0/(h1 - h0)
for i in range(3):# (int i = 0; i < 3; ++i)
p01[i] = p0[i] + K[i]*(h - a)*b
return p01
@if_njit
def volume_fast(v0, v1, v2, v3, h):
p01 = h_midpoint(v0, v1, h)
p02 = h_midpoint(v0, v2, h)
p03 = h_midpoint(v0, v3, h)
p12 = h_midpoint(v1, v2, h)
p13 = h_midpoint(v1, v3, h)
P_tet1 = np.array([v0[0], v0[1], v0[2],
p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]])
P_tet2 = np.array([v1[0], v1[1], v1[2],
p01[0], p01[1], p01[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]])
volume_1 = full_tet_volume(P_tet1)
volume_2 = full_tet_volume(P_tet2)
return volume_1 - volume_2
@if_njit
def volume_safe(v0, v1, v2, v3, h):
p02 = h_midpoint(v0, v2, h)
p03 = h_midpoint(v0, v3, h)
p12 = h_midpoint(v1, v2, h)
p13 = h_midpoint(v1, v3, h)
P_tet1 = np.array([v0[0], v0[1], v0[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2]])
P_tet2 = np.array([v0[0], v0[1], v0[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]])
P_tet3 = np.array([v0[0], v0[1], v0[2],
v1[0], v1[1], v1[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]])
volume_1 = full_tet_volume(P_tet1)
volume_2 = full_tet_volume(P_tet2)
volume_3 = full_tet_volume(P_tet3)
return volume_1 + volume_2 + volume_3
@if_njit
def sliced_tet_volume(verts, full_volume, i, h):
v0 = np.array([verts[i*12 + 0], verts[i*12 + 1], verts[i*12 + 2]])
v1 = np.array([verts[i*12 + 3], verts[i*12 + 4], verts[i*12 + 5]])
v2 = np.array([verts[i*12 + 6], verts[i*12 + 7], verts[i*12 + 8]])
v3 = np.array([verts[i*12 + 9], verts[i*12 + 10], verts[i*12 + 11]])
if (h <= v0[0]):
# VOL 0
return 0.0
elif (h <= v1[0]):
# VOL 1
p01 = h_midpoint(v0, v1, h)
p02 = h_midpoint(v0, v2, h)
p03 = h_midpoint(v0, v3, h)
P = np.array([v0[0], v0[1], v0[2],
p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]])
volume = full_tet_volume(P)
return volume
elif (h <= v2[0]):
# VOL 2
if (np.abs(v0[0] - v1[0]) < TOL):
# Degenerate case
volume = volume_safe(v0, v1, v2, v3, h)
return volume
else:
# Non-degenerate case
volume = volume_fast(v0, v1, v2, v3, h)
return volume
elif (h <= v3[0]):
# VOL 3
p03 = h_midpoint(v0, v3, h)
p13 = h_midpoint(v1, v3, h)
p23 = h_midpoint(v2, v3, h)
P_remaining = np.array([v3[0], v3[1], v3[2],
p03[0], p03[1], p03[2],
p13[0], p13[1], p13[2],
p23[0], p23[1], p23[2]])
volume_remaining = full_tet_volume(P_remaining)
if np.isnan(volume_remaining):
print('NAN in volume_remaining on " << {P_remaining} {h}')
if np.isnan(full_volume):
print('NAN in full_volume on {v0} {v1} {v2} {v3} {h}')
return full_volume - volume_remaining
else:
# VOL 4
return full_volume
@if_njit
def sliced_tet_area(verts, i, h):
v0 = np.array([verts[i*12 + 0], verts[i*12 + 1], verts[i*12 + 2]])
v1 = np.array([verts[i*12 + 3], verts[i*12 + 4], verts[i*12 + 5]])
v2 = np.array([verts[i*12 + 6], verts[i*12 + 7], verts[i*12 + 8]])
v3 = np.array([verts[i*12 + 9], verts[i*12 + 10], verts[i*12 + 11]])
if h <= v0[0]:
# AREA 0
return 0.0
elif h <= v1[0]:
# AREA 1
p01 = h_midpoint(v0, v1, h)
p02 = h_midpoint(v0, v2, h)
p03 = h_midpoint(v0, v3, h)
P = np.array([p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]])
return full_tri_area(P)
elif (h <= v2[0]):
# AREA 2
p02 = h_midpoint(v0, v2, h)
p03 = h_midpoint(v0, v3, h)
p12 = h_midpoint(v1, v2, h)
p13 = h_midpoint(v1, v3, h)
P_1 = np.array([p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2]])
P_2 = np.array([p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]])
return full_tri_area(P_1) + full_tri_area(P_2)
elif h <= v3[0]:
# AREA 3
p03 = h_midpoint(v0, v3, h)
p13 = h_midpoint(v1, v3, h)
p23 = h_midpoint(v2, v3, h)
P = np.array([p03[0], p03[1], p03[2],
p13[0], p13[1], p13[2],
p23[0], p23[1], p23[2]])
return full_tri_area(P)
else:
# AREA 4
return 0.0
def TetMeshIntegrator(verts, elements, dists, x, adaptive):
num_verts = len(verts)
num_elements = len(elements)
num_dists = len(dists)
t_verts = np.empty(num_elements * 12, dtype=float)
t_scales = np.empty(num_elements, dtype=float)
full_volumes = np.empty(num_elements, dtype=float)
tet_h_bounds = np.empty(num_elements * 2, dtype=float)
verts = verts.flatten()
elements = elements.flatten()
#pragma omp parallel for
for i in range(num_elements):
t_verts_tmp = np.empty(12, dtype=float)
v_indices = np.empty(4, dtype=np.int32)
t_dists_tmp = np.empty(4, dtype=float)
for j in range(4):
v_indices[j] = elements[i*4+j]
t_dists_tmp[j] = dists[v_indices[j]]
t_dists_tmp, ordering = get_sorted_indices(t_dists_tmp)
t_verts_tmp = copy_verts(verts, v_indices, ordering)
tet_h_bounds[i*2] = t_dists_tmp[0]
tet_h_bounds[i*2+1] = t_dists_tmp[3]
t_verts_tmp, t_scale = transform_tet(t_verts_tmp, t_dists_tmp)
#print(t_verts_tmp)
t_scales[i] = t_scale
full_volumes[i] = full_tet_volume(t_verts_tmp)
for j in range(12):
t_verts[i*12 + j] = t_verts_tmp[j]
if adaptive:
x = np.sort(np.unique(np.concatenate([x, tet_h_bounds+TOL, tet_h_bounds-TOL])))
num_x = len(x)
y_volume = np.empty(num_x, dtype=float)
y_area = np.empty(num_x, dtype=float)
#pragma omp parallel for
for j in range(num_x):
y_volume[j] = 0.0
y_area[j] = 0.0
for i in range(num_elements):
if t_scales[i] != 0:
#print(t_scales[i])
y_volume[j] += sliced_tet_volume(t_verts, full_volumes[i], i, x[j]) / t_scales[i]
y_area[j] += sliced_tet_area(t_verts, i, x[j])
return x, y_volume, y_area
###############################################################################
# #
# Written by: Hans Jacob Teglbjaerg Stephensen #
# Date: Nov 10th 2020 #
# #
# Please reference on use: https://arxiv.org/abs/2008.03927 #
# MIT license #
# #
###############################################################################
import numpy as np
numba_available = True
try:
from numba import jit, njit
except ImportError:
numba_available = False
def if_njit(function, **kwargs):
if not numba_available:
return function
else:
return njit()(function)
def if_jit(function, **kwargs):
if not numba_available:
return function
else:
return jit()(function)
TOL = 1e-10
@if_njit
def copy_verts(verts, v_indices, ordering):
t_verts = np.empty(9, dtype=verts.dtype)
for j in range(3):
vn = v_indices[ordering[j]]
for i in range(3):
t_verts[3*j + i] = verts[3*vn + i]
return t_verts
@if_njit
def get_sorted_indices(dists):
ordering = np.array([0,1,2], dtype=np.int32)
for i in range(3):
for j in range(2-i):
if (dists[j] > dists[j+1]):
dists[j], dists[j+1] = dists[j+1], dists[j]
ordering[j], ordering[j+1] = ordering[j+1], ordering[j]
return dists, ordering
@if_njit
def rotx(vy, vz, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vz = vy * st + vz * ct
vy = vy * ct - vz * st
vz = tmp_vz
return vy, vz
@if_njit
def roty(vx, vz, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vz = ct * vz - st * vx
vx = vx * ct + st * vz
vz = tmp_vz
return vx, vz
@if_njit
def rotz(vx, vy, t):
ct = np.cos(t)
st = np.sin(t)
tmp_vy = vx * st + vy * ct
vx = vx * ct - vy * st
vy = tmp_vy
return vx, vy
def get_linear_params(verts, dists):
A = np.empty((3,3))
v = np.empty(3)
for i in range(3):
v[i] = dists[i]
for j in range(3):
if j < 2:
A[i, j] = verts[i*3 + j]
else:
A[i, j] = 1.0
if np.abs(np.linalg.det(A)) < TOL:
# this will be caught later
return 0., 0.
x = np.linalg.solve(A, v)
a, b = x[:2]
return a, b
@if_njit
def get_edge_linear_params(verts, i):
#x0, y0, x1, y1, x2, y2 = verts[i*6:(i+1)*6] # should work as well
linear_params = np.zeros(6, dtype=np.float64)
x0 = verts[i*6]
y0 = verts[i*6 + 1]
x1 = verts[i*6 + 2]
y1 = verts[i*6 + 3]
x2 = verts[i*6 + 4]
y2 = verts[i*6 + 5]
if np.abs(x1 - x0) < TOL:
linear_params[0] = np.inf
else:
linear_params[0] = (y1 - y0) / (x1 - x0)
linear_params[1] = y0 - linear_params[0]*x0
if np.abs(x2 - x0) < TOL:
linear_params[2] = np.inf
else:
linear_params[2] = (y2 - y0) / (x2 - x0)
linear_params[3] = y0 - linear_params[2]*x0
if np.abs(x2 - x1) < TOL:
linear_params[4] = np.inf
else:
linear_params[4] = (y2 - y1) / (x2 - x1)
linear_params[5] = y1 - linear_params[4]*x1
return linear_params
def transform_triangle(verts, dists):
# Normalize such that the first vertex is at the origin
for i in range(3):
verts[3+i] -= verts[i]
verts[6+i] -= verts[i]
verts[i] = 0
# Rotate so it corresponds to a flat triangWle
t = -np.arctan2(verts[4], verts[3])
verts[3], verts[4] = rotz(verts[3],verts[4], t)
verts[6], verts[7] = rotz(verts[6],verts[7], t)
t = np.arctan2(verts[5], verts[3])
verts[3], verts[5] = roty(verts[3],verts[5], t)
verts[6], verts[8] = roty(verts[6],verts[8], t)
t = -np.arctan2(verts[8], verts[7])
verts[4], verts[5] = rotx(verts[4],verts[5], t)
verts[7], verts[8] = rotx(verts[7],verts[8], t)
# calculate full area of triangle
ax = verts[0] - verts[3]
ay = verts[1] - verts[4]
bx = verts[0] - verts[6]
by = verts[1] - verts[7]
full_area = np.abs(ax*by - bx*ay) / 2.0
if (full_area < TOL):
verts[0] = dists[0]
verts[1] = dists[0]
verts[2] = dists[0] # dont think we need this now we have the bounds elsewhere (Update: well, apparently we do...)
return verts, 1, full_area
if (dists[0] == dists[2]):
verts[0] = dists[0]
verts[1] = dists[0]
verts[2] = dists[0] # dont think we need this now we have the bounds elsewhere (Update: well, apparently we do...)
return verts, 1, full_area
a, b = get_linear_params(verts, dists)
gradient_norm = np.sqrt(a*a + b*b)
scale = 1.0/gradient_norm
# Change basis
t = -np.arctan2(b, a)
verts[3],verts[4] = rotz(verts[3],verts[4], t)
verts[6],verts[7] = rotz(verts[6],verts[7], t)
# Scale to h scale
verts[0] = verts[0] + dists[0]
verts[3] = verts[3] * gradient_norm + dists[0]
verts[6] = verts[6] * gradient_norm + dists[0]
return verts, scale, full_area
def sliced_triangle_area(verts, edge_linear_params, i, h):
#x0, y0, x1, y1, x2, y2 = verts[i*6:(i+1)*6] # should work as well
x0 = verts[i*6 + 0]
y0 = verts[i*6 + 1]
x1 = verts[i*6 + 2]
y1 = verts[i*6 + 3]
x2 = verts[i*6 + 4]
y2 = verts[i*6 + 5]
if (x1 == x2) and (y1 == y2):
return 0.0
if x0 == x2:
return 0.0
if x0 == x1:
if y0 == y1:
return 0.0
if x0 == x2:
return 0.0
f02 = lambda x: edge_linear_params[i*6 + 2]*x + edge_linear_params[i*6 + 3]
f12 = lambda x: edge_linear_params[i*6 + 4]*x + edge_linear_params[i*6 + 5]
height_full = x2 - x0
base_full = y1 - y0
full_area = np.abs(0.5*height_full*base_full)
right_height = x2 - h
right_base = f02(h) - f12(h)
right_area = np.abs(0.5*right_height*right_base)
area = full_area - right_area
return area
#print('linear params', edge_linear_params[i*6 + 0], edge_linear_params[i*6 + 1])
if np.abs(edge_linear_params[i*6 + 0]) == np.inf or \
np.abs(edge_linear_params[i*6 + 1]) == np.inf:
f01 = lambda x: 0
else:
f01 = lambda x: edge_linear_params[i*6 + 0]*x + edge_linear_params[i*6 + 1]
f02 = lambda x: edge_linear_params[i*6 + 2]*x + edge_linear_params[i*6 + 3]
f12 = lambda x: edge_linear_params[i*6 + 4]*x + edge_linear_params[i*6 + 5]
if h < x1:
height = h - x0
base = f01(h) - f02(h)
area = np.abs(0.5*base*height)
return np.abs(0.5*base*height)
left_height = x1 - x0
base = f01(x1) - f02(x1)
left_area = np.abs(base*left_height)
right_height = x2 - x1
mid_base = y1 - f02(x1)
full_right_area = np.abs(right_height*mid_base)
remaining_right_area = np.abs((f12(h) - f02(h)) * (x2 - h))
area = 0.5*(left_area + full_right_area - remaining_right_area)
#print('area', area)
return area
@if_njit
def sliced_triangle_circ(verts, edge_linear_params, i, h):
#x0, y0, x1, y1, x2, y2 = verts[i*6:(i+1)*6] # should work as well
x0 = verts[i*6 + 0]
y0 = verts[i*6 + 1]
x1 = verts[i*6 + 2]
y1 = verts[i*6 + 3]
x2 = verts[i*6 + 4]
y2 = verts[i*6 + 5]
if np.abs(x0 - x1) < TOL and np.abs(y0 - y1) < TOL:
# degenerate flat triangle at left side
return 0.0
if np.abs(x0 - x2) < TOL and np.abs(y0 - y2) < TOL:
# degenerate point triangle
return 0.0
if np.abs(x1 - x2) < TOL and np.abs(y1 - y2) < TOL:
# degenerate flat triangle at right side
return 0.0
if np.abs(x0 - x1) < TOL and np.abs(x1 - x2) < TOL:
# degenerate vertical triangle, the measure is not well defined here,
# but will be 0.0 in the limit from both sides for this triangle
return 0.0
if np.abs(y0 - y1) < TOL and np.abs(y1 - y2) < TOL:
# degenerate case where we have a horizontal flat triangle
return 0.0
f02 = lambda x: edge_linear_params[i*6 + 2]*x + edge_linear_params[i*6 + 3]
if (h < x1):
f01 = lambda x: edge_linear_params[i*6 + 0]*x + edge_linear_params[i*6 + 1]
return np.abs(f01(h) - f02(h))
f12 = lambda x: edge_linear_params[i*6 + 4]*x + edge_linear_params[i*6 + 5]
return np.abs(f12(h) - f02(h))
def TriMeshIntegrator(verts, faces, dists, x, adaptive):
num_verts = len(verts)
num_faces = len(faces)
num_dists = len(dists)
t_verts = np.empty(num_faces * 6, dtype=float)
t_scales = np.empty(num_faces, dtype=float)
full_areas = np.empty(num_faces, dtype=float)
triangle_h_bounds = np.empty(num_faces * 2, dtype=float)
edge_linear_params = np.empty(num_faces * 6, dtype=float)
verts = verts.flatten()
faces = faces.flatten()
#pragma omp parallel for
for i in range(num_faces):
v_indices = np.array([faces[i*3 + 0], faces[i*3 + 1], faces[i*3 + 2]])
t_dists_tmp = np.array([dists[v_indices[0]], dists[v_indices[1]], dists[v_indices[2]]])
t_dists_tmp, ordering = get_sorted_indices(t_dists_tmp)
t_verts_tmp = copy_verts(verts, v_indices, ordering)
triangle_h_bounds[i*2] = t_dists_tmp[0]
triangle_h_bounds[i*2+1] = t_dists_tmp[2]
t_verts_tmp, t_scale, full_area = transform_triangle(t_verts_tmp, t_dists_tmp)
for j in range(6):
g_idx = (j // 2) * 3 + j%2 # This skips the z-coordinate which is now 0
t_verts[i*6 + j] = t_verts_tmp[g_idx]
t_scales[i] = t_scale
full_areas[i] = full_area
edge_linear_params[i*6:(i+1)*6] = get_edge_linear_params(t_verts, i)
if adaptive:
x = np.sort(np.unique(np.concatenate([x, triangle_h_bounds+TOL, triangle_h_bounds-TOL])))
num_x = len(x)
y_area = np.empty(num_x, dtype=float)
y_circumference = np.empty(num_x, dtype=float)
#pragma omp parallel for
for j in range(num_x):
y_area[j] = 0.0
y_circumference[j] = 0.0
for i in range(num_faces):
if x[j] <= triangle_h_bounds[i*2]:
continue
elif x[j] >= triangle_h_bounds[i*2+1]:
y_area[j] += full_areas[i]
else:
y_area[j] += sliced_triangle_area(t_verts, edge_linear_params, i, x[j]) * t_scales[i]
y_circumference[j] += sliced_triangle_circ(t_verts, edge_linear_params, i, x[j])
return x, y_area, y_circumference
from .RelationalShape import measure_interior, measure_exterior
\ No newline at end of file
File added
File added
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment