媒体数据管理上机实验

媒体数据管理上机实验

  • 2020-07-07
  • 次

这门课的上机题目还都挺有意思，尤其是第三个 LSH 索引耗费了不少精力，故记录下来。

编程语言为 Python 3

1. LZW 编码算法.

如果字典数超过 10 可能会出 bug，因为编码是 1 位的（可通过增加编码位数解决）

算法流程：

  1. 对输入字符串初始化生成单字符的字典，设P为空字符串；
  2. 遍历字符串每一位C，判断P+C是否在字典中，若是则将P+C赋予P，否则在字典中添加P+C，在结果添加P对应的值，并将C赋予P；
  3. 最后一位字符在结果添加P的值；
  4. 输出编码结果与字典；
  5. 解码程序拿到编码结果与初始化的单字符字典；
  6. 将字典键值翻转，对编码的每一位求得对应字符前缀，保存到字典；
  7. 除第一次，将前缀添加到上一条的末尾，并更新前缀，记录；
  8. 输出记录下的字符。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def init(s):
    # s: 输入字符串
    # 返回: 包含所有单字符的字典
    table = {}
    tmp = []
    for i in s:
        if i not in tmp:
            tmp.append(i)
    j = 1
    for i in range(len(tmp)):
        t = str(j)
        table[tmp[i]] = t
        j = j + 1
    return table

def encoding(s):
    # s: 输入字符串
    # r: 编码结果 table: 编码表，初始化是单字符字典
    r = []
    table = init(s)
    P = ''
    for i in range(len(s)):
        C = s[i]
        t = P+C
        # 判断 P+C 是否在词典中
        if t in table.keys():
            P = t
        if t not in table.keys():
            table[t] = str(len(table)+1)
            r.append(table[P])
            P = C
        # 最后一位用单字符编码表示
        if i == len(s)-1:
            r.append(table[P])
    print("After compress: " + ''.join(r))
    print("coding table is: ")
    print(table)

    # 将单字符字典与编码结果传入 decoding 解码
    default = init(s)
    decoding(default,r)

def decoding(table,r):
    # table:包含单独字符的字典, r: LZW编码
    inverse_table = dict(zip(table.values(),table.keys()))
    num = len(table.keys())
    tablelen = num
    tem = []
    for i in range(len(r)):
        # P 前缀
        P = inverse_table[str(r[i])]
        inverse_table[str(num+1)] = P
        # 除第一个编码外，前缀第一个字符添加到上一条的末尾
        if num > tablelen:
            inverse_table[str(num)] = inverse_table[str(num)]+P[0]
        # 重新获取前缀，添加至 tem
        P = inverse_table[str(r[i])]
        tem.append(P)
        num = num+1
    print("decoding:"+''.join(tem))

测试：

1
2
3
4
5
input a string: abababcabcab
After compress: 1244374
coding table is: 
{'a': '1', 'b': '2', 'c': '3', 'ab': '4', 'ba': '5', 'aba': '6', 'abc': '7', 'ca': '8', 'abca': '9'}
decoding:abababcabcab

2. k-l 变换与矢量量化.

在 ColorHistogram.asc 数据集上实现 k-l 变换

算法流程：

  1. 将输入数据预处理（中心化）；
  2. 求得协方差矩阵C；
  3. 求C的特征值和特征向量；
  4. 取最大特征值的特征向量，求得最佳投影矩阵P；
  5. 对原始样本矩阵投影，得到降维后的新样本矩阵。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np

def kl(data):
	# 中心化
	mu = np.mean(data, axis=0)
	data = data - mu
	C = 1/data.shape[0]*np.dot(np.transpose(data),data)

	eig_val, eig_vec = np.linalg.eig(C)
	eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(len(eig_val))]
	eig_pairs.sort(reverse=True)

	P = eig_pairs[0][1]
	Y = np.transpose(np.dot(P,np.transpose(data)))
	print(Y)

data = np.loadtxt('ColorHistogram.asc', usecols = range(1, 33), unpack= False)
kl(data)

通过 scipy.cluster.vq 实现对图片的矢量量化

算法流程

  1. 读取输入图片并转化为 (height * width, channel)
  2. 计算 kmeans
  3. 矢量量化
  4. 输出图像
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from scipy.cluster.vq import kmeans, vq
from numpy import reshape
import cv2
# 矢量量化的三个级别
vqclst = [2, 10, 256]

data = cv2.imread('cs.jpg').astype(float)
(height, width, channel) = data.shape

data = reshape(data, (height * width, channel))
for k in vqclst:
print('Generating vq-%d...' % k)
# 计算kmeans
(centroids, distor) = kmeans(data, k)
# 矢量量化
    (code, distor) = vq(data, centroids)
    print('distor: %.6f' % distor.sum())
    im_vq = centroids[code, :]
cv2.imwrite('result-%d.jpg' % k, reshape(im_vq, (height, width, channel)))

3. LSH索引实现近似最近邻搜索.

查找 corel数据集 前 1000 行每个点的最近邻点，使用原始汉明编码，计算汉明距离，使用 LSH 索引，实现近似最近邻搜索。输入：corel数据集，汉明码长度K，建表数量L，哈希桶最大容量B；

算法流程：

  1. 对数据集预处理：乘100取整，获得最大值C与数据维数n；
  2. 在n*C的范围内生成k个随机数，存入表H，为要取到的汉明码位数。新建哈希索引表；
  3. 将数据点的坐标转化为0与1的汉明编码，取2中生成的随机数位，组成字符串pi；
  4. 若哈希索引表pi桶的数量大于B则不记录，否则在pi桶记录该数据点；
  5. 重复操作2，3，4，生成L张H表与哈希索引表；
  6. 对前1000行的点进行最邻近搜索，即对数据点求它的汉明编码，对每张H表和每张哈希索引表，根据H表的位数生成字符串，将哈希索引表的该字符串编码处的点作为候选项，将所有的候选项排序选择频次最高的前十个作为结果；
  7. 计算召回率 Recall、准确率 Accuracy、所用时间 timecost

结果

自己写的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
import numpy as np
import random
import heapq
import time

# 返回点的坐标的汉明编码
def Embedding(C,line):
	# C:坐标最大值，line:点的坐标
	string = []
	for element in line:
		for i in range(element):
			string.append("1")
		for j in range(C-element):
			string.append("0")
	return  string

# 生成哈希索引表
def set_table(H,C,B,data):
	# H：生成的k个随机数，C：取到的最大值，B：桶的最大数量
	# data：处理过的数据集
	table = {}
	f = {}
	for i, line in enumerate(data):
		pi = []
		# var是汉明编码
		var = ''.join(Embedding(C,line))
		# 取汉明编码的H对应位
		for j in H:
			pi.append(var[j])
		# f字典用于记录每个哈希桶包含的行数，小于B则+1
		if (''.join(pi)) in f:
			if f[''.join(pi)]<B:
				f[''.join(pi)] += 1
				table[i] = ''.join(pi)
		else:
			f[''.join(pi)]=1
			table[i] = ''.join(pi)
	return table


def set_H(Cn,k):
	H = []
	for i in range(k):
		H.append(random.randint(0,Cn-1))
	H.sort()
	return H

# 求line的LSH最邻近点
def search(Hlist,tblist,var):
	# var：一行的汉明编码
	# 返回该行的所有匹配点（有重复项）
	keys = []
	# 用Hlist的每个H求对应位组合
	for i,h in enumerate(Hlist):
		p = []
		for hi in h:
			p.append(var[hi])
		# 寻找所以哈希相同的点
		for key, value in tblist[i].items():
			if value== ''.join(p):
				keys.append(key)
	return keys


def LSH(line, Hlist, tblist):
	array = line
	var = ''.join(Embedding(C, array))
	# keys是求得的最邻近点
	keys = search(Hlist, tblist, var)
	keys.sort()
	# 按照出现的次数从大到小排序，去除重复项
	set1 = set(keys)
	dict01 = {item: keys.count(item) for item in set1}
	sorted_x = sorted(dict01.items(), key=lambda x: x[1], reverse=True)
	keys = []
	num = 0
	for i,n in sorted_x:
		keys.append(i)
		num += 1
	return keys


def distance(line1,line2):
	dist = np.linalg.norm(line1-line2)
	return dist

#暴力搜索获得最近邻
def mindistancedata(linenum,data):
	dist = []
	for i, line in enumerate(data):
		dist.append(distance(data[linenum],data[i]))
	min_index_list = map(dist.index, heapq.nsmallest(10, dist))
	return list(min_index_list)

#LSH搜索获得近似最近邻
def mindistance(linenum,data,linelist):
	min_index_list = []
	dist = {}
	for i in linelist:
		dist[i] = distance(data[linenum], data[i])
	L = sorted(dist.items(), key=lambda item: item[1])
	for i in range(10):
		min_index_list.append(L[i][0])
	return min_index_list

# 数据预处理，乘上10的平方取整
bit = 2
odata = data = np.loadtxt('ColorHistogram.asc', usecols = range(1, 33), unpack= False)
data = np.loadtxt('ColorHistogram.asc',usecols = range(1, 33), unpack= False)
data = data*(10**bit)
data = data.astype(np.int)

k = input("input K:")
k = int(k)

L = input("input L:")
L = int(L)

B = input("input B:")
B = int(B)

C = int(np.max(data)+1)
n = data.shape[1]
hamming_code = []
Cn = C*n
Hlist=[]
tblist = []

# 创建索引
for i in range(L):
	H = set_H(Cn,k)
	Hlist.append(H)
	tblist.append(set_table(H,C,B,data))
	print("创建索引哈希表： "+str(i))


while(True):
	true = []
	my = []
	linenum = 1000
	
	# 暴力搜索，存入true
	for i in range(linenum):
		trueindex = []
		trueindex = mindistancedata(i, odata)
		true.append(trueindex)
	time_start = time.time()
	tinm = 0
	mint = 0
	
	# LSH索引搜索，记录正确匹配的点数
	for i in range(linenum):
		myindex = []
		line  = data[i]
		keys = LSH(line, Hlist, tblist)
		myindex = mindistance(i,odata,keys)
		for ni in range(10):
			if true[i][ni] in myindex:
				tinm += 1
	time_end = time.time()		
	# 输出结果
	print("K=" + str(k) +",L=" + str(L) + "B=" + str(B))
	print("召回率"+str(tinm/linenum/10))
	print("准确率"+str(1-(linenum*10-tinm)*2/linenum/68040))
	print('time cost', time_end - time_start, 's')

4. SIFT特征提取与图像匹配.

算法描述：

输入：原始图像、train图像集；

输出：最相似的图片；

算法流程：

  1. 使用SIFT提取原图的特征，得到关键点 kp 及描述符 des；
  2. 获取train中图片，使用SIFT提取特征，得到关键点 kp1 及描述符 des1；
  3. 初始化flann匹配器；
  4. knnMatch 返回原图每个特征点的两个匹配；
  5. 若第一个匹配点的距离小于第二个的0.7倍，则该对特征匹配成功，计入匹配成功列表；
  6. 若匹配成功的数量大于最大匹配成功数量 ，则替换；
  7. 重复2到6的操作，直到train中图片搜索完毕。
  8. 输出最大匹配的图片名。

注意：cv2.xfeatures2d.SIFT_create() 在新版 opencv-python中删除，我的解决方案是通过 pip uninstall opencv-python 卸载原有的 opencv-python ，再通过 pip install opencv-contrib-python==3.4.2.16 安装旧版本的 opencv-contrib-python，可以运行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import cv2
import os
import time
# 读取第一张图片，计算关键点 kp 及描述符 des
img = cv2.imread('./test/0/0_3.jpg', cv2.IMREAD_COLOR)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
cv2.imshow('origin', img);
cv2.waitKey(20)
img = gray
sift = cv2.xfeatures2d.SIFT_create()
kp, des = sift.detectAndCompute(img,None)
matchnum = 0
matchid = ""

# 遍历./train/文件夹所有的.jpg文件
path = './train/'
time_start = time.time()
for i in range(10):
	dirpath = path + str(i) + '/'
	imagelist = os.listdir(dirpath)
	for imgname in imagelist:
		if (imgname.endswith(".jpg")):
			image = cv2.imread(dirpath + imgname)
			gray1 = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
			image = gray1
			kp1, des1 = sift.detectAndCompute(image, None)
			if len(kp1)<=1:
				continue

			# 生成 flann 匹配器
			FLANN_INDEX_KDTREE = 0
			indexParams = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
			searchParams = dict(checks=50)
			flann = cv2.FlannBasedMatcher(indexParams, searchParams)

			# 匹配，让 knnMatch 返回原图每个特征点的两个匹配点
			# 若第一个匹配点的距离小于第二个的0.7倍，则该对特征匹配成功，计入 good 列表
			matches = flann.knnMatch(des, des1, k=2)
			good = []
			for m, n in matches:
				if m.distance < 0.7 * n.distance:
					good.append([m])

			# 匹配成功的数量大于 matchnum ，则替换
			if len(good) > matchnum:
				matchnum = len(good)
				matchid = imgname
			print(imgname)
time_end = time.time()
print('time cost', time_end - time_start, 's')
print("flann 关键点匹配个数:"+str(matchnum))
print("最相似图片文件名:"+matchid)
感谢您的阅读，本文由 夏月 版权所有。如若转载，请注明出处：夏月（http://www.eteluna.top/2020/07/07/Mediadataexp/）
What happens when （二）TCP
BlogUPUP 折腾记录