Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
59 changes: 37 additions & 22 deletions .github/workflows/python-package.yml
Original file line number Diff line number Diff line change
Expand Up @@ -5,36 +5,51 @@ name: Python package

on:
push:
branches: [ "master" ]
branches: ["master"]
pull_request:
branches: [ "master" ]
branches: ["master"]

jobs:
build:

runs-on: ubuntu-latest
strategy:
fail-fast: false
matrix:
python-version: ["3.8", "3.9", "3.10"]

steps:
- uses: actions/checkout@v3
- name: Set up Python ${{ matrix.python-version }}
uses: actions/setup-python@v3
with:
python-version: ${{ matrix.python-version }}
- name: Install dependencies
run: |
python -m pip install --upgrade pip
python -m pip install flake8 pytest
if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
- name: Lint with flake8
run: |
# stop the build if there are Python syntax errors or undefined names
flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
# exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
- name: Test with pytest
run: |
pytest
- uses: actions/checkout@v3
- name: Set up Python ${{ matrix.python-version }}
uses: actions/setup-python@v3
with:
python-version: ${{ matrix.python-version }}
- name: Install Linux GL Dependencies
if: runner.os == 'Linux'
shell: bash
run: sudo apt-get update && sudo apt-get install libgl1-mesa-glx xvfb -y
- name: Configure headless display on Linux
if: runner.os == 'Linux'
shell: bash
run: |
export DISPLAY=:99.0
echo "DISPLAY=:99.0" >> $GITHUB_ENV
Xvfb :99 -screen 0 1024x768x24 > /dev/null 2>&1 &
- name: Give xvfb some time to start on Linux
if: runner.os == 'Linux'
shell: bash
run: sleep 3
- name: Install dependencies
run: |
python -m pip install --upgrade pip
python -m pip install flake8 pytest
if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
- name: Lint with flake8
run: |
# stop the build if there are Python syntax errors or undefined names
flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
# exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
- name: Test with pytest
run: |
pip install -e .
pytest
113 changes: 62 additions & 51 deletions meshparty/ray_tracing.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,19 @@
import trimesh.ray
from trimesh.ray import ray_pyembree
from trimesh.ray import ray_pyembree
try:
from trimesh.ray import ray_pyembree
except ImportError:
print('You must install pyembree in order to utilize SDF functionality')
from scipy.linalg import block_diag
import numpy as np
import multiwrapper.multiprocessing_utils as mu
from meshparty import trimesh_io
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 +37,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 +54,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 +75,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 +120,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 +148,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 +170,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 +187,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 +212,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 +225,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