dct: Improve DCT1D to O(nlogn) algorithm (#29)

* hashcompute_test: Add benchmark

* dct: Improve DCT1D to O(nlogn) algorithm

AS-IS:
BenchmarkPerceptionHash-8  500  2893930 ns/op  456698 B/op  4455 allocs/op

TO-BE:
BenchmarkPerceptionHash-8  2000  890306 ns/op  456382 B/op  4455 allocs/op

reference: DCT type II, unscaled. Algorithm by Byeong Gi Lee, 1984.

_examples/examples.go

@@ -29,6 +29,8 @@ func main() {
 	hash2, _ = goimagehash.PerceptionHash(img2)
 	distance, _ = hash1.Distance(hash2)
 	fmt.Printf("Distance between images: %v\n", distance)
+	fmt.Println(hash1.ToString())
+	fmt.Println(hash2.ToString())
 	fmt.Println(hash1.Bits())
 	fmt.Println(hash2.Bits())
 }

_examples/load_and_dump.go

@@ -19,9 +19,9 @@ func main() {
 	foo := bufio.NewWriter(&b)
 	img1, _ := jpeg.Decode(file1)
 	img2, _ := jpeg.Decode(file2)
-	width, height := 15, 15
+	width, height := 16, 16
-	hash1, _ := goimagehash.ExtAverageHash(img1, width, height)
+	hash1, _ := goimagehash.ExtPerceptionHash(img1, width, height)
-	hash2, _ := goimagehash.ExtAverageHash(img2, width, height)
+	hash2, _ := goimagehash.ExtPerceptionHash(img2, width, height)
 	hash1024, _ := goimagehash.ExtAverageHash(img2, 32, 32)
 	distance, _ := hash1.Distance(hash2)
 	fmt.Printf("Distance between images: %v\n", distance)

hashcompute.go

@@ -87,12 +87,16 @@ func PerceptionHash(img image.Image) (*ImageHash, error) {
 // ExtPerceptionHash function returns phash of which the size can be set larger than uint64
 // Some variable name refer to https://github.com/JohannesBuchner/imagehash/blob/master/imagehash/__init__.py
 // Support 64bits phash (width=8, height=8) and 256bits phash (width=16, height=16)
+// Important: width * height should be the power of 2
 func ExtPerceptionHash(img image.Image, width, height int) (*ExtImageHash, error) {
+	imgSize := width * height
 	if img == nil {
 		return nil, errors.New("Image object can not be nil")
 	}
+	if imgSize <= 0 || imgSize&(imgSize-1) != 0 {
+		return nil, errors.New("width * height should be power of 2")
+	}
 	var phash []uint64
-	imgSize := width * height
 	resized := resize.Resize(uint(imgSize), uint(imgSize), img, resize.Bilinear)
 	pixels := transforms.Rgb2Gray(resized)
 	dct := transforms.DCT2D(pixels, imgSize, imgSize)

hashcompute_test.go

@@ -192,10 +192,6 @@ func TestExtImageHashCompute(t *testing.T) {
 		{"_examples/sample1.jpg", "_examples/sample4.jpg", 16, 16, ExtPerceptionHash, "ExtPerceptionHash", 122},
 		{"_examples/sample2.jpg", "_examples/sample3.jpg", 16, 16, ExtPerceptionHash, "ExtPerceptionHash", 118},
 		{"_examples/sample2.jpg", "_examples/sample4.jpg", 16, 16, ExtPerceptionHash, "ExtPerceptionHash", 104},
-		{"_examples/sample1.jpg", "_examples/sample1.jpg", 17, 17, ExtPerceptionHash, "ExtPerceptionHash", 0},
-		{"_examples/sample2.jpg", "_examples/sample2.jpg", 17, 17, ExtPerceptionHash, "ExtPerceptionHash", 0},
-		{"_examples/sample3.jpg", "_examples/sample3.jpg", 17, 17, ExtPerceptionHash, "ExtPerceptionHash", 0},
-		{"_examples/sample4.jpg", "_examples/sample4.jpg", 17, 17, ExtPerceptionHash, "ExtPerceptionHash", 0},
 		{"_examples/sample1.jpg", "_examples/sample1.jpg", 8, 8, ExtDifferenceHash, "ExtDifferenceHash", 0},
 		{"_examples/sample2.jpg", "_examples/sample2.jpg", 8, 8, ExtDifferenceHash, "ExtDifferenceHash", 0},
 		{"_examples/sample3.jpg", "_examples/sample3.jpg", 8, 8, ExtDifferenceHash, "ExtDifferenceHash", 0},
@@ -281,3 +277,21 @@ func BenchmarkExtImageHashDistanceDifferent(b *testing.B) {
 		}
 	}
 }
+
+func BenchmarkPerceptionHash(b *testing.B) {
+	file1, err := os.Open("_examples/sample3.jpg")
+	if err != nil {
+		b.Errorf("%s", err)
+	}
+	defer file1.Close()
+	img1, err := jpeg.Decode(file1)
+	if err != nil {
+		b.Errorf("%s", err)
+	}
+	for i := 0; i < b.N; i++ {
+		_, err := ExtPerceptionHash(img1, 8, 8)
+		if err != nil {
+			b.Errorf("%s", err)
+		}
+	}
+}

transforms/dct.go

@@ -10,25 +10,31 @@ import (
 )
 
 // DCT1D function returns result of DCT-II.
-// Follows Matlab dct().
+// DCT type II, unscaled. Algorithm by Byeong Gi Lee, 1984.
-// Implementation reference:
-// https://unix4lyfe.org/dct-1d/
 func DCT1D(input []float64) []float64 {
-	n := len(input)
+	temp := make([]float64, len(input))
-	out := make([]float64, n)
+	forwardTransform(input, temp, len(input))
-	for i := 0; i < n; i++ {
+	return input
-		z := 0.0
+}
-		for j := 0; j < n; j++ {
-			z += input[j] * math.Cos(math.Pi*(float64(j)+0.5)*float64(i)/float64(n))
-		}
-
-		if i == 0 {
-			z *= math.Sqrt(0.5)
-		}
-		out[i] = z * math.Sqrt(2.0/float64(n))
 
+func forwardTransform(input, temp []float64, Len int) {
+	if Len == 1 {
+		return
 	}
-	return out
+
+	halfLen := Len / 2
+	for i := 0; i < halfLen; i++ {
+		x, y := input[i], input[Len-1-i]
+		temp[i] = x + y
+		temp[i+halfLen] = (x - y) / (math.Cos((float64(i)+0.5)*math.Pi/float64(Len)) * 2)
+	}
+	forwardTransform(temp, input, halfLen)
+	forwardTransform(temp[halfLen:], input, halfLen)
+	for i := 0; i < halfLen-1; i++ {
+		input[i*2+0] = temp[i]
+		input[i*2+1] = temp[i+halfLen] + temp[i+halfLen+1]
+	}
+	input[Len-2], input[Len-1] = temp[halfLen-1], temp[Len-1]
 }
 
 // DCT2D function returns a  result of DCT2D by using the seperable property.

transforms/dct_test.go

@@ -17,7 +17,7 @@ func TestDCT1D(t *testing.T) {
 		input  []float64
 		output []float64
 	}{
-		{[]float64{1.0, 1.0, 1.0, 1.0}, []float64{2.0, 0, 0, 0}},
+		{[]float64{1.0, 1.0, 1.0, 1.0}, []float64{4.0, 0, 0, 0}},
 	} {
 
 		out := DCT1D(tt.input)
@@ -50,10 +50,10 @@ func TestDCT2D(t *testing.T) {
 			{5.0, 6.0, 7.0, 8.0},
 			{9.0, 10.0, 11.0, 12.0},
 			{13.0, 14.0, 15.0, 16.0}},
-			[][]float64{{34.0, -4.46088499, 0.0, -0.31702534},
+			[][]float64{{136.0, -12.6172881195958, 0.0, -0.8966830583359305},
-				{-17.84353998, 0.0, 0.0, 0.0},
+				{-50.4691524783832, 0.0, 0.0, 0.0},
 				{0.0, 0.0, 0.0, 0.0},
-				{-1.26810134, 0.0, 0.0, 0.0}},
+				{-3.586732233343722, 0.0, 0.0, 0.0}},
 			4, 4},
 	} {
 		out := DCT2D(tt.input, tt.w, tt.h)