Region Reconstruction
Uses & Techniques

Alec Shedelbower, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Talk Overview

Definition

Applications

Techniques

Fitting Methods

Mesh Reconstruction Methods

What is Region Reconstruction?

Goal: Given a finite set of points sampled from (or near) a region, return a close approximation of the original region.

({{, , }})

Types of Reconstruction

Reconstruction is not just limited to curves in 2D.

({{"", "Curve", "Surface", "Solid"}, {"1D", , , }, {"2D", , , }, {"3D", , , }})

Applications – 3D Object Acquisition

Consumer-grade 3D scanners are becoming more common. Even newer smart phones come with built-in LiDAR support.

These scanners produce point clouds – collections of points with optional auxiliary data. Post processing is required to reconstruct the object surface.

({{, , }, {"Point Cloud", "", "Reconstructed Surface"}})

Applications – Plane Detection

Another use for 3D scanning is plane detection. In this case, we try to find the planar regions that can approximate the point cloud.

This tends to be useful for manmade / urban environments.

Applications – Function Reconstruction

Region reconstruction can be used to approximate a function from a sample points.

({{, , }})

Techniques

Fitting Methods

Assume the general shape – or model – of the sampled region is known. Models are parameterized, such as a circle parameterized by its radius and center.

({{}})

Unique Fit

Models can often be uniquely defined by only a handful of points. For example, the symbols InfiniteLineThrough and CircleThrough require only two (2) and three (3) unique points, respectively:

({{, }, {InfiniteLineThrough, CircleThrough}})

Unique Fit – Considerations

Pros

Requires only a few sample points

Cons

Restricted to specific models

Sample points must be directly sampled from model region

Real world data is often slightly randomized due to noise. How do we choose the best fit when no single parameterization can explain all the sample points?

Error Minimization

Define an error metric between sample points and a parameterized model. The best fit is the one where the total error is minimized.

Most common metric is least squares, where the error being minimized is the sum of the squared distances from each point to the model.

LeastSquaresDemo[]

Example Code

data=Table[{x,x/2+RandomReal[1/2]},{x,0,20,1}];
lineEq=Fit[data,{1,x},x]

Show[Plot[lineEq,{x,0,20},PlotStyle->{,Dashed}],ListPlot[data],ImageSize->Large]

Error Minimization – Considerations

Pros

Robust to noise

Fast for certain cases (lines, polynomials, ...)

Cons

Restricted to specific models

Not robust to outliers

Error Minimization – Outliers

Outliers are erroneous points that do not represent the sampled region. Error minimization techniques can be easily skewed by a small percentage of outliers.

LeastSquaresDemo[]

Random Sample Consensus (RANSAC)

Instead of trying to collectively reduce the error for all sample points, try to find a subset of points that agree with each other.

How do we do this? Hold an election!

Randomly select a candidate parameterization for our model using a handful of points.

Calculate the distance from each point to the candidate region. If a point is close enough to the candidate region, then the point “votes” for that candidate.

After selecting some number of random candidates, the candidate with the most votes “wins” the election and is returned as the best fit.

RANSAC – Workflow

Below shows the process of fitting a circle to a set of points with double the number of outliers than inliers.

({{({{}, {"Minimal Random Sample"}}), , ({{}, {"Candidate (Unique Fit)"}}), , ({{}, {"Tolerance"}}), , ({{}, {"Votes (50)"}})}, {({{}, {"Minimal Random Sample"}}), , ({{}, {"Candidate (Unique Fit)"}}), , ({{}, {"Tolerance"}}), , ({{}, {"Votes (29)"}})}, {({{}, {"Minimal Random Sample"}}), , ({{}, {"Candidate (Unique Fit)"}}), , ({{}, {"Tolerance"}}), , ({{}, {"Votes (7)"}})}})

RANSAC – RegionFit

RegionFit performs RANSAC to fit a region model to a given point cloud.

Basic Example

points=Join[
RandomPoint[Circle[{1,1},3],300],
RandomReal[{-3,5},{2000,2}]
];
Graphics[Point[points]]

RegionFit[points,"Circle"]

Graphics[{%,Point[points]}]

Properties

Determine how “noisy” the data was:

RegionFit[points,"Circle","DistanceVariance"]

Color points based on their distance to the best fit:

dists=RegionFit[points,"Circle","Distances"];
dists//Short

cfunc=ColorData["BlueGreenYellow"];
colors=cfunc/@Rescale[dists,MinMax[dists],{0,1}];
Graphics[{PointSize[0.01],Transpose[{colors,Point/@points}]}]

Determine which points were considered inliers/outliers:

{inliers,outliers}=RegionFit[points,"Circle",{"Inliers","Outliers"}];
Graphics[{Point[inliers],LightGray,Point[outliers]}]

Models

RegionFit currently supports the following models:

({{"Line", "Circle", "Plane", "Cylinder", "Sphere", "Cone"}, {, , , , , }})

Would like to extend RegionFit to more models in the future.

RANSAC – Considerations

Pros

Robust to noise and outliers

Configurable via tolerance and candidate count

Tolerance and ConfidenceLevel options for RegionFit.

More flexible than least squares

Cons

Restricted to specific models

Stochastic

Mesh Reconstruction Methods

Mesh reconstruction returns a discrete approximation of the original region. The examples we will cover today fall into two general categories:

Filtering – Vertices of mesh form a subset of the original point cloud, with simplicial cells (lines, triangles, tetrahedra) connecting them.

Implicit – Point cloud is used to approximate an implicit function. The resulting mesh is an approximation of a level set of this function.

Delaunay Mesh

Foundational meshing method that constructs simplicial cells from a set of points such that the resulting mesh has “nice” properties:

Maximizes minimum interior angle (avoids thin “sliver” simplices)

Boundary is convex hull

Empty circumcircle/circumsphere of simplices

({{, , }})

Delaunay Mesh

By itself, the Delaunay mesh creates a good reconstruction for convex, full-dimensional regions such as disks and cylinders.

DelaunayMesh/@{RandomPoint[Disk[],500],RandomPoint[Cylinder[],5000]}

The same is true for its boundary – the convex hull:

ConvexHullMesh/@{RandomPoint[Disk[],500],RandomPoint[Cylinder[],5000]}

However, in order to reconstruct non-convex and/or lower-dimensional regions, additional work needs to be done.

Concave Hull

Generalization of convex hull mesh to non-convex regions. Also known as “Alpha Shapes”.

Two key insights:

Insight 1: Delaunay mesh can act as a superset of mesh cells for reconstruction.

Insight 2: Cells with larger circumsphere radii are less likely to be part of the original region.

This means we can filter the Delaunay cells with large circumsphere radii to get the reconstructed region.

Concave Hull – Algorithm

Compute Delaunay mesh of the sample points.

({{, , , , }, {"Region", "", "Point Samples", "", "Delaunay"}})

Determine the circumsphere radius of each full-dimensional cell.

({{, , (Legended[, (BarLegend)])}, {"Circumspheres", "", "Radii"}})

Remove cells with a radius greater than some value α.

({{}, {"Concave Hull"}})

Note that if α is large enough, none of the cells will be filtered and the resulting region will be equivalent to the convex hull.

Concave Hull – Interactive (Radius)

ConcaveHullDemo[]

Concave Hull – Interactive (Scoop)

Imagine vertices are rigid points. Can try to “scoop” out the area using a disk with a minimum radius of α.

ConcaveHullScoopDemo[]

Concave Hull – Dimensions

Concave hull works best on full-dimensional regions.

1D

pts=RandomPoint[Ball[{0}],10];
ConcaveHullMesh[pts,1]

2D

pts=RandomPoint[Annulus[],500];
ConcaveHullMesh[pts,0.3]

3D

pts=RandomPoint[DiscretizeRegion[Torus[]],1000];
ConcaveHullMesh[pts,0.3]

Concave Hull – Lower Dimensional Regions

With some adjustments, we can reconstruct lower-dimensional regions. This is done by filtering these lower-dimensional simplices using the same circumsphere criteria.

pts=Table[{t/2,Cos[t],Sin[t]},{t,0,8 Pi,Pi/20}];
ConcaveHullMesh[pts,Automatic,1]

pts=Flatten[Table[,{θ,1/100,π,0.05},{ϕ,1/100,2 π,0.05}],1];
ConcaveHullMesh[pts,Automatic,2]

Concave Hull – Considerations

Pros

Extendable to any dimension

Good at reconstructing densely sampled full-dimensional regions

Cons

Reconstructed mesh vertices are restricted to a subset of sample points

Weak to noise or outliers

Lower dimensional reconstructions tend to be open and non-manifold.

No obvious choice for the best α value

“Good” α value may not exist due to sparse or uneven sampling

Crust Method

Another Delaunay filtering method designed to reconstruct curves in 2D. Utilizes the concept of a medial axis.

Medial Axis

Set of all points that have more than one closest point on the region’s boundary. Forms a sort of skeleton for a region.

({{, , }})

Crust Method

The crust method improves upon curve construction using the following insights:

Insight 1: Delaunay edges that cross the medial axis should be filtered out.

Insight 2: Voronoi diagram can approximate medial axis of a region from points sampled on its boundary.

Crust – Algorithm

Compute the Voronoi diagram of the sample points.

({{, , , , }, {"Region", "", "Point Samples", "", "Voronoi"}})

Compute the Delaunay mesh of the sample points augmented by the Voronoi vertices.

({{, , }, {"Points (Black) &\nVoronoi Vertices (Red)", "", "Augmented Delaunay"}})

Select the Delaunay edges that connect the original sample points.

({{, , }, {"Remaining Edges", "", "Reconstructed Region"}})

Crust – Demo

FindCurvePath can be used to find the ordering of the sample points using the Crust method.

pts=RandomPoint[Circle[],100];
Graphics[{Line[pts],PointSize[0.01],Point[pts]}]

path=FindCurvePath[pts];
path//Short

Graphics[{Line[pts,path]}]

ReconstructionMesh

ReconstructionMesh also uses the crust method when applicable.

ReconstructionMesh[pts]

ReconstructionMesh[]

Crust – Considerations

Pros

Gives good results

Less likely to give open or non-manifold curves than Concave Hull

Provides provable guarantees assuming some minimum sampling requirements

Cons

Not robust to noise

Limited to curves in 2D

Extensions exist for surfaces in 3D (Power Crust), but this is not currently implemented in the Wolfram Language.

Implicit Reconstruction Methods

Filtering methods are limited by the density of the point samples and struggle in the presence of noise.

In contrast, implicit methods avoid this by reconstructing an implicit function to describe the original region. We can then obtain a mesh by taking the level set of this function using some type of isosurface extraction (ex. marching cubes).

Poisson Reconstruction

Poisson surface reconstruction is the standard method used for reconstructing scanned objects from point clouds.

Its goal is to first approximate the indicator function (i.e. implicit function) of the sampled region. The indicator function returns 1 if a point is inside the region, and 0 otherwise.

Poisson Reconstruction – Poisson’s Equation

This indicator function is approximated using the solution to Poisson’s equation (hence the name).

In this case, the solution lets us calculate the indicator function from its gradient.

Insight: The gradient of the indicator function for a region can be approximated from oriented points sampled from its boundary.

({{, " = ∇", }, {"Oriented Points", "", "Indicator Function"}})

Poisson Reconstruction – Algorithm

Partition oriented sample points into cells (quadtree or octree).

Solve Poisson’s equation to approximate the indicator function.

Perform isosurface extraction (i.e. marching cubes) on the indicator function.

({{, , , , , , }, {"Oriented Points", "", "Quadtree (or Octree)", "", "Estimated Indicator Function\n(Poisson Solver)", "", "Isosurface Extraction"}})

Poisson Reconstruction – Normals

Like most implicit methods, Poisson surface reconstruction requires oriented points, or points with normal directions.

This normal is a unit vector perpendicular to the surface at that point. In practice, point clouds rarely include point normals, so they must first be estimated.

Poisson Reconstruction – Normal Estimation

Normal estimation can be accomplished using the plane fitting method. It uses principal component analysis (PCA) to fit a plane to each point on the surface.

Select a point from the point cloud.

Find the k-nearest neighbors of the selected point (k = 4 in this example).

Calculate the normal vector as the eigenvector with the minimum eigenvalue of the covariance matrix of the point and its neighbors.

Repeat steps 1-3 for all points.

({{, , , , , , }})

Poisson Reconstruction – Examples

Poisson surface reconstruction is currently achieved with GradientFittedMesh.

pts=Point[,VertexNormals->];
Graphics3D[pts]

GradientFittedMesh[pts]

models={ResourceData["Stanford Bunny"],ResourceData["Horse"],ResourceData["Cow"]};
coords=RandomPoint[#,10000]&/@models;
GradientFittedMesh/@coords

Normal Estimation

Can estimate normals manually using EstimatedPointNormals:

pts=RandomPoint[ResourceData["Stanford Bunny"],10000];
normals=EstimatedPointNormals[pts];
Graphics3D[Point[pts,VertexNormals->normals]]

GradientFittedMesh[Point[pts,VertexNormals->normals]]

Poisson Reconstruction – Considerations

Pros

Good for reconstructing “smooth” surfaces

Can extrapolate surface to missing or sparsely sampled areas.

Resulting mesh resolution can be adjusted

Cons

Requires oriented points (though these can be estimated)

Struggles with sharp corners or intentional holes in the surface.

Reconstruction Mesh

As we have seen, there are a variety of mesh reconstruction methods. They each have their own advantages and limitations.

As a user, you don’t want to have to worry about which method is best for your particular use case. This led to the creation of ReconstructionMesh.

ReconstructionMesh takes a point cloud and attempts to automatically determine and apply the best available reconstruction method.

Reconstruction Mesh – Examples

pts=Table[(2+2/3.5Sin[3 t])*RotationMatrix[t].{1,0},{t,0,2 Pi,0.15}];
ReconstructionMesh[pts]

ReconstructionMesh[RandomPoint[Annulus[],1000]]

ReconstructionMesh[MeshCoordinates[ResourceData["Horse"]]]

Conclusion

For more information:

Watch for upcoming Wolfram R&D Livestreams

View the documentation pages

Questions?

Library (Initialization Code)

Sampling

EquidistantRandomPoints

Clear[EquidistantRandomPoints]
Options[EquidistantRandomPoints]={RandomSeeding->Automatic};
EquidistantRandomPoints[region_,pointNum_,OptionsPattern[]]:=
Module[{seed,edim,pts,stepDist,rm,regionRange,beFar,InRegion,moveDir},
seed=OptionValue[RandomSeeding]/.Automatic->RandomInteger[10000];
SeedRandom[seed];
edim=RegionEmbeddingDimension[region];
pts=RandomPoint[region,pointNum];
rm=RegionMember[region];
regionRange=EuclideanDistance@@Transpose[RegionBounds[region]]/100;
stepDist=Subdivide[regionRange,0,500];
moveDir=RandomReal[#,pointNum]*RandomPoint[If[edim===2,Circle[],Sphere[]],pointNum]&/@stepDist;
Fold[(beFar=Boole@Thread[EuclideanDistance@@@Nearest[#,#,2]<EuclideanDistance@@@Nearest[#1+#2,#1+#2,2]];
InRegion=Boole[rm/@(#1+#2)];
#+#2*beFar*InRegion)&,N[pts],moveDir]
]
EquidistantRandomPoints[___]:=$Failed

Formatting

Format Graphic

Clear[FormatGraphic]
FormatGraphic[expr_]:=CellPrint[Cell[ToBoxes[expr],"Text",ShowCellLabel->False,ShowCellTags->False,ShowStringCharacters->False,GeneratedCell->False]]

Illustrations

Constants

$PointSizeSmall=AbsolutePointSize[2];
$PointSizeMedium=AbsolutePointSize[4];
$PointSizeLarge=AbsolutePointSize[6];

Arrow

Clear[DocumentationArrow]
DocumentationArrow[size_:40]:=Graphics[{LightGray,EdgeForm[LightGray],Rectangle[{-2,-1},{0,1}],Triangle[{{0,-2},{2,0},{0,2}}]},ImageSize->size]

Omitted Slides

Region Reconstruction – Ambiguity

Can you determine which of the dashed regions the points were originally sampled from?

({{, , }})

The answer is ambiguous, as the points could have been sampled from any of the regions. In fact, those points belong to an infinite number of regions, all of which are valid candidates for the original. This lack of a unique solution is why region reconstruction is considered an ill-posed problem.

To avoid this, reconstruction techniques must make assumptions regarding either the sample points or underlying region.

Concave Hull – Choice of α

There is no obvious choice for α in the general case, and different values will drastically change the resulting region. This means the optimal value often has to be found through trial and error.

Heuristics

There have been some proposed heuristics for finding automatically determining a “good” α value.

ConcaveHullMesh will attempt to find a good approximation for α by finding the maximum distance of any point in the set to its nearest neighbor. This generally gives good results when the point cloud provides a sufficiently dense and uniform sampling of the original region.

pts=RandomPoint[Annulus[],2000];
ConcaveHullMesh[pts]

No Good Alpha

Reconstructing lower-dimensional regions (ex. curve in 2D), there may be no “good” choice for α.

Manipulate[ConcaveHullMesh[,α,d,ImageSize->1000],{{α,0.4},0,0.4},{d,{1,2}}]