L4_Representations
L4: Mesh and Point Cloud
Hao Su
Machine Learning meets Geometry
Shape Representation:
Origin- and Application-Dependent
• Acquired real-world objects
• Modeling “by hand”
• Procedural modeling
• …
2
Rasterized form
(regular grids)
Geometric form
(irregular)
Point Cloud
Mesh
Implicit Shape
F(x) = 0
Volumetric
Multi-view
Depth Map
Other than parametric representations, we also study
these in this course:
3
Agenda
4
Point CloudMesh
Polygonal Meshes
• Representation
• Storage
• Curvature Computation
Polygonal Meshes
• Piece-wise Linear Surface Representation
6
Triangle Mesh
7
http://graphics.stanford.edu/data/3Dscanrep/stanford-bunny-cebal-ssh.jpg
http://www.stat.washington.edu/wxs/images/BUNMID.gif
Triangle Mesh
8
Plus manifold conditions
V = {v1, v2, …, vn} ⊂ ℝ
3
E = {e1, e2, …, ek} ⊆ V × V
F = {f1, f2, …, fm} ⊆ V × V × V
Bad Surfaces
9
http://igl.ethz.ch/projects/parameterization/rangemap-param/rangemap-param.pdf
Nonmanifold Edge
10
http://blog.mixamo.com/wp-content/uploads/2011/01/nonmanifold.jpg
Manifold Mesh
11
http://www.cs.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf
1.Each edge is incident to one or two faces
2.Faces incident to a vertex form a closed or open fan
This is not a fan:
Manifold Mesh
12
http://www.cs.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf
Assume meshes are manifo
ld
(for now)
1.Each edge is incident to one or two faces
2.Faces incident to a vertex form a closed or open fan
13
http://www.pmp-book.org/download/slides/Representations.pdf
“Triangle soup”
Bad Meshes
14
http://www.sciencedirect.com/science/article/pii/S0168874X06000795
Nonuniform
areas and angles
Why is Meshing an Issue?
15
How do you interpret
one value per vertex?
http://www.sciencedirect.com/science/article/pii/S0168874X06000795
Assume Storing Scalar Functions
on Surface
16
http://www.ieeta.pt/polymeco/Screenshots/PolyMeCo_OneView.jpg
Map points to real numbers
Approximation Properties
17
Ex: Taylor’s
Theorem
functions defined at vertices
(e.g., Gaussian curvature)
f :
Techniques to Improve Mesh Quality
• Cleaning
• Repairing
• Remeshing
• …
18
Polygonal Meshes
• Representation
• Storage
• Curvature Computation
Data Structures for Surfaces
• What should be stored?
– Geometry: 3D coordinates
– Topology
– Attributes
‣ Normal, color, texture
coordinates
‣ Per vertex, face, edge
20
Simple Data Structures: Triangle List
• STL format (used in CAD)
• Storage
– Face: 3 positions
• No connectivity information
21
Triangles
0 x0 y0 z0
1 x1 x1 z1
2 x2 y2 z2
3 x3 y3 z3
4 x4 y4 z4
5 x5 y5 z5
6 x6 y6 z6
… … … …
Simple Data Structures: Indexed Face Set
• Used in formats
– OBJ, OFF, WRL
• Storage
– Vertex: position
– Face: vertex indices
– Convention is to save
vertices in counter-
clockwise order for
normal computation
Vertices
v0 x0 y0 z0
v1 x1 x1 z1
v2 x2 y2 z2
v3 x3 y3 z3
v4 x4 y4 z4
v5 x5 y5 z5
v6 x6 y6 z6
… … … …
Triangles
t0 v0 v1 v2
t1 v0 v1 v3
t2 v2 v4 v3
t3 v5 v2 v6
… … … …
22
Right-Hand Rule
23
http://viz.aset.psu.edu/gho/sem_notes/3d_fundamentals/html/3d_coordinates.html
http://mathinsight.org/stokes_theorem_orientation
Normal Computation
24
Orientability
25
orientable non-orientable
Summary of Polygonal Meshes
• Polygonal meshes are piece-wise linear
approximation of smooth surfaces
• Good triangulation is important (manifold, equi-length)
• Vertices, edges, and faces are basic elements
• While real-data 3D are often point clouds, meshes are
quite often used to visualize 3D and generate ground
truth for machine learning algorithms
26
Polygonal Meshes
• Representation
• Storage
• Curvature Computation
Challenge on Meshes
28
http://upload.wikimedia.org/wikipedia/commons/f/fb/Dolphin_triangle_mesh.png
Curvature is a
second-order derivative,
but triangles are flat.
29
Assume a local at a small triangle
Assume that ’s are roughly parallel
Assume that , i.e.,
( is also aligned with ’s)
f : U → ℝ3
Tpi
Df [uv] = u ⃗ξ u + v ⃗ξ v Df = [ ⃗ξ u, ⃗ξ v]
U Tpi
e0
e2
e1
p1
p2
p0
ξu
ξv
ξu
ξv
ξu
ξv
Rusinkiewicz’s Method
30
Recall shape operator: .
If we have , we can compute principal curvatures!
How to estimate ?
DN = Df ⋅ S
∵ Df = [ ⃗ξ u, ⃗ξ v], ∴ S = DfTDN
S
S
Rusinkiewicz’s Method
e0
e2
e1
p1
p2
p0
ξu
ξv
ξu
ξv
ξu
ξv
31
e0
e2
e1
p1
p2
p0
ξu
ξv
ξu
ξv
ξu
ξv
DfT (DN [uv]) ≈ DfTΔ ⃗n ⟹ S [uv] ≈ DfTΔ ⃗n
∵ Df [uv] = Y ∈ T(ℝ3) and Df = [ ⃗ξ u, ⃗ξ v] ∴ [uv] = DfTY
∴ S[Df ]TY ≈ [Df ]TΔ ⃗n
32
DfT (DN [uv]) ≈ DfTΔ ⃗n ⟹ S [uv] ≈ DfTΔ ⃗n
e0
e2
e1
p1
p2
p0
ξu
ξv
ξu
ξv
ξu
ξv
S[Df ]Te0 = Df
T( ⃗n 2 − ⃗n 1),
S[Df ]Te1 = Df
T( ⃗n 0 − ⃗n 2),
S[Df ]Te2 = Df
T( ⃗n 1 − ⃗n 0),
So we can solve by
least square(6 equations and
4 unknowns)
S ∈ ℝ2×2
∵ Df [uv] = Y ∈ T(ℝ3) and Df = [ ⃗ξ u, ⃗ξ v] ∴ [uv] = DfTY
∴ S[Df ]TY ≈ [Df ]TΔ ⃗n
Summary of Mesh Curvature Estimation
• Rusinkiewicz’s method is an effective approach for
face curvature estimation
– Szymon Rusinkiewicz, “Estimating Curvatures and Their
Derivatives on Triangle Meshes”, 3DPVT, 2004
• Good robustness to moderate amount of noise and
free of degenerate configurations
• Can be used to compute curvatures for point cloud as
well
33
Point Cloud
• Representation
• Sampling Points on Surfaces
• Normal Computation
Ack: Sid Chaudhuri
Acquiring Point Clouds
35
• From the real world
• 3D scanning
• Data is “striped”
• Need multiple views to
compensate occlusion
• Many techniques
• Laser (LIDAR, e.g., StreetView)
• Infrared (e.g., Kinect)
• Stereo (e.g., Bundler)
• Many challenges: resolution,
occlusion, noise, registration
Acquisition Challenges
36 htp://grail.cs.washington.edu/
Occlusion
Interiors not
captured
→
Some data was
not properly
registered with
the rest
Noise Poor detail reproduction→ Low resolution further obscures detail
• From existing virtual shapes
• Why would we want to do this?
Acquiring Point Clouds
37
Light-weight Shape Representation
Point cloud:
• Simple to understand
• Compact to store
• Generally easy to build algorithms
Yet already carries rich information!
38
N = 125 N = 250 N = 500 N = 1000
Point Cloud
• Representation
• Sampling Points on Surfaces
• Normal Computation
Ack: Sid Chaudhuri
Application-based Sampling
• For storage or analysis purposes (e.g., shape
retrieval, signature extraction),
– the objective is often to preserve surface
information as much as possible
• For learning data generation purposes (e.g.,
sim2real),
– the objective is often to minimize virtual-real domain
gap
40
Application-based Sampling
• For storage or analysis purposes (e.g., shape
retrieval, signature extraction),
– the objective is often to preserve surface
information as much as possible
• For learning data generation purposes (e.g.,
sim2real),
– the objective is often to minimize virtual-real domain
gap
41
Naive Strategy: Uniform Sampling
• Independent identically distributed (i.i.d.) samples by
surface area:
• Usually the easiest to implement (as in your HW0)
• Issue: Irregularly spaced sampling
42
Farthest Point Sampling
• Goal: Sampled points are far away from each other
• NP-hard problem
• What is a greedy approximation method?
43
Iterative Furthest Point Sampling
• Step 1: Over sample the shape by any fast method
(e.g., uniformly sample N=10,000 i.i.d. samples)
44
Iterative Furthest Point Sampling
• Step 2: Iteratively select K points
45
is the initial big set of points
add a random point from to
for i=1 to K
find a point with the largest distance to
add to
U
S = {}
U S
u ∈ U S
u S
Issues Relevant to Speed
• Theoretically, naive implementation gives , but
how to improve from is an open question
• Implementation can cause large speed difference
– As this is a serial algorithm in K, engineers optimize
the efficiency in N (computing point-set distance)
– CPU: Suggest using vectorization (e.g., numpy,
scipy.spatial.distance.cdist)
– GPU: By using shared memory, the complexity can
be reduced to , where M is
the number of threads (M=512 in practice for
modern GPU).
𝒪(KN)
𝒪(KN)
𝒪(K(N/M + log M))
46Read by yourself!
Implementation Tricks
• References:
– https://github.com/maxjaritz/mvpnet/blob/master/
mvpnet/ops/cuda/fps_kernel.cu
– https://github.com/erikwijmans/Pointnet2_PyTorch/
blob/master/pointnet2_ops_lib/pointnet2_ops/_ext-
src/src/sampling_gpu.cu
• By courtesy of Jiayuan Gu, we share a GPU version
code with you (through Piazza)
47Read by yourself!
https://github.com/maxjaritz/mvpnet/blob/master/mvpnet/ops/cuda/fps_kernel.cu
https://github.com/maxjaritz/mvpnet/blob/master/mvpnet/ops/cuda/fps_kernel.cu
https://github.com/erikwijmans/Pointnet2_PyTorch/blob/master/pointnet2_ops_lib/pointnet2_ops/_ext-src/src/sampling_gpu.cu
https://github.com/erikwijmans/Pointnet2_PyTorch/blob/master/pointnet2_ops_lib/pointnet2_ops/_ext-src/src/sampling_gpu.cu
https://github.com/erikwijmans/Pointnet2_PyTorch/blob/master/pointnet2_ops_lib/pointnet2_ops/_ext-src/src/sampling_gpu.cu
An Implementation
in Numpy
def fps_downsample(points, number_of_points_to_sample):
selected_points = np.zeros((number_of_points_to_sample, 3))
dist = np.ones(points.shape[0]) * np.inf # distance to the selected set
for i in range(number_of_points_to_sample):
# pick the point with max dist
idx = np.argmax(dist)
selected_points[i] = points[idx]
dist_ = ((points – selected_points[i]) ** 2).sum(-1)
dist = np.minimum(dist, dist_)
return selected_points
48Read by yourself!
Voxel Downsampling
• Uses a regular voxel grid to downsample, taking one
point per grid
• Allows higher parallelization level
• Generates regularly spaced sampling (with noticeable
artifacts)
49 Credit to: Open3D
Issues Relevant to Speed
• Need to map each point to a bin. Often implemented
as adding elements into a hash table
• (assuming that the inserting into hash table
takes O(1))
• On GPUs, parallelization reduces complexity of
– Mapping each point to an integer value
– Assign each value to an index so that the same
value shares the same index
– Aggregate indexes and form the output (called
scattering in CUDA)
𝒪(N)
50Read by yourself!
A Dictionary-based Implementation
in Numpy
def voxel_downsample(points: np.ndarray, voxel_size: float):
“””Voxel downsample (first).
Args:
points: [N, 3]
voxel_size: scalar
Returns:
np.ndarray: [M, 3]
“””
points_downsampled = dict() # point in each voxel cell
points_voxel_coords = (points / voxel_size).astype(int) # discretize to voxel
coordinate
for point_idx, voxel_coord in enumerate(points_voxel_coords):
key = tuple(voxel_coord.tolist()) # voxel coordinate
if key not in points_downsampled:
# assign the point to a voxel cell
points_downsampled[key] = points[point_idx]
points_downsampled = np.array(list(points_downsampled.values()))
return points_downsampled
51Read by yourself!
A Unique-based Implementation
in Torch
def voxel_downsample_torch(points: torch.Tensor, voxel_size: float):
“””Voxel downsample (average).
Args:
points: [N, 3]
voxel_size: scalar
Returns:
torch.Tensor: [M, 3]
“””
points = torch.as_tensor(points, dtype=torch.float32)
points_voxel_coords = (points / voxel_size).long() # discretize
# Generate the assignment between points and voxel cells
unique_voxel_coords, points_voxel_indices, count_voxel_coords = torch.unique(
points_voxel_coords, return_inverse=True, return_counts=True, dim=0)
M = unique_voxel_coords.size(0) # the number of voxel cells
points_downsampled = points.new_zeros([M, 3])
points_downsampled.scatter_add_(
dim=0,
index=points_voxel_indices.unsqueeze(-1).expand(-1, 3),
src=points)
points_downsampled = points_downsampled / count_voxel_coords.unsqueeze(-1)
return points_downsampled
52Read by yourself!
Application-based Sampling
• For storage or analysis purposes (e.g., shape
retrieval, signature extraction),
– the objective is often to preserve surface
information as much as possible
• For learning data generation purposes (e.g.,
sim2real),
– the objective is often to minimize virtual-real domain
gap
– a good research topic (e.g., GAN? Adversarial
training? Differentiable sampling?)
53
Point Cloud
• Representation
• Sampling Points on Surfaces
• Normal Computation
Ack: Sid Chaudhuri
Estimating Normals
• Plane-fitting: find the plane that best fits the
neighborhood of a point of interest
55
Least-square Formulation
• Assume the plane equation is:
with
• Plane-fitting solves the least square problem:
where is the neighborhood of a point
that you query the normal
wT(x − c) = 0 ∥w∥ = 1
minimize ∑
i
∥wT(xi − c)∥
2
2
subject to ∥w∥2 = 1
{xi} x
56
• Doing Lagrangian multiplier and the solution is:
– Let and ,
– : the smallest eigenvector of
–
• also corresponds to the third principal component of
(yet another usage of PCA)
– Where are the first and second principal
components?
M = ∑
i
(xi − x̄)(xi − x̄)
T x̄ =
1
n ∑
i
xi
w M
c = wT x̄
w
M
57
Summary of Normal Computation
• The normal of a point cloud can be computed through
PCA over a local neighborhood
• Remark:
– The choice of neighborhood size is important
– When outlier points exist, RANSAC can improve
quality
58