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 import java.util.ArrayList; 19 import java.util.HashMap; 20 import java.util.List; 21 import java.util.Map; 22 23 public strictfp class S2CellTest extends GeometryTestCase { 24 25 public static final boolean DEBUG_MODE = true; 26 testFaces()27 public void testFaces() { 28 Map<S2Point, Integer> edgeCounts = new HashMap<S2Point, Integer>(); 29 Map<S2Point, Integer> vertexCounts = new HashMap<S2Point, Integer>(); 30 for (int face = 0; face < 6; ++face) { 31 S2CellId id = S2CellId.fromFacePosLevel(face, 0, 0); 32 S2Cell cell = new S2Cell(id); 33 assertEquals(cell.id(), id); 34 assertEquals(cell.face(), face); 35 assertEquals(cell.level(), 0); 36 // Top-level faces have alternating orientations to get RHS coordinates. 37 assertEquals(cell.orientation(), face & S2.SWAP_MASK); 38 assertTrue(!cell.isLeaf()); 39 for (int k = 0; k < 4; ++k) { 40 if (edgeCounts.containsKey(cell.getEdgeRaw(k))) { 41 edgeCounts.put(cell.getEdgeRaw(k), edgeCounts.get(cell 42 .getEdgeRaw(k)) + 1); 43 } else { 44 edgeCounts.put(cell.getEdgeRaw(k), 1); 45 } 46 47 if (vertexCounts.containsKey(cell.getVertexRaw(k))) { 48 vertexCounts.put(cell.getVertexRaw(k), vertexCounts.get(cell 49 .getVertexRaw(k)) + 1); 50 } else { 51 vertexCounts.put(cell.getVertexRaw(k), 1); 52 } 53 assertDoubleNear(cell.getVertexRaw(k).dotProd(cell.getEdgeRaw(k)), 0); 54 assertDoubleNear(cell.getVertexRaw((k + 1) & 3).dotProd( 55 cell.getEdgeRaw(k)), 0); 56 assertDoubleNear(S2Point.normalize( 57 S2Point.crossProd(cell.getVertexRaw(k), cell 58 .getVertexRaw((k + 1) & 3))).dotProd(cell.getEdge(k)), 1.0); 59 } 60 } 61 // Check that edges have multiplicity 2 and vertices have multiplicity 3. 62 for (Integer i : edgeCounts.values()) { 63 assertEquals(i.intValue(), 2); 64 } 65 for (Integer i : vertexCounts.values()) { 66 assertEquals(i.intValue(), 3); 67 } 68 } 69 70 static class LevelStats { 71 double count; 72 double minArea, maxArea, avgArea; 73 double minWidth, maxWidth, avgWidth; 74 double minEdge, maxEdge, avgEdge, maxEdgeAspect; 75 double minDiag, maxDiag, avgDiag, maxDiagAspect; 76 double minAngleSpan, maxAngleSpan, avgAngleSpan; 77 double minApproxRatio, maxApproxRatio; 78 LevelStats()79 LevelStats() { 80 count = 0; 81 minArea = 100; 82 maxArea = 0; 83 avgArea = 0; 84 minWidth = 100; 85 maxWidth = 0; 86 avgWidth = 0; 87 minEdge = 100; 88 maxEdge = 0; 89 avgEdge = 0; 90 maxEdgeAspect = 0; 91 minDiag = 100; 92 maxDiag = 0; 93 avgDiag = 0; 94 maxDiagAspect = 0; 95 minAngleSpan = 100; 96 maxAngleSpan = 0; 97 avgAngleSpan = 0; 98 minApproxRatio = 100; 99 maxApproxRatio = 0; 100 } 101 } 102 103 static List<LevelStats> levelStats = new ArrayList<LevelStats>( 104 S2CellId.MAX_LEVEL + 1); 105 106 static { 107 for (int i = 0; i < S2CellId.MAX_LEVEL + 1; ++i) { levelStats.add(new LevelStats())108 levelStats.add(new LevelStats()); 109 } 110 } 111 gatherStats(S2Cell cell)112 static void gatherStats(S2Cell cell) { 113 LevelStats s = levelStats.get(cell.level()); 114 double exactArea = cell.exactArea(); 115 double approxArea = cell.approxArea(); 116 double minEdge = 100, maxEdge = 0, avgEdge = 0; 117 double minDiag = 100, maxDiag = 0; 118 double minWidth = 100, maxWidth = 0; 119 double minAngleSpan = 100, maxAngleSpan = 0; 120 for (int i = 0; i < 4; ++i) { 121 double edge = cell.getVertexRaw(i).angle(cell.getVertexRaw((i + 1) & 3)); 122 minEdge = Math.min(edge, minEdge); 123 maxEdge = Math.max(edge, maxEdge); 124 avgEdge += 0.25 * edge; 125 S2Point mid = S2Point.add(cell.getVertexRaw(i), cell 126 .getVertexRaw((i + 1) & 3)); 127 double width = S2.M_PI_2 - mid.angle(cell.getEdgeRaw(i ^ 2)); 128 minWidth = Math.min(width, minWidth); 129 maxWidth = Math.max(width, maxWidth); 130 if (i < 2) { 131 double diag = cell.getVertexRaw(i).angle(cell.getVertexRaw(i ^ 2)); 132 minDiag = Math.min(diag, minDiag); 133 maxDiag = Math.max(diag, maxDiag); 134 double angleSpan = cell.getEdgeRaw(i).angle( 135 S2Point.neg(cell.getEdgeRaw(i ^ 2))); 136 minAngleSpan = Math.min(angleSpan, minAngleSpan); 137 maxAngleSpan = Math.max(angleSpan, maxAngleSpan); 138 } 139 } 140 s.count += 1; 141 s.minArea = Math.min(exactArea, s.minArea); 142 s.maxArea = Math.max(exactArea, s.maxArea); 143 s.avgArea += exactArea; 144 s.minWidth = Math.min(minWidth, s.minWidth); 145 s.maxWidth = Math.max(maxWidth, s.maxWidth); 146 s.avgWidth += 0.5 * (minWidth + maxWidth); 147 s.minEdge = Math.min(minEdge, s.minEdge); 148 s.maxEdge = Math.max(maxEdge, s.maxEdge); 149 s.avgEdge += avgEdge; 150 s.maxEdgeAspect = Math.max(maxEdge / minEdge, s.maxEdgeAspect); 151 s.minDiag = Math.min(minDiag, s.minDiag); 152 s.maxDiag = Math.max(maxDiag, s.maxDiag); 153 s.avgDiag += 0.5 * (minDiag + maxDiag); 154 s.maxDiagAspect = Math.max(maxDiag / minDiag, s.maxDiagAspect); 155 s.minAngleSpan = Math.min(minAngleSpan, s.minAngleSpan); 156 s.maxAngleSpan = Math.max(maxAngleSpan, s.maxAngleSpan); 157 s.avgAngleSpan += 0.5 * (minAngleSpan + maxAngleSpan); 158 double approxRatio = approxArea / exactArea; 159 s.minApproxRatio = Math.min(approxRatio, s.minApproxRatio); 160 s.maxApproxRatio = Math.max(approxRatio, s.maxApproxRatio); 161 } 162 testSubdivide(S2Cell cell)163 public void testSubdivide(S2Cell cell) { 164 gatherStats(cell); 165 if (cell.isLeaf()) { 166 return; 167 } 168 169 S2Cell[] children = new S2Cell[4]; 170 for (int i = 0; i < children.length; ++i) { 171 children[i] = new S2Cell(); 172 } 173 assertTrue(cell.subdivide(children)); 174 S2CellId childId = cell.id().childBegin(); 175 double exactArea = 0; 176 double approxArea = 0; 177 double averageArea = 0; 178 for (int i = 0; i < 4; ++i, childId = childId.next()) { 179 exactArea += children[i].exactArea(); 180 approxArea += children[i].approxArea(); 181 averageArea += children[i].averageArea(); 182 183 // Check that the child geometry is consistent with its cell id. 184 assertEquals(children[i].id(), childId); 185 assertTrue(children[i].getCenter().aequal(childId.toPoint(), 1e-15)); 186 S2Cell direct = new S2Cell(childId); 187 assertEquals(children[i].face(), direct.face()); 188 assertEquals(children[i].level(), direct.level()); 189 assertEquals(children[i].orientation(), direct.orientation()); 190 assertEquals(children[i].getCenterRaw(), direct.getCenterRaw()); 191 for (int k = 0; k < 4; ++k) { 192 assertEquals(children[i].getVertexRaw(k), direct.getVertexRaw(k)); 193 assertEquals(children[i].getEdgeRaw(k), direct.getEdgeRaw(k)); 194 } 195 196 // Test Contains() and MayIntersect(). 197 assertTrue(cell.contains(children[i])); 198 assertTrue(cell.mayIntersect(children[i])); 199 assertTrue(!children[i].contains(cell)); 200 assertTrue(cell.contains(children[i].getCenterRaw())); 201 for (int j = 0; j < 4; ++j) { 202 assertTrue(cell.contains(children[i].getVertexRaw(j))); 203 if (j != i) { 204 assertTrue(!children[i].contains(children[j].getCenterRaw())); 205 assertTrue(!children[i].mayIntersect(children[j])); 206 } 207 } 208 209 // Test GetCapBound and GetRectBound. 210 S2Cap parentCap = cell.getCapBound(); 211 S2LatLngRect parentRect = cell.getRectBound(); 212 if (cell.contains(new S2Point(0, 0, 1)) 213 || cell.contains(new S2Point(0, 0, -1))) { 214 assertTrue(parentRect.lng().isFull()); 215 } 216 S2Cap childCap = children[i].getCapBound(); 217 S2LatLngRect childRect = children[i].getRectBound(); 218 assertTrue(childCap.contains(children[i].getCenter())); 219 assertTrue(childRect.contains(children[i].getCenterRaw())); 220 assertTrue(parentCap.contains(children[i].getCenter())); 221 assertTrue(parentRect.contains(children[i].getCenterRaw())); 222 for (int j = 0; j < 4; ++j) { 223 assertTrue(childCap.contains(children[i].getVertex(j))); 224 assertTrue(childRect.contains(children[i].getVertex(j))); 225 assertTrue(childRect.contains(children[i].getVertexRaw(j))); 226 assertTrue(parentCap.contains(children[i].getVertex(j))); 227 if (!parentRect.contains(children[i].getVertex(j))) { 228 System.out.println("cell: " + cell + " i: " + i + " j: " + j); 229 System.out.println("Children " + i + ": " + children[i]); 230 System.out.println("Parent rect: " + parentRect); 231 System.out.println("Vertex raw(j) " + children[i].getVertex(j)); 232 System.out.println("Latlng of vertex: " + new S2LatLng(children[i].getVertex(j))); 233 cell.getRectBound(); 234 } 235 assertTrue(parentRect.contains(children[i].getVertex(j))); 236 if (!parentRect.contains(children[i].getVertexRaw(j))) { 237 System.out.println("cell: " + cell + " i: " + i + " j: " + j); 238 System.out.println("Children " + i + ": " + children[i]); 239 System.out.println("Parent rect: " + parentRect); 240 System.out.println("Vertex raw(j) " + children[i].getVertexRaw(j)); 241 System.out.println("Latlng of vertex: " + new S2LatLng(children[i].getVertexRaw(j))); 242 cell.getRectBound(); 243 } 244 assertTrue(parentRect.contains(children[i].getVertexRaw(j))); 245 if (j != i) { 246 // The bounding caps and rectangles should be tight enough so that 247 // they exclude at least two vertices of each adjacent cell. 248 int capCount = 0; 249 int rectCount = 0; 250 for (int k = 0; k < 4; ++k) { 251 if (childCap.contains(children[j].getVertex(k))) { 252 ++capCount; 253 } 254 if (childRect.contains(children[j].getVertexRaw(k))) { 255 ++rectCount; 256 } 257 } 258 assertTrue(capCount <= 2); 259 if (childRect.latLo().radians() > -S2.M_PI_2 260 && childRect.latHi().radians() < S2.M_PI_2) { 261 // Bounding rectangles may be too large at the poles because the 262 // pole itself has an arbitrary fixed longitude. 263 assertTrue(rectCount <= 2); 264 } 265 } 266 } 267 268 // Check all children for the first few levels, and then sample randomly. 269 // Also subdivide one corner cell, one edge cell, and one center cell 270 // so that we have a better chance of sample the minimum metric values. 271 boolean forceSubdivide = false; 272 S2Point center = S2Projections.getNorm(children[i].face()); 273 S2Point edge = S2Point.add(center, S2Projections.getUAxis(children[i].face())); 274 S2Point corner = S2Point.add(edge, S2Projections.getVAxis(children[i].face())); 275 for (int j = 0; j < 4; ++j) { 276 S2Point p = children[i].getVertexRaw(j); 277 if (p.equals(center) || p.equals(edge) || p.equals(corner)) { 278 forceSubdivide = true; 279 } 280 } 281 if (forceSubdivide || cell.level() < (DEBUG_MODE ? 5 : 6) 282 || random(DEBUG_MODE ? 10 : 4) == 0) { 283 testSubdivide(children[i]); 284 } 285 } 286 287 // Check sum of child areas equals parent area. 288 // 289 // For ExactArea(), the best relative error we can expect is about 1e-6 290 // because the precision of the unit vector coordinates is only about 1e-15 291 // and the edge length of a leaf cell is about 1e-9. 292 // 293 // For ApproxArea(), the areas are accurate to within a few percent. 294 // 295 // For AverageArea(), the areas themselves are not very accurate, but 296 // the average area of a parent is exactly 4 times the area of a child. 297 298 assertTrue(Math.abs(Math.log(exactArea / cell.exactArea())) <= Math 299 .abs(Math.log(1 + 1e-6))); 300 assertTrue(Math.abs(Math.log(approxArea / cell.approxArea())) <= Math 301 .abs(Math.log(1.03))); 302 assertTrue(Math.abs(Math.log(averageArea / cell.averageArea())) <= Math 303 .abs(Math.log(1 + 1e-15))); 304 } 305 testMinMaxAvg(String label, int level, double count, double absError, double minValue, double maxValue, double avgValue, S2.Metric minMetric, S2.Metric maxMetric, S2.Metric avgMetric)306 public void testMinMaxAvg(String label, int level, double count, 307 double absError, double minValue, double maxValue, double avgValue, 308 S2.Metric minMetric, S2.Metric maxMetric, S2.Metric avgMetric) { 309 310 // All metrics are minimums, maximums, or averages of differential 311 // quantities, and therefore will not be exact for cells at any finite 312 // level. The differential minimum is always a lower bound, and the maximum 313 // is always an upper bound, but these minimums and maximums may not be 314 // achieved for two different reasons. First, the cells at each level are 315 // sampled and we may miss the most extreme examples. Second, the actual 316 // metric for a cell is obtained by integrating the differential quantity, 317 // which is not constant across the cell. Therefore cells at low levels 318 // (bigger cells) have smaller variations. 319 // 320 // The "tolerance" below is an attempt to model both of these effects. 321 // At low levels, error is dominated by the variation of differential 322 // quantities across the cells, while at high levels error is dominated by 323 // the effects of random sampling. 324 double tolerance = (maxMetric.getValue(level) - minMetric.getValue(level)) 325 / Math.sqrt(Math.min(count, 0.5 * (1L << level))) * 10; 326 if (tolerance == 0) { 327 tolerance = absError; 328 } 329 330 double minError = minValue - minMetric.getValue(level); 331 double maxError = maxMetric.getValue(level) - maxValue; 332 double avgError = Math.abs(avgMetric.getValue(level) - avgValue); 333 System.out.printf( 334 "%-10s (%6.0f samples, tolerance %8.3g) - min (%9.3g : %9.3g) " 335 + "max (%9.3g : %9.3g), avg (%9.3g : %9.3g)\n", label, count, 336 tolerance, minError / minValue, minError / tolerance, maxError 337 / maxValue, maxError / tolerance, avgError / avgValue, avgError 338 / tolerance); 339 340 assertTrue(minMetric.getValue(level) <= minValue + absError); 341 assertTrue(minMetric.getValue(level) >= minValue - tolerance); 342 System.out.println("Level: " + maxMetric.getValue(level) + " max " + (maxValue + tolerance)); 343 assertTrue(maxMetric.getValue(level) <= maxValue + tolerance); 344 assertTrue(maxMetric.getValue(level) >= maxValue - absError); 345 assertDoubleNear(avgMetric.getValue(level), avgValue, 10 * tolerance); 346 } 347 testSubdivide()348 public void testSubdivide() { 349 for (int face = 0; face < 6; ++face) { 350 testSubdivide(S2Cell.fromFacePosLevel(face, (byte) 0, 0)); 351 } 352 353 // The maximum edge *ratio* is the ratio of the longest edge of any cell to 354 // the shortest edge of any cell at the same level (and similarly for the 355 // maximum diagonal ratio). 356 // 357 // The maximum edge *aspect* is the maximum ratio of the longest edge of a 358 // cell to the shortest edge of that same cell (and similarly for the 359 // maximum diagonal aspect). 360 361 System.out 362 .printf("Level Area Edge Diag Approx Average\n"); 363 System.out 364 .printf(" Ratio Ratio Aspect Ratio Aspect Min Max Min Max\n"); 365 for (int i = 0; i <= S2CellId.MAX_LEVEL; ++i) { 366 LevelStats s = levelStats.get(i); 367 if (s.count > 0) { 368 s.avgArea /= s.count; 369 s.avgWidth /= s.count; 370 s.avgEdge /= s.count; 371 s.avgDiag /= s.count; 372 s.avgAngleSpan /= s.count; 373 } 374 System.out.printf( 375 "%5d %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n", i, 376 s.maxArea / s.minArea, s.maxEdge / s.minEdge, s.maxEdgeAspect, 377 s.maxDiag / s.minDiag, s.maxDiagAspect, s.minApproxRatio, 378 s.maxApproxRatio, S2Cell.averageArea(i) / s.maxArea, S2Cell 379 .averageArea(i) 380 / s.minArea); 381 } 382 383 // Now check the validity of the S2 length and area metrics. 384 for (int i = 0; i <= S2CellId.MAX_LEVEL; ++i) { 385 LevelStats s = levelStats.get(i); 386 if (s.count == 0) { 387 continue; 388 } 389 390 System.out.printf( 391 "Level %2d - metric (error/actual : error/tolerance)\n", i); 392 393 // The various length calculations are only accurate to 1e-15 or so, 394 // so we need to allow for this amount of discrepancy with the theoretical 395 // minimums and maximums. The area calculation is accurate to about 1e-15 396 // times the cell width. 397 testMinMaxAvg("area", i, s.count, 1e-15 * s.minWidth, s.minArea, 398 s.maxArea, s.avgArea, S2Projections.MIN_AREA, S2Projections.MAX_AREA, 399 S2Projections.AVG_AREA); 400 testMinMaxAvg("width", i, s.count, 1e-15, s.minWidth, s.maxWidth, 401 s.avgWidth, S2Projections.MIN_WIDTH, S2Projections.MAX_WIDTH, 402 S2Projections.AVG_WIDTH); 403 testMinMaxAvg("edge", i, s.count, 1e-15, s.minEdge, s.maxEdge, 404 s.avgEdge, S2Projections.MIN_EDGE, S2Projections.MAX_EDGE, 405 S2Projections.AVG_EDGE); 406 testMinMaxAvg("diagonal", i, s.count, 1e-15, s.minDiag, s.maxDiag, 407 s.avgDiag, S2Projections.MIN_DIAG, S2Projections.MAX_DIAG, 408 S2Projections.AVG_DIAG); 409 testMinMaxAvg("angle span", i, s.count, 1e-15, s.minAngleSpan, 410 s.maxAngleSpan, s.avgAngleSpan, S2Projections.MIN_ANGLE_SPAN, 411 S2Projections.MAX_ANGLE_SPAN, S2Projections.AVG_ANGLE_SPAN); 412 413 // The aspect ratio calculations are ratios of lengths and are therefore 414 // less accurate at higher subdivision levels. 415 assertTrue(s.maxEdgeAspect <= S2Projections.MAX_EDGE_ASPECT + 1e-15 416 * (1 << i)); 417 assertTrue(s.maxDiagAspect <= S2Projections.MAX_DIAG_ASPECT + 1e-15 418 * (1 << i)); 419 } 420 } 421 422 static final int MAX_LEVEL = DEBUG_MODE ? 6 : 10; 423 expandChildren1(S2Cell cell)424 public void expandChildren1(S2Cell cell) { 425 S2Cell[] children = new S2Cell[4]; 426 assertTrue(cell.subdivide(children)); 427 if (children[0].level() < MAX_LEVEL) { 428 for (int pos = 0; pos < 4; ++pos) { 429 expandChildren1(children[pos]); 430 } 431 } 432 } 433 expandChildren2(S2Cell cell)434 public void expandChildren2(S2Cell cell) { 435 S2CellId id = cell.id().childBegin(); 436 for (int pos = 0; pos < 4; ++pos, id = id.next()) { 437 S2Cell child = new S2Cell(id); 438 if (child.level() < MAX_LEVEL) { 439 expandChildren2(child); 440 } 441 } 442 } 443 } 444