Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix typo #86

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fix typo #86

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 58 additions & 49 deletions meshparty/ray_tracing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,10 @@
import logging


def ray_trace_distance(vertex_inds, mesh, max_iter=10, rand_jitter=0.001, verbose=False, ray_inter=None):
'''
def ray_trace_distance(
vertex_inds, mesh, max_iter=10, rand_jitter=0.001, verbose=False, ray_inter=None
):
"""
Compute distance to opposite side of the mesh for specified vertex indices on the mesh.

Parameters
Expand All @@ -33,10 +35,11 @@ def ray_trace_distance(vertex_inds, mesh, max_iter=10, rand_jitter=0.001, verbos
np.array
rs, a K long array of sdf values. rays with no result after max_iters will contain zeros.

'''
"""
if not trimesh.ray.has_embree:
logging.warning(
"calculating rays without pyembree, conda install pyembree for large speedup")
"calculating rays without pyembree, conda install pyembree for large speedup"
)

if ray_inter is None:
ray_inter = ray_pyembree.RayMeshIntersector(mesh)
Expand All @@ -49,19 +52,19 @@ def ray_trace_distance(vertex_inds, mesh, max_iter=10, rand_jitter=0.001, verbos
if verbose:
print(np.sum(~good_rs))
blank_inds = np.where(~good_rs)[0]
starts = (mesh.vertices -
mesh.vertex_normals)[vertex_inds, :][~good_rs, :]
vs = -mesh.vertex_normals[vertex_inds, :] \
+ (1.2**it)*rand_jitter*np.random.rand(*
mesh.vertex_normals[vertex_inds, :].shape)
starts = (mesh.vertices - mesh.vertex_normals)[vertex_inds, :][~good_rs, :]
vs = -mesh.vertex_normals[vertex_inds, :] + (
1.2**it
) * rand_jitter * np.random.rand(*mesh.vertex_normals[vertex_inds, :].shape)
vs = vs[~good_rs, :]

rtrace = ray_inter.intersects_location(starts, vs, multiple_hits=False)

if len(rtrace[0] > 0):
# radius values
rs[blank_inds[rtrace[1]]] = np.linalg.norm(
mesh.vertices[vertex_inds, :][rtrace[1]]-rtrace[0], axis=1)
mesh.vertices[vertex_inds, :][rtrace[1]] - rtrace[0], axis=1
)
good_rs[blank_inds[rtrace[1]]] = True
it += 1
if it > max_iter:
Expand All @@ -70,44 +73,43 @@ def ray_trace_distance(vertex_inds, mesh, max_iter=10, rand_jitter=0.001, verbos


def vogel_disk_sampler(num_points, radius=1):
"""Uniform sampler of the unit disk.
"""
"""Uniform sampler of the unit disk."""
golden_angle = np.pi * (3 - np.sqrt(5))
thetas = golden_angle * np.arange(num_points)
radii = radius * np.sqrt(np.arange(num_points)/num_points)
radii = radius * np.sqrt(np.arange(num_points) / num_points)

xs = radii * np.cos(thetas)
ys = radii * np.sin(thetas)

return xs, ys


def unit_vector_sampler(num_points, widest_angle=np.pi/3):
def unit_vector_sampler(num_points, widest_angle=np.pi / 3):
if np.abs(widest_angle) > np.pi:
print('')
print("")
return None
xs, ys = vogel_disk_sampler(num_points, radius=1)
zs = 1/np.tan(widest_angle) * np.ones(num_points)
zs = 1 / np.tan(widest_angle) * np.ones(num_points)
Vs = np.vstack((xs, ys, zs)).T
return Vs / np.linalg.norm(Vs, axis=1)[:, np.newaxis]


def Rx(ang):
return np.array([[1, 0, 0],
[0, np.cos(ang), -np.sin(ang)],
[0, np.sin(ang), np.cos(ang)]])
return np.array(
[[1, 0, 0], [0, np.cos(ang), -np.sin(ang)], [0, np.sin(ang), np.cos(ang)]]
)


def Ry(ang):
return np.array([[np.cos(ang), 0, np.sin(ang)],
[0, 1, 0],
[-np.sin(ang), 0, np.cos(ang)]])
return np.array(
[[np.cos(ang), 0, np.sin(ang)], [0, 1, 0], [-np.sin(ang), 0, np.cos(ang)]]
)


def Rz(ang):
return np.array([[np.cos(ang), -np.sin(ang), 0],
[np.sin(ang), np.cos(ang), 0],
[0, 0, 1]])
return np.array(
[[np.cos(ang), -np.sin(ang), 0], [np.sin(ang), np.cos(ang), 0], [0, 0, 1]]
)


def _rotated_cone(data):
Expand All @@ -116,12 +118,12 @@ def _rotated_cone(data):
return np.dot(Rtrans, vs_raw.T).T


def oriented_vector_cones(center_vectors, num_points, widest_angle=np.pi/3, normalize=False):
"""Produces all ray cones
"""
def oriented_vector_cones(
center_vectors, num_points, widest_angle=np.pi / 3, normalize=False
):
"""Produces all ray cones"""
if normalize:
cv_norm = center_vectors / \
np.linalg.norm(center_vector, axis=1)[:, np.newaxis]
cv_norm = center_vectors / np.linalg.norm(center_vectors, axis=1)[:, np.newaxis]
else:
cv_norm = center_vectors

Expand All @@ -144,19 +146,19 @@ def _multi_angle_weighted_distance(data):


def angle_weighted_distance(ds, angles, weights):
"""Does angle-weighted distance averaging. Ignores outliers and emphasizes normal hits.
"""
"""Does angle-weighted distance averaging. Ignores outliers and emphasizes normal hits."""
if len(ds) == 0 or np.nansum(weights) == 0:
return np.nan

med_angle = np.median(angles)
std_angle = np.std(angles)
min_angle = med_angle-std_angle
max_angle = max(med_angle+std_angle, np.pi / 2)
min_angle = med_angle - std_angle
max_angle = max(med_angle + std_angle, np.pi / 2)
good_rows = np.logical_and(angles >= min_angle, angles <= max_angle)

nanaverage = np.nansum(
ds[good_rows] * weights[good_rows]) / np.nansum(weights[good_rows])
nanaverage = np.nansum(ds[good_rows] * weights[good_rows]) / np.nansum(
weights[good_rows]
)
return nanaverage


Expand All @@ -166,8 +168,8 @@ def all_angle_weighted_distances(ds, angles, weights, rep_inds, inds):
real_inds, slice_bnds = np.unique(rep_inds, return_index=True)

ind_map = []
for ii in range(len(real_inds)-1):
row = slice(slice_bnds[ii], slice_bnds[ii+1])
for ii in range(len(real_inds) - 1):
row = slice(slice_bnds[ii], slice_bnds[ii + 1])
data.append((ds[row], angles[row], weights[row]))
ind_map.append(real_inds[ii])
row = slice(slice_bnds[-1], len(ds))
Expand All @@ -183,17 +185,19 @@ def all_angle_weighted_distances(ds, angles, weights, rep_inds, inds):


def _compute_ray_vectors(mesh, mesh_inds, num_points, cone_angle):
return np.vstack(oriented_vector_cones(-mesh.vertex_normals[mesh_inds], num_points, cone_angle))
return np.vstack(
oriented_vector_cones(-mesh.vertex_normals[mesh_inds], num_points, cone_angle)
)


def shape_diameter_function(mesh_inds, mesh, num_points=30, cone_angle=np.pi/3):
def shape_diameter_function(mesh_inds, mesh, num_points=30, cone_angle=np.pi / 3):
"""Computes shape diameter function by sending a cone of rays from each specified vertex point
and doing a weighted average of where they hit the opposite side of the mesh.

Parameters
----------
mesh : trimesh.Mesh
Mesh
Mesh
mesh_inds : list of indices
Vertex indices at which to compute SDF
num_points : int, optional
Expand All @@ -206,9 +210,10 @@ def shape_diameter_function(mesh_inds, mesh, num_points=30, cone_angle=np.pi/3):

Description
"""
start = (mesh.vertices-mesh.vertex_normals)[mesh_inds, :]
rep_inds = np.concatenate([ii*np.ones(num_points, dtype=int)
for ii in range(start.shape[0])])
start = (mesh.vertices - mesh.vertex_normals)[mesh_inds, :]
rep_inds = np.concatenate(
[ii * np.ones(num_points, dtype=int) for ii in range(start.shape[0])]
)
starts = start[rep_inds]

vs = _compute_ray_vectors(mesh, mesh_inds, num_points, cone_angle)
Expand All @@ -218,12 +223,16 @@ def shape_diameter_function(mesh_inds, mesh, num_points=30, cone_angle=np.pi/3):

hit_rows = rtrace[1]
ds = np.linalg.norm(rtrace[0] - starts[hit_rows], axis=1)
angles = np.arccos(
np.sum(mesh.face_normals[rtrace[2]] * vs[hit_rows], axis=1))
with np.errstate(divide='ignore', invalid='ignore'):
weights = 1/angles
angles = np.arccos(np.sum(mesh.face_normals[rtrace[2]] * vs[hit_rows], axis=1))
with np.errstate(divide="ignore", invalid="ignore"):
weights = 1 / angles
good_rows = np.isfinite(weights)

rs = all_angle_weighted_distances(
ds[good_rows], angles[good_rows], weights[good_rows], rep_inds[hit_rows[good_rows]], mesh_inds)
ds[good_rows],
angles[good_rows],
weights[good_rows],
rep_inds[hit_rows[good_rows]],
mesh_inds,
)
return rs