Fitting planes to Point Clouds

Computing the best fit plane for point clouds using principal component analysis (PCA).

vector centroid = detail ( 0, "centroid" );
vector p = v@P - centroid;
setdetailattrib ( geoself ( ), "sumxx", p.x * p.x, "add" );
setdetailattrib ( geoself ( ), "sumxy", p.x * p.y, "add" );
setdetailattrib ( geoself ( ), "sumxz", p.x * p.z, "add" );
setdetailattrib ( geoself ( ), "sumyy", p.y * p.y, "add" );
setdetailattrib ( geoself ( ), "sumyz", p.y * p.z, "add" );
setdetailattrib ( geoself ( ), "sumzz", p.z * p.z, "add" );

matrix3 m = set ( @opinput1_sumxx, @opinput1_sumxy, @opinput1_sumxz, @opinput1_sumxy, @opinput1_sumyy, @opinput1_sumyz, @opinput1_sumxz, @opinput1_sumyz, @opinput1_sumzz );

matrix3 u = 0;
vector s = 0;
matrix3 v = 0;

svddecomp ( m, u, s, v );
v = transpose ( v );

vector x = set ( v.xx, v.xy, v.xz );
vector y = set ( v.yx, v.yy, v.yz );
vector z = set ( v.zx, v.zy, v.zz );

vector centroid = v@opinput1_centroid;
int pt = addpoint ( geoself ( ), centroid );

setpointattrib ( geoself ( ), "x", pt, x );
setpointattrib ( geoself ( ), "y", pt, y );
setpointattrib ( geoself ( ), "z", pt, z );

setpointattrib ( geoself ( ), "N", pt, z );