Best-fit-plane
through a point cloud

by Vince Cronin, after Paláncz and others (2013)

Revised 24 November 2017, based on earlier working drafts.

Input

Start with a clean slate by clearing the input variables

Clear[XYZ,ℕ]

The file "Composite-Loc-Data-V4" contains all of the aftershock locations published by Greensfelder (1968) and Ryall et al. (1968) with the following exceptions :  Greensfelder’s events 32, 70, and 91, because their respective epicenters were outside the geographic box lat 39.370 to 39.544, long -120.211 to -120.037.  In cases where the same event was reported by both sources, the mean of the locations was used in preference to double-reporting the event or arbitrarily choosing a preferred source.

;

The data is in an Excel spreadsheet in which column A is the latitude and column B is the longitude -- both in decimal degrees where negative latitude is south and negative longitude is west.  Column C is the focal depth in kilometers.  Column D is an event ID.

Import the Point Cloud Data File

To modify this code so that it can be used for your dataset, do the following.

1.  Enter the following as an input line (under the black box, below):

pointCloudDataIn = Import[];

2. Put your cursor between the square brackets in the Import[] statement and click to establish the insertion point we will need in the next step.

3. Go to the Insert menu, select File Path, navigate to the correct input data file and choose it, and the correct file path will be inserted at the cursor.  In this example, the Excel file is located on the desktop of Vince’s computer, and so Mathematica will insert the path "/Users/vince/Desktop/Composite-Loc-Data-V4.xls" between the brackets.  The path will be different for every different file, and on every different computer.

pointCloudDataIn = Import[“/Users/vince/Desktop/Composite-Loc-Data-V4.xls”];

You can make sure it is an input line by clicking on that line’s bracket on the far right edge of this window, going to the Format menu, selecting the Style submenu, and then choosing Input.

pointCloudDataIn=Import["/Users/vincecronin/Desktop/Composite-Loc-Data-V4.xls"];

This won' t work unless you follow directions and add the appropriate input line above this line.  Your input line should look something like the example:
pointCloudDataIn = Import[“/Users/vincecronin/Desktop/Composite-Loc-Data-V4.xls”];

The imported dataset will have to be modified so that it has the correct dimensions for this analysis, which is accomplished with the following code that defines the list pointCloudData:

pointCloudData=Flatten[pointCloudDataIn,1];

noRecords=Length[pointCloudData]

The origin of the coordinate system used in this analysis is just beyond the south and west corner of the dataset.  That is, the coordinate system origin latitude (csOriginLat) is smaller than the minimum latitude of any point in the point cloud, and the coordinate system origin longitude (csOriginLong) is smaller than the minimum longitude of any point in the point cloud.

Variable minLat is the minimum latitude value (associated with the farthest south point) in the point cloud.

minLat=Min[Table[pointCloudData[[i,1]],{i,1,noRecords}]]

Variable minLong is the maximum latitude value (associated with the farthest west point) in the point cloud.

minLong=Min[Table[pointCloudData[[i,2]],{i,1,noRecords}]]

csOriginLat=39.36;

csOriginLong=-120.22;

User-Defined Functions

unitVect3D[vect_]:={(vect[[1]]/Norm[vect]),(vect[[2]]/Norm[vect]),(vect[[3]]/Norm[vect])};

The function PlanePCA is adapted from Appendix A - 1 of Paláncz and others (2013)

PlanePCA[XYZ_]:=Module[{𝒩,co,Cco,nx,ny,nz,𝒹,sol𝒹},
𝒩=Length[XYZ];
co=Mean[XYZ];
Cco=Map[#-co&,XYZ];
{nx,ny,nz}=Eigenvectors[Covariance[Cco]][[3]];
sol𝒹=Solve [{nx,ny,nz}.Total[XYZ]-𝒹 𝒩==0,𝒹]//Flatten;
{-(nx/nz),-(ny/nz),(𝒹/nz)}/.sol𝒹];

Computation

Variable maxLat is the maximum latitude value (associated with the farthest north point) in the point cloud.

maxLat=Max[Table[pointCloudData[[i,1]],{i,1,noRecords}]]

Variable midLat is the latitude that is half-way between maxLat and minLat.

midLat=(maxLat+minLat)/2

Variable maxLong is the maximum latitude value (associated with the farthest east point) in the point cloud.

maxLong=Max[Table[pointCloudData[[i,2]],{i,1,noRecords}]]

Variable midLong is the latitude that is half-way between maxLong and minLong.

midLong=(maxLong+minLong)/2

Variable deg2kmLat is the number of kilometers per degree of latitude, assuming a mean Earth radius of 6371 km.

deg2kmLat=(6371*2*π)/360;

N[%]

Variable deg2kmLongMidLat is the number of kilometers per degree of longitude, computed at the middle latitude of the point cloud.

deg2kmLongMidLat=deg2kmLat*Cos[midLat Degree]

The vectors computed in the following expression are in km units.  The origin {0,0,0} is at {csOriginLong, csOriginLat, mean sea level}, where the difference between csOrginLat and eqData[[i,1]] is multiplied by deg2kmLat or around 111.195 km/degree to compute the y coordinate for any given vector, and the difference between csOriginLong and eqData[[i,2]] is multiplied by deg2kmLongMidLat or around 85.89 km/degree to compute the x coordinate for any given vector.  The result is a rectangular horizontal mapping of the epicenters.  The sign of the input focal depths is reversed so that negative depths are below sea level.

XYZ=Table[{(((pointCloudData[[i,2]]-csOriginLong)*deg2kmLongMidLat)),((pointCloudData[[i,1]]-csOriginLat)*deg2kmLat),(-1)*pointCloudData[[i,3]]},{i,1,noRecords}];

maxX=Max[Table[XYZ[[i,1]],{i,1,noRecords}]]

minX=Min[Table[XYZ[[i,1]],{i,1,noRecords}]]

midX=(maxX+minX)/2

maxY=Max[Table[XYZ[[i,2]],{i,1,noRecords}]]

minY=Min[Table[XYZ[[i,2]],{i,1,noRecords}]]

midY=(maxY+minY)/2

maxZ=Max[Table[XYZ[[i,3]],{i,1,noRecords}]]

minZ=Min[Table[XYZ[[i,3]],{i,1,noRecords}]]

ℕ=Length[XYZ]

north={0,1,0};

Where is the Centroid of the Point Cloud?

meanVector=Mean[XYZ]

latitude=csOriginLat+(meanVector[[2]]/deg2kmLat)

longitude=csOriginLong+(meanVector[[1]]/deg2kmLongMidLat)

depth=meanVector[[3]]

Check

meanVector

recomputedMeanVector={((longitude-csOriginLong)*deg2kmLongMidLat),((latitude-csOriginLat)*deg2kmLat),depth}

End Check

Solution via SVD

after Paláncz and others (2013)

c0=Mean[XYZ]

Cc0=Map[#-c0&,XYZ];

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];]

ν=Last[Transpose[v]]

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify

NSolve[%,𝓏]//Expand

p1=Plot3D[𝓏/.%[[1]],{𝓍,minX,maxX},{𝓎,minY,maxY},BoxRatios->{1,1,0.5}, ClippingStyle->None,Mesh->None,PlotRange->{minZ,maxZ}];

The equation of a plane can be written as z = γ - αx + βy, so the result above provides the value of α, β, and γ.

Solution via PCA

after Paláncz and others (2013)

solnPCA=PlanePCA[XYZ]

The equation of a plane can be written as z = γ - αx + βy, so the result above provides the value of α, β, and γ.

α=solnPCA[[1]]

β=solnPCA[[2]]

γ=solnPCA[[3]]

The value of meanX is equal to (γ + (β*meanY) - meanZ)/(-α).

The value of meanY is equal to (γ + (α*meanX) - meanZ)/(-β).

The value of meanZ is equal to γ + (α*meanX) + (β*meanY).

Find the Vectors Normal to the Plane

We already know the coordinates of one point on the plane -- the mean point in the XYZ location - vector dataset (meanVector).  We need to find the coordinates of two other points on the plane in order to find the vectors normal to the plane, which will allow us to find the strike, dip, and dip angle of the plane.

If we know the x and z coordinates of a point on the plane, we can use the equation of the plane to calculate the y coordinate.  We use the z coordinate of the mean vector (meanZ), select an x coordinate that is between xMax and xMin (midX), and compute the corresponding y value (trialY).  The result is another point that is known to be on the plane (point1) and is different from the meanVector.

meanZ=meanVector[[3]];

trialY=(γ+(α*midX)-meanZ)/(-β)

point1={midX,trialY,meanZ}

This time, we use midX as defined above, select a z coordinate (lowZ) that is below the z coordinate of the mean vector (meanZ), and compute the corresponding y value (lowY).  The result is another point that is known to be on the plane (point2) and is different from the meanVector.

lowZ=meanZ-(0.5(meanZ-minZ))

lowY=(γ+(α*midX)-lowZ)/(-β)

point2={midX,lowY,lowZ}

We now know the location vectors to three points that are on the plane -- point1, point2, and meanVector.  From these location vectors, we can determine two vectors parallel to the plane -- vectorOnPlane1 and vectorOnPlane2.

vectorOnPlane1=point1-meanVector

vectorOnPlane2=point2-meanVector

The vectors vectorOnPlane1 and vectorOnPlane2 are parallel to the plane but not to each other, so we can use their cross product and the user-defined function unitVect3D to determine unit vectors normal to the plane.

normal1=unitVect3D[Cross[vectorOnPlane1,vectorOnPlane2]]

normal2=-normal1

normal=If[normal1[[3]]≥0,normal1,normal2]

The dip angle of the best-fit plane (dipAngle) is measured in degrees and is the same as the angle between a vector pointing up along the vertical (Z) axis (vertUpVect) and the vector normal to the plane (normal) that has a positive or zero Z coordinate.

vertUpVect={0,0,1};

dipAngle=VectorAngle[vertUpVect,normal]*(180/π)

The trend of the dip vector (dipDirection) is also given in degrees of azimuth, measured clockwise relative to north.  To find the dip direction, we start by projecting the normal vector onto a horizontal plane (projectedVector) and then finding the unit vector of the resulting projection (unitProjVector).

projectedVector={normal[[1]],normal[[2]],0}

unitProjVector=unitVect3D[projectedVector]

dipDirection=If[unitProjVector[[1]]≥0,VectorAngle[north,unitProjVector]*(180/π),360-(VectorAngle[north,unitProjVector]*(180/π))]

direction=If[dipDirection≤15,"north",If[dipDirection≤75,"northeast",If[dipDirection≤105,"east",If[dipDirection≤165,"southeast",If[dipDirection≤195,"south",If[dipDirection≤255,"southwest",If[dipDirection≤285,"west",If[dipDirection≤345,"northwest","north"]]]]]]]]

The right-hand-rule strike azimuth (rhrStrike) is 90° less than the dip direction.

rhrStrike=If[dipDirection<90,(270+dipDirection),(dipDirection-90)]

Check

strikeVector={Sin[rhrStrike Degree],Cos[rhrStrike Degree],0}

Graphics

p2=ListPointPlot3D[XYZ,PlotStyle->Blue ];

Show[{p1,p2}]

Show[{p1,p2},ViewPoint->Above]

To use the ViewVector option, we need to specify the location vector to a point from which we will view the aftershock foci and best-fit plane (viewOrigin) and the vector of the point towards which we are looking (meanVector).  The viewOrigin is located a given distance (centroid2ViewOrigin) from the centroid of the point cloud, along strike of the best-fit plane in the opposite direction from the rhrStrike.   So the final view direction is toward the rhrStrike, and the point cloud is projected onto a plane that is perpendicular to the best-fit plane.

centroid2ViewOrigin=10*(Sqrt[(maxX^2)+(maxY^2)]);

viewOrigin={(-centroid2ViewOrigin*strikeVector[[1]]),(-centroid2ViewOrigin*strikeVector[[2]]),meanVector[[3]]}+meanVector;

Show[{p1,p2},ViewVector->{viewOrigin,strikeVector}]

Summary

The best-fit plane through the input list of aftershock focal locations has a right-hand-rule strike of

rhrStrike

degrees and a dip angle of

dipAngle

degrees toward

dipDirection

degrees (dip azimuth), or toward the

direction

That plane passes through a point at latitude

latitude

degrees, longitude

longitude

degrees, at a depth of

depth

km (negative depth value indicates an elevation below sea level).

References

Chen C C and Stamos I (2007) Range image segmentation for modeling and object detection in urban scenes, 3 DIM2007 .1

DalleMole V, do Rego R and Araújo A (2010) The Self - Organizing Approach for Surface Reconstruction from Unstructured Point Clouds,  in Matsopoulos G.K.(Ed.) Self - Organizing Map, In Tech, 

Diebel  J R, Thrun S and Brunig M (2006) A Bayesian method for probable surface reconstruction and decimation, ACM Transactions on Graphics (TOG). 1

Fischler M A and Bolles R C (1981) Random sample consensus : A paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, 24, pp .381 - 395

Huang C - M and Tseng Y - H. Plane fitting methods of Lidar point cloud, Dept. of Geomatics, National Cheng Kung Uni. Taiwan, tseng@mail.ncku.edu.tw

Hubert M , Rousseeuw P J and Van den Branden K (2005) ROBPCA : A new approach to robust principal component analysis, TECHNOMETRICS, Vol. 47. No. 1. pp .64 - 79.

Kohonen T (1998) The self - organizing map, Neurocomputing, Vol. 21, pp .1 - 6.

Lakaemper R and Latecki L J (2006) Extended EM for Planar Approximation  of 3 D Data, IEEE Int. Conf. on Robotics and Automation (ICRa), Orlando, Florida, May 2006.

Lukács G, Martin R and Marshall D (1998) Faithful Leat - Squares Fitting of Spheres, cylinders, Cones and Tori for Reliable Segmentation, in Burkhardt H and Neumann B (Eds.) Computer Vision, ECCV ' 98, Vol I, LNCS 1406, pp .671 - 686, Springer - Verlag Berlin Heidelberg.

Mitra N J and Nguyen (2003) SoCG' 03, June 8 - 10, San Diego, California, USA, ACM 1 - 58113 - 663 - 3/03/0006, pp. 322 - 328.

Nievergelt Y (2000)  A tutorial history of least square with applications to astronomy and geodesy, J. Comp. and Appl. Mathematics, 121, pp. 37 - 72.

Nurunnabi A, Belton D and West G. (2012) Diagnostic - Robust Statistical Analysis for Local Surface Fitting in 3 D Point Cloud Data, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol.I - 3, 2012 XXII ISPRS Congress, 25 Aug. 2012, Melbourne, Australia, pp .269 - 275

Paláncz B (2014) Fitting Data with Different Error Models, The Mathematica Journal (accepted for publication, to be published in April)

Paláncz, B, Somogyi, A, Lovas, T, and Molnár, B (2013) Plane Fitting to Point Cloud via Gröbner Basis, available via http://library.wolfram.com/infocenter/MathSource/8769/#downloads

Poppinga J, Vaskevicius N, Birk A and Pathak K (2006) Fast Plane Detection and Polygonalization in noisy 3 D Range Images, Int. Conf. on Intelligent Robots and Systems (IROS), Nice, France, IEEE Press 2008

Rose C and Smith D (2000) Symbolic Maximum Likelihood Estimation with Mathematica, The Statistician 49, pp .229 - 240

Russeeuw P J and Van Driessen K (1999)  A Fast Algorithm for the Minimum Covariance Determinant Estimator, TECHNOMETRICS, Vol. 41. No. 3. pp. 212 - 223.

Stathas D, Arabatzi O, Dogouris S, Piniotis G, Tsini D and Tsinis D (2003) New Monitoring Techniques on the Determination of Structure Deformation, Proc. 11 th FIG Symp. on Deformation Measurements, Santorini, Greece. 

Zuliani M (2012) RANSAC for Dummies, vision.ece.ucsb.edu/~zuliani

Yang M Y and Förtsner W (2010) Plane Detection in Point Cloud Data, TR - IGG - P - 2010 - 01, Techn.Report. Nr .1, Dept. of Photogrammetry Inst. of Geodesy and Geo - information, Uni., Bonn, Germany.

Yaniv Z (2010) Random Sample Consensus (RANSAC) Algorithm, A Generic Implementation, Georgetown University Medical Center, Washington, DC, USA, zivy@isis.georgetown.edu