Last active
January 13, 2019 19:29
-
-
Save chronodm/f8f8af5d2764b2a89375d5bde22c0839 to your computer and use it in GitHub Desktop.
Naive total perpendicular distance from points to plane
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 = math.sqrt(a ** 2 + b ** 2 + c ** 2) | |
return ((np.abs(a * x + b * y + c * z + d)) ** 2).sum() / length | |
# ############################################################ | |
# 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://mathinsight.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://mathinsight.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]) | |
d = np.abs(np.dot(v, n)) | |
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) | |
return total | |
# ############################################################ | |
# Main program | |
# points (2, 4, 6), (8, 10, 12), (14, 16, 18) | |
xyz = [ | |
[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(f"Expected: {expected}") | |
print(f"Actual: {actual}") | |
# => | |
# Expected: 54.52129335013458 | |
# Actual: 1160.9828348675715 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment