Compute the N-dimentional (hyper-)circumsphere

Hypercircumsphere.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute the N-dimentional (hyper-)circumsphere\n",
    "\n",
    "This function computes the N-dimentional hypersphere that circumscribes N+1 points in N dimensions.\n",
    "\n",
    "The code uses the fact that all points that are equidistant between two points lie on the hyperplane that perpendicularly bisects the line between the two points. Since the centre of the hyper-circumsphere covering a set of points is equidistant from each point, it must lie on all of the bisecting planes between pairs of points. By finding the intersection of a set of these planes we can find the centre.\n",
    "\n",
    "This code makes fairly efficient use of `numpy` such that the hard work is done in a total of four lines of code (which could easily be crammed into one if clarity was not a goal):\n",
    "```\n",
    "    directions = points[1:] - points[0]\n",
    "    midpoints = (points[1:] + points[0])/2\n",
    "    plane_k = (directions * midpoints).sum(axis=1)\n",
    "    answer = np.linalg.solve(directions, plane_k)\n",
    "```\n",
    "\n",
    "As a result this code can efficiently find a circumcircle of a triangle or the circumsphere of a tetrahedron.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def hypercircumsphere(points, return_R=False):\n",
    "    # First, ensure that we are working with a numpy array\n",
    "    if not isinstance(points, np.ndarray):\n",
    "        points = np.array(points)\n",
    "    # Also ensure that we are working with a floating point type\n",
    "    if not isinstance(points[0][0], np.floating):\n",
    "        points = np.array(points, dtype=np.float)\n",
    "        \n",
    "    # Check that for N dimensions we have N+1 points\n",
    "    shape = points.shape\n",
    "    if len(shape) != 2 or shape[0] != shape[1] + 1:\n",
    "        raise ValueError(\"Exactly N+1 points needed for N dimensions\")\n",
    "\n",
    "    # Compute the  directions of the vectors between the points.\n",
    "    # Implicitly these are the normals of the planes that are normal to them\n",
    "    directions = points[1:] - points[0]\n",
    "\n",
    "    # If the points all lie in some lower-dimensional subspace then there is\n",
    "    # no solution to this problem\n",
    "    if np.linalg.det(directions) == 0:\n",
    "        raise ValueError(f\"Points all lie in a subspace smaller than {shape[1]} dimensions\")\n",
    "        \n",
    "    # Compute the mid-points along these connections\n",
    "    midpoints = (points[1:] + points[0])/2\n",
    "\n",
    "    # A plane can be defined by its normal and a constant value for the dot\n",
    "    # product of the points and that nomrmal. We want the planes that go through\n",
    "    # the midpoint between each pair of points, so the constant can bbe found\n",
    "    # by just evaluating at the midpoint\n",
    "    plane_k = (directions * midpoints).sum(axis=1)\n",
    "    \n",
    "    # The center of the circumsphere is where the planes intersect, so we solve\n",
    "    # for the simultanious set of plane equations\n",
    "    answer = np.linalg.solve(directions, plane_k)\n",
    "\n",
    "    if return_R:\n",
    "        return (answer, ((points[0] - answer) ** 2).sum() ** 0.5)\n",
    "    else:\n",
    "        return answer    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test for 2 dimensions: checks PASSED\n",
      "Test for 3 dimensions: checks PASSED\n",
      "Test for 4 dimensions: checks PASSED\n"
     ]
    }
   ],
   "source": [
    "p1_2d = [0, 1]\n",
    "p2_2d = [1, 0]\n",
    "p3_2d = [2, 2]\n",
    "points_2d = [p1_2d, p2_2d, p3_2d]\n",
    "\n",
    "p1_3d = [1., 1., 1.]\n",
    "p2_3d = [2., 1., 1.]\n",
    "p3_3d = [1., 3., 1.]\n",
    "p4_3d = [1., 1., 4.]\n",
    "points_3d = [p1_3d, p2_3d, p3_3d, p4_3d]\n",
    "\n",
    "p1_4d = [1, 2, 3, 4]\n",
    "p2_4d = [1, 3, 3, 4]\n",
    "p3_4d = [0, 0, 1, 1]\n",
    "p4_4d = [0, 0, 0, 0]\n",
    "p5_4d = [4, 4, 4, 4]\n",
    "points_4d = [p1_4d, p2_4d, p3_4d, p4_4d, p5_4d]\n",
    "\n",
    "tests = [points_2d, points_3d, points_4d]\n",
    "\n",
    "for test in tests:\n",
    "    answer, radius = hypercircumsphere(test, return_R=True)\n",
    "    checks = ((test - answer) ** 2).sum(1) ** 0.5\n",
    "    # Check with np.allclose() due to potential for rounding errors\n",
    "    correct = np.allclose(checks, radius)\n",
    "    print(\"Test for {} dimensions: checks {}\".format(len(answer), \"PASSED\" if correct else \"FAILED\"))\n",
    "    if not correct:\n",
    "        print(f\"Radius={radius}, checks={checks}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}