point_measures_test.go

Documentation: github.com/golang/geo/s2

     1  // Copyright 2018 Google Inc. All rights reserved.
     2  //
     3  // Licensed under the Apache License, Version 2.0 (the "License");
     4  // you may not use this file except in compliance with the License.
     5  // You may obtain a copy of the License at
     6  //
     7  //     http://www.apache.org/licenses/LICENSE-2.0
     8  //
     9  // Unless required by applicable law or agreed to in writing, software
    10  // distributed under the License is distributed on an "AS IS" BASIS,
    11  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    12  // See the License for the specific language governing permissions and
    13  // limitations under the License.
    14  
    15  package s2
    16  
    17  import (
    18  	"math"
    19  	"testing"
    20  
    21  	"github.com/golang/geo/r3"
    22  	"github.com/golang/geo/s1"
    23  )
    24  
    25  var (
    26  	pz   = Point{r3.Vector{0, 0, 1}}
    27  	p000 = Point{r3.Vector{1, 0, 0}}
    28  	p045 = Point{r3.Vector{1, 1, 0}.Normalize()}
    29  	p090 = Point{r3.Vector{0, 1, 0}}
    30  	p180 = Point{r3.Vector{-1, 0, 0}}
    31  	// Degenerate triangles.
    32  	pr = Point{r3.Vector{0.257, -0.5723, 0.112}}
    33  	pq = Point{r3.Vector{-0.747, 0.401, 0.2235}}
    34  
    35  	// For testing the Girard area fall through case.
    36  	g1 = Point{r3.Vector{1, 1, 1}}
    37  	g2 = Point{g1.Add(pr.Mul(1e-15)).Normalize()}
    38  	g3 = Point{g1.Add(pq.Mul(1e-15)).Normalize()}
    39  )
    40  
    41  func TestPointMeasuresPointArea(t *testing.T) {
    42  	const eps = 1e-10
    43  	const exp1 = 0.5 * eps * eps
    44  	const exp2 = 5.8578643762690495119753e-11
    45  	tests := []struct {
    46  		a, b, c  Point
    47  		want     float64
    48  		nearness float64
    49  	}{
    50  		{p000, p090, pz, math.Pi / 2.0, 0},
    51  		{p045, pz, p180, 3.0 * math.Pi / 4.0, 0},
    52  		// Make sure that Area has good *relative* accuracy even for very small areas.
    53  		{Point{r3.Vector{eps, 0, 1}.Normalize()}, Point{r3.Vector{0, eps, 1}.Normalize()}, pz, exp1, 1e-14 * exp1},
    54  		// Make sure that it can handle degenerate triangles.
    55  		{pr, pr, pr, 0.0, 0},
    56  		{pr, pq, pr, 0.0, 1e-15},
    57  		{p000, p045, p090, 0.0, 0},
    58  		// Try a very long and skinny triangle.
    59  		{p000, Point{r3.Vector{1, 1, eps}.Normalize()}, p090, exp2, 1e-9 * exp2},
    60  	}
    61  
    62  	for d, test := range tests {
    63  		if got := PointArea(test.a, test.b, test.c); !float64Near(got, test.want, test.nearness) {
    64  			t.Errorf("%d, PointArea(%v, %v, %v), got %v want %v", d, test.a, test.b, test.c, got, test.want)
    65  		}
    66  	}
    67  
    68  	maxGirard := 0.0
    69  	for i := 0; i < 10000; i++ {
    70  		p0 := randomPoint()
    71  		d1 := randomPoint()
    72  		d2 := randomPoint()
    73  		p1 := Point{p0.Add(d1.Mul(1e-15)).Normalize()}
    74  		p2 := Point{p0.Add(d2.Mul(1e-15)).Normalize()}
    75  		// The actual displacement can be as much as 1.2e-15 due to roundoff.
    76  		// This yields a maximum triangle area of about 0.7e-30.
    77  		if got := PointArea(p0, p1, p2); got > 0.7e-30 {
    78  			t.Errorf("PointArea(%v, %v, %v) = %v, want <= %v", p1, p1, p2, got, 0.7e-30)
    79  		}
    80  		if a := GirardArea(p0, p1, p2); a > maxGirard {
    81  			maxGirard = a
    82  		}
    83  	}
    84  	// This check only passes if GirardArea uses PointCross.
    85  	if maxGirard > 1e-14 {
    86  		t.Errorf("maximum GirardArea = %v, want <= %v", maxGirard, 1e-14)
    87  	}
    88  
    89  	// This tests a case where the triangle has zero area, but PointArea()
    90  	// computes (dmin > 0) due to rounding errors.
    91  	a := PointFromLatLng(LatLngFromDegrees(-45, -170))
    92  	b := PointFromLatLng(LatLngFromDegrees(45, -170))
    93  	c := PointFromLatLng(LatLngFromDegrees(0, -170))
    94  	if area := PointArea(a, b, c); area != 0.0 {
    95  		t.Errorf("PointArea(%v, %v, %v) = %v, want 0.0", a, b, c, area)
    96  	}
    97  }
    98  
    99  func TestPointMeasuresPointAreaQuarterHemisphere(t *testing.T) {
   100  	const eps2 = 1e-14
   101  	tests := []struct {
   102  		a, b, c, d, e Point
   103  		want          float64
   104  	}{
   105  		// Triangles with near-180 degree edges that sum to a quarter-sphere.
   106  		{PointFromCoords(1, 0.1*eps2, eps2), p000, p045, p180, pz, math.Pi},
   107  		// Four other triangles that sum to a quarter-sphere.
   108  		{PointFromCoords(1, 1, eps2), p000, p045, p180, pz, math.Pi},
   109  	}
   110  	for _, test := range tests {
   111  		area := PointArea(test.a, test.b, test.c) +
   112  			PointArea(test.a, test.c, test.d) +
   113  			PointArea(test.a, test.d, test.e) +
   114  			PointArea(test.a, test.e, test.b)
   115  
   116  		if !float64Eq(area, test.want) {
   117  			t.Errorf("Adding up 4 quarter hemispheres with PointArea(), got %v want %v", area, test.want)
   118  		}
   119  	}
   120  
   121  	// Compute the area of a hemisphere using four triangles with one near-180
   122  	// degree edge and one near-degenerate edge.
   123  	for i := 0; i < 100; i++ {
   124  		lng := s1.Angle(2 * math.Pi * randomFloat64())
   125  		p2Lng := lng + s1.Angle(randomFloat64())
   126  		p0 := PointFromLatLng(LatLng{1e-20, lng}.Normalized())
   127  		p1 := PointFromLatLng(LatLng{0, lng}.Normalized())
   128  		p2 := PointFromLatLng(LatLng{0, p2Lng}.Normalized())
   129  		p3 := PointFromLatLng(LatLng{0, lng + math.Pi}.Normalized())
   130  		p4 := PointFromLatLng(LatLng{0, lng + 5.0}.Normalized())
   131  		area := PointArea(p0, p1, p2) + PointArea(p0, p2, p3) + PointArea(p0, p3, p4) + PointArea(p0, p4, p1)
   132  		if !float64Near(area, 2*math.Pi, 2e-15) {
   133  			t.Errorf("hemisphere area of %v, %v, %v, %v, %v = %v, want %v", p1, p1, p2, p3, p4, area, 2*math.Pi)
   134  		}
   135  	}
   136  }
   137  
   138  func TestPointMeasuresAngleMethods(t *testing.T) {
   139  
   140  	tests := []struct {
   141  		a, b, c       Point
   142  		wantAngle     s1.Angle
   143  		wantTurnAngle s1.Angle
   144  	}{
   145  		{p000, pz, p045, math.Pi / 4, -3 * math.Pi / 4},
   146  		{p045, pz, p180, 3 * math.Pi / 4, -math.Pi / 4},
   147  		{p000, pz, p180, math.Pi, 0},
   148  		{pz, p000, p045, math.Pi / 2, math.Pi / 2},
   149  		{pz, p000, pz, 0, -math.Pi},
   150  	}
   151  
   152  	for _, test := range tests {
   153  		if got := Angle(test.a, test.b, test.c); math.Abs(float64(got-test.wantAngle)) > epsilon {
   154  			t.Errorf("Angle(%v, %v, %v) = %v, want %v", test.a, test.b, test.c, got, test.wantAngle)
   155  		}
   156  		if got := TurnAngle(test.a, test.b, test.c); math.Abs(float64(got-test.wantTurnAngle)) > epsilon {
   157  			t.Errorf("TurnAngle(%v, %v, %v) = %v, want %v", test.a, test.b, test.c, got, test.wantTurnAngle)
   158  		}
   159  	}
   160  }
   161  
   162  // Previously these three points shows catastrophic error in their cross product
   163  // which prevented Area() from falling back to the Girard method properly. They
   164  // returned an area on the order of 1e-14 and the real area is ~1e-21, 7 orders
   165  // of magnitude relative error. Check that they return zero now.
   166  func TestPointMeasuresPointAreaRegression(t *testing.T) {
   167  	a := Point{r3.Vector{-1.705424004316021258e-01, -8.242696197922716461e-01,
   168  		5.399026611737816062e-01}}
   169  	b := Point{r3.Vector{-1.706078905422188652e-01, -8.246067119418969416e-01,
   170  		5.393669607095969987e-01}}
   171  	c := Point{r3.Vector{-1.705800600596222294e-01, -8.244634596153025408e-01,
   172  		5.395947061167500891e-01}}
   173  	if area := PointArea(a, b, c); area != 0 {
   174  		t.Errorf("PointArea(%v, %v, %v) should have been 0, got %v", a, b, c, area)
   175  	}
   176  }
   177  
   178  func BenchmarkPointArea(b *testing.B) {
   179  	for i := 0; i < b.N; i++ {
   180  		PointArea(p000, p090, pz)
   181  	}
   182  }
   183  
   184  func BenchmarkPointAreaGirardCase(b *testing.B) {
   185  	for i := 0; i < b.N; i++ {
   186  		PointArea(g1, g2, g3)
   187  	}
   188  }
   189
View as plain text