Search code examples
pythonscipycurve-fitting

How to fit a curve (3d), where points have specified normal in python?

I have following three points in 3d space with their respective normal vectors.

A, B, C are positions and N_A, N_B,and N_C are their respective normal vectors.

A = np.array([  348.92065834, -1402.3305998,     32.69313966])
N_A = np.array([-0.86925426,  0.02836434, -0.49355091])

B =  np.array([282.19332067,  82.52027998,  -5.92595371])
N_B = np.array([-0.82339849,  0.43041935,  0.3698028])

C = np.array([247.37475615,  -3.70129865, -22.10494737])
N_C = np.array([-0.83989222, 0.23796899, 0.48780305])

Three points are almost in one plane, but there is a slight directional change between the two closest point B and C. From there to the point A, I assumed there could be a curvature in X-Y coordinates as shown XY, but parabole in X-Z coordinates as shown in XZ

Supposed condition is cylinder can fit the three points if the curvature was much obvious. And considering their respective normal vectors from X and Z coordinates, B and C normal vectors face lower direction while the that of A faces upper direction. So, all in all it it could be a paraboloid. Question is how to fit them taking their normal vectors into account. If not possible, then how to fit them with some curvature from X and Y direction.

Here is a code for the plot

import numpy as np
import matplotlib.pyplot as plt
def set_axes_radius(ax, origin, radius):

    ax.set_xlim3d([origin[0] - radius, origin[0] + radius])
    ax.set_ylim3d([origin[1] - radius, origin[1] + radius])
    ax.set_zlim3d([origin[2] - radius, origin[2] + radius])


def set_axes_equal(ax, zoom=1.):
    '''
        Make axes of 3D plot have equal scale so that spheres appear as spheres,
        cubes as cubes, etc..  This is one possible solution to Matplotlib's
        ax.set_aspect("equal") and ax.axis("equal") not working for 3D.
        input:
          ax:   a matplotlib axis, e.g., as output from plt.gca().

    '''

    limits = np.array([
        ax.get_xlim3d(),
        ax.get_ylim3d(),
        ax.get_zlim3d(),
    ])

    origin = np.mean(limits, axis=1)
    radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0])) / zoom
    set_axes_radius(ax, origin, radius)
%matplotlib qt 


# positions and their respective normal vectors

A = np.array([  348.92065834, -1402.3305998,     32.69313966])
N_A = np.array([-0.86925426,  0.02836434, -0.49355091])

B =  np.array([282.19332067,  82.52027998,  -5.92595371])
N_B = np.array([-0.82339849,  0.43041935,  0.3698028])

C = np.array([247.37475615,  -3.70129865, -22.10494737])
N_C = np.array([-0.83989222, 0.23796899, 0.48780305])

# A plane is given by
# a*x + b*y + c*z + d = 0
# where (a, b, c) is the normal.
# If the point (x, y, z) lies on the plane, then solving for d yield:
# d = -(a*x + b*y + c*z)
d_A = -np.sum(N_A * A)
d_B = -np.sum(N_B * B)
d_C = -np.sum(N_C * C)

# Create a meshgrid:
delta = 200

xlim_A = A[0] - delta, A[0] + delta
ylim_A = A[1] - delta, A[1] + delta
xx_A, yy_A = np.meshgrid(np.arange(*xlim_A), np.arange(*ylim_A))

xlim_B = B[0] - delta, B[0] + delta
ylim_B = B[1] - delta, B[1] + delta
xx_B, yy_B = np.meshgrid(np.arange(*xlim_B), np.arange(*ylim_B))

xlim_C = C[0] - delta, C[0] + delta
ylim_C = C[1] - delta, C[1] + delta
xx_C, yy_C = np.meshgrid(np.arange(*xlim_C), np.arange(*ylim_C))

# Solving the equation above for z:
# z = -(a*x + b*y +d) / c
zz_A = -(N_A[0] * xx_A + N_A[1] * yy_A + d_A) / N_A[2]
zz_B = -(N_B[0] * xx_B + N_B[1] * yy_B + d_B) / N_B[2]
zz_C = -(N_C[0] * xx_C + N_C[1] * yy_C + d_C) / N_C[2]


fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(xx_A, yy_A, zz_A, alpha=0.5, color='g')
ax.plot_surface(xx_B, yy_B, zz_B, alpha=0.5, color='cyan')
ax.plot_surface(xx_C, yy_C, zz_C, alpha=0.5, color='crimson')






# Plot point.


x_A, y_A, z_A = A
x_B, y_B, z_B = B
x_C, y_C, z_C = C

Plane_A, = ax.plot(x_A, y_A, z_A, marker='o', markersize=5, color='g')
Plane_A.set_label('A position')
ax.legend()

Plane_B, = ax.plot(x_B, y_B, z_B, marker='o', markersize=5, color='cyan')
Plane_B.set_label('B position')
ax.legend()

Plane_C, = ax.plot(x_C, y_C, z_C, marker='o', markersize=5, color='crimson')
Plane_C.set_label('C position')
ax.legend()

# Plot normal.


dx_A, dy_A, dz_A = delta * N_A
ax.quiver(x_A, y_A, z_A, dx_A, dy_A, dz_A, arrow_length_ratio=0.15, linewidth=3, color='g')

dx_B, dy_B, dz_B = delta * N_B
ax.quiver(x_B, y_B, z_B, dx_B, dy_B, dz_B, arrow_length_ratio=0.15, linewidth=3, color='cyan')

dx_C, dy_C, dz_C = delta * N_C
ax.quiver(x_B, y_C, z_C, dx_C, dy_C, dz_C, arrow_length_ratio=0.15, linewidth=3, color='crimson')

# Enforce equal axis aspects so that the normal also appears to be normal.



ax.set_xlim(xmax=1500,xmin=-1500)
ax.set_ylim(ymax=400, ymin=-400)

zlim = max(A[2], B[2], C[2]) - delta, max(A[2], B[2], C[2]) + delta
ax.set_zlim(*zlim)

ax = plt.gca()

#ax.set_box_aspect([1,1,1])
set_axes_equal(ax)
ax.set_xlabel('X', fontsize=20)
ax.set_ylabel('Y', fontsize=20)
ax.set_zlabel('Z', fontsize=20)

plt.show()
Solution

  • Meta-problems: next time,

    • explain what you're actually doing (attempting regression on interplanetary electromagnetic bow shock from satellite readings) or else this is clearly an x/y problem
    • invest some time in learning fundamental vector calculus, which is inescapable for this class of problem

    Do not discard the magnitude from your field readings.

    It is physically invalid to attempt a surface fit containing the three satellite points. A surface is only applicable if the magnitude read at each satellite is exactly the same, and you've described that it isn't. The magnetic field is literally a field: a vector field with orientations defined at all points in spacetime.

    When you perform regression to establish an estimate for the vector field orientation as a function of position, you will not have one surface. You will have three different isosurfaces, one for each of your satellites, each with a different path through space where the field magnitude is constant.

    Don't perform your calculations in rectilinear space. Perform them in polar space. Your unit normals have too many degrees of freedom for regression (they have three when functionally they need to only have two).

    You need many, many, many more data points to do meaningful regression. As it stands, "you could fit anything". So the following results are numerically accurate but physically meaningless.

    Interpret the system as three independent variables (x, y, z), and two dependent variables (θ and φ, the unit normal rotation angles). This is vulnerable to gimbal lock but we disregard that here.

    After all of the above bad news, the good news is that - as interpreted - the system is highly linear. (It's no surprise given the poverty of data.) You can express the fit as a linear matrix transformation:

    import numpy as np

points = np.array((
    (348.92065834, -1402.33059980,  32.69313966),
    (282.19332067,    82.52027998,  -5.92595371),
    (247.37475615,    -3.70129865, -22.10494737),
))

rect_normals = np.array((
    (-0.86925426, 0.02836434, -0.49355091),
    (-0.82339849, 0.43041935,  0.36980280),
    (-0.83989222, 0.23796899,  0.48780305),
))
x, y, z = rect_normals.T
angles = np.stack((
    np.arccos(z),      # θ
    np.arctan2(y, x),  # φ
)).T

transform, *_ = np.linalg.lstsq(points, angles, rcond=None)

print('Transformation matrix:')
print(transform)
print()

est_angles = points @ transform
print('Measured vs. estimated angles:')
print(np.stack((angles, est_angles)))
print()

θ, φ = est_angles.T
est_normals = np.stack((
    np.sin(θ) * np.cos(φ),
    np.sin(θ) * np.sin(φ),
    np.cos(θ),
)).T
print('Measured vs. estimated unit normals:')
print(np.stack((rect_normals, est_normals)))
print()

    Transformation matrix:
[[ 0.00435349  0.00901185]
 [-0.00038691 -0.00064318]
 [ 0.00077585 -0.02867279]]

Measured vs. estimated angles:
[[[2.08696421 3.10897357]
  [1.19199956 2.65992277]
  [1.06122503 2.86549615]]

 [[2.08696421 3.10897357]
  [1.19199956 2.65992277]
  [1.06122503 2.86549615]]]

Measured vs. estimated unit normals:
[[[-0.86925426  0.02836434 -0.49355091]
  [-0.82339849  0.43041935  0.3698028 ]
  [-0.83989222  0.23796899  0.48780305]]

 [[-0.86925426  0.02836434 -0.49355091]
  [-0.82339849  0.43041935  0.3698028 ]
  [-0.83989222  0.23796899  0.48780305]]]

    So the fit is essentially perfect. The above result is related to (though not equal to) the gradient. To find your isosurfaces - one for each satellite - you need to perform a vector integral. If you do not know how to do this, it's time for research in a different community such as physics or mathematics StackExchange sites.