CS计算机代考程序代写 data structure cuda GPU algorithm L4_Representations

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

Mobius Strip Sphere LightForm Blog

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