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.
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:
defined by a point,
and a vector normal to the plane,
The distance of any point p from the plane is:
We measure the total error as the sum of the squares of the distances for all points:
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:
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.
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.