N-puzzle Problem

N-Puzzle Problem

[TOC]

Preview

N-Puzzle Problem:

在人工智能中,常常以N-Puzzle 问题的求解为例子来说明各种搜索策略。

N=K*K-1;

当N=8时,即为重排九宫格问题;当N=15时,即为十五迷问题;

一个状态:N个数码及一个空格的每一组置放形式,就形成了该问题的一个状态。

N-Puzzle问题空间:所有(N+1)!个互不相同的状态构成了一个N-Puzzle问题空间。

128M 1e7

操作:平移操作——把与空格相邻 的某一数码移入空格,称这种操作为平移操作。

Problem:对于给定的任一初始状态,使用平移操作,求解到目标状态的最优步数 (最少步骤的移法) 以及移动的每一步状态。

Goal position:

1556368035445

N-Puzzle Problem 的可解性判断

分奇偶性然后根据逆序数判断。

https://blog.csdn.net/Wood_Du/article/details/88730885

Algorithms

​ N-Puzzle问题之所以能用来为例说明各种搜索策略,在于对于N-Puzzle问题,根据所需求解的初始状态的复杂性,我们可以从爆搜开始,一步一步优化,当某种算法不足于求解时,可从优化数据结构,更换搜索算法上升一个境界,求得更加复杂的初始状态的解。

:你会搜索吗?

:我会爆搜 dfs+STL库,但是DFS不能找到最优解;找最优解我会bfs

:那么存状态呢?

:STL库会很废效率,我可以用hash判重;正好学了Zobrist Hash

:你还有什么可以进一步优化吗?

:为了减少状态的膨胀,我可以使用双向广搜,从初始状态和目标状态同时开始搜索,当某方向遇到另一个方向搜索过的状态的时候,则搜索成功。两个方向对接得到最后结果。

:这些都是无信息搜索算法,你可以使用有信息搜索吗?

:启发式搜索,那就想办法减少搜索范围,降低问题复杂度。这个时候我们就可以上A*+Hash了。启发函数可以用错位数或者曼哈顿距离。

:对于A*,需要每次找到Open表中f(f(n) = g(n)+h(n)) 最小的的元素。如果排序那么工作量会很大应该怎么办呢?

:上优先队列

:但是对于15码问题,内存不能存下状态了

: IDA*,解决内存问题。

是基于迭代加深的A*算法。并且不需要判重,不需要估价排序,用不到哈希表,空间需求量变得超级少。启发函数用曼哈顿距离。在找最优解的时候,对超过目前最优解的地方进行剪枝,可导致搜索深度急剧减少。再check一下不要回到上个状态。

:一般的15码样例都能解决,但是对于15码最复杂的一些样例,仍然需要三四十分钟才能解出答案。还有什么可以优化吗?

:那就只能去啃英文论文了。A.Felner,R.Korf的Disjoint Pattern Database Heuristics, 不相交模式数据库启发式算法。Additive Pattern Database Heuristics ——E.Korf,A.Felner,S.Hanan。这是一中新的设计更精确的可容许启发式评价函数的方法,该方法基于模式数据库。

本文我们将略去无信息搜索算法,求解最小步数,详解算法为:

从A* 开始,优化数据结构;

再到IDA*;

最后用Pattern Database Heuristics解决15码最复杂的样例。

针对不同复杂度的样例,本文后面将分为三个阶段。

每个阶段将说明解决的样例,求解的要求时间,算法以及启发函数,代码实现以及每个样例的解决时间。

对于展示代码因为代码量很大,并且是一个框架,所以只能展示核心代码,若需要查看整个项目代码,我会push到github上。

First phase:

所需解决样例以及最多时间:

对于以下两个实例,均能够在1秒之类解出

1556377705625

Algorithm:A*

​ A*算法是一种静态网中求解最短路径最有效的直接搜索方法也是也是一种启发式算法。在搜索过程中使用启发函数。估价函数与实际值越接近,最终搜索到目标格局的时间越快。

估价函数表示:

1
2
3
4
5
6
f(n) = g(n) + h(n)
//f(n) 表示从初始状态到目标状态的估测代价
//g(n) 表示从初始状态到当前状态(节点n)的代价(已经确定)
//h(n) 表示从当前状态(节点n)到目标状态的估测代价(预测)

h(n) < h'(n) //h'(n)为当前状态到目标状态的实际步数

h(n) 的好坏直接影响评估函数的好坏。

流程图:

//N-puzzle-Problem\

1556381217714

​ 从初始状态出发 =>经过一系列中间状态 =>最终到达目标状态(或者无法到达)。

​ 该算法用于经过中间状态时候的行进策略(其中的中间状态或者由题目给出,或者在前边已经推导得出)。

优点:与广度优先搜索策略和深度优先搜索策略相比,A*算法不是盲目搜索,而是有提示的搜索

缺点:该算法一般要使用大量的空间用于存储已搜索过的中间状态,防止重复搜索。随着样例变复杂将会爆内存。

Code:

完整代码传送门:???

Search为A*算法,此代码由于调用了框架中其他类的代码所以比较抽象,我是边百度java边写这个项目,没学过java的娃娃面向百度编程。所以如果你会java,相信你吨吨吨就看明白我代码了。

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
package core.astar;

import core.problem.Action;
import core.problem.Problem; //操作
import core.problem.State;

import java.util.ArrayList;

public class AStar {

public AStar(Problem problem) {
super();
this.problem = problem;
}

public Node childNode(Node parent, Action action) {
State state = problem.result(parent.getState(), action);
//getPathCost 初始到当前的实际代价 stepCost 从父级到后续任务的路径成本
int pathCost = parent.getPathCost() + problem.stepCost(parent.getState(), action);
int heuristic = problem.heuristic(state); //估计从状态到目标状态最便宜路径的成本
return new Node(state, parent, action, pathCost, heuristic);
}

public Problem getProblem() {
return problem;
}

public void setProblem(Problem problem) {
this.problem = problem;
}

public Node Search()
{
//起始节点
State initState = problem.getInitialState();
Node node = new Node(initState, null, null, 0, problem.heuristic(initState));

Fringe fringe = new Fringe();
fringe.insert(node);

Explored explored = new Explored();

while (true) {

if (fringe.isEmpty()) return null; //失败

node = fringe.pop(); //choose the lowest-cost node in frontier
if (problem.goalTest(node.getState())) return node;
explored.insert(node.getState());

for (Action action : problem.Actions(node.getState())) {
Node child = childNode(node, action);
//child.getState()不再在扩展节点的集合(Close表中)且fringe(Open表)中不存在状态为state的节点 则将节点插入fringe中
if (!explored.contains(child.getState()) && !fringe.contains(child.getState())) {
fringe.insert(child);
}
else {
/**
如果邻接结点N’在CLOSED表中,比较CLOSED表中的g(N')值和当前的g(N')值,如果当前的g(N')更小,那么删除CLOSED表中的N',把N设成N'的父节点,重新将N'插入OPEN表。
如果邻接结点N’在OPEN表中,比较OPEN表中的g(N')值和当前的g(N')值,如果当前的g(N')更小,那么删除OPEN表中的N',把N设成N'的父节点,重新将N'插入OPEN表。
***/
Node revisited = fringe.revisited(child.getState());
if (revisited != null && revisited.evaluation() > child.evaluation()) {
fringe.replace(revisited, child);//用child节点代替Fringe中的 revisited节点
}
}
}

}
}

//用动画展示问题的解路径
public void solution(Node node)
{
// Fix me
// 调用Problem的drawWorld方法,和simulateResult方法
problem.drawWorld();
}

private Problem problem;
}

Result:

项目目录:Searching_Student: E:\University\AI\Searching_student

No Initial State Steps Time Limited Time
1 8, 6, 7, 2, 5, 4, 3, 0, 1 31 0.132s 1s
2 6, 4, 7, 8, 5, 0, 3, 2, 1 31 0.133s 1s
3 8, 13, 0, 6, 1, 15, 9, 14, 3, 4, 5, 11, 7, 2, 10, 12 52 1.929s
4 2,9,5,11, 8,3,4,14, 7,10,1,12, 0,15,6,13 GC
5 4,7,0,9,12,10,11,8,14,6,15,1,2,5,3,13 GC

Second phase:

所需解决样例以及最多时间:

4阶的,能够在与老师程序同级别的时间内,解决相同 的问题实例。

1556382042355

1556382108542

对于以上后面几个初始状态,A*就已经无法解决。因为会生成太多新状态并消耗大量内存维护这些列表。直接爆内存。

Algorithm:IDA*

IDA*算法, ID(Iterative Deepening)指的是迭代加深.

迭代加深搜索算法,在搜索过程中采用估值函数,以减少不必要的搜索。

IDA*算法核心:

​ 设置每次可达的最大深度depth,若没有到达目标状态则加深最大深度。

​ 采用估值函数,剪掉f(n)大于depth的路径

IDA*算法的步骤

a、首先对初始状态进行评估, 评估值作为最小限度, 而最大限度为自己的设置.这个评估值在这个问题中可以用此状态到正确状态的每个位置的曼哈顿距离来表示.

b、从最小限度到最大限度进行遍历, 此值作为当前dfs的限度值, 这个限度不断在有效范围内递增的过程就称作迭代加深

c、进行dfs, 调整状态, 将新状态加入到新的dfs中, 直到找到了一个解(由于迭代加深, 此解为最优解). 进行回溯, 加入路径, 算法结束.

IDA* is described as follows:

  • Set threshold equal to the heuristic evaluation of the initial state.
  • Conduct a depth-first search, pruning a branch when the cost of the latest node exceeds threshold. If a solution is found during the search, return it.
  • If no solution is found during the current iteration, increment threshold by the minimum amount it was exceeded, and go back to the previous step.

  • 设置阈值等于初始状态的启发式评估。

  • 进行深度优先搜索,在最新节点的成本超过阈值时修剪分支。如果在搜索过程中找到解决方案,请将其返回。
  • 如果在当前迭代期间未找到解决方案,则将阈值增加超过最小值,然后返回上一步骤。

Code:

完整代码传送门:???

IDA*代码除了读取数据时调用了框架中的类,其他在一个类完成。

当时已经对框架有点绝望。由于我使用反向搜索使用框架会改很多类。加上不能与滑块公用一个算法,,所以写在一个类里面。

第二阶段我采用IDA*反向搜索,从目标状态到初始状态搜索,所以看代码要注意。速度会快很多,对于给定样例有的比老师快几十倍,分析过原因。不过后面理论又被自己否定,应该是由于样例的特殊性所以加快了速度。

main函数中调用:

1
2
ArrayList<Problem> problems = ProblemFactory.getProblems(1);
test2(problems);

test2():

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

public static void test2(ArrayList<Problem> problems) {


for (Problem problem : problems) {

Rstep=0;
flg=0;
result = new int[200][20];
PuzzleState state = (PuzzleState) problem.getInitialState();
side = state.getSide();
tmp = state.getStatus();
System.out.println();
xx = new int[20];
yy = new int[20];
for (int i = 0; i < side * side; i++) {
if (tmp[i] == 0) pos = i;
else {
xx[tmp[i]] = (i / side);
yy[tmp[i]] = (i % side);
}
}

if (!check(tmp)) System.out.println("No");
mat = new int[20];
for (int i = 0; i < side * side; i++) mat[i] = i + 1;
mat[side * side - 1] = 0;
pos = side * side - 1;

double time1 = 0;
Stopwatch timer1 = new Stopwatch();
// for(int i=0;i<100;i++){
// layer[i]=-1;
// }
//System.out.println(H(mat));
for (bound = H(mat); bound <= 100; bound = dfs(0, H(mat), 4))
if (flg == 1) break;
time1 = timer1.elapsedTime();
Ans_show();
System.out.printf("行走了 %s 步,执行了 %.3f 秒\n", Rstep, time1);
}

}

check()方法:检查是否可行解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
public static boolean check(int[] a) {
int tot=0;
int flag=0;
for(int i=0;i<side*side;i++){
if(a[i]==0) {
flag = i/side;
continue;
}
for(int j=0;j<i;j++){
if(a[j]>a[i]){
tot++;
}
}
}
if(side%2==1) return tot%2==0;
else{
tot+= abs(flag-(side-1));
if(tot%2==0) return true;
else return false;
}
}

H()方法:算曼哈顿距离

1
2
3
4
5
6
7
public static int H(int[] mat) {
int h = 0;
for (int i = 0; i < side*side-1; i++) {
h += abs(xx[mat[i]] - (i / side)) + abs(yy[mat[i]] - (i % side));
}
return h;
}

ok()方法:判断是否回到上一个状态

1
2
3
4
5
6
7
8
9
10
11

public static boolean ok(int x,int y){
if(x>y){
int ff=x;
x=y;
y=ff;
}
if(x==0&y==1) return false; //与操作 1 1 为1
if(x==2&&y==3) return false;
return true;
}

IDA*:

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
public static int dfs(int step, int h, int las) {
//g(n)+h(n)=f(n) 更新f(n) bound 最大限度 采用估值函数,剪掉f(n)大于depth的路径
if (step + h > bound) return step + h; //
// if(layer[h]==-1)layer[h]=step;
// if(step>layer[h]+15){
// System.out.printf("enter\n");
// return step+h;
// }
if (h == 0) { // 到达最终状态 输出g(n)即可
// System.out.println(step);
flg = 1;
Rstep = step;
return step;
}
int pos1 = pos;
int ret = 127, x = pos1 / side, y = pos1 % side;
int dx, dy, tar, ht, temp, i;
for (i = 0; i < 4; i++) { //四个方向扩展
dx = x + u[i];
dy = y + p[i];
if (dx < 0 || dy < 0 || dx > side - 1 || dy > side - 1 || !ok(i, las)) continue;
tar = (dx * side) + dy; //计算出扩展出的新节点的一维坐标 2,2 2*4+2= 10
mat[pos1] = mat[tar]; // 0的坐标等于扩展出来的点的坐标 a.mat[10]=11;
mat[tar] = 0;//相当于swap()
pos = tar;
//计算新的h值
//阶段三主要在这修改 用不相交模式数据库 h(n) 为当前节点的各点的曼哈顿距离和。
ht = h - (Math.abs(xx[mat[pos1]] - dx) + Math.abs(yy[mat[pos1]] - dy)) + Math.abs(xx[mat[pos1]] - x) + Math.abs(yy[mat[pos1]] - y);
if (step + ht <= bound)
for (int k = 0; k < side*side; k++)
result[step][k] = mat[k];
temp = dfs(step + 1, ht, i);
if (flg == 1) return temp;
if (ret > temp) ret = temp;
mat[tar] = mat[pos1];
mat[pos1] = 0;
pos = pos1;
}
return ret;
}

Result:

项目目录:SearchingAstar: E:\University\AI\SearchingAstar

No Initial State Steps Time Limited Time
1 8, 6, 7, 2, 5, 4, 3, 0, 1 31 0.000s 0.047s
2 6, 4, 7, 8, 5, 0, 3, 2, 1 31 0.003s 0.047s
3 8, 13, 0, 6, 1, 15, 9, 14, 3, 4, 5, 11, 7, 2, 10, 12 52 0.147s 0.304s
4 2,9,5,11, 8,3,4,14, 7,10,1,12, 0,15,6,13 51 2.423s 3.652s
5 4,7,0,9,12,10,11,8,14,6,15,1,2,5,3,13 56 0.554s 12.367s
6 12, 10, 3, 2, 0, 7, 14, 9, 1, 15, 5, 6, 8, 4, 13, 11 57 4.341s 75.458s
7 12, 1, 5, 6, 2, 11, 7, 9, 14, 10, 0, 4, 15, 3, 13, 8 50 0.058s 3.671s
8 4, 6, 15, 13, 12, 9, 10, 2, 8, 0, 7, 3, 14, 5, 1, 11 61 10.918s 299.286s
9 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 1, 2, 0 70 1761.019s

Third phase:

所需解决样例以及最多时间:

4阶的,能够在1分钟之内,解出以下实例

说明:对于这个样例,曾听说过老师用IDA*大概跑了四十多分钟得出结果,不过由于样例的特殊性,我们可以使用反向搜索,这样会增快速度,我跑的时间不到30分钟,不过这个时间肯定也不是我们所能接受的,我们需要的是在一分钟内

1556378270793

对于15码最复杂的状态,A*肯定不能解,会爆内存。

IDA*虽然在空间上消耗少,但是由于初始状态逆序数大,要搜很久,即大概三四十分钟才能解决。

解决办法:

a、曼哈顿距离加线性冲突——优化不了多少

b、添加不相交模式数据库(Disjoint Pattern Database Heuristics)与静态6-6-3分区

Algorithm:IDA*+Disjoint Pattern Database Heuristics

这个算法就只能硬着头皮去读E.Korf和A.Felner的论文Disjoint Pattern Database Heuristics。

就刚开始看了Tsan-sheng Hsu的一个应该是用来讲课的pdf,看完之后???开始思考三连:我是,我去,我思???大概花了一下午的时间看完那篇pdf,当时是很崩溃的,伴随着这门课的无理取闹与老师的嗯。。。就已经绝望了,此时,滑块那边也没有进展,老师张口闭口严格证明就很无语老师你是真的没搞懂滑块吧,真准备弃了。当天晚上又因为一个评估函数的优劣性被同学怼了,心情很糟糕。

第二天醒来已经中午,开始接手了N-Puzzle和滑块两个问题,开始了找牛逼网友过程。也找到了那篇论文开始啃,(英语真的毁我一生)。

参考论文:

Additive Pattern Database Heuristics ——E.Korf,A.Felner,S.Hanan

Disjoint Pattern Database Heuristics——E.Korf,A.Felner

以下内容来自论文(有些地方翻译会有问题)

不断优化的启发式

1、曼哈顿距离

​ 例如,在滑动拼图拼图中,要将拼贴从位置x移动到位置y,x和y必须相邻,并且位置y必须为空。如果我们忽略空约束,我们得到一个简化的问题,任何瓦片都可以移动到任何相邻位置,并且多个瓦片可以占据相同的位置。在这个新问题中,瓦片彼此独立,我们可以通过沿着最短路径将每个瓦片移动到其目标位置来最佳地解决任何实例,计算所做的移动的数量。

解决这个简化问题的最佳解决方案的成本恰好是曼哈顿从初始状态到目标状态的距离。由于我们删除了对移动的约束,原始问题的任何解决方案也是简化问题的解决方案,并且针对简化问题的最优解决方案的成本是对原始问题的最优解决方案的成本的下限。

不足:曼哈顿距离之所以只是实际解决方案成本的下限,就是没有考虑到滑块的相互作用。

可改进:

通过考虑其中一些相互作用,我们可以计算出更准确的可接受的启发函数。

2、Manhattan Distance+Linear Conflict

Historically, the linear-conflict heuristic was the first significant improvementover Manhattan distance [5]. It applies to tiles in their goal row or column, but reversed relative to each other. For example, assume the top row of a given state contains the tiles (2 1) in that order, but in the goal state they appear in the order (1 2). To reverse them, one of the tiles must move out of the top row, to allow the other tile to pass by, and then move back into the top row. Since these two moves are not counted in the Manhattan distance of either tile, two moves can be added to the sum of the Manhattan distances of these two tiles without violating admissibility.

(从历史上看,线性冲突启发式是曼哈顿距离的第一个显着改进[5]。它适用于目标行或列中的切片,但相对于彼此反转。例如,假设给定状态的顶行按顺序包含tile(2 1),但在目标状态下它们按顺序出现(1 2)。要反转它们,其中一个切片必须移出顶行,以允许另一个切片通过,然后移回顶行。由于这两个移动不计入任一区块的曼哈顿距离,因此可以将两个移动添加到这两个区块的曼哈顿距离的总和而不违反可接受性。)

同样的想法也可以应用于其目标列中的切片。实际上,目标位置的图块可以同时参与行和列的冲突。由于解决行冲突所需的额外移动是垂直移动,而解决列冲突所需的额外移动是水平移动,因此两组移动都可以添加到曼哈顿距离,而不会违反允许性。与曼哈顿距离相比,线性冲突启发式减少了一个数量级的节点数量,成本几乎是速度的两倍,总体加速比超过五倍。接下来的两行用于不相交的模式数据库启发式。

3、非可加性的模式数据库

1556464725925

图显示了十五个拼图拼贴的子集,称为边缘拼贴。对于给定状态,使边缘拼贴到其目标位置所需的最小移动次数(包括其他拼块的所需移动)是解决整个拼图所需的移动次数的下限。

与其他瓦片的位置无关。

因此,我们可以预先计算所有这些值,将它们存储在内存中,并在搜索过程中查找它们。由于有七个边缘拼贴和一个空白,以及十六个不同的位置,边缘拼贴和空白的可能排列总数为16!/(16-8)!= 518,918,400。

对于每个排列,我们存储将数字图块和空白移动到其目标位置所需的移动次数,这需要不到一个字节。因此,我们可以将整个模式数据库表存储在少于495兆字节的内存中。

存储表后,我们使用IDA *搜索特定问题实例的最佳解决方案。在生成每个状态时,使用数字图块和空白的位置来计算图案数据库中的索引,并且使用相应的条目,即解决数字图块和空白所需的移动的数量,作为该状态的启发式值。

使用这种模式数据库,Culberson和Schaeffer将生成的节点数量减少了346倍,并将运行时间缩短了6倍,与曼哈顿距离相比[1]。将此与另一个模式数据库相结合,并将两个数据库值的最大值作为整体启发式值,将生成的节点减少大约一千,将运行时间减少十二。

局限:

非加性模式数据库的主要限制是它们无法解决更大的问题。例如,由于二十四个拼图包含25个不同的位置,一个覆盖n个拼贴和空白的模式数据库需要25!/(25 - n - 1)!条目。仅有六个图块和空白的数据库将需要超过24亿个条目。此外,来自仅六个瓦片的数据库的值将小于所有瓦片的曼哈顿距离。对于多个数据库,允许组合它们的最佳方法是获取其值的最大值,即使这些图块组是不相交的。原因是非加性模式数据库值包括解决模式切片所需的所有移动,包括其他切片的移动

可改进:

我们希望能够在不违反可接受性的情况下对其值进行求和,以获得更准确的启发式,而不是采用最大的不同模式数据库值。这是不相交模式数据库的主要思想。

4、不相交的模式数据库(可加性的)

为了构建滑动拼图的不相交模式数据库,我们将拼贴划分为不相交的组,这样任何拼贴都不属于多个组。然后,我们预先计算每组中瓦片最小移动次数的表格,这些移动是将这些瓦片移动到其目标位置所需的。我们称这些表为一组,每组一个瓦片,一个不相交的模式数据库,或简称为adisjoint数据库。然后,给定搜索中的特定状态,对于每组图块,我们使用这些图块的位置来计算相应表格中的索引,检索解决该组中图块所需的移动次数,然后将每个组的值相加,以计算给定状态的整体启发式。该值至少与曼哈顿距离一样大,并且通常更大,因为它考虑了同一组中的切片之间的交互。

上述不相交数据库和非附加数据库之间的关键区别在于,非加法数据库包括解决图案图块所需的所有移动,包括不在图案集中的图块移动。结果,给定两个这样的数据库,即使它们的区块之间没有重叠,我们也只能将这两个值中的最大值作为可接受的启发式,因为在一个数据库中计数的移动可能会移动另一个数据库中的区块,因此这些举动将被计算两次。在不相交的数据库中,我们只计算组中切片的移动。

这两种类型的数据库之间的第二个区别是我们的不相交数据库不考虑空白位置,减小它们的大小。一个不相交的数据库包含了解决一组图块所需的最小移动量,对于所有可能的空白位置。

该搜索的状态由所讨论的图块的位置和空白的位置唯一地确定,并且仅计算感兴趣的图块的移动。由于滑动拼图拼图的操作员是可逆的,我们可以从其目标位置开始对每个拼贴执行单个搜索,并记录将拼块移动到每个其他位置需要多少移动。对所有图块执行此操作会生成一组表格,这些表格为每个图块的每个可能位置提供距其目标位置的曼哈顿距离。

对于每个图块,我们可以使用上述广度搜索来自动计算将图块从任何位置移动到其目标位置所需的最小移动次数的表格。这样的一组表格将包含从每个位置到其目标位置的每个瓦片的曼哈顿距离。然后,给定一个特定的状态,我们只需在相应的表中查找每个图块的位置并对结果值求和,从而计算曼哈顿距离的总和

因此我们可以将曼哈顿距离相加以获得可接受的启发式

作为一般规则,在对图块进行分区时,我们希望将目标状态中彼此靠近的图块组合在一起,因为这些图块将彼此交互最多。

1556466424757

1556466472252

第一行给出了曼哈顿距离启发式的结果。

第二行是通过线性冲突增强的曼哈顿距离。

第三行表示启发式,它是图4左侧所示的七和八数据库值的总和

第四行表示通过从第三行的启发式开始计算的启发式。然后,我们计算图4右侧所示的七和八数据库值的总和。最后,整体启发式是这两个总和中的最大值

第一个数据列显示了启发函数在1000个初始状态下的平均值。第二列给出了每个问题实例生成的平均节点数,以找到第一个最优解。第三列显示算法的平均速度,以每秒节点数为单位,, on a 440 MegaHertz Sun Ultra10 workstation. 。第四列表示找到第一个最佳解决方案的平均运行时间(以秒为单位)。最后一列给出了生成的平均节点数,以找到问题实例的所有最优解。

使用[12]中开发的分析结果,我们可以预测,单独使用曼哈顿距离解决二十四难题每个问题平均需要大约5万年!对于10,000个初始状态的随机样本,曼哈顿平均距离是76.078个移动,而对于我们的不相交数据库启发式,它是81.607个移动。

采用其最大值来组合来自不同模式数据库的启发式算法。这是最通用的方法,因为任何两个可接受的启发式方法的最大值始终是另一个可接受的启发式方法。我们引入了不相交模式数据库,以允许将来自不同数据库的值相加,从而产生更准确的启发式值。不相交的模式数据库将子目标集划分为不相交的组,然后将解决每个组中所有子目标的成本加在一起。这要求组不相交,并且单个运算符仅影响单个组中的子目标。例如,在滑动拼图游戏中,每个操作员仅移动一个拼贴。这与获取不同值的最大值一样有效,但更准确,并且仍然可以接受。

5、裁剪

used a technique, based on finite-state machines (FSMs), to prune duplicate nodes representing the same state arrived at via different paths in the graph

使用了一种基于有限状态机(FSM)的技术来修剪表示通过图中不同路径到达的相同状态的重复节点[16]。 FSM修剪减少了IDA *在五个问题上产生的节点数量,范围从2.4到3.6。对于这项工作,我们没有使用FSM修剪,因为该技术很复杂,结果取决于所使用的特定FSM,使得其他研究人员难以重现相同的结果

6、动态分区模式数据库启发式 Dynamically-Partitioned Database Heuristics

之前(Korf&Felner,2002)我们展示了如何将滑动拼图块静态划分为不相交的拼贴组以计算可接受的启发式,为每个状态和问题实例使用相同的分区。在这里,我们扩展该方法并显示它也适用于其他域。我们还提出了另一种加法启发式方法,我们将其称为动态分区模式数据库。在这里,我们动态地将问题划分为每个搜索状态的不相交的子问题

论文:Additive Pattern Database Heuristics

15-puzzle ;24-puzzle 动态模式数据库会比静态模式数据库慢。

35-puzzle 动态模式数据库会比静态模式数据库快。

Code:不相交模式数据库静态6-6-3分区

完整代码传送门:???

模式数据库生成器:

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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
/**
* File: PatternDatabaseGenerator.java
* Author: Brian Borowski
* Date created: June 10, 2010
* Date last modified: May 13, 2011
*/
import java.io.BufferedOutputStream;
import java.io.DataOutputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Queue;
import java.util.Set;
import java.util.StringTokenizer;

/**
* Examples:
* When generating a complete pattern database (i.e. no dummy tiles):
* java PatternDatabaseGenerator 8 1,2,3,4,5,6,7,8,0 8-puzzle.db
*
* When generating a disjoint additive pattern database (i.e. dummy tiles 'x'):
* java PatternDatabaseGenerator 15 0,2,3,4,x,x,x,x,x,x,x,x,x,x,x,0 15-puzzle-663-0.db
* java PatternDatabaseGenerator 15 1,x,x,x,5,6,x,x,9,10,x,x,13,x,x,0 15-puzzle-663-1.db
* java PatternDatabaseGenerator 15 x,x,x,x,x,x,7,8,x,x,11,12,x,14,15,0 15-puzzle-663-2.db
*/
public final class PatternDatabaseGenerator {
public static final byte KEY_NOT_FOUND = -1;

final int numOfTiles, numOfTilesMinusOne, dimension;

PrimitiveHashMap tempMap;
List<Entry> stateToCostEntries_8_puzzle;
byte[] costTable_15_puzzle;
long[] configTable_15_puzzle;

public PatternDatabaseGenerator(final int numOfTiles, final long boardConfig,
final byte dummyTile, final String filename) {

this.numOfTiles = Node.numOfTiles = numOfTiles;
this.numOfTilesMinusOne = numOfTiles - 1;
this.dimension = Node.dimension = (int)Math.sqrt(numOfTiles);

DataOutputStream dos = null;
if (filename != null) {
try {
dos = new DataOutputStream(new BufferedOutputStream(
new FileOutputStream(filename)));
} catch (final FileNotFoundException fnfe) {
System.err.println("Error: Cannot open file '" + filename + "' for output.");
System.exit(1);
}
}

if (numOfTilesMinusOne == 15) {
generateFifteenPuzzleDB(dummyTile, boardConfig);
outputFifteenPuzzleData(filename, dos);
} else {
generateEightPuzzleDB(boardConfig);
outputEightPuzzleData(filename, dos);
}
}

private void generateFifteenPuzzleDB(final byte dummyTile,
final long boardConfig) {

final boolean[] tilesInSubset = computeSubset(dummyTile, boardConfig);
final int[] tilePositions = new int[tilesInSubset.length];
int numTilesInSubset = 0;
for (int i = 0; i < tilePositions.length; ++i) {
if (tilesInSubset[i]) {
tilePositions[i] = numTilesInSubset++;
}
}
breadthFirstSearch(boardConfig, tilesInSubset);
final int tableLength = (int)Math.pow(2, numTilesInSubset * 4);
costTable_15_puzzle = new byte[tableLength];
configTable_15_puzzle = new long[tableLength];
for (int i = tableLength - 1; i >= 0; --i) {
costTable_15_puzzle[i] = KEY_NOT_FOUND;
}

System.err.println("Total states visited: " + tempMap.size());
System.err.print("Removing duplicates...");
final Set<Entry> entries = tempMap.entrySet();
final Iterator<Entry> it = entries.iterator();
while (it.hasNext()) {
final Entry entry = it.next();
final long config = entry.getKey();
final byte movesRequired = entry.getValue();

final int index = indexFor(config, true, tilesInSubset, tilePositions);
final byte moves = costTable_15_puzzle[index];
if (moves == KEY_NOT_FOUND || movesRequired < moves) {
configTable_15_puzzle[index] = config;
costTable_15_puzzle[index] = movesRequired;
}
}

int numElements = 0;
for (int i = tableLength - 1; i >= 0; --i) {
if (costTable_15_puzzle[i] != KEY_NOT_FOUND) {
++numElements;
}
}

System.err.println("done");
System.err.println("States in subset: " + numElements);
}

private void generateEightPuzzleDB(final long boardConfig) {
breadthFirstSearch(boardConfig, null);
final PrimitiveHashMap costMap_8_puzzle = new PrimitiveHashMap();

System.err.println("Total states visited: " + tempMap.size());
System.err.print("Removing duplicates...");
final Set<Entry> entries = tempMap.entrySet();
final Iterator<Entry> it = entries.iterator();
while (it.hasNext()) {
final Entry entry = it.next();
final long config = entry.getKey();
final byte movesRequired = entry.getValue();

final byte moves = costMap_8_puzzle.get(config);
if (moves == PrimitiveHashMap.KEY_NOT_FOUND || movesRequired < moves) {
costMap_8_puzzle.put(config, movesRequired);
}
it.remove();
}
System.err.println("done");

final int numElements = costMap_8_puzzle.size();
stateToCostEntries_8_puzzle = new LinkedList<Entry>(costMap_8_puzzle.entrySet());
System.err.print("Sorting entries...");
Collections.sort(stateToCostEntries_8_puzzle, new Comparator<Entry>() {
public int compare(final Entry e1, final Entry e2) {
return e1.getValue() - e2.getValue();
}
});
System.err.println("done");

System.err.println("States in subset: " + numElements);
}

private void outputFifteenPuzzleData(final String filename,
final DataOutputStream dos) {
System.err.print("Writing file...");
if (filename == null) {
// Write values to stdout. User can redirect stdout to a file, if desired.
for (int i = 0; i < configTable_15_puzzle.length; ++i) {
final long config = configTable_15_puzzle[i];
System.out.println((i + 1) + "," + config + "," +
costTable_15_puzzle[i] + "," +
Utility.longToString(config, numOfTiles));
}
} else {
int i = 0;
try {
for ( ; i < costTable_15_puzzle.length; ++i) {
dos.writeByte(costTable_15_puzzle[i]);
}
} catch (final IOException ioe) {
System.err.println("Error: Cannot write entry " + i + " to file.");
System.exit(1);
} finally {
try {
if (dos != null) {
dos.close();
}
} catch (final IOException ioe) { }
}
}
System.err.println("done");
}

private void outputEightPuzzleData(final String filename,
final DataOutputStream dos) {
System.err.print("Writing file...");
final Iterator<Entry> listIter = stateToCostEntries_8_puzzle.iterator();
if (filename == null) {
// Write values to stdout. User can redirect stdout to a file, if desired.
int i = 1;
while (listIter.hasNext()) {
final Entry entry = listIter.next();
final long config = entry.getKey();
System.out.println(i + "," + config + "," + entry.getValue() +
"," + Utility.longToString(config, numOfTiles));
++i;
}
} else {
int i = 0;
try {
while (listIter.hasNext()) {
final Entry entry = listIter.next();
dos.writeLong(((Long)entry.getKey()).longValue());
dos.writeByte(((Byte)entry.getValue()).byteValue());
++i;
}
} catch (final IOException ioe) {
System.err.println("Error: Cannot write entry " + i + " to file.");
System.exit(1);
} finally {
try {
if (dos != null) {
dos.close();
}
} catch (final IOException ioe) { }
}
}
System.err.println("done");
}

private void breadthFirstSearch(final long boardConfig,
final boolean[] tilesInSubset) {
BFSNode currentState = new BFSNode(boardConfig);
currentState.cost = 0;

tempMap = new PrimitiveHashMap();
tempMap.put(boardConfig, (byte)0);

final Queue<BFSNode> queue = new LinkedList<BFSNode>();
queue.add(currentState);

byte previous = 1;
while (true) {
final char fromDirection = currentState.direction;

if (fromDirection != 'R') {
final BFSNode left = currentState.moveLeftNode(tilesInSubset);
if (left != null) {
final byte moves = tempMap.get(left.boardConfig);
if (moves == PrimitiveHashMap.KEY_NOT_FOUND || left.cost < moves) {
tempMap.put(left.boardConfig, left.cost);
queue.add(left);
}
}
}

if (fromDirection != 'L') {
final BFSNode right = currentState.moveRightNode(tilesInSubset);
if (right != null) {
final byte moves = tempMap.get(right.boardConfig);
if (moves == PrimitiveHashMap.KEY_NOT_FOUND || right.cost < moves) {
tempMap.put(right.boardConfig, right.cost);
queue.add(right);
}
}
}

if (fromDirection != 'D') {
final BFSNode up = currentState.moveUpNode(tilesInSubset);
if (up != null) {
final byte moves = tempMap.get(up.boardConfig);
if (moves == PrimitiveHashMap.KEY_NOT_FOUND || up.cost < moves) {
tempMap.put(up.boardConfig, up.cost);
queue.add(up);
}
}
}

if (fromDirection != 'U') {
final BFSNode down = currentState.moveDownNode(tilesInSubset);
if (down != null) {
final byte moves = tempMap.get(down.boardConfig);
if (moves == PrimitiveHashMap.KEY_NOT_FOUND || down.cost < moves) {
tempMap.put(down.boardConfig, down.cost);
queue.add(down);
}
}
}

if (!queue.isEmpty()) {
currentState = queue.remove();
final byte movesRequired = currentState.cost;
if (movesRequired > previous) {
System.err.println("Generating boards with up to " + previous +
" moves; map size: " + tempMap.size());
previous = movesRequired;
}
} else {
break;
}
}
}

private int indexFor(final long boardConfig, final boolean isFifteenPuzzle,
final boolean[] tilesInSubset, final int[] tilePositions) {
if (isFifteenPuzzle) {
int hashCode = 0;
for (int i = numOfTilesMinusOne; i >= 0; --i) {
final int tile = (int)((boardConfig >> (i << 2)) & 0xF);
if (tilesInSubset[tile]) {
hashCode |= i << (tilePositions[tile] << 2);
}
}
return hashCode;
}
return (int)boardConfig;
}

private boolean[] computeSubset(final byte dummyValue, final long boardConfig) {
final boolean[] tilesInSubset = new boolean[numOfTiles];
for (int pos = numOfTiles - 1; pos >= 0; --pos) {
final byte tile = (byte)((boardConfig >> (pos << 2)) & 0xF);
if (tile != dummyValue && tile != 0) {
tilesInSubset[tile] = true;
}
}
return tilesInSubset;
}

private static byte getArray(final String arrayString, final byte[] tiles,
final int numOfTiles) {

final StringTokenizer st = new StringTokenizer(arrayString, ",");
final int numOfTokens = st.countTokens();

// Validate the number of tiles entered.
if (numOfTokens < numOfTiles) {
System.out.println("Error: Input contains only " + numOfTokens +
" of the " + numOfTiles + " tiles.");
System.exit(1);
} else if (numOfTokens > numOfTiles) {
System.out.println("Error: Input exceeds required " +
numOfTiles + " tiles.");
System.exit(1);
}

// Create an array of String representations of the tile numbers.
final String[] numStrings = new String[numOfTokens];
int i = 0;
while (st.hasMoreTokens()) {
numStrings[i++] = st.nextToken();
}

// Make sure each string is a number.
final int[] tilePositions = new int[numOfTiles];
for (i = 0; i < tiles.length; ++i) {
tiles[i] = -1;
tilePositions[i] = -1;
}
for (i = 0; i < numStrings.length; ++i) {
final String s = numStrings[i];
try {
final byte tile = Byte.parseByte(s);
tiles[i] = tile;
tilePositions[tile] = i;
} catch (final NumberFormatException nfe) {
if (!s.trim().toLowerCase().equals("x")) {
System.err.println("Error: Expected integer or 'X' at index " + (i + 1) +
", received '" + numStrings[i] + "'.");
System.exit(1);
}
}
}

byte dummyTile = -1;
// Make sure no tile number is missing from the input.
for (i = 0; i < numOfTiles; ++i) {
if (tilePositions[i] == -1) {
if (i == 0) {
System.err.println("Error: Tile 0 is missing for input.");
System.exit(1);
}
dummyTile = (byte)i;
break;
}
}

// Replace 'X' (-1) tiles with the dummy value.
for (i = 0; i < tiles.length; ++i) {
if (tiles[i] == -1) {
tiles[i] = dummyTile;
}
}

return dummyTile;
}

public static void main(final String[] args) {
if (args.length < 2 || args.length > 3) {
System.err.println(
"Usage: java PatternDatabaseGenerator <num of tiles> <tile order> [filename]");
System.exit(1);
}
int puzzleSize = 0;
try {
puzzleSize = Integer.parseInt(args[0].trim());
if ((puzzleSize != 8) && (puzzleSize != 15)) {
System.err.println("Error: Puzzle size must be either 8 or 15.");
System.exit(1);
}
++puzzleSize;
} catch (final NumberFormatException nfe) {
System.err.println("Error: Puzzle size must be either 8 or 15.");
System.exit(1);
}

final String tileOrder = args[1].trim();
final byte[] tiles = new byte[puzzleSize];
final byte dummyTile = getArray(tileOrder, tiles, puzzleSize);
System.err.println("Using dummy tile: " + dummyTile);
String filename = null;
if (args.length == 3) {
filename = args[2];
}
new PatternDatabaseGenerator(
puzzleSize, Utility.byteArrayToLong(tiles), dummyTile, filename);
}
}

IDA*:第三阶段由于使用不相交可加性的静态模式数据库,所以不能像第二阶段从目标状态找到初始状态。

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
public static int dfs(int step, int h, int las) {
if (step + h > bound) return step + h; // 从目标开始找 f(n)刚开始最小 如果有更大的则更新 f(n) 反方向找
// if(layer[h]==-1)layer[h]=step;
// if(step>layer[h]+15){
// System.out.printf("enter\n");
// return step+h;
// }
if (h == 0) { // 到达最终状态 输出g(n)即可
// System.out.println(step);
flg = 1;
Rstep = step;
return step;
}
int pos1 = pos;
int ret = 127, x = pos1 / side, y = pos1 % side;
int dx, dy, tar, ht, temp, i;
for (i = 0; i < 4; i++) { //四个方向扩展
dx = x + u[i];
dy = y + p[i];
if (dx < 0 || dy < 0 || dx > side - 1 || dy > side - 1 || !ok(i, las)) continue;
tar = (dx * side) + dy; //计算出扩展出的新节点的一维坐标 2,2 2*4+2= 10
tmp[pos1] = tmp[tar]; // 0的坐标等于扩展出来的点的坐标 a.mat[10]=11;
tmp[tar] = 0;//相当于swap()
pos = tar;
//如果换一个评估函数呢
//阶段三 使用不相交模式数据库得到评估函数 ht值为分为的块的h的和加起来
//

ht = getH(tmp);
//ht = h - (Math.abs(xx[mat[pos1]] - dx) + Math.abs(yy[mat[pos1]] - dy)) + Math.abs(xx[mat[pos1]] - x) + Math.abs(yy[mat[pos1]] - y);
if (step + ht <= bound) {
for (int k = 0; k < side * side; k++) {
result[step][k] = tmp[k];
}
//System.out.println(i);
if(i==0){
redir[step]='r';
}
else if(i==1){
redir[step]='l';
}
else if(i==2){
redir[step]='d';
}
else{
redir[step]='u';
}
}
temp = dfs(step + 1, ht, i);
if (flg == 1) return temp;
if (ret > temp) ret = temp;
tmp[tar] = tmp[pos1];
tmp[pos1] = 0;
pos = pos1;
}
return ret;
}

获取h值:

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
static final int[] tilePositions = {-1, 0, 0, 1, 2, 1, 2, 0, 1, 3, 4, 2, 3, 5, 4, 5};
static final int[] tileSubsets = {-1, 1, 0, 0, 0, 1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2};
public static int getH(int [] tmp) {

int index0 = 0, index1 = 0, index2 = 0;
for (int pos = 15; pos >= 0; --pos) {
final int tile = tmp[pos];
if (tile != 0) {
final int subsetNumber = tileSubsets[tile];
switch (subsetNumber) {
case 2:
index2 |= pos << (tilePositions[tile] << 2);
break;
case 1:
index1 |= pos << (tilePositions[tile] << 2);
break;
default:
index0 |= pos << (tilePositions[tile] << 2);
break;
}
}
}
return PuzzleConfiguration.costTable_15_puzzle_0[index0] +
PuzzleConfiguration.costTable_15_puzzle_1[index1] +
PuzzleConfiguration.costTable_15_puzzle_2[index2];
}

Result:

项目目录: SearchingAstar: E:\University\AI\SearchingAstar3\SearchingAstar

No Initial State Steps Time Limited Time
3 8, 13, 0, 6, 1, 15, 9, 14, 3, 4, 5, 11, 7, 2, 10, 12 52 1.045s
4 2,9,5,11, 8,3,4,14, 7,10,1,12, 0,15,6,13 51 0.070s
5 4,7,0,9,12,10,11,8,14,6,15,1,2,5,3,13 56 0.028s
6 12, 10, 3, 2, 0, 7, 14, 9, 1, 15, 5, 6, 8, 4, 13, 11 57 0.180s
7 12, 1, 5, 6, 2, 11, 7, 9, 14, 10, 0, 4, 15, 3, 13, 8 50 0.009s
8 4, 6, 15, 13, 12, 9, 10, 2, 8, 0, 7, 3, 14, 5, 1, 11 61 0.987s
9 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 1, 2, 0 70 4.758s 60s

Reference resources

(1)Richard E. Korf a,∗, Ariel Felner b ,Disjoint pattern database heuristics

(2)Ariel Felner,Richard E. Korf,Sarit Hanan,Additive Pattern Database Heuristics

https://blog.csdn.net/u013009575/article/details/17140915

http://www.brian-borowski.com/software/puzzle/

http://www.ise.bgu.ac.il/engineering/upload/4273/naij.pdf

http://www.cnblogs.com/beilin/p/5981483.html

https://www.cnblogs.com/fujudge/p/7398153.html?utm_source=itdadao&utm_medium=referral

https://blog.csdn.net/shen_gan/article/details/8146027

https://www.gamedev.net/articles/programming/artificial-intelligence/a-pathfinding-for-beginners-r2003/

空间域图像增强

空间域图像增强

空间域图像增强

基本概念(引用):

参考:https://www.mathworks.com/help/images/ref/imguidedfilter.html

https://blog.csdn.net/wp1603710463/article/details/50408152

https://www.zhihu.com/question/27780598

我一直怀疑我参考的是我某位学长,这跟我们学校图像方面的本科教学意料之中的一致。

​ 图像增强可分成两大类:频率域法和空间域法。前者把图像看成一种二维信号,对其进行基于二维傅里叶变换的信号增强。采用低通滤波(即只让低频信号通过)法,可去掉图中的噪声;采用高通滤波法,则可增强边缘等高频信号,使模糊的图片变得清晰。后者空间域法中具有代表性的算法有局部求平均值法和中值滤波(取局部邻域中的中间像素值)法等,它们可用于去除或减弱噪声。

​ 图像增强的方法是通过一定手段对原图像附加一些信息或变换数据,有选择地突出图像中感兴趣的特征或者抑制(掩盖)图像中某些不需要的特征,使图像与视觉响应特性相匹配。在图像增强过程中,不分析图像降质的原因,处理后的图像不一定逼近原始图像。图像增强技术根据增强处理过程所在的空间不同,可分为基于空域的算法和基于频域的算法两大类。基于空域的算法处理时直接对图像灰度级做运算,基于频域的算法是在图像的某种变换域内对图像的变换系数值进行某种修正,是一种间接增强的算法。

​ 基于空域的算法分为点运算算法和邻域去噪算法。点运算算法即灰度级校正、灰度变换和直方图修正等,目的或使图像成像均匀,或扩大图像动态范围,扩展对比度。邻域增强算法分为图像平滑和锐化两种。平滑一般用于消除图像噪声,但是也容易引起边缘的模糊。常用算法有均值滤波、中值滤波。锐化的目的在于突出物体的边缘轮廓,便于目标识别。常用算法有梯度法、算子、高通滤波、掩模匹配法、统计差值法等。

概念理解:图像增强?怎么才算增强了呢?
从应用角度看,怎么处理以下情况:
(1)图像太暗;
(2)图像太亮;
(3)有噪声点;
(4)对比度不明显。

空间域图像增强的方法:
1)针对每一个个像素处理
① 简单变换(取反,线性变换,指数变换,对数变换,幂次变换);
② 使用滤波器(算子)。
2)针对一组像素处理(直方图)
① 直方图均衡;
② 直方图匹配。
2.常用函数:
⑴ 取反(底片效果):imcomplement()
⑵ 二值化:~im2bw
⑶ 线性变换
⑷ 指数变换: exp()
⑸ 对数变换: log()

⑹ 幂次变换: power()
⑺ 查看直方图: imhist()
⑻ 使用直方图均衡:Histeq()
⑼ 使用平滑滤波器:(imfilter实现函数滤波)
⑽ 使用锐化滤波器
⑾ 旋转、缩放、剪切等: imrotate imresize imcrop

⑿ 对比度增强:imadjust()

实验:

  1. 用Matlab写一段程序,针对提供的图片IMG_2546.JPG,实现:

① (1)查看直方图

② (2)取反,再查看直方图

③ (3)使用直方图均衡,再查看直方图

④ (4) 通过旋转、切割,仅保留“爱丁堡花园”部分

  1. 针对图像100_3228.JPG,使用图像增强的方法使图像的效果好一点,并对比增强前后的直方图变化。

    3.人脸1.jpg、2.jpg、3.jpg、4.jpg进行滤波等操作实现类似美图秀秀磨皮功能,并对比磨皮前后直方图变化。

1、图片IMG_2546.JPG

(1)查看直方图

Code:

1
2
3
4
5
6
7
8
9
10
11
12
clc
clear all
close all
subplot(3,2,1);
%查看直方图
init_Img=imread('E:\University\Digital image\IMG_2546.JPG');
%imshit()直方图的显示
%imhist需要输入一个二维的输入参数,如果输入的图像是一个彩色图像的话,
%不能直接用imhist命令,需要先将图像转成灰度图。
i=rgb2gray(init_Img);
imhist(i);
title('org hist');

直方图:

(2)取反,再查看直方图

Code:

1
2
3
4
5
%取反,查看
subplot(3,2,2);
contray_i = imcomplement(i);
imhist(contray_i);
title('contray hist');

结果:

(3)使用直方图均衡,再查看直方图

Code

1
2
3
4
5
%使用直方图均衡,再查看直方图
subplot(3,2,3);
banlance_i = histeq(i);
imhist(banlance_i);
title('均衡后');

结果:

(4)通过旋转,切割,仅保留“爱丁堡花园”部分

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%通过旋转,切割,仅保留“爱丁堡花园”部分
figure(2);
subplot(1,1,1);
%axis( [xmin xmax ymin ymax] ) 设置当前坐标轴 x轴 和 y轴的限制范围
axis([50,250,50,200]);
imshow(init_Img);
grid on; %显示网格线
axis on; %显示坐标系
title('原图');
%旋转
figure(3);
subplot(1,1,1);
angle_i=imrotate(init_Img,-10,'bilinear','crop');
axis([50,200,50,200]);
imshow(angle_i);
grid on;
axis on;
title('旋转后');
%600 550 900 620 裁剪
figure(4);
subplot(1,1,1);
crop_i = imcrop(angle_i,[600,550,abs(600-900),abs(550-620)]);
imshow(crop_i);
title('裁剪后');

结果:

原图

旋转后:

裁剪后

2、图片100_3228.JPG,使用图像增强的方法使图像效果好一点,并对比增强前后的直方图变化

原图

1、使用中值滤波器 medfilt2

Code:

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
clc
clear all
close all
figure(1);
subplot(1,1,1);
init_Img=imread('E:\University\Digital image\100_3228.JPG');
init1 = imread('E:\University\Digital image\1.jpg');
init2 = imread('E:\University\Digital image\2.jpg');
init3 = imread('E:\University\Digital image\3.jpg');
init4 = imread('E:\University\Digital image\4.jpg');

imshow(init_Img);

figure(2);
subplot(1,1,1);
init_rgb = rgb2gray(init_Img);
imhist(init_rgb);
title('init hist');

%用图像增强的方法使图像的效果好一点,并对比增强前后的直方图变化。

%直方图均衡化实现
figure(3);
subplot(1,1,1);
tmp1=init_Img;
R = tmp1(:,:,1);
G = tmp1(:,:,2);
B = tmp1(:,:,3);
%medfilt2 消除噪声, 中值滤波器, 椒盐噪声
r=medfilt2(R); %medfilt2()中值滤波
g=medfilt2(G);
b=medfilt2(B);
o=histeq(r); %直方图均衡
p=histeq(g);
q=histeq(b);
%cat:用来联结数组
Photo1 = cat(3,o,p,q);
imshow(Photo1,[]);
title('均衡化后的图像');

2、使用高斯滤波器

Code:

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
%使用imfilter滤波
%https://blog.csdn.net/u010740725/article/details/51557202
figure(4);
subplot(1,1,1);
tmp2 = init_Img;
R1 = tmp2(:,:,1);
G1 = tmp2(:,:,2);
B1 = tmp2(:,:,3);
r1 = histeq(R1);
g1 = histeq(G1);
b1 = histeq(B1);
%fspecial() 生成滤波器(也叫算子)的函数
%fspecial(type,para)
%gaussian 高斯滤波器
gaussianFilter = fspecial('gaussian',[7,7],5); %滤波器size:7*7 标准差:5
%imfilter()
%函数语法:g=imfilter(f,w,filtering_mode,boundary_options,size_optinos)
%函数功能:对任意类型数组或多维图像进行滤波
%参数介绍:f是输入图像,w为滤波模板,g为滤波结果;表1-1总结了其他参数的含义。
o1 = imfilter(r1,gaussianFilter, 'symmetric','conv');
p1 = imfilter(g1,gaussianFilter, 'symmetric','conv');
q1 = imfilter(b1,gaussianFilter, 'symmetric','conv');
Photo2 = cat(3,o1,p1,q1);
imshow(Photo2,[]);
title('均衡化后的图像');

比较结果:

1
2
3
4
5
6
7
8
9
%比较效果
figure(5);
subplot(1,3,1);
imshow(init_Img);
subplot(1,3,2);
imshow(Photo1);
subplot(1,3,3);
imshow(Photo2);
title('比较');

1
2
3
4
5
6
7
8
9
10
figure(5);
subplot(1,3,1);
imhist(init_rgb);
subplot(1,3,2);
rgb1 = rgb2gray(Photo1);
imhist(rgb1);
subplot(1,3,3);
rgb2 = rgb2gray(Photo2);
imhist(rgb2);
title('直方图比较');

3、实现类似美图秀秀磨皮功能,并对比磨皮前后直方图变化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
clc
clear all
close all
init_Img=imread('E:\University\Digital image\100_3228.JPG');
init1 = imread('E:\University\Digital image\1.jpg');
init2 = imread('E:\University\Digital image\2.jpg');
init3 = imread('E:\University\Digital image\3.jpg');
init4 = imread('E:\University\Digital image\4.jpg');
figure(16);
LL = double(init4);
HH= double(imguidedfilter(uint8(LL))) -LL +135;
%HH= double(imgaussfilt(uint8(LL),2)) -LL +135;
GG = imfilter(HH,fspecial('gaussian',[3,3],100));
opacity = 50;
Dest = (LL*(100-opacity) + (LL+2*GG-256)*opacity)/100;
imshow([uint8(LL) uint8(Dest)]);
title('数学公式');

最后一张直方图对比:

代码:

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
clc
clear all
close all
subplot(3,2,1);
%查看直方图
init_Img=imread('E:\University\Digital image\IMG_2546.JPG');
%imshit()直方图的显示
%imhist需要输入一个二维的输入参数,如果输入的图像是一个彩色图像的话,
%不能直接用imhist命令,需要先将图像转成灰度图。
i=rgb2gray(init_Img);
imhist(i);
title('org hist');
%取反,查看
subplot(3,2,2);
contray_i = imcomplement(i);
imhist(contray_i);
title('contray hist');

%使用直方图均衡,再查看直方图
subplot(3,2,3);
banlance_i = histeq(i);
imhist(banlance_i);
title('均衡后');

%通过旋转,切割,仅保留“爱丁堡花园”部分
figure(2);
subplot(1,1,1);
%axis( [xmin xmax ymin ymax] ) 设置当前坐标轴 x轴 和 y轴的限制范围
axis([50,250,50,200]);
imshow(init_Img);
grid on; %显示网格线
axis on; %显示坐标系
title('原图');
%旋转
figure(3);
subplot(1,1,1);
angle_i=imrotate(init_Img,-10,'bilinear','crop');
axis([50,200,50,200]);
imshow(angle_i);
grid on;
axis on;
title('旋转后');
%600 550 900 620 裁剪
figure(4);
subplot(1,1,1);
crop_i = imcrop(angle_i,[600,550,abs(600-900),abs(550-620)]);
imshow(crop_i);
title('裁剪后');
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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
clc
clear all
close all
figure(1);
init_Img=imread('E:\University\Digital image\100_3228.JPG');
init1 = imread('E:\University\Digital image\1.jpg');
init2 = imread('E:\University\Digital image\2.jpg');
init3 = imread('E:\University\Digital image\3.jpg');
init4 = imread('E:\University\Digital image\4.jpg');
%人脸1.jpg、2.jpg、3.jpg、4.jpg进行滤波等操作实现类似美图秀秀磨皮功能,并对比磨皮前后直方图变化
%可使用平滑滤波器
%常数128,其实也不一定是个定值,如果把他调大,则处理后的图像整体偏亮,调小则图像偏暗。
% 第五步的图层的不透明度参数也是一个道理,如果不透明度值越大,则图片整体的斑点可能会偏多,
%偏小,那么图像又会过于模糊,也许取个50%是个不错的选择吧,
%或者自己根据处理的纹理图的某个指标做个算法更好吧。
imshow(init_Img);

figure(2);
subplot(1,1,1);
init_rgb = rgb2gray(init_Img);
imhist(init_rgb);
title('init hist');

%用图像增强的方法使图像的效果好一点,并对比增强前后的直方图变化。

%直方图均衡化实现
figure(3);
subplot(1,1,1);
tmp1=init_Img;
R = tmp1(:,:,1);
G = tmp1(:,:,2);
B = tmp1(:,:,3);
%https://blog.csdn.net/cy_543/article/details/41548399
%http://www.ilovematlab.cn/forum.php?mod=viewthread&tid=537123
%medfilt2 消除噪声, 中值滤波器, 椒盐噪声
r=medfilt2(R); %medfilt2()中值滤波
g=medfilt2(G);
b=medfilt2(B);
o=histeq(r); %直方图均衡
p=histeq(g);
q=histeq(b);
%cat:用来联结数组
Photo1 = cat(3,o,p,q);
imshow(Photo1,[]);
title('均衡化后的图像');


%使用imfilter滤波
%https://blog.csdn.net/u010740725/article/details/51557202
figure(4);
subplot(1,1,1);
tmp2 = init_Img;
R1 = tmp2(:,:,1);
G1 = tmp2(:,:,2);
B1 = tmp2(:,:,3);
r1 = histeq(R1);
g1 = histeq(G1);
b1 = histeq(B1);
%fspecial() 生成滤波器(也叫算子)的函数
%fspecial(type,para)
%gaussian 高斯滤波器
%https://blog.csdn.net/chaolei3/article/details/79400658
gaussianFilter = fspecial('gaussian',[7,7],5); %滤波器size:7*7 标准差:5
%imfilter()
%函数语法:g=imfilter(f,w,filtering_mode,boundary_options,size_optinos)
%函数功能:对任意类型数组或多维图像进行滤波
%参数介绍:f是输入图像,w为滤波模板,g为滤波结果;表1-1总结了其他参数的含义。
o1 = imfilter(r1,gaussianFilter, 'symmetric','conv');
p1 = imfilter(g1,gaussianFilter, 'symmetric','conv');
q1 = imfilter(b1,gaussianFilter, 'symmetric','conv');
Photo2 = cat(3,o1,p1,q1);
imshow(Photo2,[]);
title('均衡化后的图像');

%比较效果
figure(5);
subplot(1,3,1);
imshow(init_Img);
subplot(1,3,2);
imshow(Photo1);
subplot(1,3,3);
imshow(Photo2);
title('比较');

figure(5);
subplot(1,3,1);
imhist(init_rgb);
subplot(1,3,2);
rgb1 = rgb2gray(Photo1);
imhist(rgb1);
subplot(1,3,3);
rgb2 = rgb2gray(Photo2);
imhist(rgb2);
title('直方图比较');




figure(16);
LL = double(init4);
HH= double(imguidedfilter(uint8(LL))) -LL +135;
%HH= double(imgaussfilt(uint8(LL),2)) -LL +135;
GG = imfilter(HH,fspecial('gaussian',[3,3],100));
opacity = 50;
Dest = (LL*(100-opacity) + (LL+2*GG-256)*opacity)/100;
imshow([uint8(LL) uint8(Dest)]);
title('数学公式');

figure(17);
subplot(1,2,1);
initrgb = rgb2gray(uint8(LL));
imhist(initrgb);
title('原图');
subplot(1,2,2);
tranrgb = rgb2gray(uint8(Dest));
imhist(tranrgb);
title('处理后');

%利用均值滤波对图像进行平滑处理
%M = rgb2gray(init2);
M = init4;
MR = M(:,:,1);
MG = M(:,:,2);
MB = M(:,:,3);
%添加高斯噪声 均值为0 方差w为0.02
mr = imnoise(MR,'gaussian',0,0.01);
mg = imnoise(MG,'gaussian',0,0.01);
mb = imnoise(MB,'gaussian',0,0.01);
%mr = histeq(mr);
%mg = histeq(mg);
%mb = histeq(mb);
re_Img1 = cat(3,mr,mg,mb);
% ctrl R T 注释操作

M= rgb2gray(init4);
J=imnoise(M,'gaussian',0,0.02);
J = double(J);
H1 = ones(3)/9;
H2 = ones(7)/49;
G1 = conv2(J,H1,'same');
G2 = conv2(J,H2,'same');

%J= rgb2gray(init4);
%均值滤波
%C=conv2(A,B,shape); %卷积滤波
%:输入图像,B:卷积核
%https://blog.csdn.net/jinv5/article/details/52874880
%只能处理灰度图像

% MR = double(MR);
% MG = double(MG);
% MB = double(MB);
% H1 = ones(3)/9;
% H2 = ones(7)/49;
% G1r = conv2(MR,H1,'same');
% G1g = conv2(MG,H1,'same');
% G1b = conv2(MB,H1,'same');
% G1 = cat(3,G1r,G1g,G1b);

% G2r = conv2(MR,H1,'same');
% G2g = conv2(MG,H1,'same');
% G2b = conv2(MB,H1,'same');
% G2 = cat(3,G2r,G2g,G2b);
figure(6);
subplot(2,2,1);
imshow(M);
title('原图像');

subplot(2,2,2);
imshow(re_Img1,[]);
title('高斯噪声');

subplot(2,2,3);
imshow(G1,[]);
title('3*3均值滤波图像');

subplot(2,2,4);
imshow(G2,[]);
title('7*7均值滤波图像');

% figure(7);
% imshow(re_Img1,[]);
%
% figure(8);
% imshow(G1,[]);
%
% figure(9);
% imshow(G2,[]);


%直方图均衡化实现
tmpn=init4;
Rn = tmpn(:,:,1);
Gn = tmpn(:,:,2);
Bn = tmpn(:,:,3);
%https://blog.csdn.net/cy_543/article/details/41548399
%http://www.ilovematlab.cn/forum.php?mod=viewthread&tid=537123
rn=medfilt2(Rn); %medfilt2()中值滤波
gn=medfilt2(Gn);
bn=medfilt2(Bn);
%on=histeq(rn); %直方图均衡
%pn=histeq(gn);
%qn=histeq(bn);
%cat:用来联结数组
Photon = cat(3,rn,gn,bn);
%Photon = cat(3,on,pn,qn);

tmp22 = init4;
R12 = tmp22(:,:,1);
G12 = tmp22(:,:,2);
B12 = tmp22(:,:,3);
r12 = R12;
g12 = G12;
b12 = B12;
%r12 = histeq(R12);
%g12 = histeq(G12);
%b12 = histeq(B12);
%fspecial() 生成滤波器(也叫算子)的函数
%fspecial(type,para)
%gaussian 高斯滤波器
%https://blog.csdn.net/chaolei3/article/details/79400658
gaussianFilter = fspecial('gaussian',[30,30],7); %滤波器size:7*7 标准差:5
%imfilter()
%函数语法:g=imfilter(f,w,filtering_mode,boundary_options,size_optinos)
%函数功能:对任意类型数组或多维图像进行滤波
%参数介绍:f是输入图像,w为滤波模板,g为滤波结果;表1-1总结了其他参数的含义。
o12 = imfilter(r12,gaussianFilter, 'symmetric','conv');
p12 = imfilter(g12,gaussianFilter, 'symmetric','conv');
q12 = imfilter(b12,gaussianFilter, 'symmetric','conv');
Photon2 = cat(3,o12,p12,q12);

figure(10);
subplot(2,2,1);
imshow(init4);
title('原图');

subplot(2,2,2);
imshow(re_Img1);
title('高斯去噪');

subplot(2,2,3);
imshow(Photon,[]);
title('中值滤波');

subplot(2,2,4);
imshow(Photon2,[]);
title('高斯滤波');

figure(15);
%imguide
imguide = init4;
imguided = imguidedfilter(imguide);
imshow(imguided);
title('imguidedfilter');
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
Matlab图像处理函数汇总:
1、图像的变换
① fft2:fft2函数用于数字图像的二维傅立叶变换,如:i=imread('104_8.tif');
j=fft2(i);
②ifft2::ifft2函数用于数字图像的二维傅立叶反变换,如:
i=imread('104_8.tif');
j=fft2(i);
k=ifft2(j);
2、模拟噪声生成函数和预定义滤波器
① imnoise:用于对图像生成模拟噪声,如:
i=imread('104_8.tif');
j=imnoise(i,'gaussian',0,0.02);%模拟高斯噪声
② fspecial:用于产生预定义滤波器,如:
h=fspecial('sobel');%sobel水平边缘增强滤波器
h=fspecial('gaussian');%高斯低通滤波器
h=fspecial('laplacian');%拉普拉斯滤波器
h=fspecial('log');%高斯拉普拉斯(LoG)滤波器
h=fspecial('average');%均值滤波器
2、图像的增强
①直方图:imhist函数用于数字图像的直方图显示,如:
i=imread('104_8.tif');
imhist(i);
②直方图均化:histeq函数用于数字图像的直方图均化,如:
i=imread('104_8.tif');
j=histeq(i);
③对比度调整:imadjust函数用于数字图像的对比度调整,如:i=imread('104_8.tif');
j=imadjust(i,[0.3,0.7],[]);
④对数变换:log函数用于数字图像的对数变换,如:
i=imread('104_8.tif');
j=double(i);
k=log(j);
⑤基于卷积的图像滤波函数:filter2函数用于图像滤波,如:i=imread('104_8.tif');
h=[1,2,1;0,0,0;-1,-2,-1];
j=filter2(h,i);
⑥线性滤波:利用二维卷积conv2滤波, 如:
i=imread('104_8.tif');
h=[1,1,1;1,1,1;1,1,1];
h=h/9;
j=conv2(i,h);
⑦中值滤波:medfilt2函数用于图像的中值滤波,如:
i=imread('104_8.tif');
j=medfilt2(i);
⑧锐化
1)利用Sobel算子锐化图像, 如:
i=imread('104_8.tif');
h=[1,2,1;0,0,0;-1,-2,-1];%Sobel算子
j=filter2(h,i);
2)利用拉氏算子锐化图像, 如:
i=imread('104_8.tif');
j=double(i);
h=[0,1,0;1,-4,1;0,1,0];%拉氏算子
k=conv2(j,h,'same');
m=j-k;
3、图像边缘检测
①sobel算子 如:
i=imread('104_8.tif');
j = edge(i,'sobel',thresh)

②prewitt算子 如:
i=imread('104_8.tif');
j = edge(i,'prewitt',thresh)
③roberts算子 如:
i=imread('104_8.tif');
j = edge(i,'roberts',thresh)
log算子 如:
i=imread('104_8.tif');
j = edge(i,'log',thresh)
⑤canny算子 如:
i=imread('104_8.tif');
j = edge(i,'canny',thresh)
⑥Zero-Cross算子 如:
i=imread('104_8.tif');
j = edge(i,'zerocross',thresh)
4、形态学图像处理
①膨胀:是在二值化图像中“加长”或“变粗”的操作,函数imdilate执行膨胀运算,如:
a=imread('104_7.tif'); %输入二值图像
b=[0 1 0;1 1 1;0 1 0];
c=imdilate(a,b);
②腐蚀:函数imerode执行腐蚀,如:
a=imread('104_7.tif'); %输入二值图像
b=strel('disk',1);
c=imerode(a,b);
③开运算:先腐蚀后膨胀称为开运算,用imopen来实现,如:
a=imread('104_8.tif');
b=strel('square',2);
c=imopen(a,b);
④闭运算:先膨胀后腐蚀称为闭运算,用imclose来实现,如:
a=imread('104_8.tif');
b=strel('square',2);
c=imclose(a,b);

4.12

codeforces

你爱的姑娘,在桥下洗着你最喜欢的衣裳。 ——Wood

CF438D 线段树

传送门: http://codeforces.com/contest/438/problem/D
题解:

线段树维护区间求和,区间取模(更新),单点修改

维护区间最大值剪枝,记录区间最大值。

复杂度:nlogn(很神) 如果X比区间最大值还大,不用取模更新。否则暴力更新。

y mod x,如果x>y/2,那么y取模后小于y/2;

​ 如果x<y/2, 取模后也小于y/2;

所以每个数取模不会超过logn次。

Code:
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

#include<bits/stdc++.h>
#define LL long long
//区间求和 区间取模 单点更新
using namespace std;
const int N = 100010;
struct Tree{
int l;
int r;
LL sum;
LL lazy; //延迟标记
LL maxn;
};
Tree tree[N*4];
int ma[N];
int n,m;
LL max(LL a, LL b){
return a>b?a:b;
}
void up(int x){ //更新区间和
tree[x].sum = tree[x<<1].sum+tree[x<<1|1].sum;
tree[x].maxn = max(tree[x<<1].maxn,tree[x<<1|1].maxn);
}
//延迟标记 用于区间加
void down(int x){
LL t = tree[x].lazy;
if(t){
tree[x].lazy =0;
tree[x<<1].sum+=(tree[x<<1].r-tree[x<<1].l+1)*t;
tree[x<<1].lazy+=t;
tree[x<<1|1].sum+=(tree[x<<1|1].r - tree[x<<1|1].l+1)*t;
tree[x<<1|1].lazy+=t;
}
}

//build
void build(int x,int l,int r){
tree[x].l=l;
tree[x].r=r;
tree[x].lazy=0;
if(l==r){
tree[x].maxn = tree[x].sum=ma[l];
return;
}
else{
int mid = l+( (r-l)>>1);
build(x<<1,l,mid);
build(x<<1|1,mid+1,r);
up(x);
}
}

//单点更新 x k v
void update1(int x,int k,int v){
int ll = tree[x].l;
int rr = tree[x].r;
if(ll==rr){
tree[x].maxn=tree[x].sum=v;
return;
}
int mid = ll+((rr-ll)>>1);
if(k<=mid) update1(x<<1,k,v);
else update1(x<<1|1,k,v);
up(x);

}

//区间加操作 x l r v
void update(int x,int l,int r,int v){
int ll = tree[x].l;
int rr = tree[x].r;
if(l>r){
return;
}
if(ll==l&&rr==r){
tree[x].sum += (tree[x].r - tree[x].l+1)*v;
tree[x].lazy+=v;
}
else{
int mid = ll+((rr-ll)>>1);
down(x);
update(x<<1,l,min(mid,r),v);
update(x<<1|1,max(mid+1,l),r,v);
up(x);
}
}
//区间取模
void updateModify(int x,int l,int r,int v){
int ll = tree[x].l;
int rr = tree[x].r;
if(l>tree[x].r || r<tree[x].l || v>tree[x].maxn) return; //maxn主要用于在这剪枝
if(tree[x].l==tree[x].r){
tree[x].sum = tree[x].maxn = tree[x].sum % v;
return;
}
int mid = ll+((rr-ll)>>1);
if(r<=mid) updateModify(x<<1,l,r,v);
else if(l>mid) updateModify(x<<1|1,l,r,v);
else{
updateModify(x<<1,l,mid,v);
updateModify(x<<1|1,mid+1,r,v);
}
up(x);
}

//区间查询
LL quary(int x,int l,int r){
int ll = tree[x].l;
int rr = tree[x].r;
if(l>r) return 0;
if(ll==l&&rr==r){
return tree[x].sum;
}
else{
int mid = ll + ((rr-ll)>>1);
LL ans=0;
down(x);
ans+=quary(x<<1,l,min(mid,r));
ans+=quary(x<<1|1,max(mid+1,l),r);
up(x);
return ans;
}
}

int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++){
scanf("%d",&ma[i]);
}
build(1,1,n);
while(m--){
int op,a,b;
scanf("%d%d%d",&op,&a,&b);
switch(op){
case 1:
printf("%lld\n",quary(1,a,b));
break;
case 2:
LL v;
scanf("%lld",&v);
updateModify(1,a,b,v);
break;
default:
update1(1,a,b);
break;
}
}
return 0;
}

matlab知识总结

基础知识

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
%读取
imread();
%显示
imshow();
%写入
imwrite();
%转化为灰度图像
rgb2gray();
%裁剪
imresize(Img,[size,size]);
%创建窗口显示
figure();
subplot(cloum,row,num);
title('note');
%操作快捷键
%多行注释 Ctrl+R
%消除注释 Ctrl+T
%matlab按列优先存储
%从1开始

基本函数

disp();
1
2
3
%disp函数:显示文本或数组
%用法说明:disp(X) 函数显示指定的文本或数组。如果参量是数组,则显示数组的内容;如果参量是字符串,则显示字符串文本的内容
disp();
strcat();
1
2
3
4
5
6
7
8
9
先明白strcat函数的定义:
定义
strcat 即 Strings Catenate,横向连接字符串。
语法
combinedStr= strcat(s1, s2, ..., sN)
描述
将数组 s1,s2,...,sN 水平地连接成单个字符串,并保存于变量combinedStr中。如果任一参数是元胞数组,那么结果 combinedStr 是一个元胞数组,否则,combinedStr是一个字符数组。

其实它的作用很简单就是将多个字符连接成单个字符串,关键在于这个语句中括号里面的内容。
dir();
1
2
3
4
5
dir name  (或者 dir(name)) %returns a list of files and folders that match the string name;返回指定路径name所有文件及文件夹组成的列表;

listing = dir(name)    % returns attributes about name;返回指定路径name属性的结构体类型的数组。即数组中每一个元素以结构体类型存储着每个子文件的信息。注意必须赋予变量。
%批量读取照片 返回数组
img_list=dir(strcat(src_path,'*.jpg'));
length();
1
2
3
1length函数:计算向量或矩阵的长度
2、用法说明
y = length(x) 函数计算指定向量或矩阵的长度y。如果参数变量x是向量,则返回其长度;如果参数变量是非空矩阵,则length(x)与max(size(x))等价
sort();
1
2
3
4
5
6
7
    Y = sort(X,DIM,MODE)
    has two optional parameters.  
    DIM selects a dimension along which to sort.%指定对哪个维度进行排序
    MODE selects the direction of the sort. %指定是升序还是降序
       'ascend' results in ascending order
       'descend' results in descending order
    The result is in Y which has the same shape and type as X.%返回结果的shape和类型与输入一致
zero();
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
zeros函数——生成零矩阵

zeros的使用方法】

B=zeros(n):生成n×n全零阵。

B=zeros(m,n):生成m×n全零阵。

B=zeros([m n]):生成m×n全零阵。

B=zeros(d1,d2,d3……):生成d1×d2×d3×……全零阵或数组。

B=zeros([d1 d2 d3……]):生成d1×d2×d3×……全零阵或数组。

B=zeros(size(A)):生成与矩阵A相同大小的全零阵。
imshit() 直方图的显示
imcomplement() 取反
histeq() 使用直方图均衡
imadjust() 图像灰度调整
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
在MATLAB中,通过函数imadjust()进行图像灰度的调整,该函数调用格式如下:

J=imadjust( I )  对图像I进行灰度调整

J=imadjust( I,[low_in;high_in],[low_out;high_out]) [low_in;high_in]为原图像中要变换的灰度范围,[low_out;high_out]为变换后的灰度范围

J=imadjust( I,[low_in;high_in],[low_out;high_out],gamma)  该gamma参数为映射的方式,默认值为1,即线性映射。当gamma不等于1时为非线性映射

RGB2=imadjust(RGB1,......) 该函数对彩色图像的RGB1进行调整

close all;clear all;clc;
%通过imadjust()函数调整灰度图像的灰度范围
I=imread('F:/paohui.jpg');
J=imadjust(I,[0.2 0.5],[0 1]); %调整灰度范围
figure;
subplot(121),imshow(uint8(I));
subplot(122),imshow(uint8(J));
在程序中通过函数imadjust()调整灰度图像的灰度范围。原图像灰度范围为0-255,程序将小于255×0.2的灰度值设置为0,将大于255×0.5的灰度值设置为255。程序运行后输出如下:
uint8()
1
2
3
4
5
matlab中uint16的范围是0-65535,uint8的范围是0-255
matlab提供强制转换为uint8的函数即uint8(number)。

但这个函数的做法是把大于255的数全部强制置为255,而小于255的部分则保持原样不变。
若希望将0-65535的灰度级映射到0-255上,则可采用以下的办法:
mat2gray()
1
2
3
4
5
mat2gray   把一个double类的任意数组转换成值范围在[0,1]的归一化double类数组
mat2gray

功能:将矩阵转化为灰度图像。
用法:I = mat2gray(A,[amin amax]) 把一个double类的任意数组转换成取值范围为[0 1]的亮度图像。其中图像I的取值范围也在0(黑色)到1(白色)之间。参数amin和amax表示将A中小于amin的值转换为0,将A中大于amax的值转换为1。I = mat2gray(A) 将矩阵A中实际最小值和最大值分别赋给amin和amax。
im2uint8()
1
im2uint8,函数名称,把图像数据类型转换为无符号八位整型。如果输入图像是无符号八位整型的,返回的图像和源图像相同。如果源图像不是无符号八位整型的,该函数将返回和源图像相同但数据类型为uint8的图像(必要时对图像进行调整
dftfilt()
1
2
3
4
5
6
7
8
9
10
function g = dftfilt(f,H)
% 此函数可接受输入图像和一个滤波函数,可处理所有的
% 滤波细节并输出经滤波和剪切后的图像
% 将此.m文件保存在一个文件夹
% file->set path->add with subfolder
% 将你函数所在文件夹添加到搜索路径
% save就可以将其添加到你的函数库了
F=fft2(f,size(H,1),size(H,2));
g=real(ifft2(H.*F));
g=g(1:size(f,1),1:size(f,2));

千图成像

原理:

​ 将原图片切成一个一个的小块,用一个图库比对和这张照片的某一块最相似的图片然后替换掉。

以颜色为基准,找颜色。

三种算法:HSV RGB 直方图

RGB

matlab入门写的最丑的代码。仅以此代码祭奠我!

不过这代码清晰的说明了处理的过程。

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
init_path='E:\University\Digital image\3999.jpg';
A=imread(init_path);
initp=imresize(A,[600,600]);
imshow(initp);
%imshow(initp);
%imtool(initP);
out_path='E:\University\Digital image\re2.jpg';
init_all='E:\python\WorkSpace\Image\image\';
img_path_list=dir(strcat(init_all,'*.jpg'));
img_num=length(img_path_list);
%imgtest = imread(strcat(init_all, img_path_list(15).name));
%zoom_M=zeros(10,10,3,img_num);

[X,Y,~]=size(initp);
%50*50
X=X/10;
Y=Y/10;
dstImg=zeros(X,Y,3);
for x=0:X-2
for y=0:Y-2
ref_R = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,1);
ref_G = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,2);
ref_B = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,3);
% disp(ref_R);
im_R=0;
im_G=0;
im_B=0;
for i=1:10
for j=1:10
%uint8的范围是0-255
im_R =( im_R+ double(ref_R(i,j)));
im_G =( im_G+ double(ref_G(i,j)));
im_B =( im_B+ double(ref_B(i,j)));
end
end
im_R=im_R/100;
im_G=im_G/100;
im_B=im_B/100;
ImgRGBmax=255;
for t=1:img_num-1
%Imgtest=imread('E:\python\WorkSpace\Image\image\12_1.jpg');
Img1=imread(strcat(init_all,img_path_list(t).name));
Img=imresize(Img1,[10,10]);
%disp(t)
%imshow(Img);
%Img=zoom_M(:,:,:,t);
% fprintf('%d %s\n',t,strcat(init_all,img_path_list(t).name));
Img_R = Img(:,:,1);
Img_G = Img(:,:,2);
Img_B = Img(:,:,3);
img_R=0;
img_G=0;
img_B=0;
for i=1:10
for j=1:10
%uint8的范围是0-255
img_R =( img_R+ double(Img_R(i,j)));
img_G =( img_G+ double(Img_G(i,j)));
img_B =( img_B+ double(Img_B(i,j)));
end
end
img_R=img_R/100;
img_G=img_G/100;
img_B=img_B/100;
dr = abs(im_R-img_R);
dg = abs(im_G-img_G);
db = abs(im_B-img_B);
dsum = (dr+dg+db)/3;
% fprintf('%d \n',dsum);
if(dsum<ImgRGBmax)
ImgRGBmax=dsum;
tx=t;
end
end
fprintf('%d %d %s\n',ImgRGBmax,tx,strcat(init_all,img_path_list(tx).name));
re_img1=imread(strcat(init_all,img_path_list(tx).name));
re_img = imresize(re_img1,[10,10]);
%imshow(re_img);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,1)=re_img(:,:,1);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,2)=re_img(:,:,2);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,3)=re_img(:,:,3);
end
end
imwrite(dstImg,out_path);

效果图:没有,因为这傻逼代码复杂度。。。

先爬去图片将爬取的图片处理成10*10的小图片。

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
src_path='E:\python\WorkSpace\Image\image\';
dst_path='E:\python\WorkSpace\Image\simage\';
img_list=dir(strcat(src_path,'*.jpg'));
len=length(img_list);
a = 0;
out = 0;
disp(['There are ', num2str(len), ' images in total.'])
for i = 1:len
i_image = strcat(src_path, img_list(i).name);
try
I = imread(i_image);
catch
disp([i_image, ' processing failed.'])
a = a + 1;
continue
end
[X, Y, ~] = size(I);
I = imresize(I, [10, 10]);
[x, y, ~] = size(I);
if(x~=10 || y~=10)
disp([num2str(x), ', ', num2str(y), ' ¡ª¡ª', img_list(i).name])
out = out + 1;
end
try
imwrite(I, strcat(dst_path, [num2str(i), '.jpg']));
catch
disp([i_image, ' image failed in writing.'])
end
end
disp([num2str(a), ' images failed in reading.']);

千图成像:

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
init_path='E:\University\Digital image\3999.jpg';
A=imread(init_path);
initp=A;
imshow(initp);
%imshow(initp);
%imtool(initP);
out_path='E:\University\Digital image\reRGB6.jpg';
init_all='E:\python\WorkSpace\Image\simage\';
img_path_list=dir(strcat(init_all,'*.jpg'));
img_num=length(img_path_list);
%imgtest = imread(strcat(init_all, img_path_list(15).name));
zoom_M=zeros(10,10,3,img_num);
img_R = zeros(500);
img_G = zeros(500);
img_B = zeros(500);
num=1;
%预处理 先算出所有图像的RGB,并分别存入数组中。
for t=1:total
zoom_M(:,:,:,t)=imread(strcat(init_all,img_path_list(t).name));
r=zoom_M(:,:,1,t);
g=zoom_M(:,:,2,t);
b=zoom_M(:,:,3,t);
for i=1:10
for j=1:10

img_R(t) =( img_R(t)+ double(r(i,j)));
img_G(t) =( img_G(t)+ double(g(i,j)));
img_B(t) =( img_B(t)+ double(b(i,j)));
end
end
img_R(t)=img_R(t)/100;
img_G(t)=img_G(t)/100;
img_B(t)=img_B(t)/100;

end

%dstImg_M = initp;
[X, Y, ~] = size(initp);
X=X/10;
Y=Y/10;
for x=0:X-2
for y=0:Y-2
ref_R = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,1);
ref_G = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,2);
ref_B = initp(x*10+1:(x+1)*10,y*10+1:(y+1)*10,3);
% disp(ref_R);
im_R=0;
im_G=0;
im_B=0;
for i=1:10
for j=1:10

im_R =( im_R+ double(ref_R(i,j)));
im_G =( im_G+ double(ref_G(i,j)));
im_B =( im_B+ double(ref_B(i,j)));
end
end
im_R=im_R/100;
im_G=im_G/100;
im_B=im_B/100;
ImgRGBmax=255;

for t=1:img_num-1

dr = abs(im_R-img_R(t));
dg = abs(im_G-img_G(t));
db = abs(im_B-img_B(t));
dsum = (dr+dg+db)/3;
% fprintf('%d \n',dsum);
if(dsum<ImgRGBmax)
ImgRGBmax=dsum;
tx=t;
end
end

fprintf('%d %d %d %s\n',num,ImgRGBmax,tx,strcat(init_all,img_path_list(tx).name));
num=num+1;

re_img=imread(strcat(init_all,img_path_list(tx).name));
%re_img = imresize(re_img1,[10,10]);
%imshow(re_img);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,1)=re_img(:,:,1);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,2)=re_img(:,:,2);
dstImg(x*10+1:(x+1)*10,y*10+1:(y+1)*10,3)=re_img(:,:,3);
end
end
imwrite(dstImg,out_path);

效果图:

原图:

结果:

颜色分布直方图算法:

参考:https://blog.csdn.net/u011600592/article/details/73224512

https://blog.csdn.net/taoyanqi8932/article/details/52758376

http://blog.sina.com.cn/s/blog_4a540be60100vjae.html

直方图介绍

强度直方图图形化显示不同的像素值在不同的强度值上的出现频率,对于灰度图像来说强度

范围为[0~255]之间,对于RGB的彩色图像可以独立显示三种颜色的强度直方图。强度直方

图是用来寻找灰度图像二值化阈值常用而且是有效的手段之一,如果一幅灰度图像的直方图

显示为两个波峰,则二值化阈值应该是这两个波峰之间的某个灰度值。同时强度直方图是调

整图像对比度的重要依据

直方图实现方法:

对一幅灰度图像从上到下,从左到右扫描每个像素值,在每个灰度值上计算像素数目,以这

些数据为基础完成图像直方图的绘制

直方图的性质

直方图反映了图像中的灰度分布规律。它描述每个灰度级具有的像元个数,但不包含这些像元在图像中的位置信息。
任何一幅特定的图像都有唯一的直方图与之对应,但不同的图像可以有相同的直方图。

如果一幅图像有两个不相连的区域组成,并且每个区域的直方图已知,则整幅图像的直方图是该两个区域的直方图之和

基于直方图距离的图像相似度计算:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
%计算图像直方图距离
%巴氏系数计算法

M=imread('1.jpg');
N=imread('2.jpg');
I=rgb2gray(M);
J=rgb2gray(N);

[Count1,x]=imhist(I);
[Count2,x]=imhist(J);
Sum1=sum(Count1);Sum2=sum(Count2);
Sumup = sqrt(Count1.*Count2);
SumDown = sqrt(Sum1*Sum2);
Sumup = sum(Sumup);
figure(1);
subplot(2,2,1);imshow(I);
subplot(2,2,2);imshow(J);
subplot(2,2,3);imhist(I);
subplot(2,2,4);imhist(J);
HistDist=1-sqrt(1-Sumup/SumDown)

弊端分析:

(参考:http://blog.sina.com.cn/s/blog_4a540be60100vjae.html)

(1)直方图匹配。

​ 比如有图像A和图像B,分别计算两幅图像的直方图,HistA,HistB,然后计算两个直方图的归一化相关系数(巴氏距离,直方图相交距离)等等。

​ 这种思想是基于简单的数学上的向量之间的差异来进行图像相似程度的度量,这种方法是目前用的比较多的一种方法,第一,直方图能够很好的归一化,比如通常的256个bin条的。那么两幅分辨率不同的图像可以直接通过计算直方图来计算相似度很方便。而且计算量比较小。

​ 这种方法的缺点:

​ 1、直方图反映的是图像像素灰度值的概率分布,比如灰度值为200的像素有多少个,但是对于这些像素原来的位置在直方图中并没有体现,所以图像的骨架,也就是图像内部到底存在什么样的物体,形状是什么,每一块的灰度分布式什么样的这些在直方图信息中是被省略掉得。那么造成的一个问题就是,比如一个上黑下白的图像和上白下黑的图像其直方图分布是一模一样的,其相似度为100%。

​ 2、两幅图像之间的距离度量,采用的是巴氏距离或者归一化相关系数,这种用分析数学向量的方法去分析图像本身就是一个很不好的办法。

​ 3、就信息量的道理来说,采用一个数值来判断两幅图像的相似程度本身就是一个信息压缩的过程,那么两个256个元素的向量(假定直方图有256个bin条)的距离用一个数值表示那么肯定就会存在不准确性。

Code:

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
%直方图  来自 crc的代码
ref_path = 'E:\University\Digital image\3999.jpg';
% image to be refered
refM = imread(ref_path);
% image to be created
dstImgName = 'E:\University\Digital image\re4.jpg';
% zoomed image file set folder(10x10)
folder = 'E:\python\WorkSpace\Image\simage\';
img_list = dir(strcat(folder, '*.jpg'));
total = length(img_list);
zoom_M = zeros(10, 10, 3, total);

for i = 1:total
zoom_M(:, :, :, i) = imread(strcat(folder, img_list(i).name));
end

dstImg_M = refM;

[X, Y, ~] = size(refM);
X = X / 10;
Y = Y / 10;

for x = 1:X
for y = 1:Y

sim_cof = 0.0;
tx = 1;
for t = 1:total
img_RGB = zoom_M(:, :, :, t);
sim_sum = 0;
for spec = 1:3
r_g_b = img_RGB(:, :, spec);
r_g_b = reshape(r_g_b, [1, 100]);
ref_r_g_b = refM((x-1)*10+1:x*10, (y-1)*10+1:y*10, spec);
ref_r_g_b = reshape(ref_r_g_b, [1, 100]);

[count1, ~] = hist(r_g_b);
[count2, ~] = hist(double(ref_r_g_b));
sum1 = sum(count1);
sum2 = sum(count2);
sumUp = sqrt(count1.*count2);
sumDown = sqrt(sum1*sum2);
sumUp = sum(sumUp);
hist_cof = 1 - sqrt(1 - sumUp / sumDown);
sim_sum = sim_sum + hist_cof;

end

if(sim_cof < sim_sum)
sim_cof = sim_sum;
tx = t;
end
end

dstImg_M((x-1)*10+1:x*10, (y-1)*10+1:y*10, :) = zoom_M(:, :, :, tx);
fprintf('block (%d, %d) located by image %s. histDist: %f\n', x, y, img_list(tx).name, sim_cof);
end
end

imwrite(dstImg_M, dstImgName);

效果图:

考研倒计时

只有经历了一些事,您才会懂得好好珍惜眼前的时光!

3/8 生日(长达几个月的冷战;生日的时候和你分手了 ,弄丢的我还要把她追回来)

3/10 开始了呢

晚上找找时间打打篮球;

打篮球的时候还可以打开水;

少呆会b站就好了;

走路的时候带着耳机听听英语找感觉也行啊;

锻炼减肥;

晚上不想带手机回宿舍了,想着每天晚上做完数学就锻炼晚上回宿舍写博客整理算法吧;

如果有时间偶尔打打比赛就更爽了;

准备去打省赛,退役选手很久都没有写过题了。

减肥,考研,还有那该死的英语。

四五月状态很糟糕,或多或少对身边的朋友有影响。哎呀对不起,你们也不会怪我吧。

我需要光

每日图书馆的数学时光很美好;

Everything the light touches is our kingdom!

英语听力咋就听不懂呜呜呜

计划

一直在做,找老师谈过之后,不留退路。

待完成

数字图像处理 (matlab)

人工智能大作业(java)

编译原理 lex yapp

CSP考试

模板整理

前端学习 (html css js)

校赛

驾照

省赛

纷扰

2019.5.25

感情揭过。一直爱,但放手。说开了。

原因:不想变成自己最讨厌的那类人。

AIProject1

1、Zobrist Hash(By 15-Puzzle problem)

目的:使用Zobrist Hash判断N-Pullze棋局新生成的状态是否存在;

​ 如何从一个状态的Zobrist值得到其后继状态的Zobrist值;

Zobrist Hash思想:

Zobrist hash是一种专门针对棋类游戏而提出来的编码方式, Zobrist哈希通过一种特殊的置换表,也就是对棋盘上每一位置的各个可能状态赋予一个编码索引值,来实现在极低冲突率的前提下在一个整型数据上对棋盘进行编码。

以15-Puzzle为例说明:

题意:4*4棋盘,由15个方块和一个空位组成。15个棋子,每个棋子上标有1至15的某一数字;空格用0表示;空格周围的棋子可以移动到空格中。移动方块,让所有的方块按照数字的顺序排列。(以下讨论将0也算入其中)

一共4x4x16种组合;设定三维数组

1
2
3
Z[4][4][16] //每一个产生一个随机数,至少保证是32位的长度(即32bit),最好是64位。一般我们设定为32位,很多时候为了校验还会再设定一个64位的随机数。

//Z[0][0][0]——Z[0][0][15]的值表示0-15在(0,0)的每个状态的随机数;

对于一个特定的棋局,将每个单位的状态对应的随机数作异或(^)运算,所得值即为Hash值;

对于生成的初始状态,判断新生成的状态是否有与原来的生成的初始状态的Hash值是否一致。若一致,则新生成的初始状态已经存在。

对于判断是否与某历史状态是否存在同理。

如何从一个状态的Zobrist值得到其后继状态的Zobrist值?

当前状态会得到一个Zobrist值a;

比如(0,0)位置上是8,(0,1)位置上是0;

操作:将8移到(0,1);

​ 0到(0,0);

1
2
3
//原来键值里异或过Z[0][0][8]  和 Z[0][1][0];
//再次异或这个值,棋子即从键值中删除;再异或新的状态;所以后继状态的Zobrist为
a1=a^Z[0][0][8]^Z[0][1][0]^Z[0][0][0]^Z[0][1][8];

参考:https://www.cnblogs.com/IThaitian/p/3616507.html

2、N-Puzzle问题的可解性判断

末状态逆序对:0

初始状态:8-Puzzle ;

​ 3x3; 3为奇数;

上下每交换一次逆序数改变偶次数,左右交换不变

初始状态的逆序数与末状态逆序数奇偶性相同则可达;

​ 15-Puzzle;

​ 4x4; 4为偶数;

上下每交换一次逆序数改变奇次数,左右交换不变;

初始状态的逆序数+初始0到最低行的行距(0最后在最低行)与末状态奇偶性相同即可互相到达;

延申:MxN Puzzle问题同理

逆序数求法:归并排序&树状数组

归并排序原理:分治法

图来源:https://www.cnblogs.com/chengxiao/p/6194356.html

合并:

Code:

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
#include<bits/stdc++.h>
#define maxn 100
int a[maxn],b[maxn];
using namespace std;
void merge_sort(int l,int r){
if(r-l>0){
int mid=(l+r)/2;
int i=l;
int p=l;
int q=mid+1;
merge_sort(l,mid);
merge_sort(mid+1,r);
while(p<=mid||q<=r){ //左右两部分只要有一部分不为空
if(q>r||(p<=mid&&a[p]<=a[q])){
b[i++]=a[p++];//q>r 即将左边剩的复制到数组,否则比较
}
else{
b[i++]=a[q++]; //将右半数组复制到辅助数组
}
}
for(i=l;i<=r;i++){
a[i]=b[i]; //拷贝
}
}
}
int main(){
int n;
scanf("%d",&n);
for(int i=1;i<=n;i++){
scanf("%d",&a[i]);
}
merge_sort(1,n);
for(int i=1;i<n;i++){
printf("%d ",a[i]);
}
cout<<a[n]<<endl;

return 0;
}

求逆序对原理:

左右两部分合并时,当右指针数字小于左指针的时候,他会比从左指针开始到中间的元素个数个逆序对。也就是mid-l+1;

所以添加一行代码即可求出逆序对:

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
#include<bits/stdc++.h>
#define maxn 100
int a[maxn],b[maxn];
int cnt;
using namespace std;
void merge_sort(int l,int r){
if(r-l>0){
int mid=(l+r)/2;
int i=l;
int p=l;
int q=mid+1;
merge_sort(l,mid);
merge_sort(mid+1,r);
while(p<=mid||q<=r){ //左右两部分只要有一部分不为空
if(q>r||(p<=mid&&a[p]<=a[q])){
b[i++]=a[p++];//q>r 即将左边剩的复制到数组,否则比较
}
else{
b[i++]=a[q++]; //将右半数组复制到辅助数组
cnt+=mid-p+1;
}
}
for(i=l;i<=r;i++){
a[i]=b[i]; //拷贝
}
}
}
int main(){
int n;
cnt=0;
scanf("%d",&n);
for(int i=1;i<=n;i++){
scanf("%d",&a[i]);
}
merge_sort(1,n);
for(int i=1;i<n;i++){
printf("%d ",a[i]);
}
cout<<a[n]<<endl;
cout<<cnt<<endl;

return 0;
}

树状数组求逆序对:

树状数组原理:

参考https://blog.csdn.net/Small_Orange_glory/article/details/81290634

Code:

1
2
3
4
5
//单点更新
void update(int x,int y,int n){
for(int i=x;i<=n;i+=lowbit(i))    //x为更新的位置,y为更新后的数,n为数组最大值
c[i] += y;
}
1
2
3
4
5
6
7
//区间查询
int getsum(int x){
int ans = 0;
for(int i=x;i;i-=lowbit(i))
ans += c[i];
return ans;
}

求逆序对原理:

https://blog.csdn.net/ssimple_y/article/details/53744096

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
#include<bits/stdc++.h>
#define maxn 1007
using namespace std;
int aa[maxn];//离散化后的数组
int c[maxn]; //树状数组
int n;

struct Node
{
int v;
int order;
}a[maxn];

bool cmp(Node a, Node b)
{
return a.v < b.v;
}

int lowbit(int k)
{
return k&(-k); //基本的lowbit函数
}

void update(int t, int value)
{ //即一开始都为0,一个个往上加(+1),
int i;
for (i = t; i <= n; i += lowbit(i))
c[i] += value;
}

int getsum(int t)
{ //即就是求和函数,求前面和多少就是小于它的个数
int i, sum = 0;
for (i = t; i >= 1; i -= lowbit(i))
sum += c[i];
return sum;
}

int main()
{
int i;
while (scanf("%d", &n), n)
{
for (i = 1; i <= n; i++) //离散化
{
scanf("%d", &a[i].v);
a[i].order = i;
}
sort(a + 1, a + n + 1,cmp);//从1到n排序,cmp容易忘
memset(c, 0, sizeof(c));
for (i = 1; i <= n; i++)
aa[a[i].order] = i;
int ans = 0;
for (i = 1; i <= n; i++)
{
update(aa[i], 1);
ans += i - getsum(aa[i]); //减去小于的数即为大于的数即为逆序数
}
printf("%d\n", ans);
}
return 0;
}