From 9cbd808449c5d1660b013e13f84f9f234811cee2 Mon Sep 17 00:00:00 2001 From: Dustin Oprea Date: Thu, 19 Jan 2017 20:57:15 -0500 Subject: [PATCH 1/2] Added uint64 support. --- hilbert64.go | 109 +++++++++++++++++++ hilbert64_test.go | 235 +++++++++++++++++++++++++++++++++++++++++ peano64.go | 103 ++++++++++++++++++ peano64_test.go | 261 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 708 insertions(+) create mode 100644 hilbert64.go create mode 100644 hilbert64_test.go create mode 100644 peano64.go create mode 100644 peano64_test.go diff --git a/hilbert64.go b/hilbert64.go new file mode 100644 index 0000000..b124cf6 --- /dev/null +++ b/hilbert64.go @@ -0,0 +1,109 @@ +// Copyright 2015 Google Inc. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Package hilbert is for mapping values to and from space-filling curves, such as Hilbert and Peano +// curves. +package hilbert + +// Hilbert represents a 2D Hilbert space of order N for mapping to and from. +// Implements SpaceFilling interface. +type Hilbert64 struct { + N uint64 +} + +// NewHilbert returns a Hilbert space which maps integers to and from the curve. +// n must be a power of two. +func NewHilbert64(n uint64) (*Hilbert64, error) { + if n == 0 { + return nil, ErrNotPositive + } + + // Test if power of two + if (n & (n - 1)) != 0 { + return nil, ErrNotPowerOfTwo + } + + return &Hilbert64{ + N: n, + }, nil +} + +// GetDimensions returns the width and height of the 2D space. +func (s *Hilbert64) GetDimensions() (uint64, uint64) { + return s.N, s.N +} + +// Map transforms a one dimension value, t, in the range [0, n^2-1] to coordinates on the Hilbert +// curve in the two-dimension space, where x and y are within [0,n-1]. +func (s *Hilbert64) Map(t uint64) (x, y uint64, err error) { + if t >= s.N*s.N { + return 0, 0, ErrOutOfRange + } + + for i := uint64(1); i < s.N; i = i * 2 { + rx := t&2 == 2 + ry := t&1 == 1 + if rx { + ry = !ry + } + + x, y = s.rotate(i, x, y, rx, ry) + + if rx { + x = x + i + } + if ry { + y = y + i + } + + t /= 4 + } + + return +} + +// MapInverse transform coordinates on Hilbert curve from (x,y) to t. +func (s *Hilbert64) MapInverse(x, y uint64) (t uint64, err error) { + if x >= s.N || y >= s.N { + return 0, ErrOutOfRange + } + + for i := s.N / 2; i > 0; i = i / 2 { + rx := (x & i) > 0 + ry := (y & i) > 0 + + a := uint64(0) + if rx { + a = 3 + } + t += i * i * (a ^ uint64(b2i(ry))) + + x, y = s.rotate(i, x, y, rx, ry) + } + + return +} + +// rotate rotates and flips the quadrant appropriately. +func (s *Hilbert64) rotate(n, x, y uint64, rx, ry bool) (uint64, uint64) { + if !ry { + if rx { + x = n - 1 - x + y = n - 1 - y + } + + x, y = y, x + } + return x, y +} diff --git a/hilbert64_test.go b/hilbert64_test.go new file mode 100644 index 0000000..24539dc --- /dev/null +++ b/hilbert64_test.go @@ -0,0 +1,235 @@ +// Copyright 2015 Google Inc. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +package hilbert + +import ( + "math/rand" + "testing" +) + +const benchmarkN64 = uint64(32) + +// Test cases below assume N=16 +var testCases64 = []struct { + d, x, y uint64 +}{ + {0, 0, 0}, + {16, 4, 0}, + {32, 4, 4}, + {48, 3, 7}, + {64, 0, 8}, + {80, 0, 12}, + {96, 4, 12}, + {112, 7, 11}, + {128, 8, 8}, + {144, 8, 12}, + {160, 12, 12}, + {170, 15, 15}, + {176, 15, 11}, + {192, 15, 7}, + {208, 11, 7}, + {224, 11, 3}, + {240, 12, 0}, + {255, 15, 0}, +} + +func TestNewErrors64(t *testing.T) { + var newTestCases = []struct { + n uint64 + wantErr error + }{ + {0, ErrNotPositive}, + {3, ErrNotPowerOfTwo}, + {5, ErrNotPowerOfTwo}, + } + + for _, tc := range newTestCases { + s, err := NewHilbert64(tc.n) + if s != nil || err != tc.wantErr { + t.Errorf("NewHilbert64(%d) did not fail, want %q, got (%+v, %q)", tc.n, tc.wantErr, s, err) + } + } +} + +func TestMapRangeErrors64(t *testing.T) { + var mapRangeTestCases = []struct { + d uint64 + wantErr error + }{ + {0, nil}, + {255, nil}, + {256, ErrOutOfRange}, + } + + s, err := NewHilbert64(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range mapRangeTestCases { + if _, _, err = s.Map(tc.d); err != tc.wantErr { + t.Errorf("Map(%d) did not fail, want %q, got %q", tc.d, tc.wantErr, err) + } + } +} + +func TestMapInverseRangeErrors64(t *testing.T) { + var mapInverseRangeTestCases = []struct { + x, y uint64 + wantErr error + }{ + {0, 0, nil}, + {15, 15, nil}, + {16, 0, ErrOutOfRange}, + {0, 16, ErrOutOfRange}, + } + + s, err := NewHilbert64(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range mapInverseRangeTestCases { + if _, err = s.MapInverse(tc.x, tc.y); err != tc.wantErr { + t.Errorf("MapInverse(%d, %d) did not fail, want %q, got %q", tc.x, tc.y, tc.wantErr, err) + } + } +} + +func TestSmallMap64(t *testing.T) { + s, err := NewHilbert64(1) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + x, y, err := s.Map(0) + if err != nil { + t.Errorf("Map(0) returned error: %s", err) + } + if x != 0 || y != 0 { + t.Errorf("Map(0) failed, want (0, 0), got (%d, %d)", x, y) + } + + d, err := s.MapInverse(0, 0) + if err != nil { + t.Errorf("MapInverse(0,0) returned error: %s", err) + } + if d != 0 { + t.Errorf("MapInverse(0, 0) failed, want 0, got %d", d) + } +} + +func TestMap64(t *testing.T) { + s, err := NewHilbert64(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range testCases64 { + x, y, err := s.Map(tc.d) + if err != nil { + t.Errorf("Map(%d) returned error: %s", tc.d, err) + } + if x != tc.x || y != tc.y { + t.Errorf("Map(%d) failed, want (%d, %d), got (%d, %d)", tc.d, tc.x, tc.y, x, y) + } + } +} + +func TestMapInverse64(t *testing.T) { + s, err := NewHilbert64(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range testCases64 { + d, err := s.MapInverse(tc.x, tc.y) + if err != nil { + t.Errorf("MapInverse(%d, %d) returned error: %s", tc.x, tc.y, err) + } + if d != tc.d { + t.Errorf("MapInverse(%d, %d) failed, want %d, got %d", tc.x, tc.y, tc.d, d) + } + } +} + +func TestAllMapValues64(t *testing.T) { + s, err := NewHilbert64(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for d := uint64(0); d < s.N*s.N; d++ { + // Map forwards and then back + x, y, err := s.Map(d) + if err != nil { + t.Errorf("Map(%d) returned error: %s", d, err) + } + if x < 0 || x >= s.N || y < 0 || y >= s.N { + t.Errorf("Map(%d) returned x,y out of range: (%d, %d)", d, x, y) + } + + dPrime, err := s.MapInverse(x, y) + if err != nil { + t.Errorf("MapInverse(%d, %d) returned error: %s", x, y, err) + } + if d != dPrime { + t.Errorf("Failed Map(%d) -> MapInverse(%d, %d) -> %d", d, x, y, dPrime) + } + } +} + +func BenchmarkMap64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := NewHilbert64(benchmarkN64) + if err != nil { + b.Fatalf("Failed to create hibert space: %s", err) + } + for d := uint64(0); d < benchmarkN64*benchmarkN64; d++ { + s.Map(d) + } + } +} + +func BenchmarkMapRandom64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := NewHilbert64(benchmarkN64) + if err != nil { + b.Fatalf("Failed to create hibert space: %s", err) + } + for d := uint64(0); d < benchmarkN64*benchmarkN64; d++ { + // Note that we're going to coerce this to an implied 32-bit type + // as required by Intn() but that the value is tiny and will always + // be within the range. + rd := rand.Intn(int(benchmarkN64 * benchmarkN64)) // Pick a random d + s.Map(uint64(rd)) + } + } +} + +func BenchmarkMapInverse64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := NewHilbert64(benchmarkN64) + if err != nil { + b.Fatalf("Failed to create hibert space: %s", err) + } + + for x := uint64(0); x < benchmarkN64; x++ { + for y := uint64(0); y < benchmarkN64; y++ { + s.MapInverse(x, y) + } + } + } +} diff --git a/peano64.go b/peano64.go new file mode 100644 index 0000000..29746f5 --- /dev/null +++ b/peano64.go @@ -0,0 +1,103 @@ +package hilbert + +// Peano represents a 2D Peano curve of order N for mapping to and from. +// Implements SpaceFilling interface. +type Peano64 struct { + N uint64 // Always a power of three, and is the width/height of the space. +} + +// NewPeano returns a new Peano space filling curve which maps integers to and from the curve. +// n must be a power of three. +func NewPeano64(n uint64) (*Peano64, error) { + if n == 0 { + return nil, ErrNotPositive + } + + if !isPow3(float64(n)) { + return nil, ErrNotPowerOfThree + } + + return &Peano64{ + N: n, + }, nil +} + +// GetDimensions returns the width and height of the 2D space. +func (p *Peano64) GetDimensions() (uint64, uint64) { + return p.N, p.N +} + +// Map transforms a one dimension value, t, in the range [0, n^3-1] to coordinates on the Peano +// curve in the two-dimension space, where x and y are within [0,n-1]. +func (p *Peano64) Map(t uint64) (x, y uint64, err error) { + if t >= p.N*p.N { + return 0, 0, ErrOutOfRange + } + + for i := uint64(1); i < p.N; i = i * 3 { + s := t % 9 + + // rx/ry are the coordinates in the 3x3 grid + rx := uint64(s / 3) + ry := uint64(s % 3) + if rx == 1 { + ry = 2 - ry + } + + // now based on depth rotate our points + if i > 1 { + x, y = p.rotate(i, x, y, s) + } + + x += rx * i + y += ry * i + + t /= 9 + } + + return x, y, nil +} + +// rotate rotates the x and y coordinates depending on the current n depth. +func (p *Peano64) rotate(n, x, y, s uint64) (uint64, uint64) { + + if n == 1 { + // Special case + return x, y + } + + n = n - 1 + switch s { + case 0: + return x, y // normal + case 1: + return n - x, y // fliph + case 2: + return x, y // normal + case 3: + return x, n - y // flipv + case 4: + return n - x, n - y // flipv and fliph + case 5: + return x, n - y // flipv + case 6: + return x, y // normal + case 7: + return n - x, y // fliph + case 8: + return x, y // normal + } + + panic("assertion failure: this line should never be reached") +} + +// MapInverse transform coordinates on the Peano curve from (x,y) to t. +// NOT IMPLEMENTED YET +func (p *Peano64) MapInverse(x, y uint64) (t uint64, err error) { + if x >= p.N || y >= p.N { + return 0, ErrOutOfRange + } + + panic("Not finished") + return 0, nil +} diff --git a/peano64_test.go b/peano64_test.go new file mode 100644 index 0000000..d7a1a06 --- /dev/null +++ b/peano64_test.go @@ -0,0 +1,261 @@ +// Copyright 2016 Google Inc. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +package hilbert + +import ( + "math/rand" + "testing" +) + +const peanoBenchmarkN64 = uint64(81) + +// Test cases below assume N=9 +var peanoTestCases64 = []struct { + d, x, y uint64 +}{ + {0, 0, 0}, + {1, 0, 1}, + {2, 0, 2}, + {3, 1, 2}, + {4, 1, 1}, + {5, 1, 0}, + {6, 2, 0}, + {7, 2, 1}, + {8, 2, 2}, + {9, 2, 3}, + // TODO Add more +} + +func TestPeanoNewErrors64(t *testing.T) { + var newTestCases = []struct { + n uint64 + want error + }{ + {0, ErrNotPositive}, + {2, ErrNotPowerOfThree}, + {4, ErrNotPowerOfThree}, + } + + for _, tc := range newTestCases { + s, err := NewPeano64(tc.n) + if s != nil || err != tc.want { + t.Errorf("NewPeano64(%d) = (%+v, %q) did not fail want (?, %q)", tc.n, s, err, tc.want) + } + } +} + +func TestPeanoMapRangeErrors64(t *testing.T) { + var mapRangeTestCases = []struct { + d uint64 + wantErr error + }{ + {0, nil}, + {80, nil}, + {81, ErrOutOfRange}, + } + + s, err := NewPeano64(9) + if err != nil { + t.Fatalf("NewPeano64(9) failed: %s", err) + } + + for _, tc := range mapRangeTestCases { + if _, _, err = s.Map(tc.d); err != tc.wantErr { + t.Errorf("Map(%d) = %q want %q", tc.d, tc.wantErr, err) + } + } +} + +/* +func TestPeanoMapInverseRangeErrors64(t *testing.T) { + var mapInverseRangeTestCases = []struct { + x, y uint64 + wantErr error + }{ + {0, 0, nil}, + {15, 15, nil}, + {16, 0, ErrOutOfRange}, + {0, 16, ErrOutOfRange}, + } + + s, err := New(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range mapInverseRangeTestCases { + if _, err = s.MapInverse(tc.x, tc.y); err != tc.wantErr { + t.Errorf("MapInverse(%d, %d) did not fail, want %q, got %q", tc.x, tc.y, tc.wantErr, err) + } + } +} +*/ + +func TestPeanoSmallMap64(t *testing.T) { + s, err := NewPeano64(1) + if err != nil { + t.Fatalf("NewPeano(1) failed: %s", err) + } + + x, y, err := s.Map(0) + if err != nil { + t.Errorf("Map(0) returned error: %s", err) + } + if x != 0 || y != 0 { + t.Errorf("Map(0) = (%d, %d) want (0, 0)", x, y) + } + + /* + // TODO Test when MapInverse is implemented + d, err := s.MapInverse(0, 0) + if err != nil { + t.Errorf("MapInverse(0,0) returned error: %s", err) + } + if d != 0 { + t.Errorf("MapInverse(0, 0) failed, want 0, got %d", d) + } + */ +} + +func TestPeanoMap64(t *testing.T) { + s, err := NewPeano64(9) + if err != nil { + t.Fatalf("NewPeano(9) failed: %s", err) + } + + for _, tc := range peanoTestCases64 { + x, y, err := s.Map(tc.d) + if err != nil { + t.Errorf("Map(%d) returned error: %s", tc.d, err) + } + if x != tc.x || y != tc.y { + t.Errorf("Map(%d) = (%d, %d) want (%d, %d)", tc.d, x, y, tc.x, tc.y) + } + } +} + +/* +func TestPeanoMapInverse64(t *testing.T) { + s, err := New(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for _, tc := range testCases64 { + d, err := s.MapInverse(tc.x, tc.y) + if err != nil { + t.Errorf("MapInverse(%d, %d) returned error: %s", tc.x, tc.y, err) + } + if d != tc.d { + t.Errorf("MapInverse(%d, %d) failed, want %d, got %d", tc.x, tc.y, tc.d, d) + } + } +} + +func TestPeanoAllMapValues64(t *testing.T) { + s, err := New(16) + if err != nil { + t.Fatalf("Failed to create hibert space: %s", err) + } + + for d := 0; d < s.N*s.N; d++ { + // Map forwards and then back + x, y, err := s.Map(d) + if err != nil { + t.Errorf("Map(%d) returned error: %s", d, err) + } + if x >= s.N || y >= s.N { + t.Errorf("Map(%d) returned x,y out of range: (%d, %d)", d, x, y) + } + + dPrime, err := s.MapInverse(x, y) + if err != nil { + t.Errorf("MapInverse(%d, %d) returned error: %s", x, y, err) + } + if d != dPrime { + t.Errorf("Failed Map(%d) -> MapInverse(%d, %d) -> %d", d, x, y, dPrime) + } + } +} +*/ +func BenchmarkPeanoMap64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := NewPeano64(peanoBenchmarkN) + if err != nil { + b.Fatalf("NewPeano64(%d) failed: %s", peanoBenchmarkN64, err) + } + for d := uint64(0); d < peanoBenchmarkN64*peanoBenchmarkN64; d++ { + s.Map(d) + } + } +} + +func BenchmarkPeanoMapRandom64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := NewPeano64(peanoBenchmarkN64) + if err != nil { + b.Fatalf("NewPeano64(%d) failed: %s", peanoBenchmarkN64, err) + } + for d := uint64(0); d < peanoBenchmarkN64*peanoBenchmarkN64; d++ { + rd := rand.Intn(int(peanoBenchmarkN64 * peanoBenchmarkN64)) // Pick a random d + s.Map(uint64(rd)) + } + } +} + +/* +func BenchmarkPeanoMapInverse64(b *testing.B) { + for i := 0; i < b.N; i++ { + s, err := New(benchmarkN64) + if err != nil { + b.Fatalf("Failed to create hibert space: %s", err) + } + + for x := 0; x < benchmarkN64; x++ { + for y := 0; y < benchmarkN64; y++ { + s.MapInverse(x, y) + } + } + } +} +*/ + +func TestIsPow364(t *testing.T) { + testCases := []struct { + in float64 + want bool + }{ + {-1, false}, + {0, false}, + {1, true}, + {2, false}, + {3, true}, + {3.1, false}, + {4, false}, + {5, false}, + {8.9999, false}, + {9, true}, + {9.00001, false}, + {27, true}, + {59049, true}, + } + + for _, tc := range testCases { + got := isPow3(tc.in) + if got != tc.want { + t.Errorf("isPow3(%f) = %t want %t", tc.in, got, tc.want) + } + } +} From 2ad04e2ae9f865cea7e044c391e77bfcecf9e9d7 Mon Sep 17 00:00:00 2001 From: Dustin Oprea Date: Sat, 28 Jan 2017 22:38:27 -0500 Subject: [PATCH 2/2] Added 64-bit comments to document. - Made example code runnable. --- README.md | 49 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 44fd7dc..9b7b526 100644 --- a/README.md +++ b/README.md @@ -20,20 +20,41 @@ go get github.com/google/hilbert Example: ```go -import "github.com/google/hilbert" - -// Create a Hilbert curve for mapping to and from a 16 by 16 space. -s, err := hilbert.NewHilbert(16) - -// Create a Peano curve for mapping to and from a 27 by 27 space. -//s, err := hilbert.NewPeano(27) - -// Now map one dimension numbers in the range [0, N*N-1], to an x,y -// coordinate on the curve where both x and y are in the range [0, N-1]. -x, y, err := s.Map(t) - -// Also map back from (x,y) to t. -t, err := s.MapInverse(x, y) +package main + +import ( + "github.com/dsoprea/hilbert" +) + +func main() { + // Create a Hilbert curve for mapping to and from a 16 by 16 space. + //s, err := hilbert.NewHilbert64(16) + s, err := hilbert.NewHilbert(16) + if err != nil { + panic(err) + } + + // Create a Peano curve for mapping to and from a 27 by 27 space. + //s, err := hilbert.NewPeano64(27) + //s, err := hilbert.NewPeano(27) + + t := 112 + + // Now map one dimension numbers in the range [0, N*N-1], to an x,y + // coordinate on the curve where both x and y are in the range [0, N-1]. + x, y, err := s.Map(t) + if err != nil { + panic(err) + } + + // (x, y) <= (7, 11) + + // Also map back from (x,y) to t. + t, err = s.MapInverse(x, y) + if err != nil { + panic(err) + } +} ``` ## Demo