Mercurial > gemma
view pkg/octree/vertex.go @ 2130:f3aabc05f9b2
Fix constraints on waterway profiles
staging_done in the UNIQUE constraint had no effect, because the
exclusion constraint prevented two rows with equal location and
validity anyhow. Adding staging_done to the exclusion constraint
makes the UNIQUE constraint checking only a corner case of what
the exclusion constraint checks. Thus, remove the UNIQUE constraint.
Casting staging_done to int is needed because there is no appropriate
operator class for booleans. Casting to smallint or even bit would have
been better (i.e. should result in smaller index size), but that would
have required creating such a CAST, in addition.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Wed, 06 Feb 2019 15:42:32 +0100 |
parents | f4dcbe8941a1 |
children | a1e751c08c56 |
line wrap: on
line source
// This is Free Software under GNU Affero General Public License v >= 3.0 // without warranty, see README.md and license for details. // // SPDX-License-Identifier: AGPL-3.0-or-later // License-Filename: LICENSES/AGPL-3.0.txt // // Copyright (C) 2018 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package octree import ( "bytes" "encoding/binary" "io" "log" "math" "sort" ) type ( // Vertex is a 3D vertex. Vertex struct { X float64 Y float64 Z float64 } // Triangle is a triangle consisting of three vertices. Triangle [3]Vertex // Line is a line defined by first vertex on that line // and the second being the direction. Line [2]Vertex // MultiPointZ is a set of vertices. MultiPointZ []Vertex // LineStringZ is a line string formed of vertices. LineStringZ []Vertex // MultiLineStringZ is a set of line strings. MultiLineStringZ []LineStringZ // Box2D is 2D area from (X1, Y1) to (X2, Y2). Box2D struct { X1 float64 Y1 float64 X2 float64 Y2 float64 } // Plane2D is a 2D plane (a line in 2D space). Plane2D struct { A float64 B float64 C float64 } ) // Minimize adjust this vertex v to hold the minimum // values component-wise of v and w. func (v *Vertex) Minimize(w Vertex) { if w.X < v.X { v.X = w.X } if w.Y < v.Y { v.Y = w.Y } if w.Z < v.Z { v.Z = w.Z } } // Maximize adjust this vertex v to hold the maximum // values component-wise of v and w. func (v *Vertex) Maximize(w Vertex) { if w.X > v.X { v.X = w.X } if w.Y > v.Y { v.Y = w.Y } if w.Z > v.Z { v.Z = w.Z } } // Sub returns (v - w) component-wise. func (v Vertex) Sub(w Vertex) Vertex { return Vertex{ v.X - w.X, v.Y - w.Y, v.Z - w.Z, } } // Add returns (v + w) component-wise. func (v Vertex) Add(w Vertex) Vertex { return Vertex{ v.X + w.X, v.Y + w.Y, v.Z + w.Z, } } // Scale returns s*v component-wise. func (v Vertex) Scale(s float64) Vertex { return Vertex{ s * v.X, s * v.Y, s * v.Z, } } // Interpolate returns a function that return s*v2 + v1 // component-wise. func Interpolate(v1, v2 Vertex) func(Vertex) Vertex { v2 = v2.Sub(v1) return func(s Vertex) Vertex { return Vertex{ v2.X*s.X + v1.X, v2.Y*s.Y + v1.Y, v2.Z*s.Z + v1.Z, } } } // Less returns if one of v component is less than the // corresponing component in w. func (v Vertex) Less(w Vertex) bool { return v.X < w.X || v.Y < w.Y || v.Z < w.Z } // NewLine return a line of point/direction. func NewLine(p1, p2 Vertex) Line { return Line{ p2.Sub(p1), p1, } } // Eval returns the vertex for t*l[0] + l[1]. func (l Line) Eval(t float64) Vertex { return l[0].Scale(t).Add(l[1]) } // IntersectHorizontal returns the intersection point // for a given z value. func (l Line) IntersectHorizontal(h float64) Vertex { t := (h - l[1].Z) / l[0].Z return l.Eval(t) } func side(z, h float64) int { switch { case z < h: return -1 case z > h: return +1 } return 0 } // IntersectHorizontal calculates the line string that // results when cutting a triangle a a certain height. // Can be empty (on intersection), // one vertex (only touching an vertex) or // two vertices (real intersection). func (t *Triangle) IntersectHorizontal(h float64) LineStringZ { sides := [3]int{ side(t[0].Z, h), side(t[1].Z, h), side(t[2].Z, h), } var points LineStringZ for i := 0; i < 3; i++ { j := (i + 1) % 3 si := sides[i] sj := sides[j] switch { case si == 0 && sj == 0: // both on plane points = append(points, t[i], t[j]) case si == 0 && sj != 0: // first on plane points = append(points, t[i]) case si != 0 && sj == 0: // second on plane points = append(points, t[j]) case si == sj: // both on same side default: // real intersection v := NewLine(t[i], t[j]).IntersectHorizontal(h) points = append(points, v) } } return points } func linearScale(x1, y1, x2, y2 float64) func(Vertex) float64 { dx := x2 - x1 dy := y2 - y1 switch { case dx != 0: return func(v Vertex) float64 { return (v.X - x1) / dx } case dy != 0: return func(v Vertex) float64 { return (v.Y - y1) / dy } default: return func(Vertex) float64 { return 0 } } } func (ls LineStringZ) order(position func(Vertex) float64) { type posVertex struct { pos float64 v Vertex } positions := make([]posVertex, len(ls)) for i, v := range ls { positions[i] = posVertex{position(v), v} } sort.Slice(positions, func(i, j int) bool { return positions[i].pos < positions[j].pos }) for i := range positions { ls[i] = positions[i].v } } // EpsEquals returns true if v and w are equal component-wise // with the values within a epsilon range. func (v Vertex) EpsEquals(w Vertex) bool { const eps = 1e-5 return math.Abs(v.X-w.X) < eps && math.Abs(v.Y-w.Y) < eps && math.Abs(v.Z-w.Z) < eps } // JoinOnLine joins the the elements of a given multi line string // under the assumption that the segments are all on the line segment // from (x1, y1) to (x2, y2). func (mls MultiLineStringZ) JoinOnLine(x1, y1, x2, y2 float64) MultiLineStringZ { position := linearScale(x1, y1, x2, y2) type posLineString struct { pos float64 line LineStringZ } positions := make([]posLineString, 0, len(mls)) for _, ls := range mls { if len(ls) == 0 { continue } // order the atoms first ls.order(position) positions = append(positions, posLineString{position(ls[0]), ls}) } sort.Slice(positions, func(i, j int) bool { return positions[i].pos < positions[j].pos }) out := make(MultiLineStringZ, 0, len(positions)) var ignored int for i := range positions { curr := positions[i].line if l := len(out); l > 0 { last := out[l-1] if last[len(last)-1].EpsEquals(curr[0]) { out[l-1] = append(last[:len(last)-1], curr...) continue } if position(last[len(last)-1]) > position(curr[0]) { ignored++ continue } } out = append(out, curr) } log.Printf("info: ignored parts: %d\n", ignored) return out } // Write writes a Vertex as three 64 bit values in little endian order // to the given writer. func (v *Vertex) Write(w io.Writer) error { if err := binary.Write( w, binary.LittleEndian, math.Float64bits(v.X)); err != nil { return err } if err := binary.Write( w, binary.LittleEndian, math.Float64bits(v.Y)); err != nil { return err } return binary.Write( w, binary.LittleEndian, math.Float64bits(v.Z)) } // Read fills this vertex with three 64 bit values stored as // little endian from the given reader. func (v *Vertex) Read(r io.Reader) error { var buf [8]byte b := buf[:] if _, err := io.ReadFull(r, b); err != nil { return nil } v.X = math.Float64frombits(binary.LittleEndian.Uint64(b)) if _, err := io.ReadFull(r, b); err != nil { return nil } v.Y = math.Float64frombits(binary.LittleEndian.Uint64(b)) if _, err := io.ReadFull(r, b); err != nil { return nil } v.Z = math.Float64frombits(binary.LittleEndian.Uint64(b)) return nil } // AsWKB returns the WKB representation of the given multi line string. func (mls MultiLineStringZ) AsWKB() []byte { // pre-calculate size to avoid reallocations. size := 1 + 4 + 4 for _, ml := range mls { size += 1 + 4 + 4 + len(ml)*3*8 } buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbMultiLineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbLineStringZ) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) } } return buf.Bytes() } // AsWKB2D returns the WKB representation of the given multi line string // leaving the z component out. func (mls MultiLineStringZ) AsWKB2D() []byte { // pre-calculate size to avoid reallocations. size := 1 + 4 + 4 for _, ml := range mls { size += 1 + 4 + 4 + len(ml)*2*8 } buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbMultiLineString) binary.Write(buf, binary.LittleEndian, uint32(len(mls))) for _, ml := range mls { binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbLineString) binary.Write(buf, binary.LittleEndian, uint32(len(ml))) for _, p := range ml { binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) } } return buf.Bytes() } // Join joins two lines leaving the first of the second out. func (ls LineStringZ) Join(other LineStringZ) LineStringZ { nline := make(LineStringZ, len(ls)+len(other)-1) copy(nline, ls) copy(nline[len(ls):], other[1:]) return nline } // Merge merges line segments of a given multi line string // by finding common start and end vertices. func (mls MultiLineStringZ) Merge() MultiLineStringZ { var out MultiLineStringZ min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} for _, line := range mls { for _, v := range line { min.Minimize(v) } } type point struct { x int64 y int64 } const precision = 1e7 quant := func(v Vertex) point { return point{ x: int64(math.Round((v.X - min.X) * precision)), y: int64(math.Round((v.Y - min.Y) * precision)), } } heads := make(map[point]*[]LineStringZ) for _, line := range mls { if len(line) < 2 { out = append(out, line) continue } head := quant(line[0]) tail := quant(line[len(line)-1]) if head == tail { // its already a ring out = append(out, line) continue } if hs := heads[tail]; hs != nil { l := len(*hs) last := (*hs)[l-1] if l == 1 { delete(heads, tail) } else { (*hs)[l-1] = nil *hs = (*hs)[:l-1] } line = line.Join(last) if head == quant(line[len(line)-1]) { // its a ring out = append(out, line) continue } } if hs := heads[head]; hs != nil { *hs = append(*hs, line) } else { heads[head] = &[]LineStringZ{line} } } again: for head, lines := range heads { for i, line := range *lines { tail := quant(line[len(line)-1]) for hs := heads[tail]; hs != nil && len(*hs) > 0; hs = heads[tail] { l := len(*hs) last := (*hs)[l-1] (*hs)[l-1] = nil *hs = (*hs)[:l-1] line = line.Join(last) if tail = quant(line[len(line)-1]); head == tail { // its a ring out = append(out, line) // remove from current lines copy((*lines)[i:], (*lines)[i+1:]) (*lines)[len(*lines)-1] = nil *lines = (*lines)[:len(*lines)-1] goto again } // overwrite in current lines (*lines)[i] = line } } } // rings := len(out) // The rest are open line strings. for _, lines := range heads { for _, line := range *lines { out = append(out, line) } } // log.Printf("segments before/after merge: %d/%d (%d rings)\n", // len(mls), len(out), rings) return out } // Intersects checks if two Box2Ds intersect. func (a Box2D) Intersects(b Box2D) bool { return !(a.X2 < a.X1 || a.X2 < b.X1 || a.Y2 < a.Y1 || a.Y2 < b.Y1) } // Xi returns the i-th x component. func (a Box2D) Xi(i int) float64 { if i == 0 { return a.X1 } return a.X2 } // Yi returns the i-th y component. func (a Box2D) Yi(i int) float64 { if i == 0 { return a.Y1 } return a.Y2 } // NewPlane2D creates a new Plane2D from two given points. func NewPlane2D(x1, y1, x2, y2 float64) Plane2D { b := x2 - x1 a := -(y2 - y1) l := math.Sqrt(a*a + b*b) a /= l b /= l // a*x1 + b*y1 + c = 0 c := -(a*x1 + b*y1) return Plane2D{a, b, c} } // Eval determines the distance of a given point // from the plane. The sign of the result indicates // the sideness. func (p Plane2D) Eval(x, y float64) float64 { return p.A*x + p.B*y + p.C } const epsPlane = 1e-5 func sides(s int, x float64) int { if math.Signbit(x) { return s | 2 } return s | 1 } // IntersectsPlane checks if a Box2D intersects with // a given Plane2D. func (a Box2D) IntersectsPlane(p Plane2D) bool { var s int for i := 0; i < 2; i++ { x := a.Xi(i) for j := 0; j < 2; j++ { y := a.Yi(j) v := p.Eval(x, y) if math.Abs(v) < epsPlane { //log.Printf("on line") return true } if s = sides(s, v); s == 3 { //log.Printf("... on both sides (djt)") return true } } } //log.Printf("side: %d\n", s) return false } // Cross calculates the cross product of two vertices. func (v Vertex) Cross(w Vertex) Vertex { return Vertex{ v.Y*w.Z - v.Z*w.Y, v.Z*w.X - v.X*w.Z, v.X*w.Y - v.Y*w.X, } } // Intersection calcultes the 2D intersection point of // two Plane2Ds. If they do not intersect the returned // bool flags is set to false. func (p Plane2D) Intersection(o Plane2D) (float64, float64, bool) { u1 := Vertex{p.A, p.B, p.C} u2 := Vertex{o.A, o.B, o.C} plane := u1.Cross(u2) if plane.Z == 0 { return 0, 0, false } return plane.X / plane.Z, plane.Y / plane.Z, true } // VerticalLine is a 2D line segment. type VerticalLine struct { x1 float64 y1 float64 x2 float64 y2 float64 line Plane2D } // NewVerticalLine creates a new 2D line segment. func NewVerticalLine(x1, y1, x2, y2 float64) *VerticalLine { return &VerticalLine{ x1: x1, y1: y1, x2: x2, y2: y2, line: NewPlane2D(x1, y1, x2, y2), } } func onPlane(x float64) bool { return math.Abs(x) < epsPlane } func relative(v1, v2 Vertex) func(x, y float64) float64 { ls := linearScale(v1.X, v1.Y, v2.X, v2.Y) return func(x, y float64) float64 { return ls(Vertex{x, y, 0}) } } func interpolate(a, b float64) func(float64) float64 { return func(x float64) float64 { return (b-a)*x + a } } func linear(v1, v2 Vertex) func(t float64) Vertex { return func(t float64) Vertex { return Vertex{ (v2.X-v1.X)*t + v1.X, (v2.Y-v1.Y)*t + v1.Y, (v2.Z-v1.Z)*t + v1.Z, } } } func inRange(a float64) bool { return 0 <= a && a <= 1 } // Intersection intersects a line segment with a triangle. func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ { var out LineStringZ //defer func() { log.Printf("length out: %d\n", len(out)) }() edges: for i := 0; i < 3 && len(out) < 2; i++ { j := (i + 1) % 3 edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y) s1 := edge.Eval(vl.x1, vl.y1) s2 := edge.Eval(vl.x2, vl.y2) o1 := onPlane(s1) o2 := onPlane(s2) // log.Printf("s1, s2: %t %t (%f %f)\n", o1, o2, s1, s2) switch { case o1 && o2: pos := relative(t[i], t[j]) t1 := pos(vl.x1, vl.y1) t2 := pos(vl.x2, vl.y2) r1 := inRange(t1) r2 := inRange(t2) switch { case r1 && r2: lin := linear(t[i], t[j]) out = append(out, lin(t1), lin(t2)) return out case !r1 && !r2: // the hole edge out = append(out, t[i], t[j]) return out case !r1: if t1 < 0 { // below first -> clip by first lin := linear(t[i], t[j]) out = append(out, t[i], lin(t2)) } else { // above second -> clip by second lin := linear(t[i], t[j]) out = append(out, lin(t2), t[j]) } return out case !r2: if t2 < 0 { // below first -> clip by first lin := linear(t[i], t[j]) out = append(out, t[i], lin(t1)) } else { // above second -> clip by second lin := linear(t[i], t[j]) out = append(out, lin(t1), t[j]) } return out } case o1: t1 := relative(t[i], t[j])(vl.x1, vl.y1) if !inRange(t1) { continue edges } out = append(out, linear(t[i], t[j])(t1)) case o2: t2 := relative(t[i], t[j])(vl.x2, vl.y2) if !inRange(t2) { continue edges } out = append(out, linear(t[i], t[j])(t2)) default: x, y, intersects := vl.line.Intersection(edge) if !intersects { continue edges } // log.Println("Intersection -----------------------------") t1 := relative(t[i], t[j])(x, y) // log.Printf("relative pos: %f\n", t1) if !inRange(t1) { continue edges } // log.Println("valid ***************") z := interpolate(t[j].Z, t[i].Z)(t1) n := Vertex{x, y, z} if math.Signbit(s1) != math.Signbit(s2) { // log.Println("\ton different sides") // different sides -> vertex on edge. out = append(out, n) } else { // same side -> inside. Find peer if len(out) > 0 { // we already have our peer out = append(out, n) return out } for k := 0; k < 3; k++ { if k == i { continue } l := (k + 1) % 3 other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y) xo, yo, intersects := vl.line.Intersection(other) if !intersects { continue } t2 := relative(t[k], t[l])(xo, yo) if !inRange(t2) { continue } zo := interpolate(t[k].Z, t[l].Z)(t2) m := Vertex{xo, yo, zo} var xn, yn, xf, yf float64 // Which is nearer to current edge? if math.Abs(s1) < math.Abs(s2) { xn, yn = vl.x1, vl.y1 xf, yf = vl.x2, vl.y2 } else { xn, yn = vl.x2, vl.y2 xf, yf = vl.x1, vl.y1 } if onPlane(other.Eval(xn, yn)) { // triangle intersect in only point // treat as no intersection break edges } pos := relative(n, m) tzn := pos(xn, yn) tzf := pos(xf, yf) if !inRange(tzn) { // if near is out of range far is, too. return out } lin := interpolate(n.Z, m.Z) if inRange(tzf) { m.X = xf m.Y = yf m.Z = lin(tzf) } // else m is clipping n.X = xn n.Y = yn n.Z = lin(tzn) out = append(out, n, m) return out } } } } // supress single point touches. if len(out) == 1 { out = out[:0] } return out } // AsWKB returns a WKB representation of the given point cloud. func (mpz MultiPointZ) AsWKB() []byte { size := 1 + 4 + 4 + len(mpz)*(1+4+3*8) buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbMultiPointZ) binary.Write(buf, binary.LittleEndian, uint32(len(mpz))) for _, p := range mpz { binary.Write(buf, binary.LittleEndian, wkbNDR) binary.Write(buf, binary.LittleEndian, wkbPointZ) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.X)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Y)) binary.Write(buf, binary.LittleEndian, math.Float64bits(p.Z)) } return buf.Bytes() }