PetraSim – How we extract a surface from a set of XYZ points

Posted: March 3rd, 2015 in News

In the latest release of PetraSim we have improved the method we use to extract a fault surface from a set of XYZ points. The challenge is that we know the points are on the surface, but we do not know how the points are connected to define the surface. Previously, we projected the 3D points onto the XY plane, triangulated the XY coordinates, and used these connections to define the 3D surface. This works well as long as the surface is approximately horizontal. However, if the surface is vertical, then this projection fails.

For example, when the points shown in Figure 1 are projected onto the XY plane and triangulated, the surface has an error caused by overlapping triangles, Figure 2.

Figure 1 - XYZ points that define a fault surface.
Figure 1 – XYZ points that define a fault surface. Click to enlarge.
Figure 2 - Error in fault surface that occurs due to projecting points onto XY plane and triangulating in XY plane.
Figure 2 – Error in fault surface that occurs due to projecting points onto XY plane and triangulating in XY plane.

In PetraSim 2015.1 we have generalized this approach so that we find a plane that “best” fits the cloud of 3D points and project the points onto that plane. We then use that projection to triangulate the points and obtain the connectivity of the surface.

Given a plane:

fault plane eqn

defined by a point,

fault point eqn

and a vector normal to the plane,

fault normal eqn

The distance of any point p from the plane is:

fault distance eqn

We measure the total error as the sum of the squares of the distances for all points:

fault error eqn

Now the goal is to find a plane that minimizes this error and project our points onto the plane. I spent a fair amount of time deriving the equations to minimize this error with respect to a, b, and c and implemented it using Newton-Raphson iteration. The solution came back telling me to scale the vector v (for example use 2*a, 2*b, and 2*c). In retrospect this is obvious, since the error does not change with respect to the magnitude of v. Instead, we need to minimize the error with respect to the Euler angles defining v.

OK. I started the new derivation, went away to lunch and when I came back one of our programmers had implemented an algorithm where he randomly selected three points from the cloud, used those points to define a plane, and then measured the error from that plane. He then repeated this trial for a specified number of times and had the result shown below:

Figure 3 - A "best fit" of a plane to the points obtained by repeated trials.
Figure 3 – A “best fit” of a plane to the points obtained by repeated trials.

This trial and error approach is more reliable than setting up a nonlinear iteration scheme and is good enough, since we are just looking for a plane to project points onto. This is the approach we implemented. The surface that results from this projection is shown below.

Figure 4 - The new surface obtained using a projection onto the "best fit" plane and triangulating in that plane.
Figure 4 – The new surface obtained using a projection onto the “best fit” plane and triangulating in that plane.

Obviously, there are cases where this new approach can still fail, but it appears to work much more reliably. Feedback on any problems or suggested improvements to the algorithm would be appreciated.