媒体数据管理上机实验

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

编程语言为 Python 3

1. LZW 编码算法

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

算法流程：

对输入字符串初始化生成单字符的字典，设P为空字符串；
遍历字符串每一位C，判断P+C是否在字典中，若是则将P+C赋予P，否则在字典中添加P+C，在结果添加P对应的值，并将C赋予P；
最后一位字符在结果添加P的值；
输出编码结果与字典；
解码程序拿到编码结果与初始化的单字符字典；
将字典键值翻转，对编码的每一位求得对应字符前缀，保存到字典；
除第一次，将前缀添加到上一条的末尾，并更新前缀，记录；
输出记录下的字符。

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))

测试：

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 变换

算法流程：

将输入数据预处理（中心化）；
求得协方差矩阵C；
求C的特征值和特征向量；
取最大特征值的特征向量，求得最佳投影矩阵P；
对原始样本矩阵投影，得到降维后的新样本矩阵。

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 实现对图片的矢量量化

算法流程

读取输入图片并转化为 (height * width, channel)
计算 kmeans
矢量量化
输出图像

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；

算法流程：

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

自己写的代码：

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图像集；

输出：最相似的图片；

算法流程：

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

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

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)

学习相关

#笔记 #python

BlogUPUP 折腾记录上一篇

What happens when （二）TCP 下一篇