Tracked down error in slab decomposition, improved slab decomposition…

… (only in bench so far) to use a mapping of faces to edges intead of a brute force check at the end

.gitignore

@@ -24,3 +24,4 @@ _testmain.go
 *.test
 *.prof
 
+test*.off

dcel/dcel.go

@@ -180,8 +180,8 @@ func (dc *DCEL) FullEdges() ([]geom.FullEdge, [][2]*Face, error) {
 // whether its chain is the inner or outer
 // portion of a face.
 func (dc *DCEL) CorrectDirectionality(f *Face) {
-	// Inners need to be going CC
-	// Outers need to be going Clockwise
+	// Inners need to be going Clockwise
+	// Outers need to be going CC
 
 	clock, err := f.Outer.IsClockwise()
 	if err == nil && clock {
@@ -194,6 +194,21 @@ func (dc *DCEL) CorrectDirectionality(f *Face) {
 	}
 }
 
+// This is probably how we actually want to do flipping, always.
+func (dc *DCEL) CorrectDirectionalityAll() {
+
+	for _, f := range dc.Faces {
+		clock, err := f.Outer.IsClockwise()
+		if err == nil && clock {
+			f.Outer = f.Outer.Twin
+		}
+		clock, err = f.Inner.IsClockwise()
+		if err == nil && !clock {
+			f.Inner = f.Inner.Twin
+		}
+	}
+}
+
 // CorrectTwins modifies the ordering on twins
 // inside the DCEL such that dc.HalfEdges[i] is
 // the twin of dc.HalfEdges[i+1] for all even

dcel/edge.go

@@ -318,6 +318,9 @@ func (e *Edge) Eq(e2 geom.Spanning) bool {
 // which actually calls this instead of the other way
 // around because that involves less code duplciation.
 func (e *Edge) AllEdges() []*Edge {
+	if e == nil {
+		return []*Edge{}
+	}
 	edges := make([]*Edge, 1)
 	edges[0] = e
 	edge := e.Twin.Next
@@ -348,6 +351,12 @@ func (e *Edge) FindSharedPoint(e2 *Edge, d int) (float64, error) {
 	e2Low := e2.Low(d)
 	eHigh := e.High(d)
 	e2High := e2.High(d)
+	// The two edges span the same distance case
+	if geom.F64eq(eLow.Val(d), e2Low.Val(d)) &&
+		geom.F64eq(eHigh.Val(d), e2High.Val(d)) {
+		return ((eHigh.Val(d) - eLow.Val(d)) / 2) + eLow.Val(d), nil
+	}
+
 	// No point exists case:
 	// if e.High is lower than e2.Low
 	// or if e.Low is higher than e2.High

dcel/face.go

@@ -21,7 +21,7 @@ func NewFace() *Face {
 // Vertices wraps around a face and
 // finds all vertices that border it.
 func (f *Face) Vertices() []*Vertex {
-	// Outer is not populated by anything as of this writing.
+	// Inner is not populated by anything as of this writing.
 
 	pts := []*Vertex{}
 	e := f.Outer

dcel/off/load.go

@@ -291,5 +291,7 @@ func decode(dc *dcel.DCEL, edges []*dcel.Edge,
 		ei++
 	}
 
+	dc.CorrectDirectionalityAll()
+
 	return dc, nil
 }

dcel/off/off.go

@@ -1,7 +1,6 @@
-// Package off describes methods for interacting with OFF files and structures
+// package off describes methods for interacting with OFF files and structures
 // formatted as OFF files. The file loading code is modeled after Ryan Holmes'
 // C++ code, http://www.holmes3d.net/graphics/offfiles/
-
 package off
 
 import "github.com/200sc/go-compgeo/geom"

dcel/pointLoc/bench/slab/compEdge.go

@@ -1,6 +1,8 @@
 package slab
 
 import (
+	"fmt"
+
 	"github.com/200sc/go-compgeo/dcel"
 	"github.com/200sc/go-compgeo/geom"
 	"github.com/200sc/go-compgeo/search"
@@ -18,6 +20,18 @@ func (fs faces) Equals(e search.Equalable) bool {
 	return false
 }
 
+type face struct {
+	*dcel.Face
+}
+
+func (f face) Equals(e search.Equalable) bool {
+	switch f2 := e.(type) {
+	case face:
+		return f.Face == f2.Face
+	}
+	return false
+}
+
 type shellNode struct {
 	k compEdge
 	v search.Equalable
@@ -46,7 +60,10 @@ func (ce compEdge) Compare(i interface{}) search.CompareResult {
 			geom.F64eq(ce.Twin.X(), c.Twin.X()) && geom.F64eq(ce.Twin.Y(), c.Twin.Y()) {
 			return search.Equal
 		}
-		compX, _ := ce.FindSharedPoint(c.Edge, 0)
+		compX, err := ce.FindSharedPoint(c.Edge, 0)
+		if err != nil {
+			fmt.Println("Edges share no point on x axis", ce, c)
+		}
 		p1, _ := ce.PointAt(0, compX)
 		p2, _ := c.PointAt(0, compX)
 		if p1[1] < p2[1] {

dcel/pointLoc/bench/slab/slabDecomp.go

@@ -3,10 +3,11 @@
 package slab
 
 import (
+	"fmt"
+
 	compgeo "github.com/200sc/go-compgeo"
 	"github.com/200sc/go-compgeo/dcel"
 	"github.com/200sc/go-compgeo/dcel/pointLoc"
-	"github.com/200sc/go-compgeo/dcel/pointLoc/visualize"
 	"github.com/200sc/go-compgeo/geom"
 	"github.com/200sc/go-compgeo/search"
 	"github.com/200sc/go-compgeo/search/tree"
@@ -29,6 +30,37 @@ func Decompose(dc *dcel.DCEL, bstType tree.Type) (pointLoc.LocatesPoints, error)
 	t := tree.New(bstType).ToPersistent()
 	pts := dc.VerticesSorted(0)
 
+	// For each edge, we need to know which face lies beneath it in its two
+	// faces.
+	faceEdgeMap := make(map[*dcel.Edge]*dcel.Face)
+	for _, f := range dc.Faces {
+		// walk each face
+		e := f.Outer
+		if e != nil {
+			if e.Origin.X() < e.Twin.Origin.X() {
+				faceEdgeMap[e] = f
+			}
+			for e = e.Next; e != f.Outer; e = e.Next {
+				// This edge points right, in an outer face.
+				// Then this face lies beneath e.
+				if e.Origin.X() < e.Twin.Origin.X() {
+					faceEdgeMap[e] = f
+				}
+			}
+		}
+		e = f.Inner
+		if e != nil {
+			if e.Origin.X() > e.Twin.Origin.X() {
+				faceEdgeMap[e] = f
+			}
+			for e = e.Next; e != f.Inner; e = e.Next {
+				if e.Origin.X() > e.Twin.Origin.X() {
+					faceEdgeMap[e] = f
+				}
+			}
+		}
+	}
+
 	i := 0
 	for i < len(pts) {
 		p := pts[i]
@@ -63,7 +95,10 @@ func Decompose(dc *dcel.DCEL, bstType tree.Type) (pointLoc.LocatesPoints, error)
 		// Remove all edges from the PersistentBST connecting to the left
 		// of the points
 		for _, e := range le {
-			ct.Delete(shellNode{compEdge{e.Twin}, search.Nil{}})
+			err := ct.Delete(shellNode{compEdge{e.Twin}, search.Nil{}})
+			if err != nil {
+				fmt.Println(err, e.Twin)
+			}
 		}
 		// Add all edges to the PersistentBST connecting to the right
 		// of the point
@@ -74,7 +109,7 @@ func Decompose(dc *dcel.DCEL, bstType tree.Type) (pointLoc.LocatesPoints, error)
 			// locate to the edge above the query point. Returning an
 			// edge for a query represents that the query is below
 			// the edge,
-			ct.Insert(shellNode{compEdge{e}, faces{e.Face, e.Twin.Face}})
+			ct.Insert(shellNode{compEdge{e}, face{faceEdgeMap[e]}})
 		}
 
 		i++
@@ -89,6 +124,10 @@ type PointLocator struct {
 	outerFace *dcel.Face
 }
 
+func (spl *PointLocator) String() string {
+	return fmt.Sprintf("%v", spl.dp)
+}
+
 // PointLocate returns which face within this SlabPointLocator
 // the query point lands, within two dimensions.
 func (spl *PointLocator) PointLocate(vs ...float64) (*dcel.Face, error) {
@@ -98,34 +137,36 @@ func (spl *PointLocator) PointLocate(vs ...float64) (*dcel.Face, error) {
 	tree := spl.dp.AtInstant(vs[0])
 	p := geom.Point{vs[0], vs[1], 0}
 
-	e, f := tree.SearchDown(p, 0)
-	if e == nil {
-		return nil, nil
-	}
-	e2, f2 := tree.SearchUp(p, 0)
-	if geom.VerticalCompare(p, e.(compEdge)) == search.Greater {
-		return nil, nil
-	}
+	// e, f := tree.SearchDown(p, 0)
+	// if e == nil {
+	// 	return nil, nil
+	// }
+	_, f2 := tree.SearchUp(p, 0)
+	// if geom.VerticalCompare(p, e.(compEdge)) == search.Greater {
+	// 	return nil, nil
+	// }
 
-	if geom.VerticalCompare(p, e2.(compEdge)) == search.Less {
-		return nil, nil
-	}
+	// if geom.VerticalCompare(p, e2.(compEdge)) == search.Less {
+	// 	return nil, nil
+	// }
 
 	// We then do PIP on each face, and return
 	// whichever is true, if any.
-	f3 := f.(faces)
-	f4 := f2.(faces)
-	faces := []*dcel.Face{f3.f1, f3.f2, f4.f1, f4.f2}
-
-	for _, f5 := range faces {
-		if f5 != spl.outerFace {
-			visualize.HighlightColor = visualize.CheckFaceColor
-			visualize.DrawFace(f5)
-			if f5.Contains(p) {
-				return f5, nil
-			}
-		}
-	}
+	// f3 := f.(faces)
+	// f4 := f2.(faces)
+	// faces := []*dcel.Face{f3.f1, f3.f2, f4.f1, f4.f2}
+
+	// for _, f5 := range faces {
+	// 	if f5 != spl.outerFace {
+	// 		visualize.HighlightColor = visualize.CheckFaceColor
+	// 		visualize.DrawFace(f5)
+	// 		if f5.Contains(p) {
+	// 			return f5, nil
+	// 		}
+	// 	}
+	// }
+
+	// return nil, nil
 
-	return nil, nil
+	return f2.(face).Face, nil
 }

dcel/pointLoc/test/parse/benchParse.go

@@ -10,11 +10,11 @@ import (
 )
 
 func main() {
-	f, err := os.Open("benchOut2.txt")
+	f, err := os.Open("benchOut3.txt")
 	if err != nil {
 		log.Fatal(err)
 	}
-	outF, err := os.Create("benchOut2.csv")
+	outF, err := os.Create("benchOut3.csv")
 	if err != nil {
 		log.Fatal(err)
 	}
@@ -27,19 +27,19 @@ func main() {
 			break
 		}
 		actualLine := []byte{}
-			for _, b := range line {
-				if b != 0 {
-					actualLine = append(actualLine, b)
-				}
+		for _, b := range line {
+			if b != 0 {
+				actualLine = append(actualLine, byte(b))
 			}
+		}
 		strs := strings.Fields(line)
 		if strs[0] == "PASS" {
 			break
 		}
 		fmt.Println(strs[0])
 		csvStrings := make([]string, 9)
 		csvStrings[0] = strs[1]
-		for i := 0; i < 2; i++ {
+		for i := 0; i < 1; i++ {
 			line, err := bufIn.ReadBytes('\n')
 			if err != nil {
 				log.Fatal(err)

dcel/pointLoc/test/point-loc_test.go

@@ -20,9 +20,9 @@ import (
 )
 
 var (
-	inputSize   = 5
+	inputSize   = 25
 	inputRange  = 10000.0
-	testCt      = 10000
+	testCt      = 500
 	slabErrors  = 0
 	trapErrors  = 0
 	rtreeErrors = 0
@@ -44,6 +44,7 @@ func testRandomPts(t *testing.T, pl pointLoc.LocatesPoints, limit int, errs *int
 		if !assert.True(t, bruteForceContains) {
 			t.Log("Error point:", pt)
 			t.Log("Error face:", structIntersected)
+
 			(*errs)++
 		}
 	}
@@ -59,7 +60,7 @@ func TestRandomDCELSlab(t *testing.T) {
 	printErrors()
 }
 
-func TestRandomDCELSlabErrors(t *testing.T) {
+func TestDCELSlabErrors(t *testing.T) {
 	errCt := 0
 	subTestCt := 50
 	for i := 0; i < testCt; i++ {
@@ -74,6 +75,9 @@ func TestRandomDCELSlabErrors(t *testing.T) {
 		testRandomPts(t, structure, subTestCt, &queryErrors)
 		if queryErrors != 0 {
 			off.Save(dc).WriteFile("testFail" + strconv.Itoa(i) + ".off")
+			t.Log("Error index:", i)
+			s := structure.(*slab.PointLocator)
+			t.Log(s)
 			errCt++
 			continue
 		}
@@ -285,6 +289,5 @@ func BenchmarkAdditional(b *testing.B) {
 		seed = time.Now().UnixNano()
 		fmt.Println("InputSize:", i)
 		b.Run("TrapSetup", BenchmarkRandomSetupTrapezoid)
-		b.Run("Rtree", BenchmarkRandomDCELRtree)
 	}
 }

dcel/shapes.go

@@ -120,14 +120,13 @@ func FourPoint(p1, p2, p3, p4 geom.D3) *DCEL {
 	// Correcting for faces[0] = the infinite exterior
 	dc.Faces[0], dc.Faces[1] = dc.Faces[1], dc.Faces[0]
 
-	dc.CorrectDirectionality(dc.Faces[0])
-	dc.CorrectDirectionality(dc.Faces[1])
+	dc.CorrectDirectionalityAll()
 
 	return dc
 }
 
 func goodrandf64() float64 {
-	return (rand.Float64() * (8 / 10)) + .1
+	return (rand.Float64() * (8.0 / 10.0)) + .1
 }
 
 func Random2DDCEL(size float64, splits int) *DCEL {
@@ -267,6 +266,7 @@ func Random2DDCEL(size float64, splits int) *DCEL {
 	}
 	//fmt.Println("Random dcel end")
 	//fmt.Println(dc)
+	dc.CorrectDirectionalityAll()
 
 	return dc
 }

geom/eq.go

@@ -7,7 +7,7 @@ import "math"
 const (
 	// Epsilon could probably be smaller than this without
 	// causing problems, but we're being overly cautious.
-	ε = 1.0e-4
+	ε = 1.0e-3
 	// Inf is shorthand for math.MaxFloat64
 	Inf = math.MaxFloat64
 	// NegInf is shorthand for math.MaxFloat64 * -1

search/tree/bst.go

@@ -241,7 +241,6 @@ func (bst *BST) updateRoot(n *node) {
 	}
 	if n != nil {
 		bst.root = n
-		return
 	}
 	if bst.root == nil {
 		return