+ d = 0, where p is a point on plane and n is normal vector\n", "#-------------------------------------------------------------------------------\n", "P_mean = P.mean(axis=0)\n", "P_centered = P - P_mean\n", "U,s,V = linalg.svd(P_centered)\n", "\n", "# Normal vector of fitting plane is given by 3rd column in V\n", "# Note linalg.svd returns V^T, so we need to select 3rd row from V^T\n", "normal = V[2,:]\n", "d = -dot(P_mean, normal) # d = -
\n", "\n", "#-------------------------------------------------------------------------------\n", "# (2) Project points to coords X-Y in 2D plane\n", "#-------------------------------------------------------------------------------\n", "P_xy = rodrigues_rot(P_centered, normal, [0,0,1])\n", "\n", "ax[0].scatter(P_xy[:,0], P_xy[:,1], alpha=alpha_pts, label='Projected points')\n", "\n", "#-------------------------------------------------------------------------------\n", "# (3) Fit circle in new 2D coords\n", "#-------------------------------------------------------------------------------\n", "xc, yc, r = fit_circle_2d(P_xy[:,0], P_xy[:,1])\n", "\n", "#--- Generate circle points in 2D\n", "t = linspace(0, 2*pi, 100)\n", "xx = xc + r*cos(t)\n", "yy = yc + r*sin(t)\n", "\n", "ax[0].plot(xx, yy, 'k--', lw=2, label='Fitting circle')\n", "ax[0].plot(xc, yc, 'k+', ms=10)\n", "ax[0].legend()\n", "\n", "#-------------------------------------------------------------------------------\n", "# (4) Transform circle center back to 3D coords\n", "#-------------------------------------------------------------------------------\n", "C = rodrigues_rot(array([xc,yc,0]), [0,0,1], normal) + P_mean\n", "C = C.flatten()\n", "\n", "#--- Generate points for fitting circle\n", "t = linspace(0, 2*pi, 100)\n", "u = P[0] - C\n", "P_fitcircle = generate_circle_by_vectors(t, C, r, normal, u)\n", "\n", "ax[1].plot(P_fitcircle[:,0], P_fitcircle[:,1], 'k--', lw=2, label='Fitting circle')\n", "ax[2].plot(P_fitcircle[:,0], P_fitcircle[:,2], 'k--', lw=2, label='Fitting circle')\n", "ax[3].plot(P_fitcircle[:,1], P_fitcircle[:,2], 'k--', lw=2, label='Fitting circle')\n", "ax[3].legend()\n", "\n", "#--- Generate points for fitting arc\n", "u = P[0] - C\n", "v = P[-1] - C\n", "theta = angle_between(u, v, normal)\n", "\n", "t = linspace(0, theta, 100)\n", "P_fitarc = generate_circle_by_vectors(t, C, r, normal, u)\n", "\n", "ax[1].plot(P_fitarc[:,0], P_fitarc[:,1], 'k-', lw=3, label='Fitting arc')\n", "ax[2].plot(P_fitarc[:,0], P_fitarc[:,2], 'k-', lw=3, label='Fitting arc')\n", "ax[3].plot(P_fitarc[:,1], P_fitarc[:,2], 'k-', lw=3, label='Fitting arc')\n", "ax[1].plot(C[0], C[1], 'k+', ms=10)\n", "ax[2].plot(C[0], C[2], 'k+', ms=10)\n", "ax[3].plot(C[1], C[2], 'k+', ms=10)\n", "ax[3].legend()\n", "\n", "print('Fitting plane: n = %s' % array_str(normal, precision=4))\n", "print('Fitting circle: center = %s, r = %.4g' % (array_str(C, precision=4), r))\n", "print('Fitting arc: u = %s, θ = %.4g' % (array_str(u, precision=4), theta*180/pi))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib\n", "\n", "fig = figure(figsize=(15,15))\n", "ax = fig.add_subplot(1,1,1,projection='3d')\n", "ax.plot(*P_gen.T, color='y', lw=3, label='Generating circle')\n", "ax.plot(*P.T, ls='', marker='o', alpha=0.5, label='Cluster points P')\n", "\n", "#--- Plot fitting plane\n", "xx, yy = meshgrid(linspace(0,6,11), linspace(0,6,11))\n", "zz = (-normal[0]*xx - normal[1]*yy - d) / normal[2]\n", "ax.plot_surface(xx, yy, zz, rstride=2, cstride=2, color='y' ,alpha=0.2, shade=False)\n", "\n", "#--- Plot fitting circle\n", "ax.plot(*P_fitcircle.T, color='k', ls='--', lw=2, label='Fitting circle')\n", "ax.plot(*P_fitarc.T, color='k', ls='-', lw=3, label='Fitting arc')\n", "ax.set_xlabel('X')\n", "ax.set_ylabel('Y')\n", "ax.set_zlabel('Z')\n", "ax.legend()\n", "\n", "ax.set_aspect('equal', 'datalim')\n", "set_axes_equal_3d(ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4 Conclusion\n", "\n", "To fit a circle to the cluster of points might sound as an easy task, but in 3D space it gets a bit more complicated and the algorithm had to be split into multiple steps. First, using **_SVD_** decomposition we found a plane that fits to the set of 3D points. Then, we projected the 3D points onto the plane and got new planar coordinations for them. Finally, we could use the method of **_least-squares_** to fit a 2D circle into the planar points and then project the 2D fitting circle back to the 3D coordinations. As visible from the pictures, this method gives very satisfying results, both the fitting circle and the fitting arc are very close to the generating circle.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "- [Circle in 3D](http://demonstrations.wolfram.com/ParametricEquationOfACircleIn3D/)\n", "- [Singular value decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition)\n", "- [Least squares](https://en.wikipedia.org/wiki/Least_squares)\n", "- [Rodrigues rotation](https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula)\n", "\n", "\n", "- [NumPy manual](http://docs.scipy.org/doc/numpy/)\n", "- [Matplotlib manual](http://matplotlib.org/)\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.1" }, "latex_envs": { "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 } }, "nbformat": 4, "nbformat_minor": 1 }