Skip to content

Commit

Permalink
add pointcloud.relax_point_positions() auto_distance() and fixes to c…
Browse files Browse the repository at this point in the history
…ollapse_edges() and clean()
  • Loading branch information
marcomusy committed Jan 8, 2024
1 parent a84da7d commit 19a18fa
Show file tree
Hide file tree
Showing 7 changed files with 170 additions and 23 deletions.
5 changes: 5 additions & 0 deletions docs/changes.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,10 @@
- add `tetmesh.generate_random_points()` to generate random points in a tet mesh
- rename `integrate_arrays_over_domain()` to `integrate_data`
- extend `volume.operation()` to support logic operations as per #1002
- add `pointcloud.relax_point_positions()` method
- add `pointcloud.auto_distance()` method calculates the distance to the closest
point in the same cloud of points.
- fixed `mesh.collapse_edges()` after #992


## Breaking changes
Expand All @@ -33,6 +37,7 @@
- change `clone2d(scale=...)` to `clone2d(size=...)`
- remove `shapes.StreamLines()` becoming `object.compute_streamlines()`
- split `mesh.decimate()` into `mesh.decimate()`, `mesh.decimate_pro()` and `mesh.decimate_binned()` as per #992
- modified `core.clean()` after #992


### Bug Fixes
Expand Down
8 changes: 4 additions & 4 deletions examples/pyplot/embed_matplotlib.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""Include background images in the rendering scene
(e.g. generated by matplotlib)"""
(generated by matplotlib)"""
import matplotlib.pyplot as plt
from vedo import *

Expand All @@ -11,7 +11,7 @@
plt.hist(msh.celldata["chem_0"], log=True)
plt.title(r"$\mathrm{Matplotlib\ Histogram\ of\ log(chem_0)}$")

pic1 = Image(fig).clone2d("bottom-right", 0.5).alpha(0.8)
pic2 = Image(dataurl + "images/embryo.jpg").clone2d("top-right")
# pic1 = Image(fig).clone2d("top-right", 0.5).alpha(0.8)
pic2 = Image(dataurl + "images/embryo.jpg").clone2d("bottom-right")

show(msh, pic1, pic2, __doc__, bg="lightgrey", axes=1)
show(msh, fig, pic2, __doc__, bg="lightgrey", zoom=1.2, axes=1)
37 changes: 24 additions & 13 deletions vedo/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -1235,14 +1235,19 @@ def delete_cells(self, ids):
return self

def collapse_edges(self, distance, iterations=1):
"""Collapse mesh edges so that are all above distance."""
self.clean()
x0, x1, y0, y1, z0, z1 = self.bounds()
fs = min(x1 - x0, y1 - y0, z1 - z0) / 10
d2 = distance * distance
if distance > fs:
vedo.logger.error(f"distance parameter is too large, should be < {fs}, skip!")
return self
"""
Collapse mesh edges so that are all above `distance`.
Example:
```python
from vedo import *
np.random.seed(2)
grid1 = Grid().add_gaussian_noise(0.8).triangulate().lw(1)
grid1.celldata['scalar'] = grid1.cell_centers[:,1]
grid2 = grid1.clone().collapse_edges(0.1)
show(grid1, grid2, N=2, axes=1)
```
"""
for _ in range(iterations):
medges = self.edges
pts = self.vertices
Expand All @@ -1252,16 +1257,22 @@ def collapse_edges(self, distance, iterations=1):
if len(e) == 2:
id0, id1 = e
p0, p1 = pts[id0], pts[id1]
d = mag2(p1 - p0)
if d < d2 and id0 not in moved and id1 not in moved:
if (np.linalg.norm(p1-p0) < distance
and id0 not in moved
and id1 not in moved
):
p = (p0 + p1) / 2
newpts[id0] = p
newpts[id1] = p
moved += [id0, id1]

self.vertices = newpts
self.clean()
self.compute_normals()
cpd = vtk.new("CleanPolyData")
cpd.ConvertLinesToPointsOff()
cpd.ConvertPolysToLinesOff()
cpd.ConvertStripsToPolysOff()
cpd.SetInputData(self.dataset)
cpd.Update()
self._update(cpd.GetOutput())

self.pipeline = OperationNode(
"collapse_edges",
Expand Down
3 changes: 3 additions & 0 deletions vedo/plotter.py
Original file line number Diff line number Diff line change
Expand Up @@ -3025,6 +3025,9 @@ def _scan_input_return_acts(self, objs):
elif "TetgenIO" in str(type(a)):
scanned_acts.append(vedo.TetMesh(a).shrink(0.9).c("pink7").actor)

elif "matplotlib.figure.Figure" in str(type(a)):
scanned_acts.append(vedo.Image(a).clone2d("top-right", 0.6))

else:
vedo.logger.error(f"cannot understand input in show(): {type(a)}")

Expand Down
137 changes: 132 additions & 5 deletions vedo/pointcloud.py
Original file line number Diff line number Diff line change
Expand Up @@ -1042,9 +1042,9 @@ def clean(self):
"""Clean pointcloud or mesh by removing coincident points."""
cpd = vtk.new("CleanPolyData")
cpd.PointMergingOn()
cpd.ConvertLinesToPointsOn()
cpd.ConvertPolysToLinesOn()
cpd.ConvertStripsToPolysOn()
cpd.ConvertLinesToPointsOff()
cpd.ConvertPolysToLinesOff()
cpd.ConvertStripsToPolysOff()
cpd.SetInputData(self.dataset)
cpd.Update()
self._update(cpd.GetOutput())
Expand Down Expand Up @@ -1543,6 +1543,30 @@ def closest_point(

return np.array(trgp)

def auto_distance(self):
"""
Calculate the distance to the closest point in the same cloud of points.
The output is stored in a new pointdata array called "AutoDistance",
and it is also returned by the function.
"""
points = self.vertices
if not self.point_locator:
self.point_locator = vtk.new("StaticPointLocator")
self.point_locator.SetDataSet(self.dataset)
self.point_locator.BuildLocator()
qs = []
vtklist = vtk.vtkIdList()
vtkpoints = self.dataset.GetPoints()
for p in points:
self.point_locator.FindClosestNPoints(2, p, vtklist)
q = [0, 0, 0]
pid = vtklist.GetId(1)
vtkpoints.GetPoint(pid, q)
qs.append(q)
dists = np.linalg.norm(points - np.array(qs), axis=1)
self.pointdata["AutoDistance"] = dists
return dists

def hausdorff_distance(self, points):
"""
Compute the Hausdorff distance to the input point set.
Expand Down Expand Up @@ -1650,11 +1674,112 @@ def remove_outliers(self, radius, neighbors=5):
self.pipeline = utils.OperationNode("remove_outliers", parents=[self])
return self

def relax_point_positions(
self,
n=10,
iters=10,
sub_iters=10,
packing_factor=1,
max_step=0,
constraints=(),
):
"""
Smooth mesh or points with a
[Laplacian algorithm](https://vtk.org/doc/nightly/html/classvtkPointSmoothingFilter.html)
variant. This modifies the coordinates of the input points by adjusting their positions
to create a smooth distribution (and thereby form a pleasing packing of the points).
Smoothing is performed by considering the effects of neighboring points on one another
it uses a cubic cutoff function to produce repulsive forces between close points
and attractive forces that are a little further away.
In general, the larger the neighborhood size, the greater the reduction in high frequency
information. The memory and computational requirements of the algorithm may also
significantly increase.
The algorithm incrementally adjusts the point positions through an iterative process.
Basically points are moved due to the influence of neighboring points.
As points move, both the local connectivity and data attributes associated with each point
must be updated. Rather than performing these expensive operations after every iteration,
a number of sub-iterations can be specified. If so, then the neighborhood and attribute
value updates occur only every sub iteration, which can improve performance significantly.
Arguments:
n : (int)
neighborhood size to calculate the Laplacian.
iters : (int)
number of iterations.
sub_iters : (int)
number of sub-iterations, i.e. the number of times the neighborhood and attribute
value updates occur during each iteration.
packing_factor : (float)
adjust convergence speed.
max_step : (float)
Specify the maximum smoothing step size for each smoothing iteration.
This limits the the distance over which a point can move in each iteration.
As in all iterative methods, the stability of the process is sensitive to this parameter.
In general, small step size and large numbers of iterations are more stable than a larger
step size and a smaller numbers of iterations.
constraints : (dict)
dictionary of constraints.
Point constraints are used to prevent points from moving,
or to move only on a plane. This can prevent shrinking or growing point clouds.
If enabled, a local topological analysis is performed to determine whether a point
should be marked as fixed" i.e., never moves, or the point only moves on a plane,
or the point can move freely.
If all points in the neighborhood surrounding a point are in the cone defined by
`fixed_angle`, then the point is classified as fixed.
If all points in the neighborhood surrounding a point are in the cone defined by
`boundary_angle`, then the point is classified as lying on a plane.
Angles are expressed in degrees.
Example:
```py
import numpy as np
from vedo import Points, show
from vedo.pyplot import histogram
vpts1 = Points(np.random.rand(10_000, 3))
dists = vpts1.auto_distance()
h1 = histogram(dists, xlim=(0,0.08)).clone2d()
vpts2 = vpts1.clone().relax_point_positions(n=100, iters=20, sub_iters=10)
dists = vpts2.auto_distance()
h2 = histogram(dists, xlim=(0,0.08)).clone2d()
show([[vpts1, h1], [vpts2, h2]], N=2).close()
```
"""
smooth = vtk.new("PointSmoothingFilter")
smooth.SetInputData(self.dataset)
smooth.SetSmoothingModeToUniform()
smooth.SetNumberOfIterations(iters)
smooth.SetNumberOfSubIterations(sub_iters)
smooth.SetPackingFactor(packing_factor)
if self.point_locator:
smooth.SetLocator(self.point_locator)
if not max_step:
max_step = self.diagonal_size() / 100
smooth.SetMaximumStepSize(max_step)
smooth.SetNeighborhoodSize(n)
if constraints:
fixed_angle = constraints.get("fixed_angle", 45)
boundary_angle = constraints.get("boundary_angle", 110)
smooth.EnableConstraintsOn()
smooth.SetFixedAngle(fixed_angle)
smooth.SetBoundaryAngle(boundary_angle)
smooth.GenerateConstraintScalarsOn()
smooth.GenerateConstraintNormalsOn()
smooth.Update()
self._update(smooth.GetOutput())
self.metadata["PackingRadius"] = smooth.GetPackingRadius()
self.pipeline = utils.OperationNode("relax_point_positions", parents=[self])
return self

def smooth_mls_1d(self, f=0.2, radius=None):
"""
Smooth mesh or points with a `Moving Least Squares` variant.
The point data array "Variances" will contain the residue calculated for each point.
Input mesh's polydata is modified.
Arguments:
f : (float)
Expand Down Expand Up @@ -2048,6 +2173,8 @@ def cut_with_planes(self, origins, normals, invert=False):
each cutting plane goes through this point
normals : (array)
normal of each of the cutting planes
invert : (bool)
if True, cut outside instead of inside
Check out also:
`cut_with_box()`, `cut_with_cylinder()`, `cut_with_sphere()`
Expand All @@ -2063,7 +2190,7 @@ def cut_with_planes(self, origins, normals, invert=False):
planes.SetNormals(utils.numpy2vtk(normals, dtype=float))

clipper = vtk.new("ClipPolyData")
clipper.SetInputData(self.dataset) # must be True
clipper.SetInputData(self.dataset)
clipper.SetInsideOut(invert)
clipper.SetClipFunction(planes)
clipper.GenerateClippedOutputOff()
Expand Down
2 changes: 1 addition & 1 deletion vedo/version.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
_version = '2023.5.0+dev19'
_version = '2024.5.0+dev20'
1 change: 1 addition & 0 deletions vedo/vtkclasses.py
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,7 @@
"vtkRadiusOutlierRemoval",
"vtkShepardKernel",
"vtkSignedDistance",
"vtkPointSmoothingFilter",
"vtkVoronoiKernel",
]: location[name] = "vtkFiltersPoints"

Expand Down

0 comments on commit 19a18fa

Please sign in to comment.