Source file src/github.com/golang/geo/s2/cap_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  const (
    26  	tinyRad = 1e-10
    27  )
    28  
    29  var (
    30  	emptyCap   = EmptyCap()
    31  	fullCap    = FullCap()
    32  	defaultCap = EmptyCap()
    33  
    34  	zeroHeight  = 0.0
    35  	fullHeight  = 2.0
    36  	emptyHeight = -1.0
    37  
    38  	xAxisPt = Point{r3.Vector{1, 0, 0}}
    39  	yAxisPt = Point{r3.Vector{0, 1, 0}}
    40  
    41  	xAxis = CapFromPoint(xAxisPt)
    42  	yAxis = CapFromPoint(yAxisPt)
    43  	xComp = xAxis.Complement()
    44  
    45  	hemi = CapFromCenterHeight(PointFromCoords(1, 0, 1), 1)
    46  	tiny = CapFromCenterAngle(PointFromCoords(1, 2, 3), s1.Angle(tinyRad))
    47  
    48  	// A concave cap. Note that the error bounds for point containment tests
    49  	// increase with the cap angle, so we need to use a larger error bound
    50  	// here.
    51  	concaveCenter = PointFromLatLng(LatLngFromDegrees(80, 10))
    52  	concaveRadius = s1.ChordAngleFromAngle(150 * s1.Degree)
    53  	maxCapError   = concaveRadius.MaxPointError() + concaveRadius.MaxAngleError() + 3*dblEpsilon
    54  	concave       = CapFromCenterChordAngle(concaveCenter, concaveRadius)
    55  	concaveMin    = CapFromCenterChordAngle(concaveCenter, concaveRadius.Expanded(-maxCapError))
    56  	concaveMax    = CapFromCenterChordAngle(concaveCenter, concaveRadius.Expanded(maxCapError))
    57  )
    58  
    59  func TestCapBasicEmptyFullValid(t *testing.T) {
    60  	tests := []struct {
    61  		got                Cap
    62  		empty, full, valid bool
    63  	}{
    64  		{Cap{}, false, false, false},
    65  
    66  		{emptyCap, true, false, true},
    67  		{emptyCap.Complement(), false, true, true},
    68  		{fullCap, false, true, true},
    69  		{fullCap.Complement(), true, false, true},
    70  		{defaultCap, true, false, true},
    71  
    72  		{xComp, false, true, true},
    73  		{xComp.Complement(), true, false, true},
    74  
    75  		{tiny, false, false, true},
    76  		{concave, false, false, true},
    77  		{hemi, false, false, true},
    78  		{tiny, false, false, true},
    79  	}
    80  	for _, test := range tests {
    81  		if e := test.got.IsEmpty(); e != test.empty {
    82  			t.Errorf("%v.IsEmpty() = %t; want %t", test.got, e, test.empty)
    83  		}
    84  		if f := test.got.IsFull(); f != test.full {
    85  			t.Errorf("%v.IsFull() = %t; want %t", test.got, f, test.full)
    86  		}
    87  		if v := test.got.IsValid(); v != test.valid {
    88  			t.Errorf("%v.IsValid() = %t; want %t", test.got, v, test.valid)
    89  		}
    90  	}
    91  }
    92  
    93  func TestCapCenterHeightRadius(t *testing.T) {
    94  	if xAxis == xAxis.Complement().Complement() {
    95  		t.Errorf("the complement of the complement is not the original. %v == %v",
    96  			xAxis, xAxis.Complement().Complement())
    97  	}
    98  
    99  	if fullCap.Height() != fullHeight {
   100  		t.Error("full Caps should be full height")
   101  	}
   102  	if fullCap.Radius().Degrees() != 180.0 {
   103  		t.Error("radius of x-axis cap should be 180 degrees")
   104  	}
   105  
   106  	if emptyCap.center != defaultCap.center {
   107  		t.Error("empty Caps should be have the same center as the default")
   108  	}
   109  	if emptyCap.Height() != defaultCap.Height() {
   110  		t.Error("empty Caps should be have the same height as the default")
   111  	}
   112  
   113  	if yAxis.Height() != zeroHeight {
   114  		t.Error("y-axis cap should not be empty height")
   115  	}
   116  
   117  	if xAxis.Height() != zeroHeight {
   118  		t.Error("x-axis cap should not be empty height")
   119  	}
   120  	if xAxis.Radius().Radians() != zeroHeight {
   121  		t.Errorf("radius of x-axis cap got %f want %f", xAxis.Radius().Radians(), emptyHeight)
   122  	}
   123  
   124  	hc := Point{hemi.center.Mul(-1.0)}
   125  	if hc != hemi.Complement().center {
   126  		t.Error("hemi center and its complement should have the same center")
   127  	}
   128  	if hemi.Height() != 1.0 {
   129  		t.Error("hemi cap should be 1.0 in height")
   130  	}
   131  }
   132  
   133  func TestCapContains(t *testing.T) {
   134  	tests := []struct {
   135  		c1, c2 Cap
   136  		want   bool
   137  	}{
   138  		{emptyCap, emptyCap, true},
   139  		{fullCap, emptyCap, true},
   140  		{fullCap, fullCap, true},
   141  		{emptyCap, xAxis, false},
   142  		{fullCap, xAxis, true},
   143  		{xAxis, fullCap, false},
   144  		{xAxis, xAxis, true},
   145  		{xAxis, emptyCap, true},
   146  		{hemi, tiny, true},
   147  		{hemi, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/4-epsilon)), true},
   148  		{hemi, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/4+epsilon)), false},
   149  		{concave, hemi, true},
   150  		{concave, CapFromCenterHeight(Point{concave.center.Mul(-1.0)}, 0.1), false},
   151  	}
   152  	for _, test := range tests {
   153  		if got := test.c1.Contains(test.c2); got != test.want {
   154  			t.Errorf("%v.Contains(%v) = %t; want %t", test.c1, test.c2, got, test.want)
   155  		}
   156  	}
   157  }
   158  
   159  func TestCapContainsPoint(t *testing.T) {
   160  	tangent := tiny.center.Cross(r3.Vector{3, 2, 1}).Normalize()
   161  	tests := []struct {
   162  		c    Cap
   163  		p    Point
   164  		want bool
   165  	}{
   166  		{xAxis, xAxisPt, true},
   167  		{xAxis, Point{r3.Vector{1, 1e-20, 0}}, false},
   168  		{yAxis, xAxis.center, false},
   169  		{xComp, xAxis.center, true},
   170  		{xComp.Complement(), xAxis.center, false},
   171  		{tiny, Point{tiny.center.Add(tangent.Mul(tinyRad * 0.99))}, true},
   172  		{tiny, Point{tiny.center.Add(tangent.Mul(tinyRad * 1.01))}, false},
   173  		{hemi, PointFromCoords(1, 0, -(1 - epsilon)), true},
   174  		{hemi, xAxisPt, true},
   175  		{hemi.Complement(), xAxisPt, false},
   176  		{concaveMax, PointFromLatLng(LatLngFromDegrees(-70*(1-epsilon), 10)), true},
   177  		{concaveMin, PointFromLatLng(LatLngFromDegrees(-70*(1+epsilon), 10)), false},
   178  		{concaveMax, PointFromLatLng(LatLngFromDegrees(-50*(1-epsilon), -170)), true},
   179  		{concaveMin, PointFromLatLng(LatLngFromDegrees(-50*(1+epsilon), -170)), false},
   180  	}
   181  	for _, test := range tests {
   182  		if got := test.c.ContainsPoint(test.p); got != test.want {
   183  			t.Errorf("%v.ContainsPoint(%v) = %t, want %t", test.c, test.p, got, test.want)
   184  		}
   185  	}
   186  }
   187  
   188  func TestCapInteriorIntersects(t *testing.T) {
   189  	tests := []struct {
   190  		c1, c2 Cap
   191  		want   bool
   192  	}{
   193  		{emptyCap, emptyCap, false},
   194  		{emptyCap, xAxis, false},
   195  		{fullCap, emptyCap, false},
   196  		{fullCap, fullCap, true},
   197  		{fullCap, xAxis, true},
   198  		{xAxis, fullCap, false},
   199  		{xAxis, xAxis, false},
   200  		{xAxis, emptyCap, false},
   201  		{concave, hemi.Complement(), true},
   202  	}
   203  	for _, test := range tests {
   204  		if got := test.c1.InteriorIntersects(test.c2); got != test.want {
   205  			t.Errorf("%v.InteriorIntersects(%v); got %t want %t", test.c1, test.c2, got, test.want)
   206  		}
   207  	}
   208  }
   209  
   210  func TestCapInteriorContains(t *testing.T) {
   211  	if hemi.InteriorContainsPoint(Point{r3.Vector{1, 0, -(1 + epsilon)}}) {
   212  		t.Errorf("hemi (%v) should not contain point just past half way(%v)", hemi,
   213  			Point{r3.Vector{1, 0, -(1 + epsilon)}})
   214  	}
   215  }
   216  
   217  func TestCapCellUnionBoundLevel1Radius(t *testing.T) {
   218  	// Check that a cap whose radius is approximately the width of a level 1
   219  	// Cell can be covered by only 3 faces.
   220  	c := CapFromCenterAngle(PointFromCoords(1, 1, 1), s1.Angle(MinWidthMetric.Value(1)))
   221  	covering := c.CellUnionBound()
   222  	if len(covering) != 3 {
   223  		t.Errorf("a cap with radius of a level 1 cell should be covered by 3 faces, got %d", len(covering))
   224  	}
   225  }
   226  
   227  func TestCapExpanded(t *testing.T) {
   228  	cap50 := CapFromCenterAngle(xAxisPt, 50.0*s1.Degree)
   229  	cap51 := CapFromCenterAngle(xAxisPt, 51.0*s1.Degree)
   230  
   231  	if !emptyCap.Expanded(s1.Angle(fullHeight)).IsEmpty() {
   232  		t.Error("Expanding empty cap should return an empty cap")
   233  	}
   234  	if !fullCap.Expanded(s1.Angle(fullHeight)).IsFull() {
   235  		t.Error("Expanding a full cap should return an full cap")
   236  	}
   237  
   238  	if !cap50.Expanded(0).ApproxEqual(cap50) {
   239  		t.Error("Expanding a cap by 0° should be equal to the original")
   240  	}
   241  	if !cap50.Expanded(1 * s1.Degree).ApproxEqual(cap51) {
   242  		t.Error("Expanding 50° by 1° should equal the 51° cap")
   243  	}
   244  
   245  	if cap50.Expanded(129.99 * s1.Degree).IsFull() {
   246  		t.Error("Expanding 50° by 129.99° should not give a full cap")
   247  	}
   248  	if !cap50.Expanded(130.01 * s1.Degree).IsFull() {
   249  		t.Error("Expanding 50° by 130.01° should give a full cap")
   250  	}
   251  }
   252  
   253  func TestCapRadiusToHeight(t *testing.T) {
   254  	tests := []struct {
   255  		got  s1.Angle
   256  		want float64
   257  	}{
   258  		// Above/below boundary checks.
   259  		{s1.Angle(-0.5), emptyHeight},
   260  		{s1.Angle(0), 0},
   261  		{s1.Angle(math.Pi), fullHeight},
   262  		{s1.Angle(2 * math.Pi), fullHeight},
   263  		// Degree tests.
   264  		{-7.0 * s1.Degree, emptyHeight},
   265  		{0.0 * s1.Degree, 0},
   266  		{12.0 * s1.Degree, 0.0218523992661943},
   267  		{30.0 * s1.Degree, 0.1339745962155613},
   268  		{45.0 * s1.Degree, 0.2928932188134525},
   269  		{90.0 * s1.Degree, 1.0},
   270  		{179.99 * s1.Degree, 1.9999999847691292},
   271  		{180.0 * s1.Degree, fullHeight},
   272  		{270.0 * s1.Degree, fullHeight},
   273  		// Radians tests.
   274  		{-1.0 * s1.Radian, emptyHeight},
   275  		{0.0 * s1.Radian, 0},
   276  		{1.0 * s1.Radian, 0.45969769413186},
   277  		{math.Pi / 2.0 * s1.Radian, 1.0},
   278  		{2.0 * s1.Radian, 1.4161468365471424},
   279  		{3.0 * s1.Radian, 1.9899924966004454},
   280  		{math.Pi * s1.Radian, fullHeight},
   281  		{4.0 * s1.Radian, fullHeight},
   282  	}
   283  	for _, test := range tests {
   284  		// float64Eq comes from s2latlng_test.go
   285  		if got := radiusToHeight(test.got); !float64Eq(got, test.want) {
   286  			t.Errorf("radiusToHeight(%v) = %v; want %v", test.got, got, test.want)
   287  		}
   288  	}
   289  }
   290  
   291  func TestCapRectBounds(t *testing.T) {
   292  	const epsilon = 1e-13
   293  	var tests = []struct {
   294  		desc     string
   295  		have     Cap
   296  		latLoDeg float64
   297  		latHiDeg float64
   298  		lngLoDeg float64
   299  		lngHiDeg float64
   300  		isFull   bool
   301  	}{
   302  		{
   303  			"Cap that includes South Pole.",
   304  			CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(-45, 57)), s1.Degree*50),
   305  			-90, 5, -180, 180, true,
   306  		},
   307  		{
   308  			"Cap that is tangent to the North Pole.",
   309  			CapFromCenterAngle(PointFromCoords(1, 0, 1), s1.Radian*(math.Pi/4.0+1e-16)),
   310  			0, 90, -180, 180, true,
   311  		},
   312  		{
   313  			"Cap that at 45 degree center that goes from equator to the pole.",
   314  			CapFromCenterAngle(PointFromCoords(1, 0, 1), s1.Degree*(45+5e-15)),
   315  			0, 90, -180, 180, true,
   316  		},
   317  		{
   318  			"The eastern hemisphere.",
   319  			CapFromCenterAngle(Point{r3.Vector{0, 1, 0}}, s1.Radian*(math.Pi/2+2e-16)),
   320  			-90, 90, -180, 180, true,
   321  		},
   322  		{
   323  			"A cap centered on the equator.",
   324  			CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(0, 50)), s1.Degree*20),
   325  			-20, 20, 30, 70, false,
   326  		},
   327  		{
   328  			"A cap centered on the North Pole.",
   329  			CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(90, 123)), s1.Degree*10),
   330  			80, 90, -180, 180, true,
   331  		},
   332  	}
   333  
   334  	for _, test := range tests {
   335  		r := test.have.RectBound()
   336  		if !float64Near(s1.Angle(r.Lat.Lo).Degrees(), test.latLoDeg, epsilon) {
   337  			t.Errorf("%s: %v.RectBound(), Lat.Lo not close enough, got %0.20f, want %0.20f",
   338  				test.desc, test.have, s1.Angle(r.Lat.Lo).Degrees(), test.latLoDeg)
   339  		}
   340  		if !float64Near(s1.Angle(r.Lat.Hi).Degrees(), test.latHiDeg, epsilon) {
   341  			t.Errorf("%s: %v.RectBound(), Lat.Hi not close enough, got %0.20f, want %0.20f",
   342  				test.desc, test.have, s1.Angle(r.Lat.Hi).Degrees(), test.latHiDeg)
   343  		}
   344  		if !float64Near(s1.Angle(r.Lng.Lo).Degrees(), test.lngLoDeg, epsilon) {
   345  			t.Errorf("%s: %v.RectBound(), Lng.Lo not close enough, got %0.20f, want %0.20f",
   346  				test.desc, test.have, s1.Angle(r.Lng.Lo).Degrees(), test.lngLoDeg)
   347  		}
   348  		if !float64Near(s1.Angle(r.Lng.Hi).Degrees(), test.lngHiDeg, epsilon) {
   349  			t.Errorf("%s: %v.RectBound(), Lng.Hi not close enough, got %0.20f, want %0.20f",
   350  				test.desc, test.have, s1.Angle(r.Lng.Hi).Degrees(), test.lngHiDeg)
   351  		}
   352  		if got := r.Lng.IsFull(); got != test.isFull {
   353  			t.Errorf("%s: RectBound(%v).isFull() = %t, want %t", test.desc, test.have, got, test.isFull)
   354  		}
   355  	}
   356  
   357  	// Empty and full caps.
   358  	if !EmptyCap().RectBound().IsEmpty() {
   359  		t.Errorf("RectBound() on EmptyCap should be empty.")
   360  	}
   361  
   362  	if !FullCap().RectBound().IsFull() {
   363  		t.Errorf("RectBound() on FullCap should be full.")
   364  	}
   365  }
   366  
   367  func TestCapAddPoint(t *testing.T) {
   368  	const epsilon = 1e-14
   369  	tests := []struct {
   370  		have Cap
   371  		p    Point
   372  		want Cap
   373  	}{
   374  		// Cap plus its center equals itself.
   375  		{xAxis, xAxisPt, xAxis},
   376  		{yAxis, yAxisPt, yAxis},
   377  
   378  		// Cap plus opposite point equals full.
   379  		{xAxis, Point{r3.Vector{-1, 0, 0}}, fullCap},
   380  		{yAxis, Point{r3.Vector{0, -1, 0}}, fullCap},
   381  
   382  		// Cap plus orthogonal axis equals half cap.
   383  		{xAxis, Point{r3.Vector{0, 0, 1}}, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},
   384  		{xAxis, Point{r3.Vector{0, 0, -1}}, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},
   385  
   386  		// The 45 degree angled hemisphere plus some points.
   387  		{
   388  			hemi,
   389  			PointFromCoords(0, 1, -1),
   390  			CapFromCenterAngle(Point{r3.Vector{1, 0, 1}},
   391  				s1.Angle(120.0)*s1.Degree),
   392  		},
   393  		{
   394  			hemi,
   395  			PointFromCoords(0, -1, -1),
   396  			CapFromCenterAngle(Point{r3.Vector{1, 0, 1}},
   397  				s1.Angle(120.0)*s1.Degree),
   398  		},
   399  		{
   400  			hemi,
   401  			PointFromCoords(-1, -1, -1),
   402  			CapFromCenterAngle(Point{r3.Vector{1, 0, 1}},
   403  				s1.Angle(math.Acos(-math.Sqrt(2.0/3.0)))),
   404  		},
   405  		{hemi, Point{r3.Vector{0, 1, 1}}, hemi},
   406  		{hemi, Point{r3.Vector{1, 0, 0}}, hemi},
   407  	}
   408  
   409  	for _, test := range tests {
   410  		got := test.have.AddPoint(test.p)
   411  		if !got.ApproxEqual(test.want) {
   412  			t.Errorf("%v.AddPoint(%v) = %v, want %v", test.have, test.p, got, test.want)
   413  		}
   414  
   415  		if !got.ContainsPoint(test.p) {
   416  			t.Errorf("%v.AddPoint(%v) did not contain added point", test.have, test.p)
   417  		}
   418  	}
   419  }
   420  
   421  func TestCapAddCap(t *testing.T) {
   422  	tests := []struct {
   423  		have  Cap
   424  		other Cap
   425  		want  Cap
   426  	}{
   427  		// Identity cases.
   428  		{emptyCap, emptyCap, emptyCap},
   429  		{fullCap, fullCap, fullCap},
   430  
   431  		// Anything plus empty equals itself.
   432  		{fullCap, emptyCap, fullCap},
   433  		{emptyCap, fullCap, fullCap},
   434  		{xAxis, emptyCap, xAxis},
   435  		{emptyCap, xAxis, xAxis},
   436  		{yAxis, emptyCap, yAxis},
   437  		{emptyCap, yAxis, yAxis},
   438  
   439  		// Two halves make a whole.
   440  		{xAxis, xComp, fullCap},
   441  
   442  		// Two zero-height orthogonal axis caps make a half-cap.
   443  		{xAxis, yAxis, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},
   444  	}
   445  
   446  	for _, test := range tests {
   447  		got := test.have.AddCap(test.other)
   448  		if !got.ApproxEqual(test.want) {
   449  			t.Errorf("%v.AddCap(%v) = %v, want %v", test.have, test.other, got, test.want)
   450  		}
   451  	}
   452  }
   453  
   454  func TestCapContainsCell(t *testing.T) {
   455  	faceRadius := math.Atan(math.Sqrt2)
   456  	for face := 0; face < 6; face++ {
   457  		// The cell consisting of the entire face.
   458  		rootCell := CellFromCellID(CellIDFromFace(face))
   459  
   460  		// A leaf cell at the midpoint of the v=1 edge.
   461  		edgeCell := CellFromPoint(Point{faceUVToXYZ(face, 0, 1-epsilon)})
   462  
   463  		// A leaf cell at the u=1, v=1 corner
   464  		cornerCell := CellFromPoint(Point{faceUVToXYZ(face, 1-epsilon, 1-epsilon)})
   465  
   466  		// Quick check for full and empty caps.
   467  		if !fullCap.ContainsCell(rootCell) {
   468  			t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", fullCap, rootCell, false, true)
   469  		}
   470  
   471  		// Check intersections with the bounding caps of the leaf cells that are adjacent to
   472  		// cornerCell along the Hilbert curve.  Because this corner is at (u=1,v=1), the curve
   473  		// stays locally within the same cube face.
   474  		first := cornerCell.id.Advance(-3)
   475  		last := cornerCell.id.Advance(4)
   476  		for id := first; id < last; id = id.Next() {
   477  			c := CellFromCellID(id).CapBound()
   478  			if got, want := c.ContainsCell(cornerCell), id == cornerCell.id; got != want {
   479  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", c, cornerCell, got, want)
   480  			}
   481  		}
   482  
   483  		for capFace := 0; capFace < 6; capFace++ {
   484  			// A cap that barely contains all of capFace.
   485  			center := unitNorm(capFace)
   486  			covering := CapFromCenterAngle(center, s1.Angle(faceRadius+epsilon))
   487  			if got, want := covering.ContainsCell(rootCell), capFace == face; got != want {
   488  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, rootCell, got, want)
   489  			}
   490  			if got, want := covering.ContainsCell(edgeCell), center.Vector.Dot(edgeCell.id.Point().Vector) > 0.1; got != want {
   491  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, edgeCell, got, want)
   492  			}
   493  			if got, want := covering.ContainsCell(edgeCell), covering.IntersectsCell(edgeCell); got != want {
   494  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, edgeCell, got, want)
   495  			}
   496  			if got, want := covering.ContainsCell(cornerCell), capFace == face; got != want {
   497  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, cornerCell, got, want)
   498  			}
   499  
   500  			// A cap that barely intersects the edges of capFace.
   501  			bulging := CapFromCenterAngle(center, s1.Angle(math.Pi/4+epsilon))
   502  			if bulging.ContainsCell(rootCell) {
   503  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, rootCell, true, false)
   504  			}
   505  			if got, want := bulging.ContainsCell(edgeCell), capFace == face; got != want {
   506  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, edgeCell, got, want)
   507  			}
   508  			if bulging.ContainsCell(cornerCell) {
   509  				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, cornerCell, true, false)
   510  			}
   511  		}
   512  	}
   513  }
   514  
   515  func TestCapIntersectsCell(t *testing.T) {
   516  	faceRadius := math.Atan(math.Sqrt2)
   517  	for face := 0; face < 6; face++ {
   518  		// The cell consisting of the entire face.
   519  		rootCell := CellFromCellID(CellIDFromFace(face))
   520  
   521  		// A leaf cell at the midpoint of the v=1 edge.
   522  		edgeCell := CellFromPoint(Point{faceUVToXYZ(face, 0, 1-epsilon)})
   523  
   524  		// A leaf cell at the u=1, v=1 corner
   525  		cornerCell := CellFromPoint(Point{faceUVToXYZ(face, 1-epsilon, 1-epsilon)})
   526  
   527  		// Quick check for full and empty caps.
   528  		if emptyCap.IntersectsCell(rootCell) {
   529  			t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", emptyCap, rootCell, true, false)
   530  		}
   531  
   532  		// Check intersections with the bounding caps of the leaf cells that are adjacent to
   533  		// cornerCell along the Hilbert curve.  Because this corner is at (u=1,v=1), the curve
   534  		// stays locally within the same cube face.
   535  		first := cornerCell.id.Advance(-3)
   536  		last := cornerCell.id.Advance(4)
   537  		for id := first; id < last; id = id.Next() {
   538  			c := CellFromCellID(id).CapBound()
   539  			if got, want := c.IntersectsCell(cornerCell), id.immediateParent().Contains(cornerCell.id); got != want {
   540  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", c, cornerCell, got, want)
   541  			}
   542  		}
   543  
   544  		antiFace := (face + 3) % 6
   545  		for capFace := 0; capFace < 6; capFace++ {
   546  			// A cap that barely contains all of capFace.
   547  			center := unitNorm(capFace)
   548  			covering := CapFromCenterAngle(center, s1.Angle(faceRadius+epsilon))
   549  			if got, want := covering.IntersectsCell(rootCell), capFace != antiFace; got != want {
   550  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", covering, rootCell, got, want)
   551  			}
   552  			if got, want := covering.IntersectsCell(edgeCell), covering.ContainsCell(edgeCell); got != want {
   553  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", covering, edgeCell, got, want)
   554  			}
   555  			if got, want := covering.IntersectsCell(cornerCell), center.Vector.Dot(cornerCell.id.Point().Vector) > 0; got != want {
   556  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", covering, cornerCell, got, want)
   557  			}
   558  
   559  			// A cap that barely intersects the edges of capFace.
   560  			bulging := CapFromCenterAngle(center, s1.Angle(math.Pi/4+epsilon))
   561  			if got, want := bulging.IntersectsCell(rootCell), capFace != antiFace; got != want {
   562  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", bulging, rootCell, got, want)
   563  			}
   564  			if got, want := bulging.IntersectsCell(edgeCell), center.Vector.Dot(edgeCell.id.Point().Vector) > 0.1; got != want {
   565  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", bulging, edgeCell, got, want)
   566  			}
   567  			if bulging.IntersectsCell(cornerCell) {
   568  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", bulging, cornerCell, true, false)
   569  			}
   570  
   571  			// A singleton cap.
   572  			singleton := CapFromCenterAngle(center, 0)
   573  			if got, want := singleton.IntersectsCell(rootCell), capFace == face; got != want {
   574  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", singleton, rootCell, got, want)
   575  			}
   576  			if singleton.IntersectsCell(edgeCell) {
   577  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", singleton, edgeCell, true, false)
   578  			}
   579  			if singleton.IntersectsCell(cornerCell) {
   580  				t.Errorf("Cap(%v).IntersectsCell(%v) = %t; want = %t", singleton, cornerCell, true, false)
   581  			}
   582  		}
   583  	}
   584  }
   585  
   586  func TestCapCentroid(t *testing.T) {
   587  	// Empty and full caps.
   588  	if got, want := EmptyCap().Centroid(), (Point{}); !got.ApproxEqual(want) {
   589  		t.Errorf("Centroid of EmptyCap should be zero point, got %v", want)
   590  	}
   591  	if got, want := FullCap().Centroid().Norm(), 1e-15; got > want {
   592  		t.Errorf("Centroid of FullCap should have a Norm of 0, got %v", want)
   593  	}
   594  
   595  	// Random caps.
   596  	for i := 0; i < 100; i++ {
   597  		center := randomPoint()
   598  		height := randomUniformFloat64(0.0, 2.0)
   599  		c := CapFromCenterHeight(center, height)
   600  		got := c.Centroid()
   601  		want := center.Mul((1.0 - height/2.0) * c.Area())
   602  		if delta := got.Sub(want).Norm(); delta > 1e-15 {
   603  			t.Errorf("%v.Sub(%v).Norm() = %v, want %v", got, want, delta, 1e-15)
   604  		}
   605  	}
   606  }
   607  
   608  func TestCapUnion(t *testing.T) {
   609  	// Two caps which have the same center but one has a larger radius.
   610  	a := CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(50.0, 10.0)), s1.Degree*0.2)
   611  	b := CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(50.0, 10.0)), s1.Degree*0.3)
   612  	if !b.Contains(a) {
   613  		t.Errorf("%v.Contains(%v) = false, want true", b, a)
   614  	}
   615  	if got := b.ApproxEqual(a.Union(b)); !got {
   616  		t.Errorf("%v.ApproxEqual(%v) = %v, want true", b, a.Union(b), got)
   617  	}
   618  
   619  	// Two caps where one is the full cap.
   620  	if got := a.Union(FullCap()); !got.IsFull() {
   621  		t.Errorf("%v.Union(%v).IsFull() = %v, want true", a, got, got.IsFull())
   622  	}
   623  
   624  	// Two caps where one is the empty cap.
   625  	if got := a.Union(EmptyCap()); !a.ApproxEqual(got) {
   626  		t.Errorf("%v.Union(EmptyCap) = %v, want %v", a, got, a)
   627  	}
   628  
   629  	// Two caps which have different centers, one entirely encompasses the other.
   630  	c := CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(51.0, 11.0)), s1.Degree*1.5)
   631  	if !c.Contains(a) {
   632  		t.Errorf("%v.Contains(%v) = false, want true", c, a)
   633  	}
   634  	if got := a.Union(c).center; !got.ApproxEqual(c.center) {
   635  		t.Errorf("%v.Union(%v).center = %v, want %v", a, c, got, c.center)
   636  	}
   637  	if got := a.Union(c); !float64Eq(float64(got.Radius()), float64(c.Radius())) {
   638  		t.Errorf("%v.Union(%v).Radius = %v, want %v", a, c, got.Radius(), c.Radius())
   639  	}
   640  
   641  	// Two entirely disjoint caps.
   642  	d := CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(51.0, 11.0)), s1.Degree*0.1)
   643  	if d.Contains(a) {
   644  		t.Errorf("%v.Contains(%v) = true, want false", d, a)
   645  	}
   646  	if d.Intersects(a) {
   647  		t.Errorf("%v.Intersects(%v) = true, want false", d, a)
   648  	}
   649  
   650  	// Check union and reverse direction are the same.
   651  	aUnionD := a.Union(d)
   652  	if !aUnionD.ApproxEqual(d.Union(a)) {
   653  		t.Errorf("%v.Union(%v).ApproxEqual(%v.Union(%v)) = false, want true", a, d, d, a)
   654  	}
   655  	if got, want := LatLngFromPoint(aUnionD.center).Lat.Degrees(), 50.4588; !float64Near(got, want, 0.001) {
   656  		t.Errorf("%v.Center.Lat = %v, want %v", aUnionD, got, want)
   657  	}
   658  	if got, want := LatLngFromPoint(aUnionD.center).Lng.Degrees(), 10.4525; !float64Near(got, want, 0.001) {
   659  		t.Errorf("%v.Center.Lng = %v, want %v", aUnionD, got, want)
   660  	}
   661  	if got, want := aUnionD.Radius().Degrees(), 0.7425; !float64Near(got, want, 0.001) {
   662  		t.Errorf("%v.Radius = %v, want %v", aUnionD, got, want)
   663  	}
   664  
   665  	// Two partially overlapping caps.
   666  	e := CapFromCenterAngle(PointFromLatLng(LatLngFromDegrees(50.3, 10.3)), s1.Degree*0.2)
   667  	aUnionE := a.Union(e)
   668  	if e.Contains(a) {
   669  		t.Errorf("%v.Contains(%v) = false, want true", e, a)
   670  	}
   671  	if !e.Intersects(a) {
   672  		t.Errorf("%v.Intersects(%v) = false, want true", e, a)
   673  	}
   674  	if !aUnionE.ApproxEqual(e.Union(a)) {
   675  		t.Errorf("%v.Union(%v).ApproxEqual(%v.Union(%v)) = false, want true", a, e, e, a)
   676  	}
   677  	if got, want := LatLngFromPoint(aUnionE.center).Lat.Degrees(), 50.1500; !float64Near(got, want, 0.001) {
   678  		t.Errorf("%v.Center.Lat = %v, want %v", aUnionE, got, want)
   679  	}
   680  	if got, want := LatLngFromPoint(aUnionE.center).Lng.Degrees(), 10.1495; !float64Near(got, want, 0.001) {
   681  		t.Errorf("%v.Center.Lng = %v, want %v", aUnionE, got, want)
   682  	}
   683  	if got, want := aUnionE.Radius().Degrees(), 0.3781; !float64Near(got, want, 0.001) {
   684  		t.Errorf("%v.Radius = %v, want %v", aUnionE, got, want)
   685  	}
   686  
   687  	p1 := Point{r3.Vector{0, 0, 1}}
   688  	p2 := Point{r3.Vector{0, 1, 0}}
   689  	// Two very large caps, whose radius sums to in excess of 180 degrees, and
   690  	// whose centers are not antipodal.
   691  	f := CapFromCenterAngle(p1, s1.Degree*150)
   692  	g := CapFromCenterAngle(p2, s1.Degree*150)
   693  	if !f.Union(g).IsFull() {
   694  		t.Errorf("%v.Union(%v).IsFull() = false, want true", f, g)
   695  	}
   696  
   697  	// Two non-overlapping hemisphere caps with antipodal centers.
   698  	hemi := CapFromCenterHeight(p1, 1)
   699  	if !hemi.Union(hemi.Complement()).IsFull() {
   700  		t.Errorf("%v.Union(%v).Complement().IsFull() = false, want true", hemi, hemi.Complement())
   701  	}
   702  }
   703  
   704  func TestCapEqual(t *testing.T) {
   705  	tests := []struct {
   706  		a, b Cap
   707  		want bool
   708  	}{
   709  		{EmptyCap(), EmptyCap(), true},
   710  		{EmptyCap(), FullCap(), false},
   711  		{FullCap(), FullCap(), true},
   712  		{
   713  			CapFromCenterAngle(PointFromCoords(0, 0, 1), s1.Degree*150),
   714  			CapFromCenterAngle(PointFromCoords(0, 0, 1), s1.Degree*151),
   715  			false,
   716  		},
   717  		{xAxis, xAxis, true},
   718  		{xAxis, yAxis, false},
   719  		{xComp, xAxis.Complement(), true},
   720  	}
   721  
   722  	for _, test := range tests {
   723  		if got := test.a.Equal(test.b); got != test.want {
   724  			t.Errorf("%v.Equal(%v) = %t, want %t", test.a, test.b, got, test.want)
   725  		}
   726  	}
   727  }
   728  

View as plain text