• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright 2005 Google Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *     http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 package com.google.common.geometry;
17 
18 
19 /**
20  * An S2Cell is an S2Region object that represents a cell. Unlike S2CellIds, it
21  * supports efficient containment and intersection tests. However, it is also a
22  * more expensive representation.
23  *
24  */
25 
26 public final strictfp class S2Cell implements S2Region {
27 
28   private static final int MAX_CELL_SIZE = 1 << S2CellId.MAX_LEVEL;
29 
30   byte face;
31   byte level;
32   byte orientation;
33   S2CellId cellId;
34   double[][] uv = new double[2][2];
35 
36   /**
37    * Default constructor used only internally.
38    */
S2Cell()39   S2Cell() {
40   }
41 
42   /**
43    * An S2Cell always corresponds to a particular S2CellId. The other
44    * constructors are just convenience methods.
45    */
S2Cell(S2CellId id)46   public S2Cell(S2CellId id) {
47     init(id);
48   }
49 
50   // This is a static method in order to provide named parameters.
fromFacePosLevel(int face, byte pos, int level)51   public static S2Cell fromFacePosLevel(int face, byte pos, int level) {
52     return new S2Cell(S2CellId.fromFacePosLevel(face, pos, level));
53   }
54 
55   // Convenience methods.
S2Cell(S2Point p)56   public S2Cell(S2Point p) {
57     init(S2CellId.fromPoint(p));
58   }
59 
S2Cell(S2LatLng ll)60   public S2Cell(S2LatLng ll) {
61     init(S2CellId.fromLatLng(ll));
62   }
63 
64 
id()65   public S2CellId id() {
66     return cellId;
67   }
68 
face()69   public int face() {
70     return face;
71   }
72 
level()73   public byte level() {
74     return level;
75   }
76 
orientation()77   public byte orientation() {
78     return orientation;
79   }
80 
isLeaf()81   public boolean isLeaf() {
82     return level == S2CellId.MAX_LEVEL;
83   }
84 
getVertex(int k)85   public S2Point getVertex(int k) {
86     return S2Point.normalize(getVertexRaw(k));
87   }
88 
89   /**
90    * Return the k-th vertex of the cell (k = 0,1,2,3). Vertices are returned in
91    * CCW order. The points returned by GetVertexRaw are not necessarily unit
92    * length.
93    */
getVertexRaw(int k)94   public S2Point getVertexRaw(int k) {
95     // Vertices are returned in the order SW, SE, NE, NW.
96     return S2Projections.faceUvToXyz(face, uv[0][(k >> 1) ^ (k & 1)], uv[1][k >> 1]);
97   }
98 
getEdge(int k)99   public S2Point getEdge(int k) {
100     return S2Point.normalize(getEdgeRaw(k));
101   }
102 
getEdgeRaw(int k)103   public S2Point getEdgeRaw(int k) {
104     switch (k) {
105       case 0:
106         return S2Projections.getVNorm(face, uv[1][0]); // South
107       case 1:
108         return S2Projections.getUNorm(face, uv[0][1]); // East
109       case 2:
110         return S2Point.neg(S2Projections.getVNorm(face, uv[1][1])); // North
111       default:
112         return S2Point.neg(S2Projections.getUNorm(face, uv[0][0])); // West
113     }
114   }
115 
116   /**
117    * Return the inward-facing normal of the great circle passing through the
118    * edge from vertex k to vertex k+1 (mod 4). The normals returned by
119    * GetEdgeRaw are not necessarily unit length.
120    *
121    *  If this is not a leaf cell, set children[0..3] to the four children of
122    * this cell (in traversal order) and return true. Otherwise returns false.
123    * This method is equivalent to the following:
124    *
125    *  for (pos=0, id=child_begin(); id != child_end(); id = id.next(), ++pos)
126    * children[i] = S2Cell(id);
127    *
128    * except that it is more than two times faster.
129    */
subdivide(S2Cell children[])130   public boolean subdivide(S2Cell children[]) {
131     // This function is equivalent to just iterating over the child cell ids
132     // and calling the S2Cell constructor, but it is about 2.5 times faster.
133 
134     if (cellId.isLeaf()) {
135       return false;
136     }
137 
138     // Compute the cell midpoint in uv-space.
139     R2Vector uvMid = getCenterUV();
140 
141     // Create four children with the appropriate bounds.
142     S2CellId id = cellId.childBegin();
143     for (int pos = 0; pos < 4; ++pos, id = id.next()) {
144       S2Cell child = children[pos];
145       child.face = face;
146       child.level = (byte) (level + 1);
147       child.orientation = (byte) (orientation ^ S2.posToOrientation(pos));
148       child.cellId = id;
149       int ij = S2.posToIJ(orientation, pos);
150       for (int d = 0; d < 2; ++d) {
151         // The dimension 0 index (i/u) is in bit 1 of ij.
152         int m = 1 - ((ij >> (1 - d)) & 1);
153         child.uv[d][m] = uvMid.get(d);
154         child.uv[d][1 - m] = uv[d][1 - m];
155       }
156     }
157     return true;
158   }
159 
160   /**
161    * Return the direction vector corresponding to the center in (s,t)-space of
162    * the given cell. This is the point at which the cell is divided into four
163    * subcells; it is not necessarily the centroid of the cell in (u,v)-space or
164    * (x,y,z)-space. The point returned by GetCenterRaw is not necessarily unit
165    * length.
166    */
getCenter()167   public S2Point getCenter() {
168     return S2Point.normalize(getCenterRaw());
169   }
170 
getCenterRaw()171   public S2Point getCenterRaw() {
172     return cellId.toPointRaw();
173   }
174 
175   /**
176    * Return the center of the cell in (u,v) coordinates (see {@code
177    * S2Projections}). Note that the center of the cell is defined as the point
178    * at which it is recursively subdivided into four children; in general, it is
179    * not at the midpoint of the (u,v) rectangle covered by the cell
180    */
getCenterUV()181   public R2Vector getCenterUV() {
182     MutableInteger i = new MutableInteger(0);
183     MutableInteger j = new MutableInteger(0);
184     cellId.toFaceIJOrientation(i, j, null);
185     int cellSize = 1 << (S2CellId.MAX_LEVEL - level);
186 
187     // TODO(dbeaumont): Figure out a better naming of the variables here (and elsewhere).
188     int si = (i.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE;
189     double x = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * si);
190 
191     int sj = (j.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE;
192     double y = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sj);
193 
194     return new R2Vector(x, y);
195   }
196 
197   /**
198    * Return the average area for cells at the given level.
199    */
averageArea(int level)200   public static double averageArea(int level) {
201     return S2Projections.AVG_AREA.getValue(level);
202   }
203 
204   /**
205    * Return the average area of cells at this level. This is accurate to within
206    * a factor of 1.7 (for S2_QUADRATIC_PROJECTION) and is extremely cheap to
207    * compute.
208    */
averageArea()209   public double averageArea() {
210     return averageArea(level);
211   }
212 
213   /**
214    * Return the approximate area of this cell. This method is accurate to within
215    * 3% percent for all cell sizes and accurate to within 0.1% for cells at
216    * level 5 or higher (i.e. 300km square or smaller). It is moderately cheap to
217    * compute.
218    */
approxArea()219   public double approxArea() {
220 
221     // All cells at the first two levels have the same area.
222     if (level < 2) {
223       return averageArea(level);
224     }
225 
226     // First, compute the approximate area of the cell when projected
227     // perpendicular to its normal. The cross product of its diagonals gives
228     // the normal, and the length of the normal is twice the projected area.
229     double flatArea = 0.5 * S2Point.crossProd(
230         S2Point.sub(getVertex(2), getVertex(0)), S2Point.sub(getVertex(3), getVertex(1))).norm();
231 
232     // Now, compensate for the curvature of the cell surface by pretending
233     // that the cell is shaped like a spherical cap. The ratio of the
234     // area of a spherical cap to the area of its projected disc turns out
235     // to be 2 / (1 + sqrt(1 - r*r)) where "r" is the radius of the disc.
236     // For example, when r=0 the ratio is 1, and when r=1 the ratio is 2.
237     // Here we set Pi*r*r == flat_area to find the equivalent disc.
238     return flatArea * 2 / (1 + Math.sqrt(1 - Math.min(S2.M_1_PI * flatArea, 1.0)));
239   }
240 
241   /**
242    * Return the area of this cell as accurately as possible. This method is more
243    * expensive but it is accurate to 6 digits of precision even for leaf cells
244    * (whose area is approximately 1e-18).
245    */
exactArea()246   public double exactArea() {
247     S2Point v0 = getVertex(0);
248     S2Point v1 = getVertex(1);
249     S2Point v2 = getVertex(2);
250     S2Point v3 = getVertex(3);
251     return S2.area(v0, v1, v2) + S2.area(v0, v2, v3);
252   }
253 
254   // //////////////////////////////////////////////////////////////////////
255   // S2Region interface (see {@code S2Region} for details):
256 
257   @Override
clone()258   public S2Region clone() {
259     S2Cell clone = new S2Cell();
260     clone.face = this.face;
261     clone.level = this.level;
262     clone.orientation = this.orientation;
263     clone.uv = this.uv.clone();
264 
265     return clone;
266   }
267 
268   @Override
getCapBound()269   public S2Cap getCapBound() {
270     // Use the cell center in (u,v)-space as the cap axis. This vector is
271     // very close to GetCenter() and faster to compute. Neither one of these
272     // vectors yields the bounding cap with minimal surface area, but they
273     // are both pretty close.
274     //
275     // It's possible to show that the two vertices that are furthest from
276     // the (u,v)-origin never determine the maximum cap size (this is a
277     // possible future optimization).
278 
279     double u = 0.5 * (uv[0][0] + uv[0][1]);
280     double v = 0.5 * (uv[1][0] + uv[1][1]);
281     S2Cap cap = S2Cap.fromAxisHeight(S2Point.normalize(S2Projections.faceUvToXyz(face, u, v)), 0);
282     for (int k = 0; k < 4; ++k) {
283       cap = cap.addPoint(getVertex(k));
284     }
285     return cap;
286   }
287 
288   // We grow the bounds slightly to make sure that the bounding rectangle
289   // also contains the normalized versions of the vertices. Note that the
290   // maximum result magnitude is Pi, with a floating-point exponent of 1.
291   // Therefore adding or subtracting 2**-51 will always change the result.
292   private static final double MAX_ERROR = 1.0 / (1L << 51);
293 
294   // The 4 cells around the equator extend to +/-45 degrees latitude at the
295   // midpoints of their top and bottom edges. The two cells covering the
296   // poles extend down to +/-35.26 degrees at their vertices.
297   // adding kMaxError (as opposed to the C version) because of asin and atan2
298   // roundoff errors
299   private static final double POLE_MIN_LAT = Math.asin(Math.sqrt(1.0 / 3.0)) - MAX_ERROR;
300   // 35.26 degrees
301 
302 
303   @Override
getRectBound()304   public S2LatLngRect getRectBound() {
305     if (level > 0) {
306       // Except for cells at level 0, the latitude and longitude extremes are
307       // attained at the vertices. Furthermore, the latitude range is
308       // determined by one pair of diagonally opposite vertices and the
309       // longitude range is determined by the other pair.
310       //
311       // We first determine which corner (i,j) of the cell has the largest
312       // absolute latitude. To maximize latitude, we want to find the point in
313       // the cell that has the largest absolute z-coordinate and the smallest
314       // absolute x- and y-coordinates. To do this we look at each coordinate
315       // (u and v), and determine whether we want to minimize or maximize that
316       // coordinate based on the axis direction and the cell's (u,v) quadrant.
317       double u = uv[0][0] + uv[0][1];
318       double v = uv[1][0] + uv[1][1];
319       int i = S2Projections.getUAxis(face).z == 0 ? (u < 0 ? 1 : 0) : (u > 0 ? 1 : 0);
320       int j = S2Projections.getVAxis(face).z == 0 ? (v < 0 ? 1 : 0) : (v > 0 ? 1 : 0);
321 
322 
323       R1Interval lat = R1Interval.fromPointPair(getLatitude(i, j), getLatitude(1 - i, 1 - j));
324       lat = lat.expanded(MAX_ERROR).intersection(S2LatLngRect.fullLat());
325       if (lat.lo() == -S2.M_PI_2 || lat.hi() == S2.M_PI_2) {
326         return new S2LatLngRect(lat, S1Interval.full());
327       }
328       S1Interval lng = S1Interval.fromPointPair(getLongitude(i, 1 - j), getLongitude(1 - i, j));
329       return new S2LatLngRect(lat, lng.expanded(MAX_ERROR));
330     }
331 
332 
333     // The face centers are the +X, +Y, +Z, -X, -Y, -Z axes in that order.
334     // assert (S2Projections.getNorm(face).get(face % 3) == ((face < 3) ? 1 : -1));
335     switch (face) {
336       case 0:
337         return new S2LatLngRect(
338             new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-S2.M_PI_4, S2.M_PI_4));
339       case 1:
340         return new S2LatLngRect(
341             new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(S2.M_PI_4, 3 * S2.M_PI_4));
342       case 2:
343         return new S2LatLngRect(
344             new R1Interval(POLE_MIN_LAT, S2.M_PI_2), new S1Interval(-S2.M_PI, S2.M_PI));
345       case 3:
346         return new S2LatLngRect(
347             new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(3 * S2.M_PI_4, -3 * S2.M_PI_4));
348       case 4:
349         return new S2LatLngRect(
350             new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-3 * S2.M_PI_4, -S2.M_PI_4));
351       default:
352         return new S2LatLngRect(
353             new R1Interval(-S2.M_PI_2, -POLE_MIN_LAT), new S1Interval(-S2.M_PI, S2.M_PI));
354     }
355 
356   }
357 
358   @Override
mayIntersect(S2Cell cell)359   public boolean mayIntersect(S2Cell cell) {
360     return cellId.intersects(cell.cellId);
361   }
362 
contains(S2Point p)363   public boolean contains(S2Point p) {
364     // We can't just call XYZtoFaceUV, because for points that lie on the
365     // boundary between two faces (i.e. u or v is +1/-1) we need to return
366     // true for both adjacent cells.
367     R2Vector uvPoint = S2Projections.faceXyzToUv(face, p);
368     if (uvPoint == null) {
369       return false;
370     }
371     return (uvPoint.x() >= uv[0][0] && uvPoint.x() <= uv[0][1]
372         && uvPoint.y() >= uv[1][0] && uvPoint.y() <= uv[1][1]);
373   }
374 
375   // The point 'p' does not need to be normalized.
376   @Override
contains(S2Cell cell)377   public boolean contains(S2Cell cell) {
378     return cellId.contains(cell.cellId);
379   }
380 
init(S2CellId id)381   private void init(S2CellId id) {
382     cellId = id;
383     MutableInteger ij[] = new MutableInteger[2];
384     MutableInteger mOrientation = new MutableInteger(0);
385 
386     for (int d = 0; d < 2; ++d) {
387       ij[d] = new MutableInteger(0);
388     }
389 
390     face = (byte) id.toFaceIJOrientation(ij[0], ij[1], mOrientation);
391     orientation = (byte) mOrientation.intValue(); // Compress int to a byte.
392     level = (byte) id.level();
393     int cellSize = 1 << (S2CellId.MAX_LEVEL - level);
394     for (int d = 0; d < 2; ++d) {
395       // Compute the cell bounds in scaled (i,j) coordinates.
396       int sijLo = (ij[d].intValue() & -cellSize) * 2 - MAX_CELL_SIZE;
397       int sijHi = sijLo + cellSize * 2;
398       uv[d][0] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijLo);
399       uv[d][1] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijHi);
400     }
401   }
402 
403 
404   // Internal method that does the actual work in the constructors.
405 
getLatitude(int i, int j)406   private double getLatitude(int i, int j) {
407     S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]);
408     return Math.atan2(p.z, Math.sqrt(p.x * p.x + p.y * p.y));
409   }
410 
getLongitude(int i, int j)411   private double getLongitude(int i, int j) {
412     S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]);
413     return Math.atan2(p.y, p.x);
414   }
415 
416   // Return the latitude or longitude of the cell vertex given by (i,j),
417   // where "i" and "j" are either 0 or 1.
418 
419   @Override
toString()420   public String toString() {
421     return "[" + face + ", " + level + ", " + orientation + ", " + cellId + "]";
422   }
423 
424   @Override
hashCode()425   public int hashCode() {
426     int value = 17;
427     value = 37 * (37 * (37 * value + face) + orientation) + level;
428     return 37 * value + id().hashCode();
429   }
430 
431   @Override
equals(Object that)432   public boolean equals(Object that) {
433     if (that instanceof S2Cell) {
434       S2Cell thatCell = (S2Cell) that;
435       return this.face == thatCell.face && this.level == thatCell.level
436           && this.orientation == thatCell.orientation && this.cellId.equals(thatCell.cellId);
437     }
438     return false;
439   }
440 
441 }
442