point_test.go

Documentation: github.com/golang/geo/s2

     1  // Copyright 2014 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  func TestOriginPoint(t *testing.T) {
    26  	if math.Abs(OriginPoint().Norm()-1) > 1e-15 {
    27  		t.Errorf("Origin point norm = %v, want 1", OriginPoint().Norm())
    28  	}
    29  
    30  	// The point chosen below is about 66km from the north pole towards the East
    31  	// Siberian Sea. The purpose of the stToUV(2/3) calculation is to keep the
    32  	// origin as far away as possible from the longitudinal edges of large
    33  	// Cells. (The line of longitude through the chosen point is always 1/3
    34  	// or 2/3 of the way across any Cell with longitudinal edges that it
    35  	// passes through.)
    36  	p := Point{r3.Vector{-0.01, 0.01 * stToUV(2.0/3), 1}}
    37  	if !p.ApproxEqual(OriginPoint()) {
    38  		t.Errorf("Origin point should fall in the Siberian Sea, but does not.")
    39  	}
    40  
    41  	// Check that the origin is not too close to either pole.
    42  	if dist := math.Acos(OriginPoint().Z) * earthRadiusKm; dist <= 50 {
    43  		t.Errorf("Origin point is to close to the North Pole. Got %v, want >= 50km", dist)
    44  	}
    45  }
    46  
    47  func TestPointCross(t *testing.T) {
    48  	tests := []struct {
    49  		p1x, p1y, p1z, p2x, p2y, p2z, norm float64
    50  	}{
    51  		{1, 0, 0, 1, 0, 0, 1},
    52  		{1, 0, 0, 0, 1, 0, 2},
    53  		{0, 1, 0, 1, 0, 0, 2},
    54  		{1, 2, 3, -4, 5, -6, 2 * math.Sqrt(934)},
    55  	}
    56  	for _, test := range tests {
    57  		p1 := Point{r3.Vector{test.p1x, test.p1y, test.p1z}}
    58  		p2 := Point{r3.Vector{test.p2x, test.p2y, test.p2z}}
    59  		result := p1.PointCross(p2)
    60  		if !float64Eq(result.Norm(), test.norm) {
    61  			t.Errorf("|%v ⨯ %v| = %v, want %v", p1, p2, result.Norm(), test.norm)
    62  		}
    63  		if x := result.Dot(p1.Vector); !float64Eq(x, 0) {
    64  			t.Errorf("|(%v ⨯ %v) · %v| = %v, want 0", p1, p2, p1, x)
    65  		}
    66  		if x := result.Dot(p2.Vector); !float64Eq(x, 0) {
    67  			t.Errorf("|(%v ⨯ %v) · %v| = %v, want 0", p1, p2, p2, x)
    68  		}
    69  	}
    70  }
    71  
    72  func TestPointDistance(t *testing.T) {
    73  	tests := []struct {
    74  		x1, y1, z1 float64
    75  		x2, y2, z2 float64
    76  		want       float64 // radians
    77  	}{
    78  		{1, 0, 0, 1, 0, 0, 0},
    79  		{1, 0, 0, 0, 1, 0, math.Pi / 2},
    80  		{1, 0, 0, 0, 1, 1, math.Pi / 2},
    81  		{1, 0, 0, -1, 0, 0, math.Pi},
    82  		{1, 2, 3, 2, 3, -1, 1.2055891055045298},
    83  	}
    84  	for _, test := range tests {
    85  		p1 := Point{r3.Vector{test.x1, test.y1, test.z1}}
    86  		p2 := Point{r3.Vector{test.x2, test.y2, test.z2}}
    87  		if a := p1.Distance(p2).Radians(); !float64Eq(a, test.want) {
    88  			t.Errorf("%v.Distance(%v) = %v, want %v", p1, p2, a, test.want)
    89  		}
    90  		if a := p2.Distance(p1).Radians(); !float64Eq(a, test.want) {
    91  			t.Errorf("%v.Distance(%v) = %v, want %v", p2, p1, a, test.want)
    92  		}
    93  	}
    94  }
    95  
    96  func TestChordAngleBetweenPoints(t *testing.T) {
    97  	for iter := 0; iter < 10; iter++ {
    98  		m := randomFrame()
    99  		x := m.col(0)
   100  		y := m.col(1)
   101  		z := m.col(2)
   102  
   103  		if got := ChordAngleBetweenPoints(z, z).Angle(); got != 0 {
   104  			t.Errorf("ChordAngleBetweenPoints(%v, %v) = %v, want 0", z, z, got)
   105  		}
   106  		if got, want := ChordAngleBetweenPoints(Point{z.Mul(-1)}, z).Angle().Radians(), math.Pi; !float64Near(got, want, 1e-7) {
   107  			t.Errorf("ChordAngleBetweenPoints(%v, %v) = %v, want %v", z.Mul(-1), z, got, want)
   108  		}
   109  		if got, want := ChordAngleBetweenPoints(x, z).Angle().Radians(), math.Pi/2; !float64Eq(got, want) {
   110  			t.Errorf("ChordAngleBetweenPoints(%v, %v) = %v, want %v", x, z, got, want)
   111  		}
   112  		w := Point{y.Add(z.Vector).Normalize()}
   113  		if got, want := ChordAngleBetweenPoints(w, z).Angle().Radians(), math.Pi/4; !float64Eq(got, want) {
   114  			t.Errorf("ChordAngleBetweenPoints(%v, %v) = %v, want %v", w, z, got, want)
   115  		}
   116  	}
   117  }
   118  
   119  func TestPointApproxEqual(t *testing.T) {
   120  	tests := []struct {
   121  		x1, y1, z1 float64
   122  		x2, y2, z2 float64
   123  		want       bool
   124  	}{
   125  		{1, 0, 0, 1, 0, 0, true},
   126  		{1, 0, 0, 0, 1, 0, false},
   127  		{1, 0, 0, 0, 1, 1, false},
   128  		{1, 0, 0, -1, 0, 0, false},
   129  		{1, 2, 3, 2, 3, -1, false},
   130  		{1, 0, 0, 1 * (1 + epsilon), 0, 0, true},
   131  		{1, 0, 0, 1 * (1 - epsilon), 0, 0, true},
   132  		{1, 0, 0, 1 + epsilon, 0, 0, true},
   133  		{1, 0, 0, 1 - epsilon, 0, 0, true},
   134  		{1, 0, 0, 1, epsilon, 0, true},
   135  		{1, 0, 0, 1, epsilon, epsilon, false},
   136  		{1, epsilon, 0, 1, -epsilon, epsilon, false},
   137  	}
   138  	for _, test := range tests {
   139  		p1 := Point{r3.Vector{test.x1, test.y1, test.z1}}
   140  		p2 := Point{r3.Vector{test.x2, test.y2, test.z2}}
   141  		if got := p1.ApproxEqual(p2); got != test.want {
   142  			t.Errorf("%v.ApproxEqual(%v), got %v want %v", p1, p2, got, test.want)
   143  		}
   144  	}
   145  }
   146  
   147  func TestPointOrtho(t *testing.T) {
   148  	tests := []struct {
   149  		have Point
   150  		want Point
   151  	}{
   152  		// Vector's Ortho returns an axis-aligned ortho for an
   153  		// axis-aligned input. Check that this does not.
   154  		{
   155  			have: Point{r3.Vector{X: 1, Y: 0, Z: 0}},
   156  			want: Point{r3.Vector{X: 0, Y: -0.999985955295886075333556, Z: 0.005299925563068195837058}},
   157  		},
   158  		{
   159  			have: Point{r3.Vector{X: 0, Y: 1, Z: 0}},
   160  			want: Point{r3.Vector{X: 0.004569952278750987959000, Y: 0.0, Z: -0.999989557713564125585037}},
   161  		},
   162  		{
   163  			have: Point{r3.Vector{X: 0, Y: 0, Z: 1}},
   164  			want: Point{r3.Vector{X: -0.999928007775066962636856, Y: 0.011999136093300803371231, Z: 0}},
   165  		},
   166  
   167  		// Test a couple other values
   168  		{
   169  			have: Point{r3.Vector{X: 1, Y: 1, Z: 1}},
   170  			want: Point{r3.Vector{X: -0.709740689278763769998193, Y: 0.005297583276916723732386, Z: 0.704443106001847008101890}},
   171  		},
   172  		{
   173  			have: Point{r3.Vector{X: 3, Y: -2, Z: 0.4}},
   174  			want: Point{r3.Vector{X: -0.555687999915428054720223, Y: -0.831317152491703792449584, Z: 0.011074236907191168863274}},
   175  		},
   176  		{
   177  			have: Point{r3.Vector{X: 0.012, Y: 0.0053, Z: 0.00457}},
   178  			want: Point{r3.Vector{X: 0.404015523469256565558538, Y: -0.914752128609637393807930, Z: 0}},
   179  		},
   180  	}
   181  
   182  	for _, test := range tests {
   183  		got := Ortho(test.have)
   184  
   185  		if got != test.want {
   186  			t.Errorf("Ortho(%v) = %v, want %v", test.have, got, test.want)
   187  		}
   188  
   189  		// Test that the dot product with the orthogonal result is zero.
   190  		if !float64Eq(test.have.Dot(got.Vector), 0) {
   191  			t.Errorf("%v = not orthogonal to %v.Ortho()", test.have, got)
   192  		}
   193  
   194  		if !got.IsUnit() {
   195  			t.Errorf("%v should be unit length, but is not", got)
   196  		}
   197  	}
   198  }
   199  
   200  func TestPointRegularPoints(t *testing.T) {
   201  	// Conversion to/from degrees has a little more variability than the default epsilon.
   202  	const epsilon = 1e-13
   203  	center := PointFromLatLng(LatLngFromDegrees(80, 135))
   204  	radius := s1.Degree * 20
   205  	pts := regularPoints(center, radius, 4)
   206  
   207  	if len(pts) != 4 {
   208  		t.Errorf("regularPoints with 4 vertices should have 4 vertices, got %d", len(pts))
   209  	}
   210  
   211  	lls := []LatLng{
   212  		LatLngFromPoint(pts[0]),
   213  		LatLngFromPoint(pts[1]),
   214  		LatLngFromPoint(pts[2]),
   215  		LatLngFromPoint(pts[3]),
   216  	}
   217  	cll := LatLngFromPoint(center)
   218  
   219  	// Make sure that the radius is correct.
   220  	wantDist := 20.0
   221  	for i, ll := range lls {
   222  		if got := cll.Distance(ll).Degrees(); !float64Near(got, wantDist, epsilon) {
   223  			t.Errorf("Vertex %d distance from center = %v, want %v", i, got, wantDist)
   224  		}
   225  	}
   226  
   227  	// Make sure the angle between each point is correct.
   228  	wantAngle := math.Pi / 2
   229  	for i := 0; i < len(pts); i++ {
   230  		// Mod the index by 4 to wrap the values at each end.
   231  		v0, v1, v2 := pts[(4+i+1)%4], pts[(4+i)%4], pts[(4+i-1)%4]
   232  		if got := float64(v0.Sub(v1.Vector).Angle(v2.Sub(v1.Vector))); !float64Eq(got, wantAngle) {
   233  			t.Errorf("(%v-%v).Angle(%v-%v) = %v, want %v", v0, v1, v1, v2, got, wantAngle)
   234  		}
   235  	}
   236  
   237  	// Make sure that all edges of the polygon have the same length.
   238  	wantLength := 27.990890717782829
   239  	for i := 0; i < len(lls); i++ {
   240  		ll1, ll2 := lls[i], lls[(i+1)%4]
   241  		if got := ll1.Distance(ll2).Degrees(); !float64Near(got, wantLength, epsilon) {
   242  			t.Errorf("%v.Distance(%v) = %v, want %v", ll1, ll2, got, wantLength)
   243  		}
   244  	}
   245  
   246  	// Spot check an actual coordinate now that we know the points are spaced
   247  	// evenly apart at the same angles and radii.
   248  	if got, want := lls[0].Lat.Degrees(), 62.162880741097204; !float64Near(got, want, epsilon) {
   249  		t.Errorf("%v.Lat = %v, want %v", lls[0], got, want)
   250  	}
   251  	if got, want := lls[0].Lng.Degrees(), 103.11051028343407; !float64Near(got, want, epsilon) {
   252  		t.Errorf("%v.Lng = %v, want %v", lls[0], got, want)
   253  	}
   254  }
   255  
   256  func TestPointRegion(t *testing.T) {
   257  	p := Point{r3.Vector{1, 0, 0}}
   258  	r := Point{r3.Vector{1, 0, 0}}
   259  	if !r.Contains(p) {
   260  		t.Errorf("%v.Contains(%v) = false, want true", r, p)
   261  	}
   262  	if !r.ContainsPoint(p) {
   263  		t.Errorf("%v.ContainsPoint(%v) = false, want true", r, p)
   264  	}
   265  	if !r.Contains(r) {
   266  		t.Errorf("%v.Contains(%v) = false, want true", r, r)
   267  	}
   268  	if !r.ContainsPoint(r) {
   269  		t.Errorf("%v.ContainsPoint(%v) = false, want true", r, r)
   270  	}
   271  	if s := (Point{r3.Vector{1, 0, 1}}); r.Contains(s) {
   272  		t.Errorf("%v.Contains(%v) = true, want false", r, s)
   273  	}
   274  	if got, want := r.CapBound(), CapFromPoint(p); !got.ApproxEqual(want) {
   275  		t.Errorf("%v.CapBound() = %v, want %v", r, got, want)
   276  	}
   277  	if got, want := r.RectBound(), RectFromLatLng(LatLngFromPoint(p)); !rectsApproxEqual(got, want, epsilon, epsilon) {
   278  		t.Errorf("%v.RectBound() = %v, want %v", r, got, want)
   279  	}
   280  
   281  	// The leaf cell containing a point is still much larger than the point.
   282  	cell := CellFromPoint(p)
   283  	if r.ContainsCell(cell) {
   284  		t.Errorf("%v.ContainsCell(%v) = true, want false", r, cell)
   285  	}
   286  	if !r.IntersectsCell(cell) {
   287  		t.Errorf("%v.IntersectsCell(%v) = false, want true", r, cell)
   288  	}
   289  
   290  }
   291  
   292  func TestPointRotate(t *testing.T) {
   293  	for iter := 0; iter < 1000; iter++ {
   294  		axis := randomPoint()
   295  		target := randomPoint()
   296  		// Choose a distance whose logarithm is uniformly distributed.
   297  		distance := s1.Angle(math.Pi * math.Pow(1e-15, randomFloat64()))
   298  		// Sometimes choose points near the far side of the axis.
   299  		if oneIn(5) {
   300  			distance = math.Pi - distance
   301  		}
   302  		p := InterpolateAtDistance(distance, axis, target)
   303  		// Choose the rotation angle.
   304  		angle := s1.Angle(2 * math.Pi * math.Pow(1e-15, randomFloat64()))
   305  		if oneIn(3) {
   306  			angle = -angle
   307  		}
   308  		if oneIn(10) {
   309  			angle = 0
   310  		}
   311  
   312  		got := Rotate(p, axis, angle)
   313  
   314  		if !got.IsUnit() {
   315  			t.Errorf("%v should be unit length", got)
   316  		}
   317  
   318  		// got and p should be the same distance from axis.
   319  		const maxPositionError = 1e-15
   320  		if (got.Distance(axis) - p.Distance(axis)).Abs().Radians() > maxPositionError {
   321  			t.Errorf("rotated point %v should be same distance as %v, got %v, want %v", got, p, got.Distance(axis), p.Distance(axis))
   322  		}
   323  
   324  		// Check that the rotation angle is correct. We allow a fixed error in the
   325  		// *position* of the result, so we need to convert this into a rotation
   326  		// angle. The allowable error can be very large as "p" approaches "axis".
   327  		axisDistance := p.Cross(axis.Vector).Norm()
   328  		maxRotationError := 0.0
   329  		if axisDistance < maxPositionError {
   330  			maxRotationError = 2 * math.Pi
   331  		} else {
   332  			maxRotationError = math.Asin(maxPositionError / axisDistance)
   333  		}
   334  		actualRotation := TurnAngle(p, axis, got) + math.Pi
   335  		rotationError := math.Remainder((angle - actualRotation).Radians(), 2*math.Pi)
   336  		if rotationError > maxRotationError {
   337  			t.Errorf("rotational angle of %v = %v, want %v", got, actualRotation, angle)
   338  		}
   339  	}
   340  }
   341
View as plain text