golang で Douglas Peuker

内容

• Douglas Peuker アルゴリズムでの曲線間引きがどのくらい時間がかかるのか確認したかった
• ので、簡単に golang で書いてみた
• こちらを参考にさせていただいた http://www.trail-note.net/tech/thinout/
• だいたいの目安を見るための書捨てコードなのでツッコミどころが多いのはご愛嬌

ベンチマーク

$ go test -bench .
goos: darwin
goarch: amd64
pkg: xxxxxx
Benchmark10-4       	  200000	       7125 ns/op  ->  7.1 us
Benchmark100-4      	    1000	    1234958 ns/op  ->  1.2 ms
Benchmark1000-4     	      10	  129394222 ns/op  ->  0.1 ms
Benchmark10000-4    	       1	 2101666871 ns/op  ->  2.1  s
Benchmark100000-4   	       1	12936923196 ns/op  -> 12.9  s
PASS
ok  	xxxxxx	25.437s

• 1万点 → 1000点 で 秒のオーダー
• 動的に計算するのは結構キツイですね

ソースコード

db_test.go

package dp_test

import (
	"math"
	"math/rand"
	"testing"
)

type Setting struct {
	i1, i2 int
}

type Point struct {
	x       float64
	y       float64
	setting *Setting
}

type Max struct {
	dist   float64
	i      int
	i1, i2 int
}

func Benchmark10(b *testing.B) {
	bench(b, 10)
}

func Benchmark100(b *testing.B) {
	bench(b, 100)
}

func Benchmark1000(b *testing.B) {
	bench(b, 1000)
}

func Benchmark10000(b *testing.B) {
	bench(b, 10000)
}

func Benchmark100000(b *testing.B) {
	bench(b, 100000)
}

func bench(b *testing.B, num int) {
	limit := 1000

	for j := 0; j < b.N; j++ {

		// 初期化
		points := make([]Point, num)
		for i := range points {
			points[i].x = rand.Float64()
			points[i].y = rand.Float64()
		}

		if len(points) < 2 {
			return
		}

		for i := range points {
			points[i].setting = &Setting{
				i1: 0,
				i2: len(points) - 1,
			}
		}
		points[0].setting = nil
		points[len(points)-1].setting = nil

		if len(points) < limit {
			limit = len(points)
		}

		b.StartTimer()
		for num := 2; num < limit; num++ {
			var max Max
			for i, p := range points {
				if p.setting == nil {
					continue
				}
				x, y := points[i].x, points[i].y

				s := p.setting
				x1, x2 := points[s.i1].x, points[s.i2].x
				y1, y2 := points[s.i1].y, points[s.i2].y

				numerator := math.Abs((y2-y1)*x - (x2-x1)*y + y1*x2 - y2*x1)
				denominator := math.Sqrt(math.Pow((y2-y1), 2) + math.Pow((x2-x1), 2))
				dist := numerator / denominator

				if max.dist < dist {
					max = Max{
						dist: dist,
						i:    i,
						i1:   s.i1,
						i2:   s.i2,
					}
				}
			}
			points[max.i].setting = nil

			for i := max.i1 + 1; i < max.i; i++ {
				points[i].setting.i2 = max.i
			}

			for i := max.i + 1; i < max.i2; i++ {
				points[i].setting.i1 = max.i
			}
		}
		b.StopTimer()
	}
}

おまけ

ちょっと最適化したら 10万点 → 1000点 で 150msくらいでした

$ go test -bench .
goos: darwin
goarch: amd64
pkg: xxxxxx
Benchmark10-4        	  500000	      3341 ns/op ->   3us
Benchmark100-4       	   50000	     35322 ns/op ->  35us
Benchmark1000-4      	    2000	   1003985 ns/op ->   1ms
Benchmark10000-4     	     100	  13531037 ns/op ->  13ms
Benchmark100000-4    	      10	 152950556 ns/op -> 152ms
Benchmark1000000-4   	       1	2231156196 ns/op ->   2s
PASS
ok  	xxxxxx	26.530s

dp_test.go

package dp_test

import (
	"math"
	"math/rand"
	"testing"
)

type Setting struct {
	i1, i2 int
}

type Point struct {
	x        float64
	y        float64
	setting  *Setting
	distance float64
}

type Max struct {
	dist   float64
	i      int
	i1, i2 int
}

func Benchmark10(b *testing.B) {
	bench(b, 10)
}

func Benchmark100(b *testing.B) {
	bench(b, 100)
}

func Benchmark1000(b *testing.B) {
	bench(b, 1000)
}

func Benchmark10000(b *testing.B) {
	bench(b, 10000)
}

func Benchmark100000(b *testing.B) {
	bench(b, 100000)
}

func Benchmark1000000(b *testing.B) {
	bench(b, 1000000)
}

func bench(b *testing.B, num int) {
	limit := 1000

	for j := 0; j < b.N; j++ {

		// 初期化
		points := make([]Point, num)
		for i := range points {
			points[i].x = rand.Float64()
			points[i].y = rand.Float64()
			points[i].distance = math.NaN()
		}

		if len(points) < 2 {
			return
		}

		for i := range points {
			points[i].setting = &Setting{
				i1: 0,
				i2: len(points) - 1,
			}
		}
		points[0].setting = nil
		points[len(points)-1].setting = nil

		if len(points) < limit {
			limit = len(points)
		}

		b.StartTimer()
		for num := 2; num < limit; num++ {
			var max Max
			for i, p := range points {
				if !math.IsNaN(p.distance) {
					continue
				}

				if p.setting == nil {
					continue
				}

				s := p.setting
				x, y := points[i].x, points[i].y
				x1, x2 := points[s.i1].x, points[s.i2].x
				y1, y2 := points[s.i1].y, points[s.i2].y

				numerator := math.Abs((y2-y1)*x - (x2-x1)*y + y1*x2 - y2*x1)
				denominator := math.Sqrt(math.Pow((y2-y1), 2) + math.Pow((x2-x1), 2))
				dist := numerator / denominator
				points[i].distance = dist

				if max.dist < dist {
					max = Max{
						dist: dist,
						i:    i,
						i1:   s.i1,
						i2:   s.i2,
					}
				}
			}
			points[max.i].setting = nil

			for i := max.i1 + 1; i < max.i; i++ {
				points[i].setting.i2 = max.i
				points[i].distance = math.NaN()
			}

			for i := max.i + 1; i < max.i2; i++ {
				points[i].setting.i1 = max.i
				points[i].distance = math.NaN()
			}
		}
		b.StopTimer()
	}
}