Mercurial > gemma

comparison cmd/isoareas/main.go @ 4542:56f4e8cbfab7 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Sew together triangles that totally inside a single class.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 26 Sep 2019 18:11:38 +0200
parents 53b55f811666
children 2a707732331f
comparison
equal deleted inserted replaced
4541:53b55f811666 4542:56f4e8cbfab7
13 13
14 package main 14 package main
15 15
16 import ( 16 import (
17 "bufio" 17 "bufio"
18 "container/list"
18 "fmt" 19 "fmt"
19 "io" 20 "io"
20 "log" 21 "log"
21 "math" 22 "math"
22 "math/bits" 23 "math/bits"
140 result := make([][]octree.LineStringZ, len(classes)) 141 result := make([][]octree.LineStringZ, len(classes))
141 142
142 log.Println("inside classes:") 143 log.Println("inside classes:")
143 for i, c := range classes { 144 for i, c := range classes {
144 145
145 var pb polygonBuilder 146 pb := polygonBuilder{open: list.New()}
146 147
147 var isolated, inside int 148 var isolated, inside int
148 allInClass: 149 allInClass:
149 for _, idx := range c { 150 for _, idx := range c {
150 ns := neighbors(tri, idx) 151 ns := neighbors(tri, idx)
151 w := where(ns, c) 152 mask := where(ns, c)
152 switch bits.OnesCount8(w) { 153 switch bits.OnesCount8(mask) {
153 case 0: 154 case 0:
154 // Totally insides do not contribute to the geometries. 155 // Totally insides do not contribute to the geometries.
155 inside++ 156 inside++
156 continue allInClass 157 continue allInClass
157 case 3: 158 case 3:
161 pb.addTriangle( 162 pb.addTriangle(
162 tri.Points[ti[0]], 163 tri.Points[ti[0]],
163 tri.Points[ti[1]], 164 tri.Points[ti[1]],
164 tri.Points[ti[2]]) 165 tri.Points[ti[2]])
165 continue allInClass 166 continue allInClass
166 } 167 default:
167 } 168 ti := tri.Triangles[idx*3 : idx*3+3]
168 169 for j := 0; j < 3; j++ {
169 log.Printf("\t%d: inside: %d / isolated: %d\n", 170 if mask&(1<<j) == 0 {
170 i, inside, isolated) 171
172 // TODO: Look into cuts to see
173 // if there are real intersections
174 k := (j + 1) % 3
175 front := tri.Points[ti[j]]
176 back := tri.Points[ti[k]]
177
178 curr := octree.LineStringZ{front, back}
179
180 segment:
181 for e := pb.open.Front(); e != nil; e = e.Next() {
182 line := e.Value.(octree.LineStringZ)
183 first, last := line[0], line[len(line)-1]
184
185 if back.EpsEquals(first) {
186 curr = append(curr, line[1:]...)
187 pb.open.Remove(e)
188 goto segment
189 }
190
191 if back.EpsEquals(last) {
192 line.Reverse()
193 curr = append(curr, line[1:]...)
194 pb.open.Remove(e)
195 goto segment
196 }
197
198 if front.EpsEquals(last) {
199 curr = append(line, curr[1:]...)
200 pb.open.Remove(e)
201 goto segment
202 }
203
204 if front.EpsEquals(first) {
205 line.Reverse()
206 curr = append(line, curr[1:]...)
207 pb.open.Remove(e)
208 goto segment
209 }
210 } // all open
211
212 // check if closed
213 if curr[0].EpsEquals(curr[len(curr)-1]) {
214 if !curr.CCW() {
215 curr.Reverse()
216 }
217 pb.polygons = append(pb.polygons, curr)
218 } else {
219 pb.open.PushFront(curr)
220 }
221 } // for all border parts
222 }
223 }
224 }
225
226 log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d\n",
227 i, inside, isolated, pb.open.Len(), len(pb.polygons))
171 228
172 result[i] = pb.polygons 229 result[i] = pb.polygons
173 } 230 }
174 231
175 log.Println("cuts:") 232 log.Println("cuts:")
183 240
184 } 241 }
185 242
186 type polygonBuilder struct { 243 type polygonBuilder struct {
187 polygons []octree.LineStringZ 244 polygons []octree.LineStringZ
245
246 open *list.List
188 } 247 }
189 248
190 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) { 249 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
191 polygon := octree.LineStringZ{v0, v1, v2, v0} 250 polygon := octree.LineStringZ{v0, v1, v2, v0}
192 pb.polygons = append(pb.polygons, polygon) 251 pb.polygons = append(pb.polygons, polygon)