兹于2017年12月，应《多媒体技术基础》课程实验的要求，本人就基于二维哈尔小波变换算法做了较为深入的理解，用MATLAB和C++语言实现二维哈尔小波变换算法。


二维哈尔小波变换算法



一、实验目的
掌握小波及小波变换相关算法。
 
二、实验设备与环境
Windows 7，MATLAB，DEV C++
 
三、实验内容、程序清单及运行结果
实验要求
根据课本内容，实现图像小波变换算法。
具体可参考《多媒体技术基础（第三版）》，可以在这里下载PDF电子书

小波与小波变换算法概述
小波是定义在有限间隔而且其平均值为零的一种函数。
小波变换算法求有限信号的均值和差值：
例如：假设有一幅分辨率只有4个像素p0、p1、p2、p3 的一维图像，对应的像素值或者叫做图像位置的系数分别为：
[ 9   7     3     5 ]
计算它的哈尔小波变换系数。
计算步骤如下：
步骤1：求均值(averaging)。计算相邻像素对的平均值，得到一幅分辨率比较低的新图像，它的像素数目变成了2个，即新的图像的分辨率是原来的1/2，相应的像素值为：
[ 8       4 ]
步骤2：求差值(differencing)。很明显，用2个像系表示这幅图像时，图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成的原始图像，就需变存储一些图像的细节系数(detail coefficient)，以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像索对的平均值.或者使用这个像素对的差值除以2。在这个例子中，第一个细节系数是(9-8)=1，因为计算得到的平均值是8，它比9小1面比7大1，存储这个细节系数就可以恢复原始图像的前两个像素值。使用同样的方法，第二个细节系数是(3-4)=-1，存储这个细节系数就可以恢复后2个像素值。因此，原始图像就可以用下面的两个平均值和两个细节系数表示：
[ 8       4     1     -1 ]
步骤3：重复步骤1和2，把由第一步分解得到的图像进一步分解成分辨率更低的图像和细节系数。在这个例子中，分解到最后，就用一个像素的平均值6和三个细节系数2、1和-1表示整幅图像：
[ 6       2     1     -1 ]
这个分解过程见下表：
    分辨率    
平均值
     细节系数     
4
      [ 9   7     3     5 ]      
 
2
[ 8   4 ]
[ 1   -1 ]
1
[ 6 ]
[ 2 ]
由此可见，通过上述分解就把由4个像素组成的一幅图像用一个平均像素值和3个细节系数表示，这个过程就叫做哈尔小波变换( Haar wavelet transform)，也称哈尔小波分解(Haarwavelet decomposition)。这个概念可以推广到使用其他小波基的变换。
从这个例子中我们可以看到：
（1）变换过程中没有丢失信息，因为能够从所记录的数据中重构出原始图像。
（2）对这个给定的变换，我们可以从所记录的数据中重构出各种分辨率的图像。例如，在分辨率为l的图像基础上重构出分辨率为2的图像，在分辨率为2的图像基础上重构出分辨率为4的图像。
（3）通过变换之后产生的细节系数的幅度值比较小，这就为图像压缩提供了一种途径。例如，去掉一些微不足道的细节系数并不影响对重构图像的理解。

小波与小波变换算法步骤及核心程序
二维哈尔小波变换是将一个由8*8像素组成的图像块的矩阵，进行小波变换时，对矩阵中的每一行进行行变换，然后对行变换后的矩阵每一列进行列变换，最后将变换后的矩阵进行编码。
其次可通过设置阈值δ，例如δ≤5，表示把绝对值小于5的细节系数当作“0”看待，从而对编码后的矩阵进行压缩。
最后通过逆变换重构源图像矩阵的近似矩阵。
 
MATLAB代码：
% % 使用线性代数方法计算哈尔小波变换
disp '设第1次相乘的矩阵为M1，第2次为M2，第3次为M3。W=M1*M2*M3：'
W = [1/8,1/8,1/4,0,1/2,0,0,0
	1/8,1/8,1/4,0,-1/2,0,0,0
    1/8,1/8,-1/4,0,0,1/2,0,0
    1/8,1/8,-1/4,0,0,-1/2,0,0
    1/8,-1/8,0,1/4,0,0,1/2,0
    1/8,-1/8,0,1/4,0,0,-1/2,0
    1/8,-1/8,0,-1/4,0,0,0,1/2
    1/8,-1/8,0,-1/4,0,0,0,-1/2]


disp '源图像的一个图像块矩阵：'
A = [64,2,3,61,60,6,7,57
	 9,55,54,12,13,51,50,16
     17,47,46,20,21,43,42,24
     40,26,27,37,36,30,31,33
     32,34,35,29,28,38,39,25
     41,23,22,44,45,19,18,48
     49,15,14,52,53,11,10,56
     8,58,59,5,4,62,63,1]
 
disp '求均值差值后的矩阵：'
 T = W'*A*W
 
% %  设阈值为5，绝对值小于等于5的元素设为0
 for i=1:8
    for j=1:8
       if(abs(T(i,j))<=5)
           T(i,j)=0;
       end
    end
 end
 
disp '设置阈值变换后的矩阵：'
 T


disp '逆变换重构后的近似矩阵：'
 LA = ((W')^(-1))*T*W^(-1)
运行结果：






C++代码：
#include<iostream>
#include<cmath>
#define N 8 //把图像划分成8*8个图像块 
using namespace std;
/*使用线性代数方法计算哈尔小波变换*/

/*矩阵转置函数*/
bool  MatrixTranspose(const double source[N][N],double transpose[N][N]){
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			transpose[i][j] = source[j][i];
		}
	}
	return true;
}

/*矩阵乘法*/
bool MatrixMultiplication(const double operand1[N][N], const double operand2[N][N], double result[N][N]){
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			double sum=0.0;
			for(int k=0;k<N;k++){
				sum = sum + operand1[i][k]*operand2[k][j];
			}
			result[i][j] = sum;
		}
	}
	return true;
} 

int main(){
	/*设第1次相乘的矩阵为M1，第2次为M2，第3次为M3。W=M1*M2*M3*/
	double W[N][N] = {1.0/8,1.0/8,1.0/4,0,1.0/2,0,0,0,
					  1.0/8,1.0/8,1.0/4,0,-1.0/2,0,0,0,
					  1.0/8,1.0/8,-1.0/4,0,0,1.0/2,0,0,
					  1.0/8,1.0/8,-1.0/4,0,0,-1.0/2,0,0,
					  1.0/8,-1.0/8,0,1.0/4,0,0,1.0/2,0,
					  1.0/8,-1.0/8,0,1.0/4,0,0,-1.0/2,0,
					  1.0/8,-1.0/8,0,-1.0/4,0,0,0,1.0/2,
					  1.0/8,-1.0/8,0,-1.0/4,0,0,0,-1.0/2};
					  
	cout<<"设第1次相乘的矩阵为M1，第2次为M2，第3次为M3。W=M1*M2*M3，W矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<W[i][j]<<"\t";
		}
		cout<<endl;
	}
	
	
	/*源图像的一个图像块矩阵*/
	double A[N][N] = {64,2,3,61,60,6,7,57,
					  9,55,54,12,13,51,50,16,
					  17,47,46,20,21,43,42,24,
					  40,26,27,37,36,30,31,33,
					  32,34,35,29,28,38,39,25,
					  41,23,22,44,45,19,18,48,
					  49,15,14,52,53,11,10,56,
					  8,58,59,5,4,62,63,1};
	
	cout<<"源图像的一个图像块矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<A[i][j]<<"\t";
		}
		cout<<endl;
	}
	
	/*矩阵W的转置矩阵*/
	double WT[N][N];
	MatrixTranspose(W,WT);
	
	/*求均值差值后的矩阵 T = W'*A*W */
	double T[N][N];
	double temp[N][N];
	MatrixMultiplication(WT,A,temp);// temp = W'*A
	MatrixMultiplication(temp,W,T);// T = temp*W
	
	cout<<"求均值差值后的矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<T[i][j]<<"\t";
		}
		cout<<endl;
	}
	
	/*设阈值为5，绝对值小于等于5的元素设为0 */
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			if(abs(T[i][j])<=5){
				T[i][j] = 0.0;
			}
		}
	}
	
	cout<<"设置阈值变换后的矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<T[i][j]<<"\t";
		}
		cout<<endl;
	}
	
	/*本来想通过用逆变换式子直接用数学公式 A = (W')^(-1)*T*W'  但是C++没有求逆矩阵的库函数，自己编写又太长，故而使用下面这照片这种模拟人工运算方法 */
		
	double tem[N];//该数组作为每一次由均值和差值计算出的结果 
	/*列逆变换*/
	for(int j=0;j<N;j++){
		for(int k=0;k<3;k++){//每列8个元素，需要循环3次均值和差值计算
			int t = pow(2,k);//t为每次运算结束的位置 
			for(int i=0;i<t;i++){
				tem[2*i] = T[i][j]+T[i+t][j];	//运算得到结果的第一个 
				tem[2*i+1] = T[i][j]-T[i+t][j]; //运算得到结果的第二个
			}
			for(int i=0;i<t;i++){				//一次循环完成后要更新T矩阵，因为下次循环要基于上一次循环的结果，但是为什么不在上面那个循环里面不直接赋值给T，不用tem呢，因为上面循环赋值会导致上面运算后续的一些运算的值被覆盖导致错误结果 
				T[2*i][j] = tem[2*i];
				T[2*i+1][j] = tem[2*i+1];
			} 
		}
	}
	
	cout<<"列逆变换后的矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<T[i][j]<<"\t";
		}
		cout<<endl;
	}


	/*行逆变换*/
	//原理与列逆变换一样 
	for(int i=0;i<N;i++){
		for(int k=0;k<3;k++){//循环3次 
			int t = pow(2,k);
			for(int j=0;j<t;j++){
				tem[2*j] = T[i][j]+T[i][j+t];
				tem[2*j+1] = T[i][j]-T[i][j+t];
			}
			for(int j=0;j<t;j++){
				T[i][2*j] = tem[2*j];
				T[i][2*j+1] = tem[2*j+1];
			} 
		}
	}
	
	cout<<"行逆变换,得到重构后的近似矩阵:"<<endl;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			cout<<T[i][j]<<"\t";
		}
		cout<<endl;
	}
	
	return 0;
} 


运行结果：




    四、实验结论、实验体会
    本次实验的算法思路很清晰、也比较简单。但是这是相对于来说的，假如是用MATLAB编程，则只要几分钟就能完成这次的实验，因为MATLAB是强大的矩阵处理语言，有非常多的系统类库进行调用，而且语法非常简单。而用C/C++则是非常的麻烦，所有的细节都要自己编写，而且每个细节都要编好长的一段，编码时非常地痛苦，花了很长时间去实现。在进行逆变换时，要用到求解逆矩阵，而求逆矩阵C/C++却没有类库，需要大量的编码工作，所以我选择通过模拟人工计算方式进行逆变换，在这个模拟过程中需要不断地计算每一时刻计算所处的位置、值、和所要存放的位置，存放时是否会跟其它运算是否冲突造成值覆盖等等，非常地麻烦。将整个算法在脑子里一步一步地演进找出所有可能的错误，然后在VS2010进行DEBUG跟踪程序执行，最终完成了此次实验。
    这次实验让我更深刻地理解使用数学方法对算法进行精简不仅可以减轻编程人员的负担，同时也能简化计算，提高效率。


写得不是很好，望各位大神多多指正，不喜勿喷。

简介
笔者曾在《Face swapping》一文中提到过人脸检测，使用流行的“哈尔特征（Haar-like features）级联”算法逐步执行，就能在图像中获得初级的人脸轮廓。
在本文中，笔者将会详细介绍一个转化这种人脸检测算法的自制脚本：取代原本在图像中确认是否包含人脸的做法，这种新办法会通过单纯的级联数据生成一张人脸图片。
跟其他文章一样，这篇文章的全部代码也能在GitHub下载到。
哈尔级联（Haar Cascades）
在2001年，Viola和Jones两位大牛推出了革命性的对象检测算法，这种算法基于哈尔特征级联，首次支持实时监测视频数据中的人脸（及其他对象）。
其核心在于通过一张小图片（一般是20x20）与一些预先算好的级联数据（见下文描述），返回对象是否存在的信息。之后用多个大小与位置不同的窗口使用核心算法对完整图片进行检测，区分目标与非目标：
笔者尝试修改的就是这个核心算法。不过，如何操作呢？实际上是基于所谓的“哈尔特征”：
每个特征都与构成所谓“弱分类器”的阈值相关。如果用黑色区域的单个像素点总和减去白色区域的单个像素点总和，结果超出阈值的话则代表通过了弱分类器。例如，在第一张图片中，弱分类器检测到，与脸颊上方对比，眼睛周围有一片深色区域。
所有特征均由轴对齐矩形定义，采用如上图所示的基本形式之一。检测窗口使用与输出图片尺寸相同的小型网格（比如20x20）。
弱分类器组合成阶段。该阶段是否通过取决于相关的弱分类器是否通过；每个弱分类器都有相关权值，如果所有通过的弱分类器总和的权值超过了阶段阈值，那么就表示这个阶段通过了。
阶段、若分类器与相关权值通过运行训练算法来计算，可以在Viola和Jones的论文中查看更多信息。
如果所有阶段都通过了，那么算法就会返回检测到该对象的信息。
相应的Python代码：
def detect(cascade, im):
    for stage in cascade.stages:
        total = 0
        for classifier in stage.weak_classifiers:
            if numpy.sum(classifier.feature * im) >= classifier.threshold:
                total += classifier.weight
        if total < stage.threshold:
            return False
    return True
假定输入图片im已经调整为与小型级联大小一致。classifier.feature是一个与im形状相同的数组，白色特征区域记为1s，黑色特征区域记为-1s，其他区域为0。注意，Viola和Jones所描述的算法实际上使用了积分图来统计像素值，因此这种算法速度很快。
分为多个阶段的原因主要是为了效率：一般来说，第一阶段只包含少量特征，却能否定大多非人脸图像。在一个特定的级联中，一般有数百个特征，与十几个阶段。
混合整数线性规划（MILP）
为了转化上述的检测功能，笔者用混合整数线性规划表述了这个问题，然后在线性规划中使用MILP求解器。
下面是按照MILP约束描述的检测功能。首先约束能够确保弱分类器在需要时通过。我们称其为特征约束：
……还有一系列约束可确保每个阶段通过，我们称其为阶段约束：
其中：
是输入图像的像素值（在代码中以im指代）。 
是与弱分类器j相关联的特征值权值。（在代码中用classifier.feature指代）
是弱分类器j的阈值。（在代码中以classifier.threshold指代）
是弱分类器j的权值（在代码中以classifier.weight指代）。
是一系列指示正权值的分类器，即
是一系列指示负权值的分类器，即
是一个二进制指示变量，指代弱分类器j是否通过。
是选中足够大的数字，确保所在项非零，不等式正确。
是一系列指示k阶段的弱分类器。
是k阶段的阶段阈值（在代码中以stage.threshold指代）。
MILP求解器寻找的变量是和的值。其余值源自级联定义自身。
主要注意点就是将用作指示变量，即其状态如何打开或关闭其中一个特征约束。比如，取，如果特定的解决方案作为1通过，则我们最好确定特征j实际上超过了分类器的阈值，因为它对阶段约束通过产生了积极作用。在这种情况下，关联特征约束中的项为0，即约束“打开”。
相反，在分类器j为负权值时，我们只关心在为0的情况下，该特征没有通过分类器阈值。
考虑到这一点，很明显当且仅当有线性规划的解决方案取自cascade：变量取im中的相应像素值时，detect(cascade, im)为True。
在Python中的MILP
笔者选用了docplex模块以Python编写上述约束。将通过这个模块所编写的约束发送到IBM的DoCloud服务，用远程解决。注意：DoCloud需要注册，而且是收费的，尽管有一个月的免费试用期。我最初用免费的PuLP模块尝试解决这个问题，但由于底层求解器过于简单或是我机器自身的限制，在这方面进展有限。
例如，下面是定义变量的方式：
from docplex.mp.model import Model

model = Model("Inverse haar cascade", docloud_context=docloud_context)

pixel_vars = {(x, y): model.continuous_var(
                name="pixel_{}_{}".format(x, y)
                lb=0., ub=1.0)
                for y in range(cascade.height)
                for x in range(cascade.width)}
passed_vars = [self.binary_var(name="passed_{}".format(idx))
                                       for idx in range(len(cascade.features))]
（节选）片段，增加阶段约束：
for stage in cascade.stages:
    model.add_constraint(sum(c.weight * passed_vars[c.feature_idx]
                                        for c in stage.weak_classifiers) >=
                                                               stage.threshold)
通过调用model.solve()解决掉了问题。
如果成功的话，会从解决方案中提取像素变量值，并转化为图像（numpy数组），之后可以用cv2（或类似的）写入磁盘。
im = numpy.array([[pixel_vars[x, y].solution_value
                    for x in range(cascade.width)]
                  for y in range(cascade.width)])
cv2.imwrite("out.png", im * 255.)
默认docplex只会寻找可行的解决方案。但可以设置下面这样的对象：
model.set_objective("max",
    sum((c.weight * passed_vars[c.feature_idx] 
            for s in self.cascade.stages
            for c in s.weak_classifiers)))
该对象会尝试找到解决方案，让大多数结果超过阶段约束。由于可能花费过长事件来找寻真正的最大值，我们将时间设置为：
model.set_time_limit(60 * 60)
求解器会在一个小时后自动停止运算，输出当时的最佳方案（如果有的话）。 

更多细节请参见源码。注意：由于支持tilted feature，另外还因为级联数据是OpenCV格式的，因此代码可能会有点复杂，不过本质是相同的。
结果
下面是在OpenCV的haarcascade_frontalface_alt.xml级联中运行程序一小时后的输出结果：
还不错。抱歉分辨率太低，不过由于仅定义在20x20的网格中，这点难以避免。下面是同样的一幅图，也很模糊，不过更像一张脸了：
要想测试检测的极限，我们需要将阶段约束最小化，而不是最大化：
绝对不太像人脸，不过用OpenCV应当能检测出来。
下面是最好的眼睛图像（基于haarcascade_eye.xml）。
……这幅是最差的：
这幅是最好的侧脸图像（基于haarcascade_profileface.xml）：
……这幅是最差的：

  
英文：Making faces with Haar cascades and mixed integer linear programming 

  作者：Matthew Earl 

  翻译：孙薇 

  审校/责任编辑：唐小引

 


文章目录

摘要
	
一、哈尔角点检测
	
1. 原理
		
2. 公式推导
		
3. 存在问题
	
	
二、拟合外接矩形求角点
	
1. 钝角
		
2. 直角
		
3. 锐角
	
	
三、检测的例子
	
1. 获取轮廓，求外接矩形
		
2. 角点的排序
	
	
总结


 

摘要

角点在三维场景重建运动估计，目标跟踪、目标识别、图像配准与匹配等计算机视觉领域起着非常重要的作用。在现实世界中，角点对应于物体的拐角，道路的十字路口、丁字路口等。在近几个项目中，通过获取的角点，多用于透视变换。

一、哈尔角点检测

1. 原理

根据角点的定义，可以分为如下两类：角点可以是两个边缘的角点，邻域内具有两个主方向的特征点。前者很大程度上依赖于图像的分割与边缘提取，具有相当大的难度和计算量。 第二种如Harris角点的检测。如果在各个方向上移动这个widow，

区域内的灰度发生较大的变化，就认为遇到了角点。



2. 公式推导

周围灰度变化的自相关函数如下：



进行泰勒展开如下，其中Ο(u2,v2)近似为0：





矩阵M又称Harris矩阵：



在不同方向的展示如下：

    a. 图像中的直线。一个特征值大，另一个特征值小，λ1>λ2或λ2>λ1。自相关函数值在某一方向上大，在其他方向上小。

    b. 图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。

    c. 图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。



3. 存在问题

实际的烟盒不完全是直角，导致角点检测出现偏离。





 

二、拟合外接矩形求角点

相机拍摄的图像如下所示，在绿框中寻找角点的位置。





1. 钝角

先求出初始的外接矩形，如下所示。



再在该初始的外接矩形进行图像的取反，再次进行外接矩形的拟合，输出对应位置的矩形点。



2. 直角

寻找最小外接矩形，进行角点的排序，从而求得对应位置的角定。



3. 锐角

用findRect在二值化的cornerRoi图像上求外接矩形，再对该矩形的四个点进行排序，烟盒左下方的角点是对应拟合矩形的左下方角点，烟盒右下方的角点是对应拟合矩形的右下方角点。



 

三、检测例子

1. 获取轮廓，求外接矩形。

#近似边缘求角点point
def getSubCornerPoint(binImg, posFlag):
    #寻找近似的轮廓
    resContourFlag, approxContour = findApproxContour(binImg, posFlag)
    if resContourFlag ==False:
        return False, [], []
    #轮廓的外接矩形
    boxPoints = getBoundRect(approxContour)
    #角点排序
    sortBoxPoints = sortPoints(boxReverPoints)
    #获取对应位置的point
    resPoint = getPosSortPoint(sortBoxPoints, posFlag)

#寻找近似的轮廓
def findApproxContour(binImg, posFlag):
    allContours, _ = cv2.findContours(binImg, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    if len(allContours) ==0:
        print("@ error findApproxContour not find the contour pos: ", posFlag)
        return False, []
    #计算轮廓面积，将最大的那组轮廓挑出来
    contour = sorted(allContours, key=cv2.contourArea, reverse=True)[0]
    lengthEps = 0.015  #轮廓边长的近似
    contourMinArea = (binImg.shape[0])*(binImg.shape[1])*fourCornerBinImgRation
    contourArea = cv2.contourArea(contour)
    if (contourArea > contourMinArea):
        epsilon = lengthEps * cv2.arcLength(contour, True)  #精度的近似
        approx = cv2.approxPolyDP(contour, epsilon, True)   #轮廓的近似
        return True, approx
    else:
        print("@ failed findApproxContour not ok area contour pos: ", posFlag)
        return False, []

#轮廓求外接矩形  返回tl, tr, bl, br
def getBoundRect(approxContour):  
    x,y,w,h = cv2.boundingRect(approxContour) #外接矩形
    boxPoints = np.array([[x, y], [x+w, y], [x, y+h], [x+w, y+h]])
    return boxPoints

2. 角点的排序

#位置点的排序 tl, tr, bl, br
def sortPoints(pts):
    xSorted = pts[np.argsort(pts[:, 0]), :]
    leftMost = xSorted[:2,:]
    rightMost = xSorted[2:,:]
    leftMost = leftMost[np.argsort(leftMost[:, 1]),:]  #np.argsort返回坐标的index
    (tl, bl) = leftMost
    rightMost = rightMost[np.argsort(rightMost[:, 1]), :]#y坐标偏上的是top
    (tr, br) = rightMost
    return np.array([tl, tr, bl, br], dtype=int)

 

总结

自带的角点检测方法，稳定行不高。可以采用外接矩形的求对应位置的角点。

【0】README
0.1）本文旨在讲解 哈儿小波变换（分解和重构）进行数据的降维和升维；


【timestamp: 1703281610】时隔几个月再来review
 哈儿小波变换算法的具体思路：
1）分解降维：首先对所有item进行分解降维，求相邻维度的两个元素的和均值和差均值，如 array[0] 和 array[1]为一组,array[2]和array[3]为一组；分别存储在 array[i] array[i+span]， 其中span=该当前分辨率下的 分辨率长度，如当前分辨率为resolution,
 则 length=2^resolution；
2）item重构升维：当前分辨率为curResolution，则 

int span = (int) pow(2, curResolution);
for (int i = 0; i < span; i++) {
	temp[2*i] = rawitem[i] + rawitem[i+span];
        temp[2*i+1] = rawitem[i] - rawitem[i+span];
}

因为 rawitem[i] 和 rawitem[i+span]是 降维时的和均值和差均值；
3）质心重构升维： 只需要将较低分辨率下的质心扩展一倍，如
if(ClusterData.resolutionLength < ClusterData.dimension) {
	for (int j = ClusterData.resolutionLength-1; j >= 0 ; j--) {
		ClusterData.centroid[i][2*j] = ClusterData.centroid[i][j];
		ClusterData.centroid[i][2*j+1] = ClusterData.centroid[i][j];
	}
}


【1】intro to 哈儿小波
1.1）定义：哈尔小波变换分为哈尔小波分解和重构。本文利用哈尔小波分解和哈尔小波重构分别对社交网络中的时间序列进行降维和升维操作，而分解降维和重构升维互为逆过程，这里以哈尔小波分解降维步骤为例进行说明。
1.2）哈尔小波分解（重构是其逆过程）：可以被看作是一系列对离散时间函数的均值和差操作。本文计算离散函数f(t)，即时间序列元素值的每两个相邻值间的均值和差。找出离散函数f(t)={9， 7， 3， 5}的哈尔小波分解过程如下表所示：





分辨率

级别



分辨率

（维度）



均值



差值



哈尔小波

分解系数



降维后的

时序



全分辨率下的时序



2



4



{9,7,3,5}







{9,7,3,5}







{9,7,3,5}



1



2



{8,4}



{1,-1}



{8,4,1,-1}



{8,4}



{8,8,4,4}



0



1



{6}



{2}



{6,2,1,-1}



{6}



{6,6,6,6}

对上表的分析（Analysis）：
A1）intro
 to 分辨率：哈儿小波分解和重构常应用于图像的降维和升维，这里的分辨率，你与可以理解为是图像的分辨率，分辨率越低，图像质量越不清晰，分辨率越高图片质量越清晰（数据显示更加丰富）；（干货——我们看到的图片为什么不清晰）
A2）上表中最后一列（全分辨率下的时序）：我们看到数据由 [9 7 3 5] 到 [8 8 4 4] 再到 [6 6 6 6]，显然在一个图像的某个区域内，分辨率越低，数据越不丰富，变化越不明显，这也是为什么低分辨率下数据模糊的一个原因（因为低分辨率下，数据经过降维处理了）；


【2】将上述哈尔小波分解过程用倒置二叉树表示如下：


【3】看个荔枝





对上图的分析（Analysis）：
A1）以上数据[2, 5, 8, 9, 7, 4, -1, 1]在不同分辨率下经过哈儿小波变换得到的数据如下：
分辨率3：它本身；
分辨率2：3.5， 8.5, 5.5, 0, -1.5, -0.5, 1.5, -1;
分辨率1：6， 2.75， -2.5， 2.75, -1.5, -0.5, 1.5, -1;
分辨率0：4.375， 1.625， -2.5， 2.75, -1.5, -0.5, 1.5, -1;
A2）哈儿小波变换分为哈儿小波分解和哈儿小波重构，且我们可以指定最低级分辨率（并不一定是0）；


【4】将哈儿小波变换应用到Kmeans聚类算法（干货——目的是降低聚类时间复杂度）
1）Kmeans聚类进行前：先对数据进行哈儿小波分解以降维到0分辨率（或你自己指定）（干货——聚类开始前，我们数据的分辨率是0，而维度是1，即2^0=1）；
2）Kmeans聚类过程中：聚类是一个迭代的过程，每轮迭代，我们都对数据进行哈儿小波重构进行升维（以降低聚类时间复杂度）；
如第1轮聚类后，我们将数据升维到分辨率1下；第2轮迭代后，将数据升维到分辨率2下，以此类推......；当达到全分辨率下时，我们不再进行哈儿小波变换，但可以继续聚类，直到满足聚类结束标志的条件为止；


【5】哈儿小波分解和重构源代码实现 
5.1）哈儿小波分解和重构的源代码 
public class HaarTransform{
	public static int maxResolution;
	/**
	 * executing haar reconstruction transform towards given array
	 * @param rawitem raw array
	 * @param curResolution current resolution
	 */
	public static void haarReconstruct(double[] rawitem, int curResolution)	{
		double[] temp = new double[rawitem.length];
		System.arraycopy(rawitem, 0, temp, 0, temp.length);
		
		int span = (int) pow(2, curResolution);
		for (int i = 0; i < span; i++) {
			temp[2*i] = rawitem[i] + rawitem[i+span];
			temp[2*i+1] = rawitem[i] - rawitem[i+span];
		}
		System.arraycopy(temp, 0, rawitem, 0, temp.length);
	}
	
	/**
	 * @param rawitem is a array
	 * @param leastResolution is the least resolution  allowed be zero.
	 * @param maxResolution is equals to log(full dimension)
	 * @return
	 */
	public static double[] haarDecompose(double[] rawitem, int leastResolution)	{
		double[] temp = new double[rawitem.length];
		int resolutionLength = rawitem.length;  
		int maxResolution = HaarTransform.maxResolution;
		// if rawitem.length=8, maxResolution = 3; 
		
		// you know least resolution equals to 0.
		for (int i = maxResolution; i > leastResolution; i--) {
			System.out.println("第"+i+"级分辨率下");
			AlgTools.printOneDimArray(rawitem);
			for (int j = 0; j < resolutionLength / 2; j++) {
				temp[j] = (rawitem[2*j] + rawitem[2*j+1]) / 2;
				temp[resolutionLength/2+j] = (rawitem[2*j] - rawitem[2*j+1]) / 2;
			}
			resolutionLength /= 2;
			System.arraycopy(temp, 0, rawitem, 0, temp.length);
		}
		return temp;
	}
}
5.2）哈儿小波分解和重构的测试用例

public class HaarTransformTest {
	public static void main(String[] args) {
		double[] data = {2, 5, 8, 9, 7, 4, -1, 1};
		int resolutionNum = (int) (log(data.length) / log(2));
		System.out.println("分辨率总级别个数为： " + resolutionNum);
		
		HaarTransform.maxResolution = resolutionNum;
		int leastResolution = 0;
		HaarTransform.haarDecompose(data, leastResolution);
		System.out.println("最低级分辨率"+leastResolution+"下");
		AlgTools.printOneDimArray(data);
		// 哈儿小波分解over.
		
		for (int i = leastResolution; i < resolutionNum; i++) {
			HaarTransform.haarReconstruct(data, i);
			System.out.println("第"+ (i+1) + "级分辨率下");
			AlgTools.printOneDimArray(data);
		}// 哈儿小波重构over.
	}
}

对以上代码的分析（Analysis）：
A1）因为总共维度是8，所以最大分辨率==log8=3；
A2）上述哈儿小波分解过程中对应的变量变化如下（为了便于更加理解分解过程）：


循环次数
round


分辨率级别 i 


分辨率级别的对应维度
resolutionLength


j(max)


下一分辨率级别
的维度sum


1


3


8


3(0~3)


4


2


2


4


1(0~1)


2


3


1


2


0(0~0)


1


4


0


1












 
 
 
 



5.2）打印结果（共8维）（也可以指定哈尔小波分解的最低级分辨率，如1,2，但最小是0）
分辨率总级别个数为： 3
第3级分辨率下

  2.000  5.000  8.000  9.000  7.000  4.000 -1.000  1.000
第2级分辨率下

  3.500  8.500  5.500  0.000 -1.500 -0.500  1.500 -1.000
第1级分辨率下

  6.000  2.750 -2.500  2.750 -1.500 -0.500  1.500 -1.000
最低级分辨率0下

  4.375  1.625 -2.500  2.750 -1.500 -0.500  1.500 -1.000
第1级分辨率下

  6.000  2.750 -2.500  2.750 -1.500 -0.500  1.500 -1.000
第2级分辨率下

  3.500  8.500  5.500  0.000 -1.500 -0.500  1.500 -1.000
第3级分辨率下

  2.000  5.000  8.000  9.000  7.000  4.000 -1.000  1.000