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.logging.Logger; 19 20 public strictfp class S2Test extends GeometryTestCase { 21 22 private static Logger logger = Logger.getLogger(S2Test.class.getName()); 23 24 // Test helper methods for testing the traversal order. swapAxes(int ij)25 private static int swapAxes(int ij) { 26 return ((ij >> 1) & 1) + ((ij & 1) << 1); 27 } 28 invertBits(int ij)29 private static int invertBits(int ij) { 30 return ij ^ 3; 31 } 32 testTraversalOrder()33 public void testTraversalOrder() { 34 for (int r = 0; r < 4; ++r) { 35 for (int i = 0; i < 4; ++i) { 36 // Check consistency with respect to swapping axes. 37 assertEquals(S2.ijToPos(r, i), 38 S2.ijToPos(r ^ S2.SWAP_MASK, swapAxes(i))); 39 assertEquals(S2.posToIJ(r, i), 40 swapAxes(S2.posToIJ(r ^ S2.SWAP_MASK, i))); 41 42 // Check consistency with respect to reversing axis directions. 43 assertEquals(S2.ijToPos(r, i), 44 S2.ijToPos(r ^ S2.INVERT_MASK, invertBits(i))); 45 assertEquals(S2.posToIJ(r, i), 46 invertBits(S2.posToIJ(r ^ S2.INVERT_MASK, i))); 47 48 // Check that the two tables are inverses of each other. 49 assertEquals(S2.ijToPos(r, S2.posToIJ(r, i)), i); 50 assertEquals(S2.posToIJ(r, S2.ijToPos(r, i)), i); 51 } 52 } 53 } 54 testSTUV()55 public void testSTUV() { 56 // Check boundary conditions. 57 for (double x = -1; x <= 1; ++x) { 58 assertEquals(S2Projections.stToUV(x), x); 59 assertEquals(S2Projections.uvToST(x), x); 60 } 61 // Check that UVtoST and STtoUV are inverses. 62 for (double x = -1; x <= 1; x += 0.0001) { 63 assertDoubleNear(S2Projections.uvToST(S2Projections.stToUV(x)), x); 64 assertDoubleNear(S2Projections.stToUV(S2Projections.uvToST(x)), x); 65 } 66 } 67 testFaceUVtoXYZ()68 public void testFaceUVtoXYZ() { 69 // Check that each face appears exactly once. 70 S2Point sum = new S2Point(); 71 for (int face = 0; face < 6; ++face) { 72 S2Point center = S2Projections.faceUvToXyz(face, 0, 0); 73 assertEquals(S2Projections.getNorm(face), center); 74 assertEquals(Math.abs(center.get(center.largestAbsComponent())), 1.0); 75 sum = S2Point.add(sum, S2Point.fabs(center)); 76 } 77 assertEquals(sum, new S2Point(2, 2, 2)); 78 79 // Check that each face has a right-handed coordinate system. 80 for (int face = 0; face < 6; ++face) { 81 assertEquals( 82 S2Point.crossProd(S2Projections.getUAxis(face), S2Projections.getVAxis(face)).dotProd( 83 S2Projections.faceUvToXyz(face, 0, 0)), 1.0); 84 } 85 86 // Check that the Hilbert curves on each face combine to form a 87 // continuous curve over the entire cube. 88 for (int face = 0; face < 6; ++face) { 89 // The Hilbert curve on each face starts at (-1,-1) and terminates 90 // at either (1,-1) (if axes not swapped) or (-1,1) (if swapped). 91 int sign = ((face & S2.SWAP_MASK) != 0) ? -1 : 1; 92 assertEquals(S2Projections.faceUvToXyz(face, sign, -sign), 93 S2Projections.faceUvToXyz((face + 1) % 6, -1, -1)); 94 } 95 } 96 testUVNorms()97 public void testUVNorms() { 98 // Check that GetUNorm and GetVNorm compute right-handed normals for 99 // an edge in the increasing U or V direction. 100 for (int face = 0; face < 6; ++face) { 101 for (double x = -1; x <= 1; x += 1 / 1024.) { 102 assertDoubleNear( 103 S2Point.crossProd( 104 S2Projections.faceUvToXyz(face, x, -1), S2Projections.faceUvToXyz(face, x, 1)) 105 .angle(S2Projections.getUNorm(face, x)), 0); 106 assertDoubleNear( 107 S2Point.crossProd( 108 S2Projections.faceUvToXyz(face, -1, x), S2Projections.faceUvToXyz(face, 1, x)) 109 .angle(S2Projections.getVNorm(face, x)), 0); 110 } 111 } 112 } 113 testUVAxes()114 public void testUVAxes() { 115 // Check that axes are consistent with FaceUVtoXYZ. 116 for (int face = 0; face < 6; ++face) { 117 assertEquals(S2Projections.getUAxis(face), S2Point.sub( 118 S2Projections.faceUvToXyz(face, 1, 0), S2Projections.faceUvToXyz(face, 0, 0))); 119 assertEquals(S2Projections.getVAxis(face), S2Point.sub( 120 S2Projections.faceUvToXyz(face, 0, 1), S2Projections.faceUvToXyz(face, 0, 0))); 121 } 122 } 123 testAngleArea()124 public void testAngleArea() { 125 S2Point pz = new S2Point(0, 0, 1); 126 S2Point p000 = new S2Point(1, 0, 0); 127 S2Point p045 = new S2Point(1, 1, 0); 128 S2Point p090 = new S2Point(0, 1, 0); 129 S2Point p180 = new S2Point(-1, 0, 0); 130 assertDoubleNear(S2.angle(p000, pz, p045), S2.M_PI_4); 131 assertDoubleNear(S2.angle(p045, pz, p180), 3 * S2.M_PI_4); 132 assertDoubleNear(S2.angle(p000, pz, p180), S2.M_PI); 133 assertDoubleNear(S2.angle(pz, p000, pz), 0); 134 assertDoubleNear(S2.angle(pz, p000, p045), S2.M_PI_2); 135 136 assertDoubleNear(S2.area(p000, p090, pz), S2.M_PI_2); 137 assertDoubleNear(S2.area(p045, pz, p180), 3 * S2.M_PI_4); 138 139 // Make sure that area() has good *relative* accuracy even for 140 // very small areas. 141 final double eps = 1e-10; 142 S2Point pepsx = new S2Point(eps, 0, 1); 143 S2Point pepsy = new S2Point(0, eps, 1); 144 double expected1 = 0.5 * eps * eps; 145 assertDoubleNear(S2.area(pepsx, pepsy, pz), expected1, 1e-14 * expected1); 146 147 // Make sure that it can handle degenerate triangles. 148 S2Point pr = new S2Point(0.257, -0.5723, 0.112); 149 S2Point pq = new S2Point(-0.747, 0.401, 0.2235); 150 assertEquals(S2.area(pr, pr, pr), 0.0); 151 // TODO: The following test is not exact in optimized mode because the 152 // compiler chooses to mix 64-bit and 80-bit intermediate results. 153 assertDoubleNear(S2.area(pr, pq, pr), 0); 154 assertEquals(S2.area(p000, p045, p090), 0.0); 155 156 double maxGirard = 0; 157 for (int i = 0; i < 10000; ++i) { 158 S2Point p0 = randomPoint(); 159 S2Point d1 = randomPoint(); 160 S2Point d2 = randomPoint(); 161 S2Point p1 = S2Point.add(p0, S2Point.mul(d1, 1e-15)); 162 S2Point p2 = S2Point.add(p0, S2Point.mul(d2, 1e-15)); 163 // The actual displacement can be as much as 1.2e-15 due to roundoff. 164 // This yields a maximum triangle area of about 0.7e-30. 165 assertTrue(S2.area(p0, p1, p2) < 0.7e-30); 166 maxGirard = Math.max(maxGirard, S2.girardArea(p0, p1, p2)); 167 } 168 logger.info("Worst case Girard for triangle area 1e-30: " + maxGirard); 169 170 // Try a very long and skinny triangle. 171 S2Point p045eps = new S2Point(1, 1, eps); 172 double expected2 = 5.8578643762690495119753e-11; // Mathematica. 173 assertDoubleNear(S2.area(p000, p045eps, p090), expected2, 1e-9 * expected2); 174 175 // Triangles with near-180 degree edges that sum to a quarter-sphere. 176 final double eps2 = 1e-10; 177 S2Point p000eps2 = new S2Point(1, 0.1 * eps2, eps2); 178 double quarterArea1 = 179 S2.area(p000eps2, p000, p090) + S2.area(p000eps2, p090, p180) + S2.area(p000eps2, p180, pz) 180 + S2.area(p000eps2, pz, p000); 181 assertDoubleNear(quarterArea1, S2.M_PI); 182 183 // Four other triangles that sum to a quarter-sphere. 184 S2Point p045eps2 = new S2Point(1, 1, eps2); 185 double quarterArea2 = 186 S2.area(p045eps2, p000, p090) + S2.area(p045eps2, p090, p180) + S2.area(p045eps2, p180, pz) 187 + S2.area(p045eps2, pz, p000); 188 assertDoubleNear(quarterArea2, S2.M_PI); 189 } 190 191 public void testCCW() { 192 S2Point a = new S2Point(0.72571927877036835, 0.46058825605889098, 0.51106749730504852); 193 S2Point b = new S2Point(0.7257192746638208, 0.46058826573818168, 0.51106749441312738); 194 S2Point c = new S2Point(0.72571927671709457, 0.46058826089853633, 0.51106749585908795); 195 assertTrue(S2.robustCCW(a, b, c) != 0); 196 } 197 198 // Note: obviously, I could have defined a bundle of metrics like this in the 199 // S2 class itself rather than just for testing. However, it's not clear that 200 // this is useful other than for testing purposes, and I find 201 // S2.kMinWidth.GetMaxLevel(width) to be slightly more readable than 202 // than S2.kWidth.min().GetMaxLevel(width). Also, there is no fundamental 203 // reason that we need to analyze the minimum, maximum, and average values of 204 // every metric; it would be perfectly reasonable to just define one of these. 205 206 class MetricBundle { 207 public MetricBundle(S2.Metric min, S2.Metric max, S2.Metric avg) { 208 this.min_ = min; 209 this.max_ = max; 210 this.avg_ = avg; 211 } 212 213 S2.Metric min_; 214 S2.Metric max_; 215 S2.Metric avg_; 216 } 217 218 public void testMinMaxAvg(MetricBundle bundle) { 219 assertTrue(bundle.min_.deriv() < bundle.avg_.deriv()); 220 assertTrue(bundle.avg_.deriv() < bundle.max_.deriv()); 221 } 222 223 public void testLessOrEqual(MetricBundle a, MetricBundle b) { 224 assertTrue(a.min_.deriv() <= b.min_.deriv()); 225 assertTrue(a.max_.deriv() <= b.max_.deriv()); 226 assertTrue(a.avg_.deriv() <= b.avg_.deriv()); 227 } 228 229 public void testMetrics() { 230 231 MetricBundle angleSpan = new MetricBundle( 232 S2Projections.MIN_ANGLE_SPAN, S2Projections.MAX_ANGLE_SPAN, S2Projections.AVG_ANGLE_SPAN); 233 MetricBundle width = 234 new MetricBundle(S2Projections.MIN_WIDTH, S2Projections.MAX_WIDTH, S2Projections.AVG_WIDTH); 235 MetricBundle edge = 236 new MetricBundle(S2Projections.MIN_EDGE, S2Projections.MAX_EDGE, S2Projections.AVG_EDGE); 237 MetricBundle diag = 238 new MetricBundle(S2Projections.MIN_DIAG, S2Projections.MAX_DIAG, S2Projections.AVG_DIAG); 239 MetricBundle area = 240 new MetricBundle(S2Projections.MIN_AREA, S2Projections.MAX_AREA, S2Projections.AVG_AREA); 241 242 // First, check that min <= avg <= max for each metric. 243 testMinMaxAvg(angleSpan); 244 testMinMaxAvg(width); 245 testMinMaxAvg(edge); 246 testMinMaxAvg(diag); 247 testMinMaxAvg(area); 248 249 // Check that the maximum aspect ratio of an individual cell is consistent 250 // with the global minimums and maximums. 251 assertTrue(S2Projections.MAX_EDGE_ASPECT >= 1.0); 252 assertTrue(S2Projections.MAX_EDGE_ASPECT 253 < S2Projections.MAX_EDGE.deriv() / S2Projections.MIN_EDGE.deriv()); 254 assertTrue(S2Projections.MAX_DIAG_ASPECT >= 1); 255 assertTrue(S2Projections.MAX_DIAG_ASPECT 256 < S2Projections.MAX_DIAG.deriv() / S2Projections.MIN_DIAG.deriv()); 257 258 // Check various conditions that are provable mathematically. 259 testLessOrEqual(width, angleSpan); 260 testLessOrEqual(width, edge); 261 testLessOrEqual(edge, diag); 262 263 assertTrue(S2Projections.MIN_AREA.deriv() 264 >= S2Projections.MIN_WIDTH.deriv() * S2Projections.MIN_EDGE.deriv() - 1e-15); 265 assertTrue(S2Projections.MAX_AREA.deriv() 266 < S2Projections.MAX_WIDTH.deriv() * S2Projections.MAX_EDGE.deriv() + 1e-15); 267 268 // GetMinLevelForLength() and friends have built-in assertions, we just need 269 // to call these functions to test them. 270 // 271 // We don't actually check that the metrics are correct here, e.g. that 272 // GetMinWidth(10) is a lower bound on the width of cells at level 10. 273 // It is easier to check these properties in s2cell_unittest, since 274 // S2Cell has methods to compute the cell vertices, etc. 275 276 for (int level = -2; level <= S2CellId.MAX_LEVEL + 3; ++level) { 277 double dWidth = (2 * S2Projections.MIN_WIDTH.deriv()) * Math.pow(2, -level); 278 if (level >= S2CellId.MAX_LEVEL + 3) { 279 dWidth = 0; 280 } 281 282 // Check boundary cases (exactly equal to a threshold value). 283 int expectedLevel = Math.max(0, Math.min(S2CellId.MAX_LEVEL, level)); 284 assertEquals(S2Projections.MIN_WIDTH.getMinLevel(dWidth), expectedLevel); 285 assertEquals(S2Projections.MIN_WIDTH.getMaxLevel(dWidth), expectedLevel); 286 assertEquals(S2Projections.MIN_WIDTH.getClosestLevel(dWidth), expectedLevel); 287 288 // Also check non-boundary cases. 289 assertEquals(S2Projections.MIN_WIDTH.getMinLevel(1.2 * dWidth), expectedLevel); 290 assertEquals(S2Projections.MIN_WIDTH.getMaxLevel(0.8 * dWidth), expectedLevel); 291 assertEquals(S2Projections.MIN_WIDTH.getClosestLevel(1.2 * dWidth), expectedLevel); 292 assertEquals(S2Projections.MIN_WIDTH.getClosestLevel(0.8 * dWidth), expectedLevel); 293 294 // Same thing for area1. 295 double area1 = (4 * S2Projections.MIN_AREA.deriv()) * Math.pow(4, -level); 296 if (level <= -3) { 297 area1 = 0; 298 } 299 assertEquals(S2Projections.MIN_AREA.getMinLevel(area1), expectedLevel); 300 assertEquals(S2Projections.MIN_AREA.getMaxLevel(area1), expectedLevel); 301 assertEquals(S2Projections.MIN_AREA.getClosestLevel(area1), expectedLevel); 302 assertEquals(S2Projections.MIN_AREA.getMinLevel(1.2 * area1), expectedLevel); 303 assertEquals(S2Projections.MIN_AREA.getMaxLevel(0.8 * area1), expectedLevel); 304 assertEquals(S2Projections.MIN_AREA.getClosestLevel(1.2 * area1), expectedLevel); 305 assertEquals(S2Projections.MIN_AREA.getClosestLevel(0.8 * area1), expectedLevel); 306 } 307 } 308 309 public void testExp() { 310 for (int i = 0; i < 10; ++i) { 311 assertEquals(i + 1, S2.exp(Math.pow(2, i))); 312 } 313 314 for (int i = 0; i < 10; ++i) { 315 assertEquals(i + 1, S2.exp(-Math.pow(2, i))); 316 } 317 318 assertEquals(0, S2.exp(0)); 319 assertEquals(2, S2.exp(3)); 320 assertEquals(3, S2.exp(5)); 321 } 322 } 323