import numpy as np
import math
# ############################################################
# unutbu's model() (https://stackoverflow.com/a/35118683/27358)
def model(params, xyz):
a, b, c, d = params
x, y, z = xyz
length = (a ** 2 + b ** 2 + c ** 2)
return ((a * x + b * y + c * z + d) ** 2 / length).sum()
# ############################################################
# naive implementation
def z_from(x, y, abcd):
"""Returns the z coordinate from a pair of (x, y) coordinates
and a list of coefficients in the form ax + by + cz + d = 0.
"""
a, b, c, d = abcd
return (a * x + b * y + d) / -c
# see https://m...content-available-to-author-only...t.org/distance_point_plane
def unit_normal_vector(abcd):
a, b, c, d = abcd
n = np.array([a, b, c])
len_n = (a ** 2 + b ** 2 + c ** 2) ** .5
return n / len_n
# see https://m...content-available-to-author-only...t.org/distance_point_plane
def dist_point_plane(point, abcd):
x0, y0, z0 = np.array([0, 0, z_from(0, 0, abcd)])
x1, y1, z1 = point
n = unit_normal_vector(abcd)
v = np.array([x1 - x0, y1 - y0, z1 - z0])
return d
def total_dist_plane_points(abcd, xyz):
total = 0
for i, x in enumerate(xyz[0]):
y, z = xyz[1][i], xyz[2][i]
point = [x, y, z]
total = total + dist_point_plane(point, abcd)**2
return total
# ############################################################
# Main program
# points (2, 4, 6), (8, 10, 12), (14, 16, 18)
xyz = np.array([
[2, 8, 14],
[4, 10, 16],
[6, 12, 18]
])
# x + 2y + 3z + 4 = 0
abcd = np.array([1, 2, 3, 4])
expected = total_dist_plane_points(abcd, xyz)
actual = model(abcd, xyz)
print("Expected: {}".format(expected))
print("Actual: {}".format(actual))
aW1wb3J0IG51bXB5IGFzIG5wCmltcG9ydCBtYXRoCgoKIyAjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKIyB1bnV0YnUncyBtb2RlbCgpIChodHRwczovL3N0YWNrb3ZlcmZsb3cuY29tL2EvMzUxMTg2ODMvMjczNTgpCgpkZWYgbW9kZWwocGFyYW1zLCB4eXopOgogICAgYSwgYiwgYywgZCA9IHBhcmFtcwogICAgeCwgeSwgeiA9IHh5egogICAgbGVuZ3RoID0gKGEgKiogMiArIGIgKiogMiArIGMgKiogMikKICAgIHJldHVybiAoKGEgKiB4ICsgYiAqIHkgKyBjICogeiArIGQpICoqIDIgLyBsZW5ndGgpLnN1bSgpCgoKIyAjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKIyBuYWl2ZSBpbXBsZW1lbnRhdGlvbgoKZGVmIHpfZnJvbSh4LCB5LCBhYmNkKToKICAgICIiIlJldHVybnMgdGhlIHogY29vcmRpbmF0ZSBmcm9tIGEgcGFpciBvZiAoeCwgeSkgY29vcmRpbmF0ZXMKICAgIGFuZCBhIGxpc3Qgb2YgY29lZmZpY2llbnRzIGluIHRoZSBmb3JtIGF4ICsgYnkgKyBjeiArIGQgPSAwLgogICAgIiIiCiAgICBhLCBiLCBjLCBkID0gYWJjZAogICAgcmV0dXJuIChhICogeCArIGIgKiB5ICsgZCkgLyAtYwoKCiMgc2VlIGh0dHBzOi8vbS4uLmNvbnRlbnQtYXZhaWxhYmxlLXRvLWF1dGhvci1vbmx5Li4udC5vcmcvZGlzdGFuY2VfcG9pbnRfcGxhbmUKZGVmIHVuaXRfbm9ybWFsX3ZlY3RvcihhYmNkKToKICAgIGEsIGIsIGMsIGQgPSBhYmNkCiAgICBuID0gbnAuYXJyYXkoW2EsIGIsIGNdKQogICAgbGVuX24gPSAoYSAqKiAyICsgYiAqKiAyICsgYyAqKiAyKSAqKiAuNQogICAgcmV0dXJuIG4gLyBsZW5fbgoKCiMgc2VlIGh0dHBzOi8vbS4uLmNvbnRlbnQtYXZhaWxhYmxlLXRvLWF1dGhvci1vbmx5Li4udC5vcmcvZGlzdGFuY2VfcG9pbnRfcGxhbmUKZGVmIGRpc3RfcG9pbnRfcGxhbmUocG9pbnQsIGFiY2QpOgogICAgeDAsIHkwLCB6MCA9IG5wLmFycmF5KFswLCAwLCB6X2Zyb20oMCwgMCwgYWJjZCldKQogICAgeDEsIHkxLCB6MSA9IHBvaW50CiAgICBuID0gdW5pdF9ub3JtYWxfdmVjdG9yKGFiY2QpCiAgICB2ID0gbnAuYXJyYXkoW3gxIC0geDAsIHkxIC0geTAsIHoxIC0gejBdKQogICAgZCA9IG5wLmFicyhucC5kb3QodiwgbikpCiAgICByZXR1cm4gZAoKCmRlZiB0b3RhbF9kaXN0X3BsYW5lX3BvaW50cyhhYmNkLCB4eXopOgogICAgdG90YWwgPSAwCiAgICBmb3IgaSwgeCBpbiBlbnVtZXJhdGUoeHl6WzBdKToKICAgICAgICB5LCB6ID0geHl6WzFdW2ldLCB4eXpbMl1baV0KICAgICAgICBwb2ludCA9IFt4LCB5LCB6XQogICAgICAgIHRvdGFsID0gdG90YWwgKyBkaXN0X3BvaW50X3BsYW5lKHBvaW50LCBhYmNkKSoqMgogICAgcmV0dXJuIHRvdGFsCgoKIyAjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKIyBNYWluIHByb2dyYW0KCiMgcG9pbnRzICgyLCA0LCA2KSwgKDgsIDEwLCAxMiksICgxNCwgMTYsIDE4KQp4eXogPSBucC5hcnJheShbCiAgICBbMiwgOCwgMTRdLAogICAgWzQsIDEwLCAxNl0sCiAgICBbNiwgMTIsIDE4XQpdKQojIHggKyAyeSArIDN6ICsgNCA9IDAKYWJjZCA9IG5wLmFycmF5KFsxLCAyLCAzLCA0XSkKCmV4cGVjdGVkID0gdG90YWxfZGlzdF9wbGFuZV9wb2ludHMoYWJjZCwgeHl6KQphY3R1YWwgPSBtb2RlbChhYmNkLCB4eXopCgpwcmludCgiRXhwZWN0ZWQ6IHt9Ii5mb3JtYXQoZXhwZWN0ZWQpKQpwcmludCgiQWN0dWFsOiAgIHt9Ii5mb3JtYXQoYWN0dWFsKSkK