1
2
3
4
5 package plotter
6
7 import (
8 "flag"
9 "fmt"
10 "math"
11 "reflect"
12 "sort"
13 "testing"
14
15 "golang.org/x/exp/rand"
16
17 "gonum.org/v1/gonum/mat"
18 "gonum.org/v1/plot"
19 "gonum.org/v1/plot/palette"
20 "gonum.org/v1/plot/vg"
21 )
22
23 var visualDebug = flag.Bool("visual", false, "output images for benchmarks and test data")
24
25 type unitGrid struct{ mat.Matrix }
26
27 func (g unitGrid) Dims() (c, r int) { r, c = g.Matrix.Dims(); return c, r }
28 func (g unitGrid) Z(c, r int) float64 { return g.Matrix.At(r, c) }
29 func (g unitGrid) X(c int) float64 {
30 _, n := g.Matrix.Dims()
31 if c < 0 || c >= n {
32 panic("index out of range")
33 }
34 return float64(c)
35 }
36 func (g unitGrid) Y(r int) float64 {
37 m, _ := g.Matrix.Dims()
38 if r < 0 || r >= m {
39 panic("index out of range")
40 }
41 return float64(r)
42 }
43
44 func TestHeatMapWithContour(t *testing.T) {
45 if !*visualDebug {
46 return
47 }
48 m := unitGrid{mat.NewDense(3, 4, []float64{
49 2, 1, 4, 3,
50 6, 7, 2, 5,
51 9, 10, 11, 12,
52 })}
53 h := NewHeatMap(m, palette.Heat(12, 1))
54
55 levels := []float64{1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5}
56 c := NewContour(m, levels, palette.Rainbow(10, palette.Blue, palette.Red, 1, 1, 1))
57 c.LineStyles[0].Width *= 5
58
59 plt := plot.New()
60 plt.Add(h)
61 plt.Add(c)
62 plt.Add(NewGlyphBoxes())
63
64 plt.X.Padding = 0
65 plt.Y.Padding = 0
66 plt.X.Max = 3.5
67 plt.Y.Max = 2.5
68
69 err := plt.Save(7*vg.Centimeter, 7*vg.Centimeter, "heat.svg")
70 if err != nil {
71 t.Fatalf("could not save plot: %+v", err)
72 }
73 }
74
75 func TestComplexContours(t *testing.T) {
76 rnd := rand.New(rand.NewSource(1))
77
78 if !*visualDebug {
79 return
80 }
81 for _, n := range []float64{0, 1, 2, 4, 8, 16, 32} {
82 data := make([]float64, 6400)
83 for i := range data {
84 r := float64(i/80) - 40
85 c := float64(i%80) - 40
86
87 data[i] = rnd.NormFloat64()*n + math.Hypot(r, c)
88 }
89
90 m := unitGrid{mat.NewDense(80, 80, data)}
91
92 levels := []float64{-1, 3, 7, 9, 13, 15, 19, 23, 27, 31}
93 c := NewContour(m, levels, palette.Rainbow(10, palette.Blue, palette.Red, 1, 1, 1))
94
95 plt := plot.New()
96 plt.X.Padding = 0
97 plt.Y.Padding = 0
98 plt.X.Max = 79.5
99 plt.Y.Max = 79.5
100
101 plt.Add(c)
102
103 err := plt.Save(7*vg.Centimeter, 7*vg.Centimeter, fmt.Sprintf("complex_contour-%v.svg", n))
104 if err != nil {
105 t.Fatalf("could not save plot: %+v", err)
106 }
107 }
108 }
109
110 func unity(f float64) vg.Length { return vg.Length(f) }
111
112 func BenchmarkComplexContour0(b *testing.B) { complexContourBench(0, b) }
113 func BenchmarkComplexContour1(b *testing.B) { complexContourBench(1, b) }
114 func BenchmarkComplexContour2(b *testing.B) { complexContourBench(2, b) }
115 func BenchmarkComplexContour4(b *testing.B) { complexContourBench(4, b) }
116 func BenchmarkComplexContour8(b *testing.B) { complexContourBench(8, b) }
117 func BenchmarkComplexContour16(b *testing.B) { complexContourBench(16, b) }
118 func BenchmarkComplexContour32(b *testing.B) { complexContourBench(32, b) }
119
120 var cp map[float64][]vg.Path
121
122 func complexContourBench(noise float64, b *testing.B) {
123 rnd := rand.New(rand.NewSource(1))
124
125 data := make([]float64, 6400)
126 for i := range data {
127 r := float64(i/80) - 40
128 c := float64(i%80) - 40
129
130 data[i] = rnd.NormFloat64()*noise + math.Hypot(r, c)
131 }
132
133 m := unitGrid{mat.NewDense(80, 80, data)}
134
135 levels := []float64{-1, 3, 7, 9, 13, 15, 19, 23, 27, 31}
136
137 var p map[float64][]vg.Path
138
139 b.ResetTimer()
140 for i := 0; i < b.N; i++ {
141 p = contourPaths(m, levels, unity, unity)
142 }
143
144 cp = p
145 }
146
147 func TestContourPaths(t *testing.T) {
148 m := unitGrid{mat.NewDense(3, 4, []float64{
149 2, 1, 4, 3,
150 6, 7, 2, 5,
151 9, 10, 11, 12,
152 })}
153
154 levels := []float64{1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5}
155
156 var (
157 wantClosed = 2
158 gotClosed int
159 )
160
161 got := contourPaths(m, levels, unity, unity)
162 for l, p := range got {
163 sort.Sort(byLength(p))
164 for i, c := range p {
165 if isLoop(c) && isLoop(wantContours[l][i]) {
166 if !circularPermutations(c[1:], wantContours[l][i][1:]) {
167 t.Errorf("unexpected path:\n\tgot:%+v\n\twant:%+v", c, wantContours[l][i])
168 }
169 } else if !reflect.DeepEqual(c, wantContours[l][i]) && !reflect.DeepEqual(c, reverseOfPath(wantContours[l][i])) {
170 t.Errorf("unexpected path:\n\tgot:%+v\n\twant:%+v", c, wantContours[l][i])
171 }
172
173 if isLoop(c) {
174 gotClosed++
175 }
176 }
177 }
178
179 if gotClosed != wantClosed {
180 t.Errorf("unexpected number of loops: got:%d want:%d", gotClosed, wantClosed)
181 }
182 }
183
184 type byLength []vg.Path
185
186 func (p byLength) Len() int { return len(p) }
187 func (p byLength) Less(i, j int) bool { return len(p[i]) < len(p[j]) }
188 func (p byLength) Swap(i, j int) { p[i], p[j] = p[j], p[i] }
189
190 func reverseOfPath(p vg.Path) vg.Path {
191 rp := make(vg.Path, 0, len(p))
192 for i := len(p) - 1; i >= 0; i-- {
193 rp = append(rp, p[i])
194 }
195 rp[0].Type = vg.MoveComp
196 rp[len(rp)-1].Type = vg.LineComp
197
198 return rp
199 }
200
201 func circularPermutations(a, b vg.Path) bool {
202 if len(a) != len(b) {
203 return false
204 }
205
206 var off int
207
208 var forward bool
209 for i, pc := range b {
210 if reflect.DeepEqual(a[0], pc) {
211 off = i
212 forward = true
213 break
214 }
215 }
216 for i, pc := range a {
217 if !reflect.DeepEqual(b[(off+i)%len(a)], pc) {
218 forward = false
219 break
220 }
221 }
222
223 var reverse bool
224 for i, pc := range b {
225 if reflect.DeepEqual(a[0], pc) {
226 off = i
227 reverse = true
228 break
229 }
230 }
231 for i, pc := range a {
232 if !reflect.DeepEqual(b[(off-i+len(a))%len(a)], pc) {
233 reverse = false
234 break
235 }
236 }
237
238 return forward || reverse
239 }
240
241
242 var wantContours = map[float64][]vg.Path{
243 1.5: {
244 {
245 {Type: vg.MoveComp, Pos: vg.Point{X: 1.1666666666666667, Y: 0}},
246 {Type: vg.LineComp, Pos: vg.Point{X: 1.1, Y: 0.1}},
247 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.08333333333333333}},
248 {Type: vg.LineComp, Pos: vg.Point{X: 0.9166666666666666, Y: 0.08333333333333333}},
249 {Type: vg.LineComp, Pos: vg.Point{X: 0.5, Y: 0}},
250 },
251 },
252 2.5: {
253 {
254 {Type: vg.MoveComp, Pos: vg.Point{X: 1.5, Y: 0}},
255 {Type: vg.LineComp, Pos: vg.Point{X: 1.3, Y: 0.3}},
256 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.25}},
257 {Type: vg.LineComp, Pos: vg.Point{X: 0.75, Y: 0.25}},
258 {Type: vg.LineComp, Pos: vg.Point{X: 0.125, Y: 0.125}},
259 {Type: vg.LineComp, Pos: vg.Point{X: 0, Y: 0.125}},
260 },
261 {
262 {Type: vg.MoveComp, Pos: vg.Point{X: 2, Y: 1.0555555555555556}},
263 {Type: vg.LineComp, Pos: vg.Point{X: 1.9545454545454546, Y: 1.0454545454545454}},
264 {Type: vg.LineComp, Pos: vg.Point{X: 1.9, Y: 1}},
265 {Type: vg.LineComp, Pos: vg.Point{X: 1.8333333333333333, Y: 0.8333333333333334}},
266 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 0.75}},
267 {Type: vg.LineComp, Pos: vg.Point{X: 2.1666666666666665, Y: 0.8333333333333334}},
268 {Type: vg.LineComp, Pos: vg.Point{X: 2.1666666666666665, Y: 1}},
269 {Type: vg.LineComp, Pos: vg.Point{X: 2.0454545454545454, Y: 1.0454545454545454}},
270 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.0555555555555556}},
271 },
272 },
273 3.5: {
274 {
275 {Type: vg.MoveComp, Pos: vg.Point{X: 3, Y: 0.25}},
276 {Type: vg.LineComp, Pos: vg.Point{X: 2.5, Y: 0.5}},
277 {Type: vg.LineComp, Pos: vg.Point{X: 2.5, Y: 0}},
278 },
279 {
280 {Type: vg.MoveComp, Pos: vg.Point{X: 1.8333333333333333, Y: 0}},
281 {Type: vg.LineComp, Pos: vg.Point{X: 1.5, Y: 0.5}},
282 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.4166666666666667}},
283 {Type: vg.LineComp, Pos: vg.Point{X: 0.5833333333333334, Y: 0.4166666666666667}},
284 {Type: vg.LineComp, Pos: vg.Point{X: 0.375, Y: 0.375}},
285 {Type: vg.LineComp, Pos: vg.Point{X: 0, Y: 0.375}},
286 },
287 {
288 {Type: vg.MoveComp, Pos: vg.Point{X: 1.5, Y: 0.5}},
289 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 0.25}},
290 {Type: vg.LineComp, Pos: vg.Point{X: 2.5, Y: 0.5}},
291 {Type: vg.LineComp, Pos: vg.Point{X: 2.5, Y: 1}},
292 {Type: vg.LineComp, Pos: vg.Point{X: 2.1363636363636362, Y: 1.1363636363636365}},
293 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.1666666666666667}},
294 {Type: vg.LineComp, Pos: vg.Point{X: 1.8636363636363635, Y: 1.1363636363636365}},
295 {Type: vg.LineComp, Pos: vg.Point{X: 1.7, Y: 1}},
296 {Type: vg.LineComp, Pos: vg.Point{X: 1.5, Y: 0.5}},
297 },
298 },
299 4.5: {
300 {
301 {Type: vg.MoveComp, Pos: vg.Point{X: 0, Y: 0.625}},
302 {Type: vg.LineComp, Pos: vg.Point{X: 0.375, Y: 0.625}},
303 {Type: vg.LineComp, Pos: vg.Point{X: 0.5833333333333334, Y: 0.5833333333333334}},
304 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.5833333333333334}},
305 {Type: vg.LineComp, Pos: vg.Point{X: 1.3571428571428572, Y: 0.6428571428571429}},
306 {Type: vg.LineComp, Pos: vg.Point{X: 1.5, Y: 1}},
307 {Type: vg.LineComp, Pos: vg.Point{X: 1.7727272727272727, Y: 1.2272727272727273}},
308 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.2777777777777777}},
309 {Type: vg.LineComp, Pos: vg.Point{X: 2.227272727272727, Y: 1.2272727272727273}},
310 {Type: vg.LineComp, Pos: vg.Point{X: 2.8333333333333335, Y: 1}},
311 {Type: vg.LineComp, Pos: vg.Point{X: 2.8333333333333335, Y: 0.8333333333333334}},
312 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 0.75}},
313 },
314 },
315 5.5: {
316 {
317 {Type: vg.MoveComp, Pos: vg.Point{X: 0, Y: 0.875}},
318 {Type: vg.LineComp, Pos: vg.Point{X: 0.125, Y: 0.875}},
319 {Type: vg.LineComp, Pos: vg.Point{X: 0.75, Y: 0.75}},
320 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.75}},
321 {Type: vg.LineComp, Pos: vg.Point{X: 1.2142857142857142, Y: 0.7857142857142857}},
322 {Type: vg.LineComp, Pos: vg.Point{X: 1.3, Y: 1}},
323 {Type: vg.LineComp, Pos: vg.Point{X: 1.6818181818181819, Y: 1.3181818181818181}},
324 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.3888888888888888}},
325 {Type: vg.LineComp, Pos: vg.Point{X: 2.3181818181818183, Y: 1.3181818181818181}},
326 {Type: vg.LineComp, Pos: vg.Point{X: 2.9, Y: 1.1}},
327 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.0714285714285714}},
328 },
329 },
330 6.5: {
331 {
332 {Type: vg.MoveComp, Pos: vg.Point{X: 0, Y: 1.1666666666666667}},
333 {Type: vg.LineComp, Pos: vg.Point{X: 0.125, Y: 1.125}},
334 {Type: vg.LineComp, Pos: vg.Point{X: 0.5, Y: 1}},
335 {Type: vg.LineComp, Pos: vg.Point{X: 0.9166666666666666, Y: 0.9166666666666666}},
336 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 0.9166666666666666}},
337 {Type: vg.LineComp, Pos: vg.Point{X: 1.0714285714285714, Y: 0.9285714285714286}},
338 {Type: vg.LineComp, Pos: vg.Point{X: 1.1, Y: 1}},
339 {Type: vg.LineComp, Pos: vg.Point{X: 1.5909090909090908, Y: 1.4090909090909092}},
340 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.5}},
341 {Type: vg.LineComp, Pos: vg.Point{X: 2.409090909090909, Y: 1.4090909090909092}},
342 {Type: vg.LineComp, Pos: vg.Point{X: 2.7, Y: 1.3}},
343 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.2142857142857142}},
344 },
345 },
346 7.5: {
347 {
348 {Type: vg.MoveComp, Pos: vg.Point{X: 0, Y: 1.5}},
349 {Type: vg.LineComp, Pos: vg.Point{X: 0.375, Y: 1.375}},
350 {Type: vg.LineComp, Pos: vg.Point{X: 0.75, Y: 1.25}},
351 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 1.1666666666666667}},
352 {Type: vg.LineComp, Pos: vg.Point{X: 1.5, Y: 1.5}},
353 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.6111111111111112}},
354 {Type: vg.LineComp, Pos: vg.Point{X: 2.5, Y: 1.5}},
355 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.3571428571428572}},
356 },
357 },
358 8.5: {
359 {
360 {Type: vg.MoveComp, Pos: vg.Point{X: 0, Y: 1.8333333333333333}},
361 {Type: vg.LineComp, Pos: vg.Point{X: 0.25, Y: 1.75}},
362 {Type: vg.LineComp, Pos: vg.Point{X: 0.625, Y: 1.625}},
363 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 1.5}},
364 {Type: vg.LineComp, Pos: vg.Point{X: 1.3, Y: 1.7}},
365 {Type: vg.LineComp, Pos: vg.Point{X: 1.6428571428571428, Y: 1.6428571428571428}},
366 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.7222222222222223}},
367 {Type: vg.LineComp, Pos: vg.Point{X: 2.357142857142857, Y: 1.6428571428571428}},
368 {Type: vg.LineComp, Pos: vg.Point{X: 2.611111111111111, Y: 1.6111111111111112}},
369 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.5}},
370 },
371 },
372 9.5: {
373 {
374 {Type: vg.MoveComp, Pos: vg.Point{X: 0.5, Y: 2}},
375 {Type: vg.LineComp, Pos: vg.Point{X: 0.875, Y: 1.875}},
376 {Type: vg.LineComp, Pos: vg.Point{X: 1, Y: 1.8333333333333333}},
377 {Type: vg.LineComp, Pos: vg.Point{X: 1.1, Y: 1.9}},
378 {Type: vg.LineComp, Pos: vg.Point{X: 1.7857142857142858, Y: 1.7857142857142858}},
379 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.8333333333333333}},
380 {Type: vg.LineComp, Pos: vg.Point{X: 2.2142857142857144, Y: 1.7857142857142858}},
381 {Type: vg.LineComp, Pos: vg.Point{X: 2.7222222222222223, Y: 1.7222222222222223}},
382 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.6428571428571428}},
383 },
384 },
385 10.5: {
386 {
387 {Type: vg.MoveComp, Pos: vg.Point{X: 1.5, Y: 2}},
388 {Type: vg.LineComp, Pos: vg.Point{X: 1.9285714285714286, Y: 1.9285714285714286}},
389 {Type: vg.LineComp, Pos: vg.Point{X: 2, Y: 1.9444444444444444}},
390 {Type: vg.LineComp, Pos: vg.Point{X: 2.0714285714285716, Y: 1.9285714285714286}},
391 {Type: vg.LineComp, Pos: vg.Point{X: 2.8333333333333335, Y: 1.8333333333333333}},
392 {Type: vg.LineComp, Pos: vg.Point{X: 3, Y: 1.7857142857142858}},
393 },
394 },
395 }
396
397 var loopTests = []struct {
398 c *contour
399
400 want []*contour
401 }{
402 {
403 c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {7, 7}, {8, 8}, {9, 9}}},
404 want: []*contour{
405 {backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {7, 7}, {8, 8}, {9, 9}}},
406 },
407 },
408 {
409 c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {4, 4}, {7, 7}, {8, 8}, {9, 9}}},
410 want: []*contour{
411 {backward: path{{4, 4}}, forward: path{{5, 5}, {6, 6}, {4, 4}}},
412 {backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {7, 7}, {8, 8}, {9, 9}}},
413 },
414 },
415 {
416 c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {3, 3}, {7, 7}, {1, 1}, {9, 9}}},
417 want: []*contour{
418 {backward: path{{0, 0}}, forward: path{{1, 1}, {9, 9}}},
419 {backward: path{{3, 3}}, forward: path{{4, 4}, {5, 5}, {3, 3}}},
420 {backward: path{{1, 1}}, forward: path{{2, 2}, {3, 3}, {7, 7}, {1, 1}}},
421 },
422 },
423 {
424 c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {2, 2}, {7, 7}, {2, 2}, {9, 9}}},
425 want: []*contour{
426 {backward: path{{2, 2}}, forward: path{{7, 7}, {2, 2}}},
427 {backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {9, 9}}},
428 {backward: path{{2, 2}}, forward: path{{3, 3}, {4, 4}, {5, 5}, {2, 2}}},
429 },
430 },
431 {
432
433 c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {3, 3}, {8, 8}, {9, 9}, {5, 5}, {10, 10}}},
434 want: []*contour{
435 {backward: path{{5, 5}}, forward: path{{10, 10}}},
436 {backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}}},
437 {backward: path{{3, 3}}, forward: path{{4, 4}, {5, 5}, {6, 6}, {3, 3}}},
438 {backward: path{{3, 3}}, forward: path{{8, 8}, {9, 9}, {5, 5}, {6, 6}, {3, 3}}},
439 },
440 },
441 }
442
443 func (c testContour) String() string {
444 var s string
445 for i, p := range c {
446 if i != 0 {
447 s += ", "
448 }
449 s += fmt.Sprintf("%v", append(p.backward.reverse(), p.forward...))
450 p.backward.reverse()
451 }
452 return s
453 }
454
455 func TestExciseLoops(t *testing.T) {
456 for _, quick := range []bool{true, false} {
457 for i, test := range loopTests {
458 gotSet := make(contourSet)
459 c := &contour{
460 backward: append(path(nil), test.c.backward...),
461 forward: append(path(nil), test.c.forward...),
462 }
463 gotSet[c] = struct{}{}
464 c.exciseLoops(gotSet, quick)
465 var got []*contour
466 for c := range gotSet {
467 got = append(got, c)
468 }
469 sort.Sort(testContour(got))
470 if !reflect.DeepEqual(got, test.want) {
471 t.Errorf("unexpected loop excision result for %d quick=%t:\n\tgot:%v\n\twant:%v",
472 i, quick, testContour(got), testContour(test.want))
473 }
474 }
475 }
476 }
477
478 type testContour []*contour
479
480 func (c testContour) Len() int { return len(c) }
481 func (c testContour) Less(i, j int) bool { return len(c[i].forward) < len(c[j].forward) }
482 func (c testContour) Swap(i, j int) { c[i], c[j] = c[j], c[i] }
483
View as plain text