Plane fitting to 4 (or more) XYZ points

Question

I have 4 points, which are very near to be at the one plane - it is the 1,4-Dihydropyridine cycle.

I need to calculate distance from C3 and N1 to the plane, which is made of C1-C2-C4-C5. Calculating distance is OK, but fitting plane is quite difficult to me.

1,4-DHP cycle:

1,4-DHP cycle, another view:

from array import *
from numpy import *
from scipy import *

# coordinates (XYZ) of C1, C2, C4 and C5
x = [0.274791784, -1.001679346, -1.851320839, 0.365840754]
y = [-1.155674199, -1.215133985, 0.053119249, 1.162878076]
z = [1.216239624, 0.764265677, 0.956099579, 1.198231236]

# plane equation Ax + By + Cz = D
# non-fitted plane
abcd = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

# creating distance variable
distance =  zeros(4, float)

# calculating distance from point to plane
for i in range(4):
    distance[i] = (x[i]*abcd[0]+y[i]*abcd[1]+z[i]*abcd[2]+abcd[3])/sqrt(abcd[0]**2 + abcd[1]**2 + abcd[2]**2)
    
print distance

# calculating squares
squares = distance**2

print squares

How to make sum(squares) minimized? I have tried least squares, but it is too had for me.

Answer 1

The fact that you are fitting to a plane is only slightly relevant here. What you are trying to do is minimize a particular function starting from a guess. For that use scipy.optimize. Note that there is no guarantee that this is the globally optimal solution, only locally optimal. A different initial condition may converge to a different result, this works well if you start close to the local minima you are seeking.

I've taken the liberty to clean up your code by taking advantage of numpy's broadcasting:

import numpy as np

# coordinates (XYZ) of C1, C2, C4 and C5
XYZ = np.array([
        [0.274791784, -1.001679346, -1.851320839, 0.365840754],
        [-1.155674199, -1.215133985, 0.053119249, 1.162878076],
        [1.216239624, 0.764265677, 0.956099579, 1.198231236]])

# Inital guess of the plane
p0 = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

def f_min(X,p):
    plane_xyz = p[0:3]
    distance = (plane_xyz*X.T).sum(axis=1) + p[3]
    return distance / np.linalg.norm(plane_xyz)

def residuals(params, signal, X):
    return f_min(X, params)

from scipy.optimize import leastsq
sol = leastsq(residuals, p0, args=(None, XYZ))[0]

print("Solution: ", sol)
print("Old Error: ", (f_min(XYZ, p0)**2).sum())
print("New Error: ", (f_min(XYZ, sol)**2).sum())

This gives:

Solution:  [  14.74286241    5.84070802 -101.4155017   114.6745077 ]
Old Error:  0.441513295404
New Error:  0.0453564286112