Region Reconstructions: Uses and Techniques

Hello everyone. My name is Alex Schellbauer. I work on the Discrete Computation team with Charles and jbum. And today I'm going to be talking about region reconstruction. So this is something we've been adding new features for over the past few years and today I'm going to kind of COVID what's going on with those. So kind of overview for today's talk. I'm going to start with some definitions and applications for this, then go into the different techniques we can use to accomplish it. So first up, what is region reconstruction? So the goal is given a finite set of points sampled from or near a region, we want to return a close approximation of that original region, which kind of like again connected dots, where you have a bunch of finite points and want to find the curve that goes through them all. However, unlike connected dots, we are not limited to just 2D curves. We can do curves in 1D, 2D, 3D, as well as surfaces and solids. So regional reconstruction covers them all. As far as applications, probably the most common application here is 3D object acquisition. So with kind of 3D scanners being more affordable now, and even some smartphones having built in wider support, being able to kind of scan real world objects is kind of more and more common. And these scanners produce point clouds, which are basically collections of points with optional metadata. So we can't really do much with these point clouds by themselves. However, we usually want to reconstruct the actual surface using region reconstruction that we can then do stuff like render or apply physics to. So for example, here I have a model Dragon that I scanned in with a desktop scanner. I use Mathematica to reconstruct it here on the right. Some applications here there's like plane detection, which is kind of similar, but instead of trying to do a more general mesh, you just look for planes. This is good for like urban environments or kind of like indoor man made interiors. There's also kind of function reconstruction where if you just want to reconstruct a function from a bunch of point samples, you can do that as well. Yeah, there's more, but that's all I'll cover for now. Right now we're going to go into techniques. So kind of the first category techniques here we have the fitting methods for these. We kind of assume the general shape or model of the sampled region is already known. This model is parameterized, for example, like a circle being parameterized by its center and its radius. And then sort of kind of the first method that uses this category is a unique fit, which is pretty Simple. It's basically just a lot of these models can be uniquely determined from only a few points. For example, a line can be determined by just two points and a circle by just three. And we have even functions for doing this. Infinite line through and circle through which will produce a line or circle respectively. It's kind of a very basic but it's good to know method for region reconstruction. So it requires only a few sample points. That's a good pro. However, cons it is restricted to a few specific models and the sample points must be directly sampled from the model region. This is a problem because a lot of times real world data is slightly randomized due to noise. So the question is, how can we kind of find a best fit when not all the points agree on what it should be? We're going to see kind of two answers to this. The first one here is error minimization. If you're taking a statistics class, this is probably pretty familiar to you. Basically define some error metric between our points in the model and we'll want to find the model that has sort of the minimum error. The most common metric is least squares. We have kind of a visualization here. For example, we have this kind of blue line that we can control and the error is the distance from each point to that line squared, kind of visualized by these red squares here. So the area of all these squares summed up will be the total error that we minimize at the best fit, which you can see is about here. Screen down with a few functions. In Mathematica, one of the probably easiest ones is this fit function. We can just run it here, get our equation of our line and then pull out it for some points. So for this one, this one's actually robust to the noise. Unlike the previous one, it's fast for certain cases like lines and polynomials. However, it's still restricted to certain models and is now robust to outliers. To go into more detail on that. So, for example, that previous example here, if we add some outliers, so points that don't really represent the surface at all, our best fit gets skewed very dramatically. So now it doesn't really fit anything well. So this is kind of a problem when you're working with real world data, especially scanner data, because that often has a lot of outliers. So kind of the solution to this is a method called random sample consensus, or ransack for short here. Instead of trying to collect reduce the error for all sample points, we just try to find a certain subset that agree with each Other. And we do this by holding a basic election. So what we're going to do is create a bunch of random candidates and we're going to calculate a distance from each point to these candidates. And if a point is in a certain tolerance of that candidate, it will vote for it. The very end, we'll tally up all the votes from all the points and the candidate with the most votes wins and is kind of considered the best fit. Kind of visualize this. Let's say you're trying to fit a circle. We have all these sort of random noise points around it. What we start by doing is for our first candidate, we choose three random points. We'll then do that unique fit to find a circle through them. We'll kind of expand that fit by some tolerance value and look at all the points in that kind of tolerance area and then count them as votes. So we can do that multiple times. Here we see the first candidate here is very lucky. It got three points on the actual circle, so it had a lot of votes, whereas last one only had one point on a circle. So it had only seven votes. So in this election, this first candidate would win and be returned as the actual circle. So in actual code, this is done with the new function region fit, I guess, newer function. So this performs ransack for us. So we have, for example, a very noisy kind of point set with a lot of outliers. We can pass those to region fit and ask for a circle. And we'll find the best fit circle, which pretty nicely matches our actual circle there. We can also get more properties about this. So for example, we can ask for the distance variance, which kind of tells us how noisy all the data was. We can get the actual distances themselves from each point to the best fit candidate. And we can kind of color points based on their distance. And we can also just get the inliers and outliers. So we can see all these darker points voted for our best fit candidate and all the other ones did not. And then as far as models, here are the currently supported models for region fit. We have line, circle, plane, cylinder, sphere and cone. We would like to add more in the future, but for now this is what we have considerations for this. So one, it's robust to noise and outliers, it's configurable via tolerance and candid count settings, and it's a little more flexible in the least squares, kind of easy to extend to new models. However, it's still restricted to specific models and it's also stochastic. So since there's some random sampling involved. You're not guaranteed to get the right answer. Just you can only work with like probabilities of getting the right answer. All right with that. That's all the kind of fitting methods. So now we'll go into the mesh reconstruction methods. So here I kind of we're trying to fit an actual mesh now, like a mesh region, to our point set. They're kind of two main strategies for this. The first one is called filtering and the other one's implicit. So filtering is more like those connected dot examples I showed earlier, where we keep the same input points, but we just try to connect them with lines, triangles or tetrahedra. The implicit examples, however, try to reconstruct the actual implicit function for the region and then they do some sort of level set in order to get the actual mesh from this. So kind of start with the filtering ones. Kind of the basis of most filtering methods is a Stalini mesh. It's kind of a foundational meshing method that kind of constructs these simplicial cells between our input points. That's like these lines, triangles and then 3D tetrahedra. And it does in such a way that they have nice properties, some main ones here. So it maximizes the minimum interior angle, which helps avoid kind of narrow sliver triangles. And also the boundary is a convex hole, which is kind of convenient in terms of using this. So the line mesh is actually a decent candidate for reconstruction method. If the kind of region in question is convex and full dimensional, like a disk or a cylinder here. And kind of its cousin, the convex hull mesh, it kind of gives the same result, except it doesn't include the interior cells. However, if you want to do some more general stuff like some non convex regions or lower dimensional regions, like a curve in 3D, we're going to have to do a little more work here. That brings us to the concave hole. So this is kind of a generalization of that convex hull mesh to non convex regions. You also hear sometimes called the alpha shape. So kind of two key insights this brings us. The first one is the line. You can act as a superset of mesh cells for reconstruction and that its cells with the larger circumsperior radii are less likely to be part of the original region. What that means is that we can, instead of having to worry about connecting points ourselves and just run a Delau mesh over all the points. And then we can selectively filter out the cells that are too big, that don't do a good Job of representing the region. So it's kind of visualize that process. Let's say we start with this annulus here as our ground truth region. We take some point samples and then we can run the line mesh on that. We can then look at each of the triangles in our mesh and do the circumsphere of it. And each of these circumspheres or circles in this case has a radius that we can then kind of color code triangles based on these radii. And you can see the ones in the center here where there aren't many points, a very large radii and are shown as red. We can then choose some cutoff value, some threshold we'll call alpha. And this will be kind of our cutoff where we don't include any triangles with a radius above this. And that gives us our final concave fold down here does a pretty decent job of representing our original region. Kind of click interactive demo of this. Here I have a smiley face that I sampled from. You can see that the eyes and the mouth didn't have any points sampled from it. So these triangles are kind of stretching across the entire surface. Right now our alpha is very high, so it includes all the triangles. But I start decreasing it when those first de mouth triangles start to get filtered out, followed by the eyes. And then we have a pretty nice shape that we want. However, if we continue to do this, we'd slowly lose all the cells as the alpha is less than all of the radii. Another kind of common way this is sometimes described is you're kind of like scooping the triangles with kind of a disk that you kind of move around and sort of destroy these triangles. And the disk has a minimum radii, which is the alpha radius, which is alpha. And it cannot sort of pass through points. So for example, I can get these like wider regions, but if I try to go inside this shape, it kind of gets stuck there and won't move. So we do this and it produces a good, you know, nice hole of this hole. All right. Concave pole is pretty flexible. It works on any dimension. So we can have 1D, 2D and 3D. These are all full dimensional. However, we can tweak it a bit to make it lower dimensional. For example, we can have a curve in 3D or a surface with triangles in 3D. We can do this by just instead of filtering out the polymenal cells of the Delaunay mesh, we can filter out these lower dimensional ones and do a little processing to help remove redundant Ones as well. All right, with all that said, so concave full considerations. So it's very kind of extendable to any dimension. It's good ever constructing those densely sampled full dimensional regions like annulus or cylinders, spheres, stuff like that, or I guess balls. The cons, however, is that the reconstructed meshes, they're still restricted to those input samples. So it's so weak to noise and outliers. If you give it some bad points, it can't magically make those points better. Right. And then next up. So lower dimensional reconstructions tend to be open and non manifold, so they have some issues where they get a little messy and there's not always an obvious or a good choice for alpha. So kind of relieve this. We have this thing called the crest method. This is another deline filtering method designed to reconstruct curves in 2D. This utilizes the concept of a medial axis. I won't go into too much detail here, but the basic idea is the medial axis is kind of like a skeleton for these shapes where it contains all the points that have more than one closest point on the region's boundary. It's usually inside, but for some shapes it can also be on the outside if they're concave. Obviously this kind of star shape. The idea is that we don't want to include any edges that cross this medial axis skeleton. And then we can also kind of approximate that skeleton using the Veroni diagram. Does a pretty good job approximating even with some sample points. See this work? Let's say you want to reconstruct this kind of curve around the outer boundary of that star shape. We first take our point samples and then calculate the Voronoi diagram of it. Then take all the vertices of this mesh and add it to our original point set. You see, the vertices are color coded in red here. We're going to pass all of these to the Delaunay meshing algorithm and get this mesh. As a result, you can kind of see that it kind of separates the sides of our curve so that they cannot sort of short circuit and reach the other side. Acts like a barrier there. So now all we need to do is just select the edges that connect the original points from our point cloud. And once we do that, we're left with this nice curve that is our reconstructed region. In order to use this in Wolfram language, we have this fine curve path between a bunch of kind of random points, can rearrange them and give the correct ordering to make them a nice curve. We can also use this reconstruction mesh function. Give it the points, and then we can even use the whole example we saw before and see that does a good job of reconstructing that hole without having any short circuiting or any other issues. All right, so overall across the pros it gives good results and it's less likely to kind of create those non manifold or open kind of shapes. And it does provide provable guarantees assuming some minimal sampling requirements are met. However, like all the other filtering methods, it is not robust to noise and it is limited to curves in 2D for now. All right, so the last category and method we'll be talking about today is implicit reconstruction methods. Unlike filtering methods that are limited by density of point samples and kind of struggle with noise, implicit reconstruction methods try to avoid this by reconstructing the implicit function that describes the actual region. So we can obtain a mesh from this function by then using the level set function with some sort of ISO service extraction, like marching cubes. So kind of the main type of this method is called Poisson reconstruction or Poisson surface reconstruction. It's kind of a standard method, especially for like 3D object exposition. So its goal is to first approximate the indicator function, which here is basically to be a function that returns 1 when it's inside a region and 0 when it's outside the region. So kind of the Poisson and Poisson reconstruction is the fact that we can use the solution to Poisson's equation to get this indicator function from its gradient doesn't help us too much initially because we also do not know its gradient. But the key insight is that if we have a bunch of points sampled on the boundary of this shape, and we also have the orientations or normals of those points, that can act as an estimate for the gradient of this entire indicator function. Here you see on the left we have a bunch of points. They all have sort of the directions of the surface at that point, and that can act as our gradient. Once we have that, we can kind of step through the entire process. So once we have this oriented points, we can kind of partition them using some sort of quadtree or octree. We can then run a Poisson solver on that sort of data structure and kind of have a discrete approximation of our implicit function. And then we can then use whatever isosurface extraction you want, like marching cubes to get the actual curve. This does require normals. As I kind of mentioned, normals are not normally sort of measured directly with a scanner, so they most often be estimated. Luckily, this isn't too hard. There's some methods like plane fitting, that makes it pretty easy where the kind of idea here is given a point, we look at its nearest neighbors. In this case it's four nearest neighbors. We can then fit a plane to those points, weigh the normal at that point and assign that plane normal to that point. We just do this for every point and then it lets us kind of assign a point to all of our points or assign normal to all of our points in order to actually do this in the Wolfram language. So service reconstruction is done with gradient fitted mesh. So over here I can have some oriented points of Beethoven, pass them to gradient fitted mesh, and we get a kind of smooth Beethoven. We do this with a lot of our other kind of standard models here. So it's going to load the models, take a bunch of random points, and then we get the outputs here and then we can even do normal estimation directly using estimated point normals, which we can use on those points, and then pass those normals directly to the function ourselves without having to do what it calculated automatically. And yeah, so that's Poisson reconstruction. So it's good for those smooth surfaces. It can extrapolate missing or sparsely sampled areas and you can kind of adjust the resulting resolution of the mesh. However, it does require oriented points though. We can and re estimate those ourselves. And it does struggle with sharp corners or intentional holes on the surface. All right, so hopefully what you got from all that is that there are a lot of different variations and methods for how to do this and it can be a little confusing to know which ones to use. They all have their own bandages limitations. Because of this, we have added this function called reconstruction mesh, which acts as a high level function that will kind of just take your point cloud, whatever one you throw at it, figure out what dimension the region should be, and figure out the best method to use. Just quick examples. We have a curve in 2D, automatically knows to use that Crust algorithm. We can have our annulus and 2D and those to use concave hole, and we can have a surface in 3D and those use gradient fitted mesh. All right, and that concludes that part of my talk. For more information, we'll have an upcoming live stream on this talk probably soon. As always, the documentation pages are available if you want to kind of get more detail about how to use these functions and their options. All right, and with that, I'll use the remaining couple minutes here to answer any questions.

Alec Paritoj is wondering, what's the method or technique you use to actually handle noise? So if you have a Set of sample points with noise. So what do we use behind the scene?

Depends on the method for the fitting. We have the ransack method which uses that's region fit and uses ransack to kind of do that election where it finds the candidates and teams which one has the most votes for more general meshing right now, Poisson is our best one at that. And during this kind of computation of approximating this implicit function, there is some sort of noise reduction there. So you have kind of a noisy surface that can sort of take care of some of that during this stage before the final kind of mesh is extracted. Thank you.

So a question for min. So what's the benchmark and comparison between the different methods? Which one is the fastest and slowest?

It depends on the input. Yeah, we currently we are working on a benchmark in general using the ABC data set. Now I would say probably gradient fitted mesh takes a little extra work. Presumably we're just talking about like 3D meshes. This one's probably a little slower. That's because it's having to usually do some normal estimation itself. And then the solver can take a little time and then also it has to do the extraction at the end, whereas the concave hole, it just runs Delaunay and I can just filter those from there. So in general I'd say gradient fit is slower, but it's usually worth it for it. And concave hole is a little simpler, so it tends to run faster.

So Paritage is asking if the triangle after the end of the gradient fitted mesh will be, you know, self intersecting. If it's possible. They shouldn't be.

I don't believe they can be. So we're doing like this.

The output is a boundary mesh. So when you solve the pd, you actually get the actual boundary. So there's no self intersection during the reconstruction by this time. I like. Thank you.

Thank you everyone.