An Area Preserving Projection from the Regular Octahedron to the Sphere

In this paper, we propose an area preserving bijective map from the regular octahedron to the unit sphere S, both centered at the origin. The construction scheme consists of two steps. First, each face Fi of the octahedron is mapped to a curved planar triangle Ti of the same area. Afterwards, each Ti is mapped onto the sphere using the inverse Lambert azimuthal equal area projection with respect to a certain point of S. The proposed map is then used to construct uniform and refinable grids on a sphere, starting from any triangular uniform and refinable grid on the triangular faces of the octahedron. Mathematics Subject Classification (2000). Primary 86A30; Secondary 85-08, 86-08.


Introduction
In many applications, especially in geosciences and astronomy, but also in computer vision, one is interested in simple, refinable grids on the sphere. In particular, one requires partitions of S 2 into regions of equal area and small diameter in order to avoid the distortions that often occur in statistical computations and function approximations using non-equal area partitions.
There exist already some constructions of equal area partitions of S 2 . Based on the construction by Zhou [12], Leopardi [4] derives a recursive zonal equal area partition of the unit sphere S d ⊂ R d+1 that consists of polar cups and rectilinear regions. This construction has the disadvantage that one has to deal with different kinds of areas. The existence of partitions of S 2 into regions of equal area and small diameter has been already used by Alexander [1], who derived lower bounds for the maximum sum of distances between This work has been supported by the German Research Foundation, grant PL 170/14-1.
points on the sphere. Tegmark [11] considered an icosahedron-based pixelizing of the sphere. Similarly, the algorithm for binning globally distributed measurements on the sphere by Teanby [10] is based on repeated subdivision of the icosahedron. In [8], a subdivision method is proposed using a spherical triangulation that is obtained by projecting the faces of an icosahedron to the sphere. The HEALPix (Hierarchical Equal Area iso-Latitude Pixelization) introduced in [3], has been frequently used in recent years for application of local pixel set operations, multiresolution applications and fast spherical harmonic transforms.
However, one simple method to construct grids on the sphere is to transfer existing planar grids by a suitable projection. For a survey of known spherical projections from the sphere or parts of the sphere to the plane we refer to [2,9]. Therefore we are especially interested in the construction of equal area partitions being obtained by application of an area-preserving bijective map from suitable planar domains to (parts of) the sphere.
In [5], one of the authors already suggested a new area preserving projection method based on a mapping of the square onto a disc in a first step, followed by a lifting to the sphere by the inverse Lambert projection. This idea can also be generalized to construct uniform and refinable grids on elliptic domains and on some surfaces of revolution, see [6]. In [7], the authors have recently constructed an equal area projection from the cube to the sphere.
In this paper we construct an area preserving map from a regular octahedron to the unit sphere S 2 . Thus, any grid on the octahedron can be transported to the sphere. Further, since arbitrary grids on the triangle can now simply be transported to the sphere, we believe that this construction may achieve an essential impact for different applications in geosciences. Since we give explicit formulas both for the map from the octahedron to the sphere, and from the sphere to the octahedron, the method is easy to implement.
The construction scheme consists of two steps. In the first step, we construct in Section 3 an area preserving bijection U from each face F i of the octahedron, onto a planar domain F i , bounded by a curved triangle T i . In the second step, we combine U with an inverse Lambert azimuthal projection, in order to map each face F i of the octahedron onto a subset F i of the sphere, A closed form of the inverse area preserving map U −1 is presented in Section 6. Further, we present some examples of the obtained spherical grids.

Preliminaries
Consider the unit sphere S 2 centered at the origin O and the regular octahedron K of the same area, centered at O and with vertices on the coordinate axes. Since the area of the sphere is 4π, the area of each face of K is π/2 and (2.1) We cut the sphere with the coordinate planes x = 0, y = 0, z = 0 and obtain the spherical triangles in Figure 1. Each face F i , i = 1, . . . , 8, of the octahedron K is thus situated in one of the following domains More precisely, F i ⊂ I i for i = 1, . . . , 8 and the portion of the sphere situated in I i will be denoted by F i , i.e., Figure 1.
We focus on the portion F 1 of the sphere and consider its center point The point G 1 is in fact the intersection of the line OG 1 with S 2 , where G 1 is the weight center of the face F 1 .
We consider the Lambert projection L G1 of the domain F 1 ⊂ S 2 with respect to the point G 1 . Thus, F 1 will be mapped on a curved triangle that is situated in the plane x + y + z = √ 3 perpendicular to the line OG 1 . For calculating this projection, we apply a transformation to the sphere such that G 1 is mapped onto the South Pole S = (0, 0, −1). Such a transformation R consists in two successive rotations, R 1 and R 2 . The first rotation R 1 is taken around the z-axis with angle −π/4 and maps G 1 onto P 1 = The second rotation R 2 is around the y-axis and maps P 1 onto S = (0, 0, −1), so its angle is β + π 2 , with β = arcsin 1 √ 3 . The corresponding rotation matrices of R 1 , R 2 and R have the form The matrices R 1 , R 2 and R are orthogonal, therefore R −1 = R T . Now, in order to obtain the Lambert projection L G1 with respect to G 1 ∈ S 2 , we use the azimuthal Lambert projection L S with respect to the South Pole S, as follows: It is well known that the image of a point (x, y, z) ∈ S 2 by L S onto the tangential plane at the South Pole is For the rotated Lambert projection on Z L = −1 we obtain Let us denote by l 1 , l 2 , l 3 the edges of the spherical triangle F 1 , situated on the planes z = 0, y = 0 and x = 0, respectively, see Figure 1. These edges are mapped by L S • R onto the planar curves in the plane Z L = −1, given by the following parametric equations We denote the curved triangle determined by these three curves by T 1 = T . It can be simply transferred to the plane x + y + z = √ 3 using the rotation map R −1 (and hence finishing the projection L G1 ). But we prefer to compare the curved triangle T with the rotated triangular face T = RT 1 of the octahedron in the plane Z L = −1. A simple calculation shows that the vertices of the equilateral triangle T = RT 1 are P 1 = (2α, 0, −1), and side length in (2.1). Figure 2 (left) shows the curved triangle T and the equilateral triangle T of the same area π/2. With the help of the projection L G1 (resp. L S • R), we have simplified now the problem of finding an area preserving map from the octahedron to the sphere S 2 to a problem of finding a two-dimensional map from an equilateral triangle T , i.e., a face of the octahedron, to the curved triangle T . In our further considerations in Section 3, we will consider only two dimensions while Z L = −1 is fixed.

Mapping a triangle onto a curved triangle
In this section we derive an area preserving bijection U : R 2 → R 2 , which maps the equilateral triangle T with vertices P 1 = (2α, 0), P 2 = (−α, /2), and P 3 = (−α, − /2) onto the curved triangle T that has been constructed in Section 2. Here, we say that the map U is area-preserving, if it has the property where A(D) denotes the area of D.
We consider the three half-lines h 1 = OP 1 , h 2 = OP 2 , h 3 = OP 3 , that form the angles 0, 2π/3, 4π/3 around O, see Figure 2(left), and determine three disjoint regions of R 2 defined by that are bounded by these half lines. We focus for the moment on the region Q + 1 defined by and from y P = ϕ(m)x P we have From this equality we can determine t P as , and further, if we replace t P in (3.1) and (3.2), we obtain the coordinates of P in the form .

Some simple calculations yield the Euklidean distances
where α is given in (2.5). In order to simplify the determination of the map U , we suppose that ON OP = OM OQ . From the above calculations we then obtain ON = x M α OP and hence .
Thus, the map U is now completely described by means of the function ϕ, and we obtain that U maps the point (x, y) ∈ Q + 1 onto the point (X, Y ) given by .

(3.4)
Next we have to ensure the area preserving property of U by a suitable determination of ϕ. For this purpose, we define the function ϕ such that 1 the Jacobian of U is 1. After simplification, the Jacobian writes as For solving the equation J(U ) = 1 we again substitute m := y x , and thus, in the considered case 0 ≤ y ≤ − √ 3x we have m ∈ [− √ 3, 0]. Hence, with the simplified notation ϕ = ϕ(m), we get The condition ϕ(0) = 0 yields C = 0. Next, in order to determine ϕ we use the formula arctan a − arctan b = arctan a−b 1+ab ∀a, b ∈ R, ab > −1, and we further obtain To simplify this term, we introduce the notation η = tan α 2 m 2 = tan πm 12 √ 3 , and equality (3.6) yields We must have ϕ > η and with this condition, (3.7) is equivalent with which gives To simplify the formulas, we denote δ = α 2 m 2 = πm 12 √

= πy
12 √ 3x and with η(1 + η 2 ) −1/2 = sin δ and (1 + η 2 ) −1/2 = cos δ we further calculate If we take into account the condition ϕ(− √ 3) = − √ 3 and the equality the only convenient solution is 3x . Finally, in order to give the explicit presentation of the area preserving map U , we need to replace ϕ in the equations (3.3)-(3.4). Some straightforward calculations show that Hence, for (x, y) ∈ Q + 1 , the formulas for the desired area preserving map U (x, y) = (X, Y ) are given by Similar arguments for the region Q − 1 = {(x, y) ∈ R 2 , x < 0, √ 3x ≤ y ≤ 0} show that these formulas also hold for (x, y) ∈ Q − 1 , so we can conclude that the image U (x, y) = (X, Y ) of a point (x, y) in the region is given by formulas (3.8)- (3.9). For the regions Q 2 and Q 3 we may use the formulas (3.8)-(3.9) after performing some rotations. More precisely, we need the rotation matrix We denote by U 1 , U 2 , U 3 the restrictions of U to Q 1 , Q 2 , Q 3 , respectively. Then, for (x, y) ∈ Q 2 we have x y = A −1 · x y ∈ Q 1 , therefore, to (x , y ) we can apply U defined by the formulas (3.8)-(3.9) and obtain'(X , Y ) ∈ Q 1 , Finally, if we apply to (X , Y ) a rotation of 2π 3 we obtain (X , Y ) ∈ Q 2 defined by In conclusion, the area preserving map U which maps equilateral triangles onto curved triangles can be written as follows: • For (x, y) . Figure 3 shows two examples of planar grids.

Mapping
The area preserving map U : R 2 → R 2 derived in Section 3 depends on two variables. We want to use it for mapping F i onto F i , where both F i and F i are included in the same plane of R 3 . We denote by P −1 the plane z = −1 and define U : P −1 → P −1 by This map allows us to write the formulas for the functions U i : F i → F i , for i = 1, . . . , 8, which map triangles (i.e., the faces of the octahedron) onto curved triangles. As shown in Section 2, we observe that in each of the cases (x, y, z) ∈ F i , i = 1, . . . , 8, where R is given in (2.2). Therefore, U (R · (|x|, |y|, |z|) T ) ∈ P −1 and So, we define U 1 : F 1 → F 1 , U 1 (x, y, z) = R −1 · U (R·(|x|, |y|, |z|) T ) =: (X, Y, Z), for (x, y, z) ∈ F 1 , (4.1) and further, with the notations in (4.1), the required applications U i : F i → F i are defined as follows: Finally, we denote by U :

Mapping the curved triangles onto the sphere
The complete mapping from the regular octahedron K to the sphere S 2 is now described in two steps. In the first step, each face F i of K will be mapped onto a planar domain F i , bounded by a curved triangle T i , using the transformation U constructed in the previous section. In the second step, each F i will be mapped onto F i ⊆ S 2 by the inverse Lambert azimuthal projection, with respect to G i = U (W i ), where W i is the weight center of F i . Obviously Applying the inverse Lambert projection L −1 S , the point (X, Y, −1), situated in the tangent plane to S 2 at the South Pole S, maps onto (x L , y L , z L ) ∈ S 2 given by Thus, the application L −1 Gi • U maps the face F i of the octahedron onto and for obtaining the other projections L −1 Gi we proceed as in the previous section. Again, R · (|x|, |y|, |z|) T ∈ P −1 , therefore we use formulas (5.1)-(5.3) and we have (x L , y L , z L ) = L −1 S (R · (|x|, |y|, |z|) T ). Further, (X, Y, Z) T := R −1 · (x L , y L , z L ) T ∈ F 1 , and then, with these notations, for the other cases we obtain the following expressions for L −1 Gi : F i → F i : G8 (x, y, z) = (−X, −Y, −Z). Figure 4 shows some grids of the sphere, where a regular partition of the faces of the octahedron into 4, 16, 25 and 64 equal area triangles has been applied.

The inverse map
To make the area preserving map U : K → ∪ 8 i=1 F i applicable in practice, we also need to derive a closed simple form for the inverse mapping U −1 . In view of the above considerations we can restrict to the calculation of the inverse U −1 of the area-preserving map proposed in Section 3, since the mapping of the portions F i toF i is obtained by the azimuthal Lambert projection, and the area preserving mapping of the curved triangleF i to the faces F i of the octahedron is completely understood, if we find U −1 : T → T with the two-dimensional triangles considered in Section 3, see Figure 2. Further, we can restrict the calculations again to the region Q 1 , and the complete map U −1 is then derived by rotation.

(6.2)
For simplicity, we introduce the notation w := Y /X. Since (X, Y ) ∈ Q 1 , we have w ∈ [− √ 3, √ 3] and δ = arccos and b = w √ 1+w 2 , and we obtain 3) For the calculation of x from X and Y , we use the equalities in (6.1) and (6.2) yielding Finally, from (6.3) we find y = |y|(sign Y ) with The obtained formulas for x and y form an explicit representation of the inverse map U −1 in the region Q 1 . For the other two regions, the map is directly obtained using the rotation matrix A in Section 3.