5

It's strange that I haven't found a nicely explained example of this on the Internet, since it seems like a fairly standard example. However, I'm trying to intersect a hypersphere and a hyperplane in $\mathbb{R}^n$. The equations are:

$$||x - C||^2 = R^2$$

for the hypersphere, where $C$ is the center and $R$ is the radius, and

$$\vec{v} \cdot \vec{n} = D$$

where $\vec{n} \in \mathbb{R}^n$ is a unit vector. From what I can tell, I do the following:

First I need a point $p$ on the hyperplane, which I can obtain. I then compute the distance between the center of the hypersphere and the hyperplane, which is given by:

$$ \rho = (C - p)\cdot \vec{n}$$

The intersection is nonempty if $-R < \rho < R$, if I am correct, and the intersection is an $n-1$ dimensional hypersphere. Then, the radii and center of the hypersphere are given by:

$$ r = \sqrt{R^2 - \rho^2}$$ $$\bar{c} = C - \rho \vec{n}$$

First of all, is that correct? Second, I want to find a point on the intersecting hypersphere. I've used generalized spherical coordinates to find such a point, but it is in the original space (and thus has $n$ coordinates). Something about that seems very wrong to me, since I appear to have an extra dimension. How do I go about finding points on the intersecting hypersphere? Even worse, when I try to use a point I find using this method, it does not seem to lie on my hyperplane. I'm not sure what went wrong.

Michael Stachowsky
  • 3,133
  • It’s hard to say, not having seen the rest of your work, but my first guess would be that when you try to find a point on the intersection, you’re using the equation of an $(n-1)$-dimensional hypershpere instead of $n-2$, i.e., that you’ve lost the constraint that the point be on the hyperplane. Try working out a concrete example in three dimensions using your method. Doing so might make your error obvious to you. – amd Jun 09 '16 at 16:52
  • I believe it has something to do with the fact that, when I compute a coordinate on the new hypersphere, I have to express it in terms of the basis vectors of the hyperplane, don't I? Even so, was my method to find the center and the radius correct? If so, I think I can work out my error – Michael Stachowsky Jun 09 '16 at 16:59
  • 1
    Here’s one problem for sure: with the sign convention you’re using for the equation of the hyperplane, you should have $\bar{c} = C - \rho\vec{n}$ instead. If you compute $\bar c\cdot n$ for your expression, you get $D+2\rho$ instead of $D$. – amd Jun 09 '16 at 17:19
  • Yes, you're right. I've also now been able to solve the problem I was having. Let's see if I'm right: $\bar{c}$ is in $\mathbb{R}^n$. I first need to project it into the hyperplane, then use the generalized spherical coordinates to produce a point on the now reduced dimensional hyperplane. These coordinates are multiplied by the basis vectors of the hyperplane itself to move me back into the original space. I tried it using a 2-sphere in $\mathbb{R}^3$ and it gave me a circle as expected. Too bad I can't visualize higher dimensions – Michael Stachowsky Jun 09 '16 at 18:03
  • Sounds about right. Rather than projecting $\bar c$ onto the hyperplane, it might be simpler to take it as your new origin, then translate at the end. You’ll need an orthonormal basis for $\vec n^\perp$ either way, which is kind of a nuisance. – amd Jun 09 '16 at 18:48

2 Answers2

4

Your method is sound. You just had a sign error in your original formulation. The center of the $(n-1)$-sphere is at $\bar c=C-\rho\vec n$ (you want to move toward the hyperplane). Given an orthonormal basis $(w_1,\dots,w_{n-1})$ for $\operatorname{span}\{\vec n\}^\perp$ you can then parameterize this $(n-1)$-sphere using generalized spherical coordinates, as you propose: $$\begin{align} \bar c &+ r w_1\cos\phi_1 \\ &+r w_2\sin\phi_1\cos\phi_2 \\ &+r w_3\sin\phi_1\sin\phi_2\cos\phi_3 \\ &+\cdots \\ &+r w_{n-2}\sin\phi_1\sin\phi_2\cdots\cos\phi_{n-2} \\ &+r w_{n-1}\sin\phi_1\sin\phi_2\cdots\sin\phi_{n-2}. \end{align}$$

amd
  • 55,082
4

I found the original question and amd's answer both .. sneaky, in the typical way to mathematicians. Correct, but .. sneaky, for lack of a better term.

A hypersphere of radius $R$ centered at $\vec{c}$ is defined by points $\vec{p}$ that fulfill $$\lVert \vec{p} - \vec{c} \rVert^2 = R^2$$

A hyperplane is defined by its unit normal vector $\hat{n}$, and its minimum signed distance $D$ from origin; its points $\vec{p}$ fulfill $$\vec{p} \cdot \hat{n} = D$$

The minimum distance between the hyperplane and the center of the hypersphere is $$\begin{array}{rl} & d = \left( \vec{c} - D\hat{n} \right ) \cdot \hat{n}\\ \iff & d = \vec{c}\cdot\hat{n} - D\end{array}$$ where the first right side is as in amd's answer (but using the point on the hyperplane closest to origin, $D\hat{n}$, as the measurement point), and the second is the minimum signed distance from the hyperplane to the origin subtracted from the length of the hypersphere center vector $\vec{C}$ projected to the hyperplane surface normal; or, in other terms, the length of the vector between the center of the hypersphere and the point on the hyperplane nearest to origin, projected to/measured along the hyperplane surface normal.

The intersection is nonempty if and only if $-R \le d \le +R$. It does correspond to a hypersphere of one less dimension, centered at $\vec{o} = \vec{c} - d\hat{n}$, $$\begin{array}{rl} &\vec{o} = \vec{c} - d\hat{n}\\ \iff&\vec{o} = \vec{c} - \left(\vec{c}\cdot\hat{n} - D\right)\hat{n}\\ \iff&\vec{o} = \vec{c} + D\hat{n} - \left(\vec{c}\cdot\hat{n}\right)\hat{n}\end{array}$$ with radius $$\begin{array}{ll} &r = \sqrt{R^2 - d^2}\\ \iff&r = \sqrt{R^2 - \left(\vec{c}\cdot\hat{n} - D\right)^2}\end{array}$$ and perpendicular to $\hat{n}$.

The part I call sneaky is related to the three last words in the previous paragraph: perpendicular to $\hat{n}$.

The intersection is the part of the hyperplane bounded by the original hypersphere. That is, for points $\vec{p}$ on the intersection, both $$\begin{cases} \left\lVert\vec{p} - \vec{c}\right\rVert^2 = R^2\\ \vec{p} \cdot \hat{n} = D\end{cases}$$ apply. If you rotate the coordinate system so that the unit normal for the hyperplane $\hat{n}$ is parallel to an axis, the coordinate for that axis can be omitted, and the intersection in the rotated coordinate system simplifies to only a hypersphere of one less dimension, centered at $\vec{o} = \vec{c} - d\hat{n}$ in the original coordinates, with radius $r = \sqrt{R^2 - d^2}$.

Original asker obviously forgot about the rotation needed. In the original coordinate system, both the original hypersphere and the hyperplane equations apply to the points; only in the rotated coordinates is the one-less-dimensional hypersphere enough, and that because the rotation rotates the hyperplane perpendicular to one axis, allowing that axis to be eliminated from consideration. Essentially, it implicitly enforces the hyperplane requirement!

Now, amd did mention this explicitly, just in mathematical terms: a new orthonormal basis spanning the space orthogonal to the hyperplane unit normal vector $\hat{n}$, is just a mathematical description for the rotation needed.

I found this sneaky because it hides the crucial, core point in the answer -- and the point the OP had missed --, as concise jargon. I felt it was like the tiny little over-powered hand cannon in movie Men in Black.

Mathematicians do this all the time. I suspect this is because they want to keep their arts secret to us mere mortals with weaker math-fu.

Implementation-wise, we haven't yet (as of this writing) mentioned any practical ways of constructing the rotation matrix (or the new orthonormal basis vectors, which basically amount to the same thing).

Again, the rotation matrix needed is an orthonormal $d$-by-$d$ square matrix, which rotates the hyperplane surface normal parallel to one of the axes. Most useful would be, if it would rotate it parallel to the positive final axis, as that would make it easier to use in practice.

Although in specific dimensions (2 and 3, in particular) it is possible to construct the rotation directly (via vector cross product and axis-angle representation for the rotation), there is at least two approaches that should work in any number of dimensions (greater than one).

First, and my preferred one, is based on Givens rotations. Start with the original hyperplane surface normal vector, and an identity matrix. Picking the dimension with the largest component in magnitude, excluding the last dimension, apply a Givens rotation to rotate that component to the final dimension, to the current rotation matrix. (This only modifies two columns or two rows in the rotation matrix, depending on how you implement it; and two components in the current unit normal vector, turning the earlier component zero, and increasing the final component in magnitude.) Repeat until the final component is one, and all other components zero, in the normal vector, and you end up with the rotation matrix needed.

Another method is to start by copying the hyperplane unit normal to the rotation matrix unit vector corresponding to the final dimension, and then orthogonalize each new basis vector against previously added new basis vectors. Mathematically, this looks appealing, but in my experience, it is numerically terribly unstable with higher number of dimensions; subtracting the vector sum of dot products with the existing basis vectors seems to be too sensitive to rounding errors to really work in practice.

I intended to include here some example code in C on how to implement the Givens rotations, but my initial attempt didn't rotate the hyperplane unit normal correctly to $(0,\dots,0,1)$, although the rotation matrix was orthonormal. Such issues are common with code doing such rotations/changing to a new orthonormal basis, since it is hard for us code-oriented math-weaklings to keep track of the order of the matrix multiplications (as the order of the rotations is important) and when to use the transpose or not (because for orthonormal matrices, the transpose is the same as the inverse, and this rotation matrix is an orthonormal matrix). The algorithm is not hard per se; there are just an annoying amount of details and complexity to keep in mind. It really needs some care from a math-oriented person, verifying the details (order and transposes taken).

Community
  • 1
Nominal Animal
  • 9,551
  • The OP never asked for an algorithm. I don’t see what you think is “sneaky”: finding a basis for the orthogonal complement of a subspace and generating an orthonormal basis given any basis (via the Gram-Schmidt process, for instance) are both well-known procedures in linear algebra that most might take for granted as known. Your point that G-S can be computationally unstable is well-taken, though, and were I trying to implement this computation, I’d probably do it via a rotation, too. – amd Jun 11 '16 at 06:51
  • @amd: I was afraid you might misunderstand me. Your answer, in my opinion, is right. The part that I call sneaky is how the key point, "Given an orthonormal basis .. for span .." is easy to miss. "sneaky" is a poor word for it, though, as it is perfectly valid definition -- just using extremely concise terminology. No slight of any shape or form was intended at all; I was actually impressed how you could express in one paragraph what takes me a page to wade through. – Nominal Animal Jun 11 '16 at 13:13
  • @amd: I only had good intentions, but it seems it backfired awfully. I suspect it might be better to delete this answer altogether. I do admire your math skills (and thank you for naming the Gram-Schmidt process; I racked my brain but both it and my search-fu failed me), and I'd rather not have you (or others with skills necessary to describe the solutions for us others with weaker math-fu) discouraged from helping us. I'd rather be banned from here, actually. – Nominal Animal Jun 11 '16 at 13:17
  • @amd: I'd delete this answer, but I cannot, as long as it is accepted. – Nominal Animal Jun 13 '16 at 20:58
  • Why delete? There’s useful information in your answer and the OP liked what you wrote. I certainly wasn’t offended, but only somewhat puzzled by your choice of words. – amd Jun 14 '16 at 19:33
  • @amd: I thought you were offended, and that was definitely not my intent. I used "sneaky" in the "difficult to catch due to constantly outwitting the adversaries" sense only. Your answer is concise, complete, and looks very simple on the surface. The new orthonormal basis that rotates the hyperplane normal to an axis, is central to the answer, and the point where OP stumbled. I only noticed this, because I see the same in other math answers and journal articles all the time: core points in the solution are passed with just a few words. – Nominal Animal Jun 14 '16 at 20:53
  • @amd: Yet, mathematicians do good work. It is just that us others, who only use math when solving problems, the path towards solutions is way more important than a mathematical description for the requirements of the solution, or a proof that such a solution is correct (or exists). Typically, we are limited to a specific range of reals for all variables, for example, and don't need to worry about singularities or such; but floating-point numerical stability is a huge concern. We have different viewpoints. I see my answer as basically a translation of yours, to us non-mathematicians. – Nominal Animal Jun 14 '16 at 20:56
  • Just to make sure I understand, to use the equations you provided for the center and radius of the intersection hypersphere, you have to first rotate the hyperplane equation so that the normal is axis-aligned with some axis--is that correct, or do your equations already reflect that rotation? – bob Nov 22 '21 at 19:56