Skip to content
Snippets Groups Projects
Commit 05c31d0e authored by Hans Stephensen's avatar Hans Stephensen
Browse files

fixed and issue which would result high negative area in triangle calculations...

fixed and issue which would result high negative area in triangle calculations when triangles were degenerate.
parent fd359efb
Branches 3d_watershed
No related tags found
No related merge requests found
import numpy as np import numpy as np
import overlap import overlap
import logging
from polygon_triangulate import polygon_triangulate, PlotTriangulation from polygon_triangulate import polygon_triangulate, PlotTriangulation
from tetrahedralize_triangles import tetrahedralize from tetrahedralize_triangles import tetrahedralize
logging.basicConfig(filename='warnings.log', level=logging.DEBUG)
def measure_volume(V_tet, T, origin, r): def measure_volume(V_tet, T, origin, r):
sphere = overlap.Sphere(origin, r) sphere = overlap.Sphere(origin, r)
volume = 0.0 volume = 0.0
volumes = []
for tet in T: for tet in T:
vertices = V_tet[tet] vertices = V_tet[tet]
overlap_tet = overlap.Tetrahedron(vertices) overlap_tet = overlap.Tetrahedron(vertices)
overlap_volume = overlap.overlap(sphere, overlap_tet) overlap_volume = overlap.overlap(sphere, overlap_tet)
volumes.append(overlap_volume)
if overlap_volume < 0:
print('Encountered negative volume overlap')
logging.warning('Encountered negative volume overlap')
volume += overlap_volume volume += overlap_volume
#import matplotlib.pyplot as plt
#plt.title(f'radius {r} total {volume}')
#plt.plot(np.arange(len(volumes)), volumes)
#plt.show()
return volume return volume
def same_point(x, y): def same_point(x, y):
...@@ -28,6 +42,16 @@ def remove_duplicate_contour_points(P): ...@@ -28,6 +42,16 @@ def remove_duplicate_contour_points(P):
return remove_duplicate_contour_points(np.delete(P, i, axis=0)) return remove_duplicate_contour_points(np.delete(P, i, axis=0))
return P return P
def remove_line_triangles(P, F):
n = len(F)
for i in range(n-1, -1, -1):
PF = P[F[i]]
v1 = PF[0] - PF[1]
v2 = PF[1] - PF[2]
if np.isclose(np.abs(np.dot(v1/np.linalg.norm(v1),v2/np.linalg.norm(v2))),1):
F = np.delete(F, i, axis=0)
return F
def measure_cumulative_volume(z, P, thickness, origin, R): def measure_cumulative_volume(z, P, thickness, origin, R):
''' '''
Integration of the intersecting volume of a thick polygon and a sphere of increasing radius. Integration of the intersecting volume of a thick polygon and a sphere of increasing radius.
...@@ -60,6 +84,13 @@ def measure_cumulative_volume(z, P, thickness, origin, R): ...@@ -60,6 +84,13 @@ def measure_cumulative_volume(z, P, thickness, origin, R):
z1 = z + thickness / 2 z1 = z + thickness / 2
_, V_tri, F = polygon_triangulate(P) _, V_tri, F = polygon_triangulate(P)
F = remove_line_triangles(V_tri, F)
#import matplotlib.pyplot as plt
#print(V_tri[:5,:])
#PlotTriangulation(V_tri, F)
#plt.show()
V_tet, T = tetrahedralize(V_tri, F, z0, z1) V_tet, T = tetrahedralize(V_tri, F, z0, z1)
mu = np.zeros(len(R)) mu = np.zeros(len(R))
...@@ -123,4 +154,4 @@ def measure_cumulative_points(z, X, thickness, origin, R): ...@@ -123,4 +154,4 @@ def measure_cumulative_points(z, X, thickness, origin, R):
s = h - d0 s = h - d0
mu[i] = np.sum(np.clip(s/thickness, 0, 1)) mu[i] = np.sum(np.clip(s/thickness, 0, 1))
return mu return mu
\ No newline at end of file
...@@ -78,17 +78,16 @@ def write_csv(path, R, dP, dV): ...@@ -78,17 +78,16 @@ def write_csv(path, R, dP, dV):
for i0,i1,a,b,c,d,e,v in zip(R[:-1], R[1:], dP, dV, Dens, intdP, intdV, intDens): for i0,i1,a,b,c,d,e,v in zip(R[:-1], R[1:], dP, dV, Dens, intdP, intdV, intDens):
f.write(f'{i0},{a},{b},{c},,{i0},{i1},{d},{e},{v}\n') f.write(f'{i0},{a},{b},{c},,{i0},{i1},{d},{e},{v}\n')
def run(base_dir): def run(base_dir, output_dir):
R = np.arange(0, 3001, 125) R = np.arange(0, 3001, 125)
output_dir = 'output'
os.makedirs(output_dir, exist_ok=True) os.makedirs(output_dir, exist_ok=True)
for experiment in tqdm(os.listdir(base_dir), desc='experiment'): for experiment in tqdm(os.listdir(base_dir), desc='experiment'):
for reading in os.listdir(os.path.join(base_dir, experiment)): for reading in os.listdir(os.path.join(base_dir, experiment)):
#if reading != '7566':
# continue
data_dir = os.path.join(base_dir, experiment, reading) data_dir = os.path.join(base_dir, experiment, reading)
origin, thickness = parse_params(data_dir) origin, thickness = parse_params(data_dir)
...@@ -116,4 +115,5 @@ def run(base_dir): ...@@ -116,4 +115,5 @@ def run(base_dir):
if __name__ == '__main__': if __name__ == '__main__':
base_dir = 'data' base_dir = 'data'
run(base_dir) output_dir = 'output'
\ No newline at end of file run(base_dir, output_dir)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment