Circles and spheresIn what follows are various notes and algorithms dealing with circles and spheres.
Spheres, equations and terminologyWritten by Paul BourkeApril 1992 OpenGL/GLUT source code demonstrating the Great Circle
The most basic definition of the surface of a sphere is "the set of points an equal distance (called the radius) from a single point called the center". Or as a function of 3 space coordinates (x,y,z), all the points satisfying the following lie on a sphere of radius r centered at the origin
For a sphere centered at a point (x_{o},y_{o},z_{o}) the equation is simply
If the expression on the left is less than r^{2} then the point (x,y,z) is on the interior of the sphere, if greater than r^{2} it is on the exterior of the sphere. A sphere may be defined parametrically in terms of (u,v)
y = y_{o} + r cos(phi) sin(theta) z = z_{o} + r sin(phi) Where 0 <= theta < 2 pi, and pi/2 <= phi <= pi/2. The convention in common usage is for lines of constant theta to run from one pole (phi = pi/2 for the south pole) to the other pole (phi = pi/2 for the north pole) and are usually referred to as lines of longitude. Lines of constant phi are often referred to as lines of latitude, for example the equator is at phi = 0. Lines through a sphere
A line can intersect a sphere at one point in which case it is called a tangent. It can not intersect the sphere at all or it can intersect the sphere at two points, the entry and exit points. For the mathematics for the intersection point(s) of a line (or line segment) and a sphere see this. Antipodal points
A line that passes through the center of a sphere has two intersection points, these are called antipodal points. Planes through a sphere
A plane can intersect a sphere at one point in which case it is called a tangent plane. Otherwise if a plane intersects a sphere the "cut" is a circle. Lines of latitude are examples of planes that intersect the Earth sphere.
Lune
A lune is the area between two great circles who share antipodal points. If the angle between the planes defining the great circle is A, then the area of a lune on a sphere of radius r is
Ellipsoid
An ellipsoid squashed along each (x,y,z) axis by a,b,c is defined as
Or parametrically
y = y_{o} + b r cos(theta) sin(phi) z = z_{o} + c r sin(theta)
Equation of a Circle from 3 Points (2 dimensions)Written by Paul BourkeJanuary 1990
This note describes a technique for determining the attributes of a circle (centre and radius) given three points P_{1}, P_{2}, and P_{3} on a plane.
Calculating Centre
Two lines can be formed through 2 pairs of the three points, the first passes
through the first two points P_{1}
and P_{2}. Line b passes through the
next two points P_{2} and P_{3}.
where m is the slope of the line given by
The centre of the circle is the intersection of the two lines perpendicular to and passing through the midpoints of the lines P_{1}P_{2} and P_{2} P_{3}. The perpendicular of a line with slope m has slope 1/m, thus equations of the lines perpendicular to lines a and b and passing through the midpoints of P_{1}P_{2} and P_{2}P_{3} are These two lines intersect at the centre, solving for x gives
Calculate the y value of the centre by substituting the x value into one of the equations of the perpendiculars. Alternatively one can also rearrange the equations of the perpendiculars and solve for y. Radius
The radius is easy, for example the point P_{1} lies on the circle and we know the centre.... Notes:
C++ code implemented as MFC (MS Foundation Class) supplied by Jae Hun Ryu. Circle.cpp, Circle.h.
Equation of a Sphere from 4 Points on the SurfaceWritten by Paul BourkeJune 2002
Given 4 points in 3 dimensional space [ (x_{1},y_{1},z_{1}) (x_{2},y_{2},z_{2}) (x_{3},y_{3},z_{3}) (x_{4},y_{4},z_{4}) ] the equation of the sphere with those points on the surface is found by solving the following determinant.
There are conditions on the 4 points, they are listed below and correspond to the determinant above being undefined (no solutions, multiple solutions, or infinite solutions).
If the determinant is found using the expansion by minors using the top row then the equation of the sphere can be written as follows.
Or more simply in term of the minors M_{1j}
Equating the terms from these two equations allows one to solve for the center and radius of the sphere, namely:
Note that these can't be solved for M_{11} equal to zero. This corresponds to no quadratic terms (x^{2}, y^{2}, z^{2}) in which case we aren't dealing with a sphere and the points are either coplanar or three are collinear. C source codeThis piece of simple C code tests the solution as described above. It creates a known sphere (center and radius) and creates 4 random points on that sphere. It then proceeds to find the original center and radius using those four random points.
Intersection of a Line and a Sphere (or circle)Written by Paul BourkeNovember 1992
Points P (x,y) on a line defined by two points P_{1} (x_{1},y_{1},z_{1}) and P_{2} (x_{2},y_{2},z_{2}) is described by
or in each coordinate y = y_{1} + u (y_{2}  y_{1}) z = z_{1} + u (z_{2}  z_{1})
A sphere centered at P_{3} (x_{3},y_{3},z_{3}) with radius r is described by Substituting the equation of the line into the sphere gives a quadratic equation of the form
where: b = 2[ (x_{2}  x_{1}) (x_{1}  x_{3}) + (y_{2}  y_{1}) (y_{1}  y_{3}) + (z_{2}  z_{1}) (z_{1}  z_{3}) ] c = x_{3}^{2} + y_{3}^{2} + z_{3}^{2} + x_{1}^{2} + y_{1}^{2} + z_{1}^{2}  2[x_{3} x_{1} + y_{3} y_{1} + z_{3} z_{1}]  r^{2} The solutions to this quadratic are described by
The exact behaviour is determined by the expression within the square root
To apply this to two dimensions, that is, the intersection of a line and a circle simply remove the z component from the above mathematics. Line SegmentFor a line segment between P_{1} and P_{2} there are 5 cases to consider.
When dealing with a line segment it may be more efficient to first determine whether the line actually intersects the sphere or circle. This is achieved by noting that the closest point on the line through P_{1}P_{2} to the point P_{3} is along a perpendicular from P_{3} to the line. In other words if P is the closest point on the line then Substituting the equation of the line into this Solving the above for u =  (x_{2}  x_{1})(x_{2}  x_{1}) + (y_{2}  y_{1})(y_{2}  y_{1}) + (z_{2}  z_{1})(z_{2}  z_{1}) If u is not between 0 and 1 then the closest point is not between P_{1} and P_{2} Given u, the intersection point can be found, it must also be less than the radius r. If these two tests succeed then the earlier calculation of the actual intersection point can be applied. Source code
C code example by author.
Intersection of two circlesWritten by Paul BourkeApril 1997
The following note describes how to find the intersection point(s) between two circles on a plane, the following notation is used. The aim is to find the two points P_{3} = (x_{3}, y_{3}) if they exist. First calculate the distance d between the center of the circles. d = P_{1}  P_{0}.
Using d = a + b we can solve for a,
It can be readily shown that this reduces to r_{0} when the two circles touch at one point, ie: d = r_{0} ± r_{1} Solve for h by substituting a into the first equation, h^{2} = r_{0}^{2}  a^{2} SoAnd finally, P_{3} = (x_{3},y_{3}) in terms of P_{0} = (x_{0},y_{0}), P_{1} = (x_{1},y_{1}) and P_{2} = (x_{2},y_{2}), is y_{3} = y_{2} + h ( x_{1}  x_{0} ) / d Source code contributions
Python version by Matt Woodhead.
Intersection of two spheresWritten by Paul BourkeNovember 1995 Consider two spheres on the x axis, one centered at the origin, separated by a distance d, and of radius r_{1} and r_{2}. The equations of the two spheres are
(x  d)^{2} + y^{2} + z^{2} = r_{2}^{2} Subtracting the first equation from the second, expanding the powers, and solving for x gives
The intersection of the two spheres is a circle perpendicular to the x axis, at a position given by x above. Substituting this into the equation of the first sphere gives where
Distributing Points on a SphereWritten by Paul BourkeJune 1996 The following describes two (inefficient) methods of evenly distributing points on a sphere. They do however allow for an arbitrary number of points to be distributed unlike many other algorithms which only work for a restricted set of points. The first approach is to randomly distribute the required number of points on a sphere of the desired radius. The iteration involves finding the closest two points and then moving them apart slightly. The following shows the results for 100 and 400 points, the disks have a radius of the minimum distance. C source code: source1.c Physical method
If we place the same electric charge on each particle (except perhaps the particle in the center) then each particle will repel every other particle. This system will tend to a stable configuration where each particle is equidistant from the center (due to spring forces) and each particle maximally separated from its closest neighbours (electric repulsive forces). It is important to model this with viscous damping as well as with spring damping to avoid oscillatory motion. An example using 31 particles randomly distributed in a cube is shown in the animation above. A midpoint ODE solver was used to solve the equations of motion, it took only 200 steps to reach a stable (minimum energy) configuration. Uniform DistributionA simple way to randomly (uniform) distribute points on sphere is called the "hypercube rejection method". To apply this to a unit cube at the origin, choose coordinates (x,y,z) each uniformly distributed on the interval [1,1]. If the length of this vector is greater than 1 then reject it, otherwise normalise it and use it as a sample. Contribution by Jonathan D. LettvinC++ source code: diffuse.cpp Contribution by Max DowneyJava implementation: java.tar.gz Contribution by Orion ElenzilOrion Elenzil proposes that by choosing uniformly distributed polar coordinates theta (0 <= theta < 360) and phi (0 <= phi <= pi/2) but the using the sqrt(phi) results in points uniformly distributed on the surface of a hemisphere. If the poles lie along the z axis then the position on a unit hemisphere sphere is y = cos(sqrt(phi)) sin(theta) z = sin(sqrt(phi)) A whole sphere is obtained by simply randomising the sign of z.
Area of multiple intersecting circlesWritten by Paul BourkeJanuary 2000 Introduction
The following is a straightforward but good example of a range of techniques called "MonteCarlo" methods. It will be used here to numerically find the area of intersection of a number of circles on a plane. Consider a single circle with radius r, the area is pi r^{2}. The minimal square enclosing that circle has sides 2 r and therefore an area of 4 r^{2}. If one was to choose random numbers from a uniform distribution within the bounding rectangle then the ratio of those falling within the circle to the total number will be the ratio of the area of the circle to the rectangle. In other words, countinside/totalcount = pi / 4, this ratio of pi / 4 would be approached closer as the totalcount increases.. This could be used as a way of estimate pi, albeit a very inefficient way!
C source that numerically estimates the intersection area of any number of circles on a plane is given here: area.c. The basic idea is to choose a random point within the bounding square of one of the circles and check to see if the point is within all the other circles. The successful count is scaled by 4 r^{2} / totalcount to give the area of the intersecting piece.
Creating a plane/disk perpendicular to a line segmentWritten by Paul BourkeFebruary 1997 Contribution by Dan Wills in MEL (Maya Embedded Language): source2.mel. There are a number of 3D geometric construction techniques that require a coordinate system perpendicular to a line segment, some examples are:
A straightforward method will be described which facilitates each of these. The key is deriving a pair of orthonormal vectors on the plane perpendicular to a line segment P_{1}, P_{2}. Procedure
1. Choose any point P randomly which doesn't lie on the line
through P_{1} and P_{2}
Points on the plane through P_{1} and perpendicular to
n = P_{2}  P_{1} can be found from linear combinations
of the unit vectors R and S, for example, a point Q might be
Q_{y} = P_{1y} + alpha R_{y} + beta S_{y} Q_{z} = P_{1z} + alpha R_{z} + beta S_{z} for all alpha and beta. A disk of radius r, centered at P_{1}, with normal n = P_{2}  P_{1} is described as follows
Q_{y} = P_{1y} + r cos(theta) R_{y} + r sin(theta) S_{y} Q_{z} = P_{1z} + r cos(theta) R_{z} + r sin(theta) S_{z} for 0 <= theta <= 2 pi Example The following is a simple example of a disk and the C source stub that generated it.
Sphere GenerationWritten by Paul BourkeMay 1992 Polar CoordinatesThe following illustrate methods for generating a facet approximation to a sphere. The most straightforward method uses polar to Cartesian coordinates, if theta and phi as shown in the diagram below are varied the resulting vector describes points on the surface of a sphere.
To create a facet approximation, theta and phi are stepped in small angles between their respective bounds. So if we take the angle step size to be dtheta and dphi, the four vertices of any facet correspond to
(theta+dtheta,phi) (theta+dtheta,phi+dphi) (theta,phi+dphi)
The sphere can be generated at any resolution, the following shows a progression from 45 degrees through to 5 degree angle increments. The number of facets being (180 / dtheta) (360 / dphi), the 5 degree example on the right contains almost 2600 facets. Source code
Surface refinement
The following illustrates the sphere after 5 iterations, the number of facets increases on each iteration by 4 so this representation has 1024 facets.
Perhaps unexpectedly, all the facets are not the same size, those nearer the vertices of the original tetrahedron are smaller. The above example resulted in a triangular faceted model, if a cube is used as the starting form then a representation with rectangular facets can be derived.
As in the tetrahedron example the facets are split into 4 and thus the number of facets increases by a factor of 4 on each iteration. The representation on the far right consists of 6144 facets.
The nonuniformity of the facets most disappears if one uses an octahedron as the starting shape. Bisecting the triangular facets results in sphere approximations with 8, 32, 128, 512, 2048, .... facets as the iteration count increases.
Source code
Uniform Distribution A simple way to randomly (uniform) distribute points on sphere is called the "hypercube rejection method". To apply this to a unit cube at the origin, choose coordinates (x,y,z) each uniformly distributed on the interval [1,1]. If the length of this vector is greater than 1 then reject it, otherwise normalise it and use it as a sample. Spherical triangle
The same technique can be used to form and represent a spherical triangle, that is, the triangle formed by three points on the surface of a sphere, bordered by three great circle segments.
The algorithm and the conventions used in the sample source code provided is illustrated below. The three vertices of the triangle are each defined by two angles, longitude and latitude, on each iteration the number of triangles increases by a factor of 4.
The following images show the cylinders with either 4 vertex faces or entirely 3 vertex facets. Note that since the 4 vertex polygons are coplanar, splitting them into two 3 vertex facets doesn't improve the resolution. End caps are normally optional, whether they are needed or not is application dependent.
q[0] = P1 + r1 * cos(theta1) * A + r1 * sin(theta1) * B
Note P1,P2,A, and B are all vectors in 3 space. r1 and r2 are the radii at the two ends. Written as some pseudo C code the facets might be created as follows create A and B for (i=0;i<NFACETS;i++) { n = 0; theta1 = i * TWOPI / N; theta2 = (i+1) * TWOPI / N; q[n].x = P1.x + r1 * cos(theta1) * A.x + r1 * sin(theta1) * B.x q[n].y = P1.y + r1 * cos(theta1) * A.y + r1 * sin(theta1) * B.y q[n].z = P1.z + r1 * cos(theta1) * A.z + r1 * sin(theta1) * B.z n++; q[n].x = P2.x + r2 * cos(theta1) * A.x + r2 * sin(theta1) * B.x q[n].y = P2.y + r2 * cos(theta1) * A.y + r2 * sin(theta1) * B.y q[n].z = P2.z + r2 * cos(theta1) * A.z + r2 * sin(theta1) * B.z n++; if (r2 != 0) { q[n].x = P2.x + r2 * cos(theta2) * A.x + r2 * sin(theta2) * B.x q[n].y = P2.y + r2 * cos(theta2) * A.y + r2 * sin(theta2) * B.y q[n].z = P2.z + r2 * cos(theta2) * A.z + r2 * sin(theta2) * B.z n++; } if (r1 != 0) { q[n].x = P1.x + r1 * cos(theta2) * A.x + r1 * sin(theta2) * B.x q[n].y = P1.y + r1 * cos(theta2) * A.y + r1 * sin(theta2) * B.y q[n].z = P1.z + r1 * cos(theta2) * A.z + r1 * sin(theta2) * B.z n++; } do something with q[0..n] }Note
Modelling with SpheresWritten by Paul BourkeAugust 1991 Many computer modelling and visualisation problems lend themselves to placing markers at points in 3 space. It may be that such markers are a natural consequence of the object being studied (for example: chaotic attractors) or it may be that forming other higher level primitives such as tubes or planar facets may be problematic given the description of the object being modelled. A simple and directionally symmetric marker is the sphere, a point is discounted here, even though it can be considered to be a sphere of zero radius, because most rendering packages do not support such ideal nonreal entities. (A ray from a raytracer will never intersect a point which occupies no volume, in the same way, lines can generally not be rendered) Another reason for wanting to model using spheres as markers is that many rendering packages handle spheres very efficiently. The computationally expensive part of raytracing geometric primitives is testing the intersection of a ray with the primitive. Such a test for a sphere is the most efficient of all primitives, one only needs to determine whether the closest position of the center of the sphere to the ray is less than the radius of the sphere. Curve ExampleThe first example will be modelling a curve in space. The figures below show the same curve represented with an increased number of points, a sphere at each point.
Chaotic ExampleModelling chaotic attractors is a natural candidate for modelling with spheres because the points are not generated sequentially.
Surface ExampleSurfaces can also be modelled with spheres although this can obviously be very inefficient. The following is an example from a project to visualise the Steiner surface.
Modelling RopeEach strand of the rope is modelled as a series of spheres, each tracing a sinusoidal route through space.
Sea ShellsSome biological forms lend themselves naturally to being modelled with spherical building blocks as it adds an existing surface texture. Some sea shells for example have a rippled effect Representing attractors
Rounded boxWritten by Paul BourkeJuly 1989
Creating box shapes is very common in computer modelling applications. The boxes used to form walls, table tops, steps, etc generally have perfectly sharp edges. Such sharpness does not normally occur in real life because of wear and for safety reasons. There are many ways of introducing curvature and ideally this would be done in the rendering phase. One modelling technique is to turn edges into cylinders and the corners into spheres. The planar facets that made up the original object are trimmed back until they are tangent to the sphere and/or cylinder surface. To illustrate this consider the following which shows the corner of a box converted into a corner with curvature. Over the whole box, each of the 6 facets reduce in size, each of the 12 edges become cylinders, and each of the 8 vertices become spheres. One problem with this technique as described here is that the resulting object does not normally have the desired effect internally. If this is important then the cylinders and spheres described above need to be turned into the appropriate cylindrical and spherical wedges/sections. In the following example a cube with sides of length 2 and increasing edge radii is used to illustrate the effect. radius = 0Just a cube Degenerate case where the four spheres are coincident, the cylinders don't exist, and there are no facets.
Pipes for Rendering EnginesWritten by Paul BourkeJune 1996 Most rendering engines support simple geometric primitives such as planes, spheres, cylinders, cones, etc. Many times a pipe is needed, by pipe I am referring to a tube like structure which passes through 3D space. The actual path is irrelevant but might be an arc or a Bezier/Spline curve defined by control points in space. The standard method of geometrically representing this structure, as illustrated here, uses combinations of cylinders and spheres. Basically the curve is split into a straight line approximation to the desired level or resolution. Each straight line segment is represented by a cylinder. Since this would lead to gaps at the intersection of cylinders, spheres of the same radius are placed at the intersection points. Optionally disks can be placed at the end points to seal the pipe.
Note 1 If the radius of the pipe is to change along the path then the cylinders need to be replaced with a cone sections, namely a cylinder with different radii at each end. The radius of each cylinder is the same at an intersection point so an appropriate sphere still fills the gaps. Note 2 This method is only suitable if the pipe is to be viewed from the outside. As an example, the following pipes are arc paths, 20 straight line sections per pipe.
