构造方向可控金字塔1. 方向可控金字塔简介
Willian T. Freeman等1991年提出方向可控金字塔。和高斯金字塔或拉普拉斯金字塔相比，方向金字塔除了分解尺度子带以外，还会分解方向子带。为什么需要用方向子带呢？因为在图像匹配等应用场合中，尺度变化外加旋转操作计算量很大。微分计算获得方向和多尺度的金字塔同时具有线性和平移不变性，所以考虑将尺度和方向结合起来操作。
2. 拉普拉斯金字塔，双尺度小波变换和方向可控金字塔比较（1）紧框架
紧框架是指变换矩阵的逆矩阵等于变换矩阵的转置。在构造图像金字塔时需用不同的带通滤波器对图像作分解，而这些带通滤波器以矩阵的形式出现，而在图像重构需用不同的带通滤波器的逆矩阵对图像重建。因为矩阵求逆比较耗时，所以紧框架的性质会大大加快图像重建的计算速度。拉普拉斯金字塔不具有紧框架性质；双尺度小波变换和尺度金字塔具有紧框架性质。
（2）过完备
拉普拉斯金字塔过完备性为4/3，双尺度小波变换过完备性为1，方向可控金字塔过完备性为4k/3。过完备是指信息的冗余量。方向金字塔包含方向信息时，描述同一幅图像占用的资源最多。
（3）混叠
拉普拉斯金字塔有可能发生混叠；双尺度小波变换一定会混叠；尺度金字塔不会混叠。
（4）旋转方向子带
拉普拉斯金字塔没有方向子带；小波变换在六角晶格上有方向子带；方向可控金字塔有方向子带。
综上，方向可控金字塔是占用资源和方向个数成正比的，重建速度快，不会混叠的，同时有尺度和方向子带的图像表示。（P.S. 双尺度小波变换还没用过，支线加上不懂先跳过）
3. 构造方向可控金字塔滤波器
用方向可控金字塔分解图像最重要的是构造方向子带的滤波器。不过，OpenCV没有给我们提供这样的库函数，所以只能自己写构造过程了。幸运的是，K. R. Castleman等的文章“方向可控滤波器的简单设计”给我们提供了滤波器的构造方法。根据它笔者实现了方向可控金字塔滤波器的构造。
（1）实现代码
import cv2import csvimport mathimport numpy as npfrom matplotlib import cmfrom matplotlib import pyplot as pltfrom mpl_toolkits.mplot3d.axes3d import Axes3D################################################################################print 'parameter set'f1 = 0f2 = math.pi * 5 / 8fn = math.pi################################################################################print 'range of u and v'umin, umax, du = -fn, fn, 0.25vmin, vmax, dv = -fn, fn, 0.25ulist = np.arange(umin, umax, du)vlist = np.arange(vmin, vmax, dv)'''ulist = []vlist = []for i in range(int((umax – umin) / du)):ulist.append(umin + du * i)for j in range(int((vmax – vmin) / dv)):vlist.append(vmin + dv * j)'''################################################################################# lowpass filter functiondef lowpass(f, a, b):if f < a or f == a:s = 1elif f > a and f < b:s = math.sqrt(0.5 * (1 + math.cos(math.pi * (f – a) / (b – a))))else:s = 0return s# highpass filter functiondef highpass(f, a, b):if f < a or f == a:s = 0elif f > a and f < b:s = math.sqrt(0.5 * (1 – math.cos(math.pi * (f – a) / (b – a))))else:s = 1return s# build lowpass or highpass filterdef passFilter(ulist, vlist, a, b, case):P = []for u in ulist:U = []for v in vlist:f = math.sqrt(u * u + v * v)if case.lower() == 'l':s = lowpass(f, a, b)elif case.lower() == 'h':s = highpass(f, a, b)U.append(s)P.append(U)return P################################################################################print 'Lowpass Filter L0(u,v)'LP0 = passFilter(ulist, vlist, f2, math.pi, 'l')print 'Lowpass Filter L1(u,v)'LP1 = passFilter(ulist, vlist, f1, fn / 2, 'l')LP = []LP.append(LP0)LP.append(LP1)################################################################################print 'Highpass Filter H0(u,v)'HP0 = passFilter(ulist, vlist, f2, fn, 'h')print 'Highpass Filter H1(u,v)'HP1 = passFilter(ulist, vlist, f1, fn / 2, 'h')HP = []HP.append(HP0)HP.append(HP1)################################################################################print 'Steerable Bandpass Filter'# k = 2, sm(m = 0,…,k – 1)def steerableBandpass(ulist, vlist, HP):B0 = []B1 = []for u in range(len(ulist)):U0 = []U1 = []for v in range(len(vlist)):s0 = HP[u][v] * math.cos(math.atan2(ulist[v], ulist[u]))s1 = HP[u][v] * math.cos(math.atan2(ulist[v], ulist[u]) + math.pi / 2)U0.append(s0)U1.append(s1)B0.append(U0)B1.append(U1)return B0, B1B0, B1 = steerableBandpass(ulist, vlist, HP1)################################################################################print 'display filters'fig = plt.figure()ax = fig.add_subplot(231, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V, LP0, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False) ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('L0')ax = fig.add_subplot(232, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V, LP1, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('L1')ax = fig.add_subplot(233, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V,HP0, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False) ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('H0')ax = fig.add_subplot(234, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V, HP1, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('H1')ax = fig.add_subplot(235, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V, B0, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False) ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('B1')ax = fig.add_subplot(236, projection = '3d')U, V = np.meshgrid(ulist, vlist) surf = ax.plot_surface(U, V, B1, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)ax.set_xlabel('U')ax.set_ylabel('V')ax.set_zlabel('B2')plt.show() ################################################################################print 'Goodbye!'（2）实验结果
而更像是听见了天地间冥冥中的呼唤，
