title: plannar fit with odrpack
body: <p>I am trying to use scipy.odr to get a best fit plane for some x, y, z points.</p>
  
  <p>I define the plane equation implicitly as ax + by + cz + d = 0 and I perform a least squares (with scipy.linalg.lstsq) to provide the odr with an initial estimation.</p>
  
  <p>The components of the beta vector (where beta = [a, b, c, d]) returned by the odr are of a magnitude between 1e167 and 1e172... Is such a result trustworthy? I find the numbers to be absurd...</p>
  
  <p>Note that the points come from 3D scanning of a relatively flat face which is almost parallel to the xz plane (nearly vertical).</p>
  
  <p>Here is the pprint() of the odr result:</p>
  
  <pre><code>'
  Beta: [  3.14570111e-170   3.21821458e-169   4.49232028e-172   4.49374557e-167]
  Beta Std Error: [ 0.  0.  0.  0.]
  Beta Covariance: [[  6.37459471e-10  -8.57690019e-09  -2.18092934e-11  -1.13009384e-06]
   [ -8.57690019e-09   5.11732570e-07   1.30123070e-09   6.74263262e-05]
   [ -2.18092934e-11   1.30123070e-09   5.22674068e-12   1.70799469e-07]
   [ -1.13009384e-06   6.74263262e-05   1.70799469e-07   8.88444676e-03]]
  Residual Variance: 0.0
  Inverse Condition #: 0.0010484041422201213
  Reason(s) for Halting:
    Sum of squares convergence
  None
  '
  </code></pre>
  
  <p>The code I am using :</p>
  
  <pre><code>import numpy as np
  import scipy.linalg
  from scipy import odr
  import pickle
  
  def planar_fit(points):
      # best-fit linear plane
      a = np.c_[points[:, 0], points[:, 1], np.ones(points.shape[0])]
      c, _, _, _ = scipy.linalg.lstsq(a, points[:, 2])  # coefficients
      # The coefficients are returned as an array beta=[a, b, c, d] from the implicit form 'a*x + b*y + c*z + d = 0'.
      beta = np.r_[c[0], c[1], -1, c[2]] / c[2]
      return beta
  
  
  def odr_planar_fit(points):
      def f_3(beta, xyz):
          """ implicit definition of the plane"""
          return beta[0] * xyz[0] + beta[1] * xyz[1] + beta[2] * xyz[2] + beta[3]
  
      # # Coordinates of the 2D points
      x = points[:, 0]
      y = points[:, 1]
      z = points[:, 2]
  
      # Use least squares for initial estimate.
      beta0 = planar_fit(points)
  
      # Create the data object for the odr. The equation is given in the implicit form 'a*x + b*y + c*z + d = 0' and
      # beta=[a, b, c, d] (beta is the vector to be fitted). The positional argument y=1 means that the dimensionality
      # of the fitting is 1.
      lsc_data = odr.Data(np.row_stack([x, y, z]), y=1)
      # Create the odr model
      lsc_model = odr.Model(f_3, implicit=True)
      # Create the odr object based on the data, the model and the first estimation vector.
      lsc_odr = odr.ODR(lsc_data, lsc_model, beta0)
      # run the regression.
      lsc_out = lsc_odr.run()
  
      return lsc_out, beta0
  
  def main():
      #import from pickle.
      with open('./points.pkl', 'rb') as f:
          points = np.array(pickle.load(f))
  
      # Perform the ODR
      odr_out, lstsq = odr_planar_fit(points)
      print(lstsq)
      print(odr_out.pprint())
  
  main()
  </code></pre>
  
  <p><a href="https://drive.google.com/file/d/1_sYm0kzBvYPZ3AFMO8aMOPqMYbAmlXg9/view?usp=sharing" rel="nofollow noreferrer">The pickle containing my points.</a></p>