I am glad to be here. First off, please excuse my almost saddening lack of knowledge. Math (in general) isn't my strong point. I might ask some really basic questions, you have been warned :)

The background: I am building a (so-to-speak) galaxy simulator. I have generated an eliptical galaxy with a large number of "stars"; currently staying at 10.000 but the final number will be around 2 million. They are represented as dots in a 3D space, generated using spherical coordinates which are also transformed and available as cartesian coordinates.

The work: There are two things I'm working on right now:

  1. Let Sphere A with radius 10000 containing 10000 dots with known spherical/cartesian coordinates. Within sphere A let Sphere B with radius R, centered around coordinates X, Y, Z. Find out all dots located both within Sphere A and Sphere B. This one I solved by using a cube with size 2R, centered at X, Y, Z, then finding all dots located between x-r and x+r, y-r and y+r, z-r and z+r. Then I calculated the distance between the center of the cube and each of the dots. If the distance is greater than R, then they are ignored, leaving only the dots which are located within the sphere inscribed in the cube. I am, however, interested to find whether there is a more elegant solution which wouldn't involve the "cube" middle step.

  2. Let Sphere A with radius 10000 containing 10000 dots with known spherical/cartesian coordinates. Within sphere A let a dot with known cartesian coordinates X, Y, Z. From this dot we generate a spherical sector B with radius R and cone angle phi (which is between 30 and 90 degrees). Find out all known dots which are located both within sphere A and spherical sector B. This one I couldn't solve. It is further complicated by the fact that it doesn't necessarily start from the center of sphere A, and its angle relative to sphere A could point pretty much anywhere.

To make an analogy, imagine you're located on Earth and you look at the sky through a telescope. You could point the telescope anywhere. You need to find only stars which are within your field of view, given the fact that you know all star coordinates relative to the center of the galaxy.

How I approached this second problem:

I generated a "virtual" sphere (S1) centered on coordinates X, Y, Z and found all dots located within this sphere using solution outlined at point 1. Then I converted all the cartesian coordinates to "local" coordinates. But then I kind of got stuck... I can't figure a reliable way to limit the volume to only cover the spherical sector.

Please let me know if that makes sense... I'll try to give more details if needed.

Additional information

Thank you for your answer, however...

The way to find the stars in some region is to loop through the stars and for each one evaluate a formula and do a test to see if it is in the region or not.

I can't evaluate two million points every time I want to find what's in a region that might find 5 of them. Combined with the fact that such a query needs to finish really fast (under 500 ms worst case scenario, preferably within 100 ms), it becomes obvious that evaluating 2 million entries every time is really out of the question.

My "sphere within a sphere" PL/SQL procedure gives an answer in about 60 ms for "maximum" thresholds of r (which is about 5% of maximum radius) so I'm good there.

As for the spherical sector, I understand most of the formulas presented, however they all gravitate around comparing all those 2 million points, which puts a lot of strain on the queries. Isn't there a way to reduce the number of checks?

The spherical sector (which is a bit more than a cone) will have the following inputs: - cone axis (theta, phi, r) - cone angle width (in degrees): 30, 45, 60, 75 or 90 degrees.

That is, the user will enter, for example: Find_stars(10, 10, 1000, 45) which means "Find all stars located in the spherical sector with axis at theta=10, phi=10, length of 1000 and width of 45 degrees", and he will be presented with a list of entries showing all stars located within those bounds.


The way to find the stars in some region is to loop through the stars and for each one evaluate a formula and do a test to see if it is in the region or not. There's not much point in using a cube first; on modern computers multiplications are cheaper than "if" statements!

So the important step is to find the formula for the shapes you are interested in.

The interior of a sphere centred on the origin with radius $r$ satisfies this inequality:

$$x^2+y^2+z^2 \leq r^2$$

This is just Pythagoras's theorem, applied to find the distance from the origin to the point $(x, y, z)$.

If you want the sphere centred elsewhere or with a different radius, you can just substitute for $x$, $y$, $z$:

$$x = x' - x_0$$ $$y = y' - y_0$$ $$z = z' - z_0$$

This gives:

$$(x'-x_0)^2 + (y'-y_0)^2 + (z'-z_0)^2 \leq r^2$$

In this formula, $(x_0, y_0, z_0)$ is the centre of the sphere.

The formula for the interior of a double cone with its point at the origin and its axis along the z-axis and a half-angle of 45 degrees is:

$$x^2 + y^2 \leq z^2$$

Again, this is just Pythagoras's theorem. The left-hand side is the squared distance from $(x, y, z)$ to the point on the z-axis at $(0, 0, z)$.

Note that this gives a double cone; if we want just one of the cones we need in addition something like:

$$z \geq 0$$

You can change the half-angle of the cone to $\phi$ by scaling $z$:

$$x^2 + y^2 \leq (z \tan(\phi))^2$$

If you want to change the axis of the cone, you can do it as follows. First, rewrite the equation as follows:

$$x^2 + y^2 + z^2 \leq z^2 (\tan^2(\phi) + 1)$$

Now the left-hand side looks like the formula for a sphere, and does not depend on the orientation of the cone. The right-hand side does. We need to replace $z$ on the right-hand side with a quantity $d$ that measures the distance along the desired axis of the cone.

Suppose you want the cone to point towards $(u, v, w)$. Then a measure of the distance in that direction is the dot-product:

$$d = \frac{u}{l}x + \frac{v}{l}y + \frac{w}{l}z$$

Where $l$ is the length of the vector $(u, v, w)$, which can be computed (again) using Pythagoras's theorem:

$$l = \sqrt{u^2 + v^2 + w^2}$$

So, putting this all together:

$$(x^2 + y^2 + z^2)(u^2 + v^2 + w^2) \leq (wx + vy + wz)^2 (\tan^2(\phi) + 1)$$ $$ux + vy + wz \geq 0$$

Finally, you can move the point of the cone around using a coordinate transform, just as for the sphere.

  • $\begingroup$ Thank you; I have added more information to my original post to explain why checking all stars is unfeasible. We're getting somewhere, but we're probably not there yet. $\endgroup$ – war4peace Dec 5 '13 at 23:46
  • $\begingroup$ Okay, if you have some clever index that can use an approximate bounding box to reduce the number of starts that need to be considered, great! That means you won't have to test so many. But you have not given enough information for me to help with the indexing. Probably you should use whatever the database gives you. $\endgroup$ – apt1002 Dec 6 '13 at 3:16
  • $\begingroup$ If you want to limit the distance as well as the viewing cone, then you just need to test the stars against a sphere and a cone. $\endgroup$ – apt1002 Dec 6 '13 at 3:17
  • $\begingroup$ You can convert between your $\theta$ and $\phi$ and my $(u, v, w)$ using the following formulae: $$u = \cos(\phi)cos(\theta)$$ $$v = \sin(\phi)\cos(\theta)$$ $$w = sin(\theta)$$ This has the pleasant side-effect of guaranteeing $l = 1$. Note, I have used the geographers' convention that $\theta = 0$ means the equator; there is another convention in common use where it means the north pole. $\endgroup$ – apt1002 Dec 6 '13 at 3:21

One way to deal with a problem like this involving very large numbers of irregularly distributed points is to organize the points within a regular intermediate structure before making any queries.

For your two problems, you could use your Cartesian coordinates to build a box around the galaxy and subdivide that box into a three-dimensional grid containing many cubical cells. The boundaries between cells are formed by regularly spaced planes perpendicular to each axis; between two adjacent parallel planes you have a "slab" of cells arranged in rows and columns.

It is a simple task to sort two million points into an array of such cells. You can organize the array of cells in memory so that any cell is easily found by the coordinates of one of its corners. (The corner with minimum $x,$ $y,$ and $z$ coordinates is a convenient index.) To find the cell in which each point resides, you add suitable offsets to its $x,$ $y,$ and $z$ coordinates, divide each of the results by the width of a cell, and truncate away the fractional parts after division. You then have three integers that identify the cell in a three-dimensional array where the point belongs.

For a sphere of radius $R$ centered at $(X,Y,Z),$ you know that lower and upper bounds on the $z$ coordinates of any point in the sphere are $Z-R$ and $Z+R.$ Take each of those coordinates and apply the offset-divide-truncate algorithm to determine the range of slabs of cells that could possibly contain points within the sphere.

Then you visit each of these slabs one at a time. You compute the radius of the largest circle on the sphere contained in that slab. For slabs "below" the center of the sphere that circle is the intersection of the upper boundary of the slab and the sphere; you can find the radius by taking the $z$ coordinate of the upper surface and applying the formula $r = \sqrt{R^2 - (Z - z)^2}.$ For slabs "above" the center of the sphere you make a similar calculation using the lower surface of the slab. For the slab containing the center of the sphere the radius is simply $r = R.$

Now that you have a radius $r,$ you find the minimum and maximum $y$ coordinates of the intersection of the sphere with the slab: $Y - r$ and $Y + r.$ The offset-divide-truncate algorithm identifies the range of columns of cells within the slab that could possibly intersect the sphere. Then within each column you find the first and last cell that could possibly intersect the sphere: for columns to the "left" of the $y$ coordinate $Y$ you take the rightmost $y$ coordinate and apply the formula $(X - x) = \pm\sqrt{r^2 - (Y - y)^2}$ to obtain minimum and maximum $x$ coordinates, from which you can find the row numbers of the cells within that column that could possibly intersect the sphere.

If you are careful, you can determine not only which cells in each column have any intersection with the sphere, you can determine which cells are entirely with in the sphere.

For each column you then visit each cell that intersects the sphere. If the cell is entirely in the sphere you simply collect all the points in that cell; every one is in the sphere. If the cell intersects the sphere but is not completely contained, you apply the distance formula to each point to decide whether it is in the sphere.

By the time you have visited every intersecting cell in every row in every slab, you will have found all the points in the sphere.

Alternatively, if the mathematics of finding the intersecting columns in each slab and the intersecting cells in each column are more complicated than you care to implement, you can just take all the cells whose $x,$ $y$, and $z$ coordinates overlap the ranges $[X-R, X+R],$ $[Y-R, Y+R],$ and $[Z-R, Z+R],$ and visit every one of those cells. This adds some computation because you are inspecting some cells that cannot possibly contain any points in the sphere and because you are checking the distance to each individual point (not just the ones in cells that are not completely contained), but as long as the sphere is smaller than the galaxy you are not visiting all two million points.

This is similar to your intermediate cube algorithm in that it involves an intermediate data structure which involves a cube. But unlike your custom-made cube, the cubes in this structure were determined just once during the entire run of your program. When the program receives the center and radius of a sphere, it does not have to examine every one of the two million points even to determine if it is in the bounding cube of the sphere; instead, you simply walk through a subset of the cells in the precomputed array of cubes and visit only the points in those cells.

The computation for a sector follows a similar idea, though it is slightly more complicated because you are dealing with both a sphere around a point and a cone with a vertex at that point. Consider the angle between the axis of the cone and the lateral surface of the cone (which may be $\phi$ or $\phi/2$ depending on how you define "cone angle"). If the maximum value of that angle is $45$ degrees then it is always possible to choose an axis such that the cone's intersection with any slab boundary perpendicular to that axis is always an ellipse. You find the formula of the cone as a quadratic equation in $x,$ $y,$ and $z.$ Then if (for example) the slabs for this particular cone are perpendicular to the $z$ axis, you take the $z$ coordinate of each slab boundary and plug it into the equation, obtaining an equation of an ellipse with variables $x$ and $y.$ You can then use that equation to determine the minimum and maximum $y$ coordinate of the ellipse, find the columns intersecting the ellipse, and find the cells in each column intersecting the ellipse. Unlike the calculation for the sphere, for each slab you have to consider the ellipse at both its upper and lower boundary plane, and you also have to consider the two lines through the coordinates $(X,Y)$ (from the vertex of the cone) tangent to both ellipses. Altogether the calculations could be rather complicated to work out and you might prefer to just figure out the minimum and maximum $x$ and $y$ coordinates of the intersection of the cone with the slab (which you can get from the two ellipses) and visit all cells in the resulting rectangular block. You will also consider the sphere of radius $R$ around $(X,Y,Z)$ as before, which provides another boundary of the sector.

If the maximum angle between the axis and the lateral side of the cone is greater than $45$ degrees then you must also deal with cases where the intersection of the cone with a slab boundary is a hyperbola. You will also need to use the circular intersection of the sphere with the slab boundary in order to get upper and lower bounds on the coordinates.

Alternatively you can do some mathematics to determine the minimum and maximum of each of the coordinates of the sector, and simply iterate through all cells that overlap those ranges of coordinates. To find the minimum and maximum coordinates you need merely consider the vertex of the cone, the circle where the cone intersects the sphere, and the points where any axis intersects the sphere inside the cone. For a narrow cone in a diagonal direction this alternative can mean that most of the cells you visit do not intersect the cone at all, but if $R$ is smaller than the diameter of the galaxy it can at least let you avoid visiting most of the two million points. But this has the advantage that the mathematics for a wide cone is not more complicated than for a narrow cone.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.