We can start with z-coordinate fit from Henrik Schachner’s answer:
In[]:=
pts={{0,0,0},{1,0,1},{1,1,2},{2,2,2},{2,4,3},{3,6,4}};model=ax+by+c;params=FindFit[pts,model,{a,b,c},{x,y}];
We can construct an infinite plane from it, for instance by sampling three random points on the plane:
In[]:=
plane1=InfinitePlane[{x,y,z}/.N@FindInstance[z==model/.params,{x,y,z},Reals,3]]
Out[]=
InfinitePlane[{{-13.4,-5.5,-10.2741},{-3.8,-10.6,-5.22532},{12.4,-13.1,5.03354}}]
I propose another fitting function, which attempts to minimize true distance to the plane. Common distance function to an implicit plane whose normal {a,b,c} is forced to be unit length is used, least squares fit computed with Minimize, and InfinitePlane object created by sampling three random points:
In[]:=
plane2=InfinitePlane[{x,y,z}/.N@FindInstance[{a,b,c,d}.{x,y,z,1}==0/.Last@Minimize[{Sum[({a,b,c}.p+d)^2,{p,pts}],{a,b,c}∈Sphere[]&&d>=0},{a,b,c,d}],{x,y,z},Reals,3]]
Out[]=
InfinitePlane[{{-13.4,-5.5,-15.8885},{-3.8,-10.6,-5.34056},{12.4,-13.1,13.0537}}]
We can see that the norm of distances measured plainly on the z direction are smaller for plane1. This is computed by solving the z value of the corresponding plane at {x,y} coordinates of each point.
In[]:=
Norm[Abs[#3-z]/.First@Solve[RegionMember[plane1,{x,y,z}]&&{x,y}=={#1,#2},{x,y,z}]&@@@pts]
Out[]=
0.838172
In[]:=
Norm[Abs[#3-z]/.First@Solve[RegionMember[plane2,{x,y,z}]&&{x,y}=={#1,#2},{x,y,z}]&@@@pts]
Out[]=
0.952435
If we compute the norm of the Euclidean distances of the points from each plane, we see that it is actually smaller for plane2:
In[]:=
Norm[RegionDistance[plane1,pts]]
Out[]=
0.675476
In[]:=
Norm[RegionDistance[plane2,pts]]
Out[]=
0.623551
All of these objects can be rendered directly in Graphics3D, but the planes are a bit... featureless. We can see that they are clearly different, though:
In[]:=
Graphics3D[{Red,Sphere[pts,.2],Opacity[0.5],Green,plane1,Blue,plane2}]
Out[]=
With some symbolic gymnastics they can be plotted the same way as in the earlier answer:
In[]:=
Show[Plot3D[Evaluate[First[SolveValues[RegionMember[#,{x,y,z}],z]]&/@{plane1,plane2}],{x,-1,4},{y,-1,7},PlotStyle->Opacity[0.5],BoxRatios->Automatic],Graphics3D[{Red,Sphere[pts,0.1]}],ViewPoint->{-1,-3,2}]
Out[]=
We can also visualise distances to corresponding planes by finding the closest point for each point and drawing a Tube between the point and the plane:
In[]:=
Show[Plot3D[Evaluate[First[SolveValues[RegionMember[plane1,{x,y,z}],z]]],{x,-1,4},{y,-1,7},PlotStyle->Opacity[0.5],BoxRatios->Automatic],Graphics3D[{Red,Sphere[pts,0.1],Tube[{#,RegionNearest[plane1,#]}&/@pts]}],ViewPoint->{1,-3,2}]
Out[]=
In[]:=
Show[Plot3D[Evaluate[First[SolveValues[RegionMember[plane2,{x,y,z}],z]]],{x,-1,4},{y,-1,7},PlotStyle->Opacity[0.5],BoxRatios->Automatic],Graphics3D[{Red,Sphere[pts,0.1],Tube[{#,RegionNearest[plane2,#]}&/@pts]}],ViewPoint->{1,-3,2}]
Out[]=
One thing neither of these fits handle are outliers. RegionFit with ”Plane” argument may be worth checking out if that is necessary.