1
$\begingroup$

I've seen this cool function, which generates a given distribution of points on the unit sphere, known as Dimroth-Watson distribution:

(* special case; κ = 0 is the uniform distribution *) dimrothWatsonRandom[μ_?VectorQ, κ_ /; κ == 0] := Normalize[RandomVariate[NormalDistribution[], 3]] dimrothWatsonRandom[μ_?VectorQ, κ_ /; NumericQ[κ] && Positive[κ]] := Module[{c, u, v, w}, c = Exp[-κ/2] Csch[κ/2]/2; While[ {u, v} = RandomReal[1, 2]; w = Log[1 + u/c]/κ; v > Exp[κ w (w - 1)]]; RotationTransform[{{0, 0, 1}, Normalize[μ]}][Append[ Sqrt[1 - w^2] Normalize[RandomVariate[NormalDistribution[], 2]], RandomChoice[{-1, 1}] w]]] dimrothWatsonRandom[μ_?VectorQ, κ_ /; NumericQ[κ] && Negative[κ]] := Module[{c, d, u, v, w}, c = Sqrt[-κ]; d = ArcTan[c]; While[ {u, v} = RandomReal[1, 2]; w = Tan[d u]/c; v > (1 - κ w^2) Exp[κ w^2]]; RotationTransform[{{0, 0, 1}, Normalize[μ]}][Append[ Sqrt[1 - w^2] Normalize[RandomVariate[NormalDistribution[], 2]], RandomChoice[{-1, 1}] w]]]

Let us consider a specific case, with $k>0$, so that the points become 'bipolar'. Also, let's just add a Sphere[] to visualize the points, the center, and the direction of the vector:

k = 4; dirVector = {1, 0, 0}; center = {0, 0, 0}; nPoints = 1000; myPoints = Table[dimrothWatsonRandom[dirVector, k], {nPoints}]; (*Visualize it*) Graphics3D[{Opacity[0.125], Sphere[], Opacity[0.5], Red, PointSize[0.02], Point[center], Opacity[0.25], Blue, PointSize[0.01], Point[myPoints], Opacity[0.5], Red, PointSize[0.02], Point[dirVector]}]

enter image description here

Question is: can we map the points from the surface of this sphere onto the surface of an Ellipsoid[]?

So far, my attempt is to create a small ellipsoid (so it's inside the unit sphere), then find the nearest points to its surface:

reg = Ellipsoid[{0, 0, 0}, {0.5, 0.25, 0.25}]; npl = Table[RegionNearest[reg, p], {p, myPoints}]; Legended[Graphics3D[{reg, {Thin, Gray, Line[Transpose[{myPoints, npl}]]}, {Red, Point[myPoints]}, {PointSize[Medium], Blue, Point[npl]}}, Lighting -> "Neutral", Boxed -> False], PointLegend[{Red, Blue}, {"start", "nearest"}]]

enter image description here

I'm hoping there's a more elegant way to do this.

Edit 1: In the example I show the circle and ellipsoid have the same orientation because of the dirVector I used. Suppose that the ellipsoid has some arbitrary orientation. Thus, the transfer of points should take the direction of dirVector and align it with the major axis of the ellipsoid. As flinty pointed out in the comments, this would require some rotation/translation transformation, which I don't know how'd be applied, so it would be great if somebody comes up with a solution for this too.

Edit 2: Could something like this be done here, preserving the relative positions of the points?

Improve this question
$\endgroup$
4
  • $\begingroup$ Can you just scale two of the coordinates down, but leave the other coordinate (assuming an axis aligned ellipse) i.e just squash the circle points in two directions? That's just a linear transformation which is kind of trivial. If it weren't axis aligned then you'd need another matrix to rotate the coordinate frame first, squash, then un-rotate. $\endgroup$ Commented Aug 17, 2021 at 20:59
  • $\begingroup$ @flinty I'm not sure I understand about squashing only two directions. This squashing would need some fine tuning on more squashing towards the center and less towards the ends, no? Also, what you say about the axis-aligned transfer is a great point I hadn't considered. I'll add a comment on that. Ideally, this would keep the direction vector aligned to the major axis of the ellipsoid. $\endgroup$ Commented Aug 17, 2021 at 22:35
  • $\begingroup$ Are the points projected on to the ellipsoid not the nearest using the spherical radius? Also could you take a distribution of points on an ellipsoid and Manipulate[] the eccentricity to a sphere? $\endgroup$ Commented Aug 17, 2021 at 23:24
  • $\begingroup$ @TumbiSapichu I meant like this: Graphics3D[{reg, Blue, Point[reg[[2]]*# & /@ myPoints], Green, Point[npl] }] - but that will change the density of points - note the discrepancy between blue and green points. For the rotation this is really easy, just use a RotationTransform[{u,{1,0,0}}] that moves the vector u to {1,0,0} and apply to all points on the sphere, do stuff with the axis aligned ellipse, then rotate back with RotationTransform[{{1,0,0},u}] $\endgroup$ Commented Aug 18, 2021 at 11:52

2 Answers 2

Reset to default
4
$\begingroup$

Let's approach the problem from the geometrical perspective. Let the ellipsoid semi-axes be $(a, b, c)$. Let the point on the sphere be $\vec{r}=(x,y,z)$, the corresponding point on the ellipsoid $\vec{r}' = (x',y',z')$, and the normal vector to the surface of ellipsoid at the point $\vec{r}'$ be $\vec{p}$ so that the following is true for some $\xi \geq 0$ (ellipsoid is smaller than sphere):

$$ \vec{r} = \vec{r}' + \xi \ \vec{p}.$$

The normal vector on the ellipsoid is, as shown here: $$ \vec{p} = \left(\frac{x'}{a^2}, \frac{y'}{b^2}, \frac{z'}{c^2}\right).$$ Now we have to simultaneously solve the following four equations:

$$\begin{align} x' &= \frac{x}{1+\xi/a^2} \\ y' &= \frac{y}{1+\xi/b^2} \\ z' &= \frac{z}{1+\xi/c^2} \\ 1 &= \left(\frac{x'}{a}\right)^2 + \left(\frac{y'}{b}\right)^2 + \left(\frac{z'}{c}\right)^2 \end{align}$$

Solving this for a generic point is not useful, because we cannot know which root is the correct one. So we solve the system for each point separately.

a = 1/2; b = 1/4; c = 1/4; reg = Ellipsoid[{0, 0, 0}, {a, b, c}]; npl = {xp, yp, zp} /. FindRoot[ {(xp/a)^2 + (yp/b)^2 + (zp/c)^2 == 1, xp == #[[1]]/(1 + \[Xi]/a^2), yp == #[[2]]/(1 + \[Xi]/b^2), zp == #[[3]]/(1 + \[Xi]/c^2)}, {{xp, #[[1]]}, {yp, #[[2]]}, {zp, #[[3]]}, {\[Xi], .1}}] & /@ myPoints; Legended[Graphics3D[{reg, {Thin, Gray, Line[Transpose[{myPoints, npl}]]}, {Red, Point[myPoints]}, {PointSize[Medium], Blue, Point[npl]}}, Lighting -> "Neutral", Boxed -> False], PointLegend[{Red, Blue}, {"start", "nearest"}]]

Mathematica graphics

Let's compare the results to RegionNearest:

npl1 = Table[RegionNearest[reg, p], {p, myPoints}]; Norm /@ (npl - npl1) (* {4.2708*10^-16, 5.55116*10^-17, 5.76388*10^-17, ...} *)

Results seem to coincide numerically. However, let's take a look at how RegionNearest actually does the job. As you might have noticed, I have defined the semi-axes with exact numbers $(1/2, 1/4, 1/4)$. Now let's take a random point with exact coordinates $(1,1,2)$ and project it to the ellipsoid:

sol1 = RegionNearest[reg, {1, 1, 2}] // FullSimplify (* {Root[-16 + 24 # + 135 #^2 - 96 #^3 + 36 #^4& , 2, 0], Root[-1 - 6 # + 135 #^2 + 480 #^3 + 720 #^4& , 2, 0], Root[-4 - 12 # + 135 #^2 + 240 #^3 + 180 #^4& , 2, 0]} *)

We can see that RegionNearest with Ellipsoid isn't doing any discretization, but it solves the above system of equations. The same result with our method can be obtained if we use Solve instead of FindRoot:

sol2 = {xp, yp, zp} /. First@Solve[{(xp/a)^2 + (yp/b)^2 + (zp/c)^2 == 1, xp == x/(1 + \[Xi]/a^2), yp == y/(1 + \[Xi]/b^2), zp == z/(1 + \[Xi]/c^2), \[Xi] > 0} /. {x -> 1, y -> 1, z -> 2}, {xp, yp, zp, \[Xi]}, Reals] // FullSimplify (* {Root[-16 + 24 # + 135 #^2 - 96 #^3 + 36 #^4& , 2, 0], Root[-1 - 6 # + 135 #^2 + 480 #^3 + 720 #^4& , 2, 0], Root[-4 - 12 # + 135 #^2 + 240 #^3 + 180 #^4& , 2, 0]} *) sol1 == sol2 (* True *)

Have we gained anything by solving the system with FindRoot?

First@RepeatedTiming@(RegionNearest[reg, #] & /@ myPoints) (* 3.37717 *) First@RepeatedTiming@(FindRoot[ (* ... *) ] & /@ myPoints) (* 0.780355 *)

It runs approximately 4 times faster.

Improve this answer
$\endgroup$
2
$\begingroup$

It's probably worth pointing out that RegionNearestFunciton[..][<list of points>] is efficiently implemented:

reg = Ellipsoid[{0, 0, 0}, {1/2, 1/4, 1/4}]; rnf = RegionNearest[reg]; pts = RandomReal[{-5, 5}, {1000, 3}]; res1 = rnf[pts]; // RepeatedTiming (* {0.556633, Null} *) res2 = rnf /@ pts; // RepeatedTiming (* {13.3916, Null} *) res1 == res2 (* True *)

Also interesting:

I don't know why rnf /@ pts is so much slower than:

RegionNearest[reg, #] & /@ pts; // RepeatedTiming (* {2.78101, Null} *)

Maybe something to do with autocompile?:

rnf[#] & /@ pts; // RepeatedTiming (* {0.541443, Null} *)

Building the RegionNearestFunction first saves negligible time, when the argument is a list of points:

RegionNearest[reg, pts]; // RepeatedTiming (* {0.565914, Null} *)

Thus we see that rnf[#] & /@ pts, rnf[pts], as well as RegionNearest[reg, pts] are equivalently fast, and RegionNearest[reg, #] & /@ pts is slower.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.