普通视图

发现新文章,点击刷新页面。
昨天以前gyro永不抽风!

深入理解 Zyzzyva 协议

2023年2月18日 23:30

Zyzzyva

动机

PBFT 的 $O(n^2)$ 开销太大,而且观察到大多数时候并会经历 Primary Fault 甚至大多数时候不会有 Fault。Zyzzyva 旨在提高没有 Fault 或者没有 Primary Fault 的情况下的协议效率。

Zyzzyva 流程

一些记号的说明

  • Checkpoint, CC, Speculative History 见上图
  • Order-Req 以及 Spec-Response阶段
    • $d = H(m)$,即消息 $m$ 的 Digest
    • $h_n = H(h_{n-1}, d)$,描述截止到 $n$ 的历史记录,即
      $$
      h_n = H(H(H(\cdots (H(d_1), d_2)\cdots), d_{n-2}), d_{n-1})
      $$
    • 注意到光靠 $h_n$ 就能判断历史记录了。正常的 Client 可以靠这个鉴别是否 Diverge。

Normal Case Operation

当有 Fault 存在

View Change

View Change 的大体过程和 PBFT 相同:

  • PBFT: Last Stable Checkpoint + Preprepare Messages
  • Zyzzyva: Last Stable Checkpoint / Committed Certificate + Speculative History

区别在于对待 Speculative History 的处理上。

Zyzzyva 这里要分两种情况处理:

  • 如果没有 Speculative History,那么就是普通的 PBFT View Change,没什么好讨论的
  • 如果有 Speculative History,我们需要判断这个 History Diverge 了还是是属于 Fast Path 的正常 History。判断方法如下:如果在 View Change 的 $2f+1$ 个 Replica 当中,有 $\geq f+1$ 个 Replica 有 CC 之后的这样一条匹配消息,那么就认为这个消息是真实消息,即属于 Fast Path。我们会在之后的章节给予证明。

一些 Intuition

对于 Zyzzyva,简要而言可以理解如下

  • 如果没有 Fault,那就非常快,Client 知道搞定就行了。这样的记录就是 Speculative 的
  • 如果存在 Fault,那也要至少要求 $2f+1$ 个 Reply,这个时候 Client 会给所有 Replica Multicast Committed Certificate。注意:有了 Committed Certificate,之前的所有 Speculative 记录就自动转正了
  • Zyzzyva 的 View Change 和 PBFT 相比也有少许不同。最大的特点就是:如果 Suspect Primary Faulty,不会暂停消息的处理! Replica 会继续处理消息,但是同时广播 I-Hate-The-Primary。只有当 I-Hate-The-Primary 的数量达到 $f+1$ 的时候,才会进入真正的 View Change。相比之下,PBFT 只要 trigger timeout 就立刻罢工,这就是 PBFT 的 In-Dark 问题。

Safety 是如何保证的?

View 内部的 Safety?

注意到 Zyzzyva 没有 Pre-Prepare - Prepare - Commit 三个节点,没有 View 内的 Consensus,完全靠的是 CC 和 Speculative History。所以我们只需要分析这两个部分。

Committed Certificate 的 Safety

Committed Certificate 的 $2f+1$ 个 Replica 和 View Change 时的 $2f+1$ 个 Replica 构成了 $f+1$ 个 Quorum,保证了 View Change 的时候至少有一个 non-faulty Replica。这个和 PBFT 的过程相同,原理也相同。

Speculative History 的 Safety 以及 View Change 时的 $f+1$ 的 Design Choice

我们要论证的大问题是,有没有可能在前一个 View Commit 了一个消息,客户端接受了,但是在新的 View 中不见了?

首先我们可以证明:如果一个 Commit 属于 Fast Path ($3f+1$),那么在 View Change 时,我们一定能在那个 $2f+1$ 的 Quorum 内找到 $f+1$ 个 Matching。

证明:反证法。如果一个 Commit 属于 Fast Path,同时找不到 $f+1$ 个 Matching,那么和至多有 $f$ 个节点是 Fault 矛盾。

接下来我们想要说明:如果一个 Commit 能够在 $2f+1$ 的 Quorum 中找到 $f+1$ Matching,那么他不可能属于 Slow Path(客户端成功执行).

证明:反证法。如果一个 Commit 属于 Slow Path 且客户端成功执行,那么他必定存在于某个 CC 当中,与 Speculative History 矛盾。

给出了上面两种情形,我们接下来说明:如果一个 Commit 能够在 $2f+1$ 的 Quorum 中找到 $f+1$ Matching,那么把他包含进 New View 中时不会 Violate Safety Condition。

证明:根据前面几个结论,我们可以得到这个 Commit 无非只有两种情形:

  1. Commit 属于 Fast Track,成功执行:不会违反 Safety Condition
  2. Commit 属于 Slow Track,Commit Certificate 只有部分(或者没有)传回,而且根据前面的论证,$2f+1$ 的 View Change 节点内不可能有 CC。在这种情况下 View Change 会促成 Commit,即在 New View 继续处理这个消息。如果客户端请求失败了进行重传,那么 Cache 会继续这个过程。
  3. Commit 尚未执行完毕,即 Client 收到 $< 2f+1$ 个 Speculative Response:回顾本章开头的大问题,全然问题无。

如果 Diverge 了怎么办?

注意到 CC 前的内容不可能 Diverge,只有可能 Speculative History 有 Diverge。文章中没有讲解,但是如果有 Diverge 就听从 New View 安排,不对的东西直接 Rollback (Discard),然后处理新的 Order-Request.

Two-Phase Protocol 的 View Change Liveness 问题

问题构造

与 PBFT 的对比

为什么 PBFT 没有这个问题?注意到 PBFT 的 View Change 原则是如果 Timeout 就自动发起 View Change,不会卡住。

解决方案——I Hate The Primary

解决方案就是 Intuition 章节提到的 I-Hate-The-Primary 消息。这里的主要思想是:只有我知道了足够多的人都怀疑 Primary Fault 的情况下,我才会进入 View Change,以防卡住的问题。这里我们可以证明为什么有 $f+1$ 个 I-Hate-The-Primary 就不会再有 liveness 问题。

注意到保证 liveness 的条件是没有 Network Partition。这就意味着 $f+1$ 个 I-Hate-The-Primary 一定可以被这 $f+1$ 个节点互相收到,也就是说会有 $f+1$ 个节点进入 View Change 节点并广播 View Change 消息。因为 Zyzzyva 和 PBFT 相同,在收到 $f+1$ 个其他节点的 View Change 之后,Replica 会立刻进入 View Change。而因为这里没有 Network Partition,所以上述操作必然可以完成,也就是说必然会有 $\geq 2f+1$ 的节点进入 View Change,完成整个操作。

有没有可能进行 View Change Attack?

参考最后一章 Design Choice 的 $f+1$ 设计问题。

Design Choice 问题汇总与 Edge Case 分析

  • Q: 为什么 Client 需要至少 $2f+1$ Matching?
  • Q: 为什么 View Change 至少需要 $2f+1$ Matching?
  • A: 这两个问题 Commited Certificate 的 Safety 章节已证明。
  • Q: 为什么 View Change 时 New Primary 计算,需要 $f+1$ 个 Ordered Request 存在才能将该记录加入 New View,而不像 PBFT 只需要 $1$ 个?
  • A: Speculative History 的 Safety 章节已证明。
  • Q: 为什么 I-Hate-The-Primary 需要 $f+1$ 个?
  • A: 我们在前面的章节已经证明了 $f+1$ 是足够的。这里我们再论证为什么 $f$ 是不行的。原因非常简单,如果需要 $f$ 个就可以进入 View Change 状态的话,假设有 $f$ 个 faulty Replica,全部不停地发 I-Hate-The-Primary,就会不断 trigger view change,liveness 便不能再保证了。

深入理解 PBFT 协议

2023年2月12日 04:08

PBFT - Practical Byzantine Fault Tolerance

论文链接:https://pmg.csail.mit.edu/papers/osdi99.pdf

BFT

BFT (Byzantine Fault Tolerant) 拜占庭容错,指的是在一个分布式系统当中,在有节点可能出现任意行为(包括但不限于崩溃、发送错误数据、出现恶意行为等)的情况下,仍可以正常运行即达成共识。

Wikipedia Link.

PBFT 的重要概念

PBFT 可以保证在一个大小为 $3f+1$ 的分布式网络下,容忍至多 $f$ 个节点的拜占庭错误。

  • Primary (Leader): PBFT 是一个 leader-based 的协议,即在任意时刻,都有一个 leader 主导协议的走向。
  • 视图 View: 在 Primary 不发生改变的一段时间我们称之为 View。注意到:View 的编号是单调递增的。给定一个视图编号 $v$,我们就可以唯一确定这个视图的 Primary 编号 $p = v \bmod |\mathcal R|$,其中 $\mathcal R$ 是所有节点的集合,即所有 Replica 的集合。
  • Replica: 可以理解为服务器节点。
  • 目标: 达成 Consensus。即在每一个 Replica 上,所有的 transaction 内容和顺序都保证一致。
  • View Change: 视图改变。正因为协议由一个 Leader 主导,那在 leader 出现 malicious behaviour 的时候,协议可以保证触发 View Change,即更换 leader,来保证 liveness。具体内容后文会详细说明。

PBFT 的大致流程

Normal Case Operation 即 Primary 不出问题的情况下

  1. [REQUEST] 客户端向 Primary 发送 $\langle Request, o, t, c \rangle _{\sigma_c}$ 请求
  • $o$: Operation
  • $c$: Client
  • $t$: Timestamp - 这是用来保证执行的 exactly-once semantics
  1. [PREPREPARE] Primary 将请求 Multicast 给所有 Replica: $\langle \langle Preprepare, v, n, d \rangle _{\sigma_p}, m \rangle$
  • $v$: view
  • $n$: Primary 给这个请求分配的 Sequence Number 【关于 Sequence Number: 一个正常的 Replica 只会按顺序 (Sequence Number) 处理消息。即若其收到了 127 号消息但始终没有收到 125、126 号消息,127 号消息便一直不会处理。】
  • $d$: $d = D(m)$ 消息 $m$ 的 digest (摘要)
  1. [PREPARE] 每个 Replica $i$ 收到请求之后,如果选择接受请求(原文有更详细的描述),就会向所有 Replica 广播 Prepare: $\langle Prepare, v, n, d, i \rangle _ {\sigma_i}$
  2. [COMMIT] 接下来,如果每个 Replica $i$ 收到了其他 Replica (包括自己在内) 一共 $2f+1$ 个 Prepare 请求,这个时候认为其进入了 Prepared 状态,并广播 Commit 消息: $\langle Commit, v, n, d, i \rangle _ {\sigma_i}$

Checkpoint Protocol

见原文,不赘述。

View Change Protocol

  • 在怀疑 Primary 作恶的情况下,Replica 会广播 View Change 消息: $\langle ViewChange, v+1, n, \mathcal C, \mathcal P, i \rangle _{\sigma_i}$
    • $\mathcal C$ 是用来证明上一个 Checkpoint 的有效性的
    • $\mathcal P$ : 对于所有编号在 $n$ 之前,上一个 Checkpoint 之后的消息 $m$,$\mathcal P_m \in \mathcal P$,其中 $\mathcal P_m$ 包含了 $m$ 的 Preprepare 消息和 $2f+1$ Matching。
    • 注意到之所以包含了这么多消息,主要就是为了防止消息伪造,全部演算一遍数字签名,出现纰漏就直接忽略,拒不执行。
  • 如果一个 Replica 收到了 $f+1$ 个 View Change 请求,那他也会广播 View Change 请求,加速 View Change 进程。
  • New Primary 收到 $2f$ 个 View Change 请求之后(不包括自己),广播 $\langle NewView, v+1, \mathcal V, \mathcal O \rangle _{\sigma_i}$,具体的计算细节原文写的很详细了。

Fault 的分类

见:https://dinhtta.github.io/pbft/

Network Failure Model

Byzantine fault tolerant protocols tolerate at most $f$ Byzantine failures. In a distributed system, however, there are two types of failures: node and network. It is important to distinguish them, especially since they determine guarantees about the protocol’s safety and liveness.

Node failure means a node (or server, peer, entity, etc.) behaves arbitrarily. This captures the strongest adversary model. A network failure refers to network partitions which can last for an unbounded amount of time. It can be quantified as the number of nodes being isolated from the rest of the network, although they can still communicate with each other.

Given this distinction, what are being counted toward $f$? In PBFT

  • $f$ refers to node failures. The protocol guarantees safety upto $f$ node failures. Safety is independent of network failure. That is, even if the network is severely partition, namely more than $f$ nodes are isolated (but across all partitions there are still fewer than $f$ node failures), the protocol is still safe.
  • However, liveness is only guaranteed with fewer than $f$ total failures, i.e. counting both node and network failures. This means at least $2f+1$ nodes must be reachable.

Note that safety being independent of network failures is a strong guarantee, and it is common in most variants of PBFT. Not until Liu et al. recent work, namely XFT, is it relaxed in a way that separates out node and network failures.

为什么达成 Prepared 状态需要 $2f+1$ 个 Matching

Prepared 保证了什么?

要回答上面的大问题首先要理解我们想要 Preprepare + Prepare Phase 保证什么。论文 Section 4.2 已经回答了这个问题,即

$$
prepared(m, v, n, i) \Rightarrow \lnot prepared(m', v, n, j), \forall i,j,m,m', s.t. D(m') \neq D(m)
$$

即我们要保证在同一个 view 内,每一个 sequence number 只对应唯一一个 prepared message。

反例构造

首先我们可以构造一个如果只有 $2f$ 个 matching 就能打破上述原则的反例。

证明

接下来我们尝试用 Quorum 的思想来证明。使用反证法:如果每个节点都有 $2f+1$ 个 matching,但是最后仍然出现了两个 Replica $a \neq b$,使得

$$
prepared(m, v, n, a) \land prepared(m', v, n, b) \land D(m) \neq D(m')
$$

根据题设,有 $2f+1$ 个节点向 $a$ 发送了 $\langle Prepare, v, n, d, i \rangle _{\sigma_i}$,有 $2f+1$ 个节点向 $b$ 发送了 $\langle Prepare, v, n, d', i \rangle _{\sigma_i}$。根据容斥原理,至少有

$$
2 \times (2f + 1) - (3f+1) = f + 1
$$

个节点向 $a, b$ 分别发送了不同的消息。这样的行为毫无疑问是 malicious 的,所以至少存在 $f+1$ 个恶意节点,与至多有 $f$ 个恶意节点的协议假设矛盾。证毕。

Commit Phase 的必要性

Commit Phase 能够保证什么?

首先,Commit Phase 是 View Change 的副产物。后面的章节会详细解释 View Change,在这个阶段我们只需要知道 View Change 用来解决 Primary 作恶时的 Liveness 问题。而 Commit Phase 就是用来解决 View Change 过程中 Safety 的问题,即在 View 切换的前后:

  • 一个 Sequence Number 对应的 Message 必定在所有正常节点之间达成 Consensus
  • 不会出现一部分节点 Commit 了,但是另一部分节点完全不知道这件事,并继续执行下一个 Commit 或者重复使用 Sequence Number 的情况。
  • 更加明确一点就是原文中的表达:

$$
Committed(n, m, v, i) \Rightarrow \lnot Committed(n, m', v', j), \forall D(m) \neq D(m'), v' > v
$$

反例构造

现在我们假设 Commit Phase 不存在,即在节点达成 Prepared 的时候就直接给 Client 发送 Reply。考虑下图的情况:

可以发现反例非常容易构造,我们甚至构造出了 $f$ 个节点使用了 $n+1$ 的 Sequence Number,但是剩余节点以及 New Primary 没有的情况。

理解

从上面这个例子已经可以看出,Commit Phase 主要就是需要保障一个 Sequence Number 在跨 View 的时候不能被 Malicious Primary 或者 Network Partition 影响,使得被重复使用,破坏 Safety。

增加一个 Commit Phase 的主要功能,就是 只有在知道有一定数量的其他 Replica 已经处于 Prepared 状态之后再进行 Commit,这个数量一定要保证能够在 View Change 的时候至少有一个节点能够证明这个请求已经 Prepared。

因为注意到上图:在节点 Commit 的时候,很多节点没有达到 Prepared 状态,无法在 Network Partition 的时候提供证据向新的 Primary 证明这个请求已经被 Commit。综上所述,这就是 Commit Phase 存在的意义。

为什么达成 Committed 状态需要 $2f+1$ 个 Prepared

估算

上一节提到,在 Commit 前,我们至少要保证有足够多的 Replica 已经处于 Prepared 阶段,使得在 Network Partition 出现的时候仍然能够向 Primary 提供信息,告诉 Primary 一个请求已经被 Commit。那这个数量到底应该是多少呢?通过简单的估算其实就可以得到。

已知 New Primary 需要 $2f$ 个除了自己以外的 View Change 请求。最悲观地考虑,出现了最大限度的 Network Partition,即有 $f$ 个节点无法通信。在这个极限情况下的极限情况是这 $f$ 个 Replica 都处于 Prepared 状态却无法与 Primary 通信。而其余 $2f+1$ 个节点中又有 $f$ 个 malicious node,可以假装不知道自己 Prepared 了。

所以,在那 $2f+1$ 个节点中至少需要 $f+1$ 个节点处于 Prepared 状态才能保证 New Primary 收到 Commit 记录,再算上 Network Partition 的 $f$ 个 Replica,一共得有 $2f+1$ 个节点处于 Prepared 状态才能保证信息能够在 Network Partition 和 View Change 的情况下送达 New Primary。

反例构造

上述文字选取的是极端情况,有一条不成立即为反例。所以 $2f$ 个 Prepared 却仍然失败的反例和没有 Commit Phase 一样非常好构造,这里不再赘述。

为什么 New Primary 需要 $2f$ 个 View Change

理解了 Commit Phase 的 $2f+1$ 便不难发现,这里 $2f$ 这个数目的选取是和 Committed 数目的选取是相关的。

两个 $2f+1$ (实际上 View Change 算上 New Primary 也是 $2f+1$) 在总共 $3f+1$ 的池子里,Intersection (Quorum) 的大小总是不小于 $f+1$ 个 Replica 的。这样就保证了总有一个 Non-faulty Replica 位于最终的 Quorum 当中(faulty replica 可以隐瞒不报,装作不知道,所以一定需要 non-faulty replica),即既是 Prepared,又没有被 Network Partition 和 Malicious Replica 影响到。

为什么 Reply 需要 $f+1$ 个 Matching

这个就非常 Intuitive 了,因为这样必定能有一个 Non-faulty Replica 回复。

Checkpoint (Garbage Collection) 的必要性

注意到 View Change 的时候需要附上所有的记录来证明自己是正常节点,这所带来的开销随着时间的推移显然是不可接受的,所以我们需要想办法回收一部分 log。

Checkpoint 为什么在数字签名模式下只需要 $f+1$?

因为 $f+1$ 保证了至少有一个 Non-faulty Node 进行了 Checkpoint,也就是在没有 Network Partition 的情况下必定可以从网络中访问到这个 Checkpoint。

但是这里其实有一个问题:如果 Network Partition 了呢?注意到 PBFT 只在没有 Network Fault 的情况下 Guarantee Liveness,所以在那种情况下,需要请求 Checkpoint 的节点就 In-dark 了。如果这样的情况足够多,Liveness 也就无法保证了。

View Change 和 Liveness 之间的关系

View Change 主要是为了解决没有 Network Fault 的情况下的 Liveness 问题。当一个 Replica 的请求迟迟无法完成的时候就会 Trigger View Change。因为请求无法完成无非几种情况:

  • Primary 作恶,比方说只发给一部分 Replica Preprepare 消息
  • Network Partition

View Change 在第二种情况下其实并不能解决问题,但是可以解决第一种情况的问题。这就是 View Change 的 Intuition。

Primary 作恶带来的 In-Dark 问题

注意到 Primary 可以故意不给某些节点发请求,即 Primary 作恶。在这种情况下,这些节点既不会触发 View Change,也不会做任何事情,也就是字面意义上的 In-Dark。

由于前面的论证,PBFT 已经保证了 Safety,同时在没有 Network fault 的情况下也保证了 Liveness,所以这样的 In-Dark 并不会有什么特别大的影响。

只有当 In-Dark Replica 的数量多到使得 Consensus 无法达成了,才会 trigger View Change,使得协议能够推进下去。

Network Recovery 的流程

  • 首先向网络请求 Checkpoint
  • 如果能顺利 Catch Up,那就无缝衔接
  • 反之,迎来的结果必然是请求积压无法处理,timeout 之后便会 Request View Change,然后在 View Change 的过程中完成 Catch Up

Misc

  • Q: 如果 Primary 无中生有虚构信息或者篡改 Client 信息怎么办?
  • A: 注意到 $\langle Request, o, t, c \rangle _{\sigma_c}$ 的角标 $\sigma_c$,客户端签名保证了 Primary 不能篡改消息也不能伪造消息。
  • Q: 如果 Client 不知道当前的 Primary 怎么办?
  • A: 如果客户端一段时间内没有收到 Reply,触发 Timeout,重新 Multicast 请求。
  • Q: 如果因为 Network Partition 正常节点没有把回复送达到 Client,Client 认为请求失败,重传怎么办?这不就等于会执行两次了吗?
  • A: 注意到 $\langle Request, o, t, c \rangle _{\sigma_c}$ 中的 timestamp 可以解决这个问题

其他小的点:

  • View Change 后面 New View 这里应该也是有 Time out 机制的,否则 New Primary 迟迟不发 New View 的话 Liveness 就寄了
  • Client 发请求之后应该也有 Timeout 机制然后 Multicast,用来规避:(1)Primary 作恶,故意不发 Preprocess(2)发送的 Replica 不是 Primary

西瓜书 绪论 习题

2022年5月17日 19:44

注:本文为《机器学习》(周志华)第一章的作业笔记。

1.1

表 1.1 的编号为 1, 4 的两个样例:

色泽 根蒂 敲声 好瓜
青绿 蜷缩 浊响
乌黑 稍蜷 沉闷

因为只有两个样本,我们只要能够保证:

$$
s_好 \in S, s_坏 \notin S
$$

1.2

首先我们可以将决策看作一棵树,而每一个取式其实就是一棵子树或者叶子节点。会过来看合取式的析合范式,其实就是在一整棵树中取 $k$ 个互相不包含的子树或者叶子节点。

这里不能直接确定答案,因为题目并没有说决策树能不能给出形如下面的节点:

$$
(色泽\in\{青绿,乌黑\}; 根蒂=*; 敲声=*)
$$

我们先假设题目不允许上述的情况出现,那我们可以断言,包含 $n$ 个并列条件的决策,整棵树的深度为 $n + 1$,因为每一个父子关系其实就是析构 $*$ 和具体内容。

这个问题其实不能简单估算,但是我们可以用树形动态规划来解决。设状态 $f[u][p]$ 代表以节点 $u$ 为根的子决策树有 $p$ 个合取式的析合范式情况数。同时我们注意到书上的旁注:会有冗余的等价情况,那么我们还需要去除这样的等价情况,即在本题中我们需要去除

$$
\bigcup_{v \in u.sons} v = u
$$

的情况。

下面是递推方程:

$$
f[u][p] = \sum_{\sum_{1 \leq i \leq |u.sons|} p_i = p, p_i \geq 0} f[v_i][p_i]
$$

减去冗余情况:

$$
f[u][|u.sons|] := f[u][|u.sons|] - 1
$$

1.3

举例:线性回归。

1.4

$$
E_{ote} (\mathfrak L_a | X, f) = \sum_h \sum_{\boldsymbol x \in \mathcal X - X} P(\boldsymbol x) \ell (h(\boldsymbol x), f(\boldsymbol x)) P(h | X, \mathfrak L_a)
$$

$$
\begin{aligned}
\sum_f E_{ote} &= \sum_f \sum_h \sum_{\boldsymbol x \in \mathcal X - X} P(\boldsymbol x) \ell (h(\boldsymbol x), f(\boldsymbol x)) P(h | X, \mathfrak L_a) \\
&=\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x)\sum_h P(h | X, \mathfrak L_a) \sum_f \ell (h(\boldsymbol x), f(\boldsymbol x))
\end{aligned}
$$

注意到最后一部分为与 $\mathfrak L$ 无关的常数。不妨设

$$
L = \sum_f \ell (h(\boldsymbol x), f(\boldsymbol x))
$$

故原式等于

$$
L\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x)\sum_h P(h | X, \mathfrak L_a) = L\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x) \cdot 1
$$

Q.E.D.

1.5

略。

PyTorch Data Augmentation 数据增广

2022年5月2日 23:15

Data Augmentation

这个词其实很容易误导人,因为 Augmentation 是增广的意思。数据增广,自然而然就会让人理解为:“增加了新的数据”。然而在实现的过程中,不一定是这样。

我们会给定一些随机化的参数,在每个 epoch 内生成随机增广后的图片。

比方说,每个 epoch 都按照一定的 random rule 来 crop 图片、调整色调,便是图像增广。

PyTorch 实现

这里拿 Fashion-MNIST 举例,同样我们会给出例子来证明我们的做法是行之有效的。

首先引入包:

%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
from matplotlib import pyplot as plt
import torch
import torchvision
from torch.utils import data

加载 FashionMNIST 数据集,并且禁用 shuffle,这样保证我们每个 epoch 图片顺序都是相同的,方便我们验证:

def load_data_fashion_mnist(batch_size, dataloader_worker_count, resize=None):
    trans = [torchvision.transforms.ToTensor()]
    if resize:
        trans.insert(0, torchvision.transforms.Resize(resize))
    trans.append(torchvision.transforms.RandomCrop((28, 28)))
    trans = torchvision.transforms.Compose(trans)
    mnist_train = torchvision.datasets.FashionMNIST(
        root="../data", train=True, transform=trans, download=True)
    mnist_test = torchvision.datasets.FashionMNIST(
        root="../data", train=False, transform=trans, download=True)
    return (data.DataLoader(mnist_train, batch_size, shuffle=False, num_workers=dataloader_worker_count),
            data.DataLoader(mnist_test, batch_size, shuffle=False, num_workers=dataloader_worker_count))

辅助渲染的函数:

def get_fashion_mnist_labels(labels):
    '''
    返回Fashion-MNIST数据集的文本标签
    '''
    text_labels = ['t-shirt', 'trouser', 'pullover', 'dress', 'coat', 'sandal', 'shirt', 'sneaker', 'bag', 'ankle boot']
    return [text_labels[int(i)] for i in labels]

def show_images(imgs, num_rows, num_cols, titles=None, scale=1.2):
    '''
    绘制图像列表
    '''
    figsize = (num_cols * scale, num_rows * scale)
    _, axes = plt.subplots(num_rows, num_cols, figsize=figsize)
    axes = axes.flatten()
    for i, (ax, img) in enumerate(zip(axes, imgs)):
        if torch.is_tensor(img):
            # 图片张量
            ax.imshow(img.numpy())
        else:
            # PIL图片
            ax.imshow(img)
        ax.axes.get_xaxis().set_visible(False)
        ax.axes.get_yaxis().set_visible(False)
        if titles:
            ax.set_title(titles[i])
    return axes

测试三个 epoch,看每次 epoch 第一个 batch 的数据:

train_iter, test_iter = load_data_fashion_mnist(18, 4, resize=40)

for _ in range(3):
    for i, (X, y) in enumerate(train_iter):
        if i == 0:
            show_images(X.reshape(18, 28, 28), 2, 9, titles=get_fashion_mnist_labels(y))

输出:

可以看出有明显的不同。

本文中的 notebook:https://github.com/JeffersonQin/deep-learning/blob/master/tests/pytorch-data-augmentation.ipynb

P3572 [POI2014] PTA-Little Bird - DP 单调队列

2022年4月8日 23:50

题面

https://www.luogu.com.cn/problem/P3572

题解

首先,不难列出方程:$f[i]$ 代表跳到 $i$ 的最小疲劳值

$$
f[i] = \min_{i - k \leq j < i} \{f[j] + (a[i] \geq a[j])\}
$$

时间复杂度为 $O(qn^2)$,不理想。但是不难发现,这道题可以用单调队列进行优化。为什么?

如下状态一定是更优的:随着 index 递增

  • $f$ 不下降(第一维比较)
  • $f$ 若相同,同时有 $a$ 不上升(第二维比较)

所以我们就可以用单调队列来维护 index,时间复杂度为 $O(qn)$

注意:

  • 初值的设定,详见代码
  • 写的顺序和单调队列的模板稍微有一点不同,因为 $i - k \leq j < i$,右边的小于号不是小于等于号,即 $f[i]$ 不能被自己更新,所以调整了一下逻辑块的顺序。

代码

#include <iostream>

using namespace std;

const int maxn = 1e6+5;

int qq, n, k, head, tail, a[maxn], f[maxn], q[maxn];

int main() {
    cin >> n;
    for (int i = 1; i <= n; i ++) cin >> a[i];
    cin >> qq;
    while (qq --> 0) {
        cin >> k;
        head = 1; tail = 1; q[1] = 1;
        for (int i = 2; i <= n; i ++) {
            while (head <= tail && q[head] < i - k) head ++;
            f[i] = f[q[head]] + (a[i] >= a[q[head]]);
            while (head <= tail && (
                f[i] < f[q[tail]] || (
                    f[i] == f[q[tail]] && a[i] >= a[q[tail]]
                )
            )) tail --;
            q[++ tail] = i;
        }
        cout << f[n] << endl;
    }
    return 0;
}

POJ 2823 滑动窗口 单调队列 - 致逝去的青春

2022年4月8日 23:03

Preface

三年前的四月,第一次接触 OI。初三下,小教室里听着 fcx 讲着滑动窗口单调队列,放学后去本部蹭了第一节斜率优化的课,随即拉开了我长达两年的竞赛摆烂之旅。

记忆犹新。第一个注册的 OJ 账户就是 POJ,第一道 AC 的题便是 POJ 2823 滑动窗口。

多年过去,今天再看这个基础中的基础,茫然。只能说这是究极摆烂的下场罢。总之就好好写一下题解。

题面

https://www.luogu.com.cn/problem/P1886

题解

接下来只讨论求最小值的情况。设 $f[i]$ 代表终止于位置 $i$ 的滑动窗口的最小值:

$$
f[i] = \min_{i - k < j \leq i} a[j]
$$

如果暴力做的话时间复杂度为 $O(n^2)$。这里我们需要使用单调队列来优化。

对于最小值:单调队列内的元素值递增,同时 index 也是递增的。为什么这样做?思考一个情况:index 递增的情况下,如果一个 $a[i] > a[j], i < j$,那么在 $i,j$ 都有效的范围内,$a[i]$ 永远不可能作为答案。

这里有一句经典的话:如果一个选手年龄比你小,但却比你强,那你就可以退役了。我们把 index 比作出生年份,$a[]$ 比作能力,那么这句话就维护了一个求最大值的单调队列。

接下来是代码中需要注意的一些问题:

  • head, tail 不是左闭右开,而是 $[head, tail]$。这样做的原因是访问 headtail 的元素就可以简单的:q[head], q[tail]
  • 队列里存的是 index (pos)

注意:第二个 while 可以修改为 if,因为每次 i ++q 内的 pos 递增,写成 while 是习惯。

代码

#include <iostream>

using namespace std;

const int maxn = 1e6+5;
int n, k, a[maxn], q[maxn], head, tail;

int main() {
    cin >> n >> k;
    for (int i = 1; i <= n; i ++) cin >> a[i];
    // min
    head = 1, tail = 0;
    for (int i = 1; i <= n; i ++) {
        while (head <= tail && a[i] <= a[q[tail]]) tail --;
        q[++ tail] = i;
        while (head < tail && q[head] <= i - k) head ++;
        if (i >= k) cout << a[q[head]] << " ";
    } cout << endl;
    // max
    head = 1, tail = 0;
    for (int i = 1; i <= n; i ++) {
        while (head <= tail && a[i] >= a[q[tail]]) tail --;
        q[++ tail] = i;
        while (head < tail && q[head] <= i - k) head ++;
        if (i >= k) cout << a[q[head]] << " ";
    } cout << endl;
    return 0;
}

VSCode绕开腾讯云COS防盗链 | Markdown魔改过程

2021年1月30日 12:05

最新更新

现在已经将其改成了blogging-tool,增加了一些主题特定支持的snippets和预览功能,但是fake-referrer仍然是核心功能之一。

项目页:https://github.com/JeffersonQin/vscode-blogging-tool

Preface

我还记得一个月前?半个月前?反正就是刚放寒假那会,给一个腾讯云COS的插件改了bug,修复了Windows下的路劲问题,现在才姑且能用。我的想法其实是想把腾讯云的COS作为图床(毕竟费用还可以)(确信)。在设置了防盗链之后,我发现VSCode是无法正常渲染图片的(毕竟就算是基于Chromium也怎么可能会设置referrer嘛...),所以就想着魔改一把市面上已有的插件来实现此功能。

所以最终魔改后的定位:可以图片请求伪造referrermarkdown preview插件

锁定魔改的插件

第一步便是要锁定需要魔改的插件,但是整个过程十分曲折迂回。一开始安装的插件是Markdown All in One: https://marketplace.visualstudio.com/items?itemName=yzhang.markdown-all-in-one

但是读代码的时候一直觉得很奇怪,那就是即便有个文件叫做preview.ts,仍然找不到WebView的相关代码。一度以为可能和html有关系,便找到了这里:

最后却发现这个只是将文件导出成html(大悲:

仔细一想,我禁用了插件,但却发现VS Code仍然可以进行Markdown的preview,这才发现原来是VS Code本身的功能(我当时直接怒摔键盘(划掉。

所以在这之后,锁定了现在的这个插件:Markdown Preview Enhancedhttps://marketplace.visualstudio.com/items?itemName=shd101wyy.markdown-preview-enhanced

Clone代码,运行

使用GitHub Desktop,Clone到了桌面(过程中开着ssr否则速度吃不消),首先尝试直接通过命令行移动:

$ wsl
$ mv vscode-markdown-preview-enhanced/ ~

但是发现:

mv: setting attribute 'system.wsl_case_sensitive' for 'system.wsl_case_sensitive': Invalid argument

好像是和wsl大小写不敏感有关,于是乎不搞这些花里胡哨的,准备通过UI界面直接移动文件夹。进入wsl子系统的用户目录:

explorer .

在Windows文件资源管理器内打开这个目录,然后直接拖动文件夹,移动完毕。

用VSCode直接打开这个文件夹:

code .

发现由于在linux环境内,换行符的问题导致所有文件git管理都显示modified,为了消除不必要的困扰以及为了以后一目了然知道哪里是修改过的代码,先commit一次:

git commit -a -m "Linux Update"

首先尝试npm install,第一次尝试使用淘宝源:

cnpm install

表面上好像并没有什么问题,但是下一步进行编译(因为项目使用的是TypeScript):

npm run compile

发现报错:

初步判断可能是因为包老,所以执行了一下Update:

cnpm update

权限不够,报错,所以加上sudo:

sudo cnpm update

仍然报错,版本冲突:

经搜索引擎搜索后无果,遂放弃。

直接选择删除outnode_modules文件夹,重新使用npm安装依赖,不使用cnpm

proxychains4 npm install

由于速度太慢,所以使用了proxychains4,有需要的可以看之前的文章。

可以看到安装成功,但不知道问什么出现了淘宝的地址,大概是装了cnpm的原因吧(我也不太懂,毕竟只是为了改个功能并非深挖nodejs,如果有了解的欢迎评论)

开始编译:

编译没有问题,所以开始调试:

虽然是在WSL内,但是还是成功地跑起来了,这里感叹一下巨硬牛逼(!

搞清大致的结构

快速浏览了一下源码,开始思考切入点:

  • 什么时候Preview会更新
  • 更新是如何触发的
  • 更新的代码写在哪里
  • 如何拦截更新

凭着之前写插件的经验,找到了extension.ts内:

测试了一下,Console确实有输出:

看到后面的一句:

contentProvider.update(event.document.uri);

contentProvider在之前定义了,是MarkdownPreviewEnhancedView的实例,暂且不去管它,event.document.uri估计是当前正在编辑的窗口的估计指针的类似的东西(纯属口胡,如有误解欢迎指正)

按下Control键单击update定位到函数,并增加两条调试语句:

调试后发现:

  • 在没有开启preview的情况下指输出here 1
  • 开启preview的情况同时输出here 2

得到第一个if只是排除了没有预览窗口的情况的结论。开始阅读第二个if的代码。

大概可以推测出thiswaiting的含义,估计只是在编辑之后有一段延时再进行渲染,所以主要着眼于:

this.updateMarkdown(sourceUri);

毕竟怎么看都是这句是关键。

继续如法炮制,通过打log证实了我的猜想:

看回代码:

第一个和第二个if估计只是业务逻辑上的一些判空,只有后面做这给的注释presentation mode我也不知道是啥,所以我姑且就没有理他,通过log的输出定位到了here 5的地方。现在看之后的代码:

// not presentation mode
vscode.workspace.openTextDocument(sourceUri).then((document) => {
    const text = document.getText();
    this.previewPostMessage(sourceUri, {
        command: "startParsingMarkdown",
    });

    const preview = this.getPreview(sourceUri);

    engine
        .parseMD(text, {
            isForPreview: true,
            useRelativeFilePath: false,
            hideFrontMatter: false,
            triggeredBySave,
            vscodePreviewPanel: preview,
        })
        .then(({ markdown, html, tocHTML, JSAndCssFiles, yamlConfig }) => {
            // check JSAndCssFiles
            if (
                JSON.stringify(JSAndCssFiles) !==
                    JSON.stringify(this.jsAndCssFilesMaps[sourceUri.fsPath]) ||
                yamlConfig["isPresentationMode"]
            ) {
                this.jsAndCssFilesMaps[sourceUri.fsPath] = JSAndCssFiles;
                // restart iframe
                this.refreshPreview(sourceUri);
            } else {
                this.previewPostMessage(sourceUri, {
                    command: "updateHTML",
                    html,
                    tocHTML,
                    totalLineCount: document.lineCount,
                    sourceUri: sourceUri.toString(),
                    id: yamlConfig.id || "",
                    class: yamlConfig.class || "",
                });
            }
        });
});

看到前两行的textpreview,我也不知道这是在干什么,所以我们需要大胆猜测。胆大心细做好备份,正是API猜测大法的方法论,世界观和方法论是统一的(dbq 胡言乱语ing)。但是通过作者的变量命名,我们发现:startParsingMarkdownstart这个词好像别有深意!所以我们大胆猜测,这大概就是一个signal,和我们要干的事情可能关系不大。为了初步印证我的猜想,我们不妨查看一下previewPostMessage到底写了啥:

里面关键调用了VS CodeWebViewpostMessage方法。扯一句题外话,关于WebView,大部分VSC插件都是通过WebView来实现界面的。可以通过Shift + Ctrl (Command) + P输入WebView来呼出调试器。想要确定一个Panel是不是WebView也很简单:呼得出就是,呼不出就不是。

好,所以好像并没有啥关系,我们不妨输出一下textpreview

所以text就是文档内容,preview倒没看出来啥,啥都不是?(逃

继续看engine那行代码,一开始的.parseMD好像并不是我们要找的,只是在解析MD构建语法树?后面的.then()貌似是在干正事。虽然目测估计正常情况下是执行else的代码,但还是打一下log确定一下执行的位置:

好,我们离成功只差最后一步了(!因为只有一行代码了(确信

this.previewPostMessage(sourceUri, {
    command: "updateHTML",
    html,
    tocHTML,
    totalLineCount: document.lineCount,
    sourceUri: sourceUri.toString(),
    id: yamlConfig.id || "",
    class: yamlConfig.class || "",
});

我现在非常的迷惑,看到这么一堆参数,但我觉得重要的只有两个:htmltocHTML,盲猜后面的是目录的html?いみふめいです。好,输出不就行了。(顺便整了点花活)

测试文件:

Console输出:

tocHTML的输出也应证了我的目录猜想。但是到这里还没完,为了保险起见,我们还需要打开WebView的调试器进行进一步的确认:

和先前的结果确实吻合,为了确认这里就是html的控制语句,我们不妨把htmltocHTML都替换成空字符串,看一下渲染效果:

内容确实消失了。

明确魔改需求

在搞清楚代码结构以及魔改方法之后,我们就需要明确我们的魔改需求:可以将腾讯云COS中的图片显示出来(之前显示不出来是因为腾讯云我开了防盗链,禁止没有Referrer的请求,而VSCode对图片的请求中并没有Referrer),所以初步的想法有两种:

  • 第一种:给请求加上Referrer,或者改变Refer-Policy
  • 第二种:发现有这个网站的图片直接下载到本地,并在渲染前更改渲染的html换为本地的。

对于第一种情况:我们不妨直接在编写markdown的时候使用img标签,里面加入referrerpolicy="origin"的属性。但是很遗憾,在WebView调试其中我发现,Referrer仍然为空,即便有了这个属性:

注:这里图片请求成功了是因为我暂时关掉了防盗链。

之后,思考到会不会是<meta>标签内有这样的信息,但是在查看了主程序的htmlWebViewhtml后都无果,遂放弃了这条思路:

对于第二种方法,我们可以绕开浏览器,通过本地的网络请求来进行实现(powershell, curl, etc.)

开始魔改

首先明确这个插件将图片渲染成html的格式:

发现格式为:

<img src="<link_addr>" alt="<alt_name>">

本来想着直接替换的,但是转念一想想到了更好的处理方法:使用现成的库来解析html不就得了。

使用搜索引擎搜索了一番之后,我锁定了cheerio这个包,正准备安装的时候,看了一下node_modules里发现作者貌似已经用了:

于是我们现在是要筛选出<img>标签中,src属性以特定字符串开头的:

import * as cheerio from "cheerio";

...

const $ = cheerio.load(html);
$("img").each((_index, element) => {
    if (element.attribs["src"].startsWith("<prefix>")) {
        element.attribs["src"] = "需要改成的东西";
    }
});

html = $.html();

反正这里还有非常多的东西需要修改,但是已经初具雏形。接下来要实现的非常重要的功能便是下载功能。由于主力生产在Windows上,所以第一反应是powershell,然而经过几次的尝试,并未成功(后来发现其实只是忘记开始了而已)

所以最后选择的思路是使用cURL。本来我以为cURL只有Linux上有,在Windows上要采取别的实现方式,但是突然之间我发现我多虑了:

接下来就是思考如何通过cURL的参数来绕过referrer的限制:

最后采取的方法:

curl -o <file_dir> -e "<referrer>" <link>

其中file_dir为保存的位置,referrer字面意思,link为下载链接。在Windows里跑了一下,发现也能运行,于是乎我就放心大胆地上了。由于原理和之前失败的powershell一样,都是通过child_process,所以这次仔细查阅了文档,得到了下面的代码:(至于console.log大家忽略就好了(悲,其实那一块就是管道的重定向,问题不是很大)

更新:突然想起来,好像也可能是因为npm忘记装child_process了(我自裁):

npm install child_process

如果觉得太慢,可以前面加上proxychains4

import { spawn } from 'child_process'

...

const curl = spawn('curl', ['-o', path.join(fileDir, path.basename(link)), '-e', '"<referrer>"', link]);
curl.stdout.on('data', (data) => {
    console.log(`stdout: ${data}`);
});
curl.stderr.on('data', (data) => {
    console.log(`stderr: ${data}`);
});
curl.on('close', (code) => {
    console.log(`child process exited with code ${code}`);
});

第一次跑的时候由于没有像上面那样设置好了path所以出现了write permission denied的情况。我百思不得其解,于是乎当时就直接跑了一下pwd命令(使用同样的方法):

const pwd = spawn('pwd');
pwd.stdout.on('data', (data) => {
    console.log(`stdout: ${data}`);
});
pwd.stderr.on('data', (data) => {
    console.log(`stderr: ${data}`);
});
pwd.on('close', (code) => {
    console.log(`child process exited with code ${code}`);
});

得到的是Program Files里面的路径,于是乎就解释地通了。

现在我的想法是,针对每一个文件,将里面地图片下载在同目录的同名文件夹内。经过搜索和先前的经验,现给出如下代码:

// Get file directory and file name
let fileDir = vscode.window.activeTextEditor.document.fileName;
let fileName = path.basename(fileDir);

fileDir = fileDir.substr(0, fileDir.length - 3);
fileName = fileName.substr(0, fileName.length - 3);

if (!fs.existsSync(fileDir)) {
    fs.mkdirSync(fileDir);
}

注意,随后实在判空,如果文件夹不存在就直接创建新文件夹。做完这部之后其实剩下的只需要拼接起来即可,但是没想到这仍然有问题。

我们会发现,我们最后html<img>标签的src属性应该是下载下来的文件路径。我一开始看的是VSCode官方是如何实现的:

我发现他是直接用相对路径的。天真无知的我欣喜若狂,以为事情就这样结束了。我如法炮制:一开始给出的路径是:

`./${fileName}/path.basename(link)`

但是运行之后却发现渲染不出来。我当时注意到,文件名之间有空格:

我当时便认定大概是空格的问题,为了验证这个想法,我随便新建了一个文件,尝试在图片路径中加入空格,果然:

无法顺利渲染。可能是夜晚使人脑子清醒,我想起了浏览器当中会自动把空格转换成%20,于是顺着这个思路去尝试了一番:

可以正确显示图片,随后我便认定了是空格的原因。搜索后得到了encodeURI的处理方法,我将前面的表达式更新为:

encodeURI(`./${fileName}/path.basename(link)`)

但是,调试后仍然发现不行。为了查明问题所在,我决定使用这个插件查看正常情况下图片是如何显示的:

真相逐渐浮出水面,是WebView里面文件协议的问题。以前开发LaTeX插件的时候我记得教程里提过这个问题,便去搜索,顺利找到了原文出处:

按照前面的方法试了一下,仍然不行,并且发现协议为vscode-resource://,而非vscode-webview-resource://。综合了一下教程发布的时间,得到了上述协议已经失效的(大概)的结论。

梅开二度,继续搜索,去VSCodeGithub上面的issues底下找到了答案:https://github.com/microsoft/vscode/issues/102959

注意到,前文使用的vscode-resource://协议为old way,发布时间为2020年,印证了我的猜想。使用了新方法尝试之后仍然不行。思考后决定把之前加上的encodeURI删去,才顺利显示,最终这部分的代码:

// Get file directory and file name
let fileDir = vscode.window.activeTextEditor.document.fileName;
let fileName = path.basename(fileDir);

fileDir = fileDir.substr(0, fileDir.length - 3);
fileName = fileName.substr(0, fileName.length - 3);

if (!fs.existsSync(fileDir)) {
    fs.mkdirSync(fileDir);
}
// replace the image attributes
const $ = cheerio.load(html);
$("img").each((_index, element) => {
    let link = element.attribs["src"];
    if (link.startsWith("<prefix>")) {
        element.attribs["src"] = this.getPreview(sourceUri).webview.asWebviewUri(vscode.Uri.file(path.join(fileDir, path.basename(link)))).toString();
        if (!fs.existsSync(path.join(fileDir, path.basename(link)))) {
            const curl = spawn('curl', ['-o', path.join(fileDir, path.basename(link)), '-e', '"<referrer>"', link]);
            curl.stdout.on('data', (data) => {
                console.log(`stdout: ${data}`);
            });
            curl.stderr.on('data', (data) => {
                console.log(`stderr: ${data}`);
            });
            curl.on('close', (code) => {
                console.log(`child process exited with code ${code}`);
            });
        }
    }
});

html = $.html();

至此,魔改告一段落,主要功能基本全部实现:

也许之后可能也会做一个Configuration的界面?大概吧,还是看心情。

打包安装

在一切准备完毕之后,由于里面涉及到太多个人链接,所以懒得做设置界面和发布了,直接通过vsix打包。如果没有安装过vsix的话:

npm install vsix

安装完毕之后执行打包命令:

vsix package

发现报错:

所以带上--no-yarn的参数再次打包:

再次报错,通过报错信息定位到preview-content-provider.ts文件内,发现问题所在,删除该行:

估计是之前手抖不小心加的(裂。

最后再次打包,成功:

由于我需要在Windows内安装(等会Linux内也装一把),反正把文件移到别处,然后再在如下位置安装:

真香,添加设置界面

转念一想,好像添加设置没什么难的,于是就做了。打开package.json,找到configuration字段,在末尾添加:

"markdown-preview-enhanced.fake-referrer": {
    "description": "The referrer used to download images from restricted servers.",
    "default": "",
    "type": "string"
},
"markdown-preview-enhanced.restricted-prefixes": {
    "description": "If the image link has one of the following prefixes, it will be downloaded using referrer previously configured.",
    "default": [],
    "type": "array"
}

看到作者有专门给config写过config.ts,所以直接到里面去改:

// referrer config
public readonly fakeReferrer: string;
public readonly restrictedPrefixes: string[];

...

this.fakeReferrer = config.get<string>("fake-referrer");
this.restrictedPrefixes = config.get<string[]>("restricted-prefixes");

然后去改一下代码:

$("img").each((_index, element) => {
    let link = element.attribs["src"];
    this.config.restrictedPrefixes.forEach((prefix) => {
        if (link.startsWith(prefix)) {
            element.attribs["src"] = this.getPreview(sourceUri).webview.asWebviewUri(vscode.Uri.file(path.join(fileDir, path.basename(link)))).toString();
            if (!fs.existsSync(path.join(fileDir, path.basename(link)))) {
                const curl = spawn('curl', ['-o', path.join(fileDir, path.basename(link)), '-e', `"${this.config.fakeReferrer}"`, link]);
                curl.stdout.on('data', (data) => {
                    console.log(`stdout: ${data}`);
                });
                curl.stderr.on('data', (data) => {
                    console.log(`stderr: ${data}`);
                });
                curl.on('close', (code) => {
                    console.log(`child process exited with code ${code}`);
                });
            }
        }
    });
});

设置界面:

准备再次发布

首先把这个项目开源了问题也不大,所以就找到原作者的GitHubfork一把:然后再clone到本地,中途用一下proxychains4:

package.json里做一些些修改:

再对LICENSE和其他设置谁做修改就可以进行发布了:

vsce publish --no-yarn

Marketplace的链接:https://marketplace.visualstudio.com/items?itemName=JeffersonQin.blogging-tool

一件小事

用这个插件的时候,右键想要导出pdf发现提示:

搜索,安装:https://www.princexml.com/

但是仍然出现相同提示,转念一想,可能是PATH的原因。打开安装目录,找到路径,然后将C:\Program Files (x86)\Prince\engine\bin添加进PATH,成功解决问题。

结束语

这种项目以后谁爱做谁做去,虽然不是特别难,但是和原作者斗智斗勇真的血腥(逃)。而且完全不熟悉前端(可能只是因为我菜),靠着搜索引擎和感觉来写代码(这真的是好文明嘛?)。不过感觉,与其最后只记录答案是如何得到的,不如像现在这样把失败的过程记录下来,反思思考,写下自己脑中的推断过程,可能才会学到更多的东西,帮助到更多的人(大概吧)

GAN(对抗生成网络)的基本原理以及数学证明

2020年8月14日 19:40

Preface

今天在PD Lib和DL斗智斗勇时,突然想起了自己非常想学的GAN,机缘巧合下便百度了,得到了以下两篇文章:

于是便对GAN有了初步的了解(以前肯定是心不在焉才没有理解的(划掉)),随后又在五楼生命科学的书架上找到了相关资料,遂学了一波。

GAN概述

2014 年,Ian Goodfellow 和他在蒙特利尔大学的同事发表了一篇震撼学界的论文。没错,我说的就是《Generative Adversarial Nets》,这标志着生成对抗网络(GAN)的诞生,而这是通过对计算图和博弈论的创新性结合。他们的研究展示,给定充分的建模能力,两个博弈模型能够通过简单的反向传播(backpropagation)来协同训练。

这两个模型的角色定位十分鲜明。给定真实数据集 R,G 是生成器(generator),它的任务是生成能以假乱真的假数据;而 D 是判别器 (discriminator),它从真实数据集或者 G 那里获取数据, 然后做出判别真假的标记。Ian Goodfellow 的比喻是,G 就像一个赝品作坊,想要让做出来的东西尽可能接近真品,蒙混过关。而 D 就是文物鉴定专家,要能区分出真品和高仿(但在这个例子中,造假者 G 看不到原始数据,而只有 D 的鉴定结果——前者是在盲干)。

理想情况下,D 和 G 都会随着不断训练,做得越来越好——直到 G 基本上成为了一个“赝品制造大师”,而 D 因无法正确区分两种数据分布输给 G。

实践中,Ian Goodfellow 展示的这项技术在本质上是:G 能够对原始数据集进行一种无监督学习,找到以更低维度的方式(lower-dimensional manner)来表示数据的某种方法。而无监督学习之所以重要,就好像 Yann LeCun 的那句话:“无监督学习是蛋糕的糕体”。这句话中的蛋糕,指的是无数学者、开发者苦苦追寻的“真正的 AI”。
——pytorch实现GAN

GAN - Generative Adversarial Nets, 生成对抗网络,简单来讲其有两个组成部分:

  • D (Discriminator) - 判别器,判断输入时捏造的还是真实的
  • G (Generator) - 生成器,从随机噪声中生成我们想要的数据

随着训练的进行,我们要提高D的辨析能力,但同时也要G的能力来骗过D,因为我们的最终目的是要让G来生成可以骗过D的信息。总结来说,通过对这两个模型的训练,我们就可以找到随机噪声与有意义数据的映射,达到创作的目的。

GAN的流程和目标函数

GAN的目标函数

GAN的目标函数如下:

$$
V(D, G) = \mathbb E_{x \sim P_{\text{data}} } [\log D(x)] + \mathbb E_{\boldsymbol {z} \sim P_z}[\log(1 - D(G(\boldsymbol {z})))]
$$

其中,$D$为Discriminator的模型函数,$G$为Generator的模型函数,随机变量$x$服从原来正确的数据集的分布$P_\text {data}$,随机变量(这里可能是高维随机变量,取决于模型具体实现)$\boldsymbol {z}$服从分布$P_z$(生成噪音),$\mathbb E$代表期望。

GAN的流程

$$
G^* = \arg \min _G \max _D V(D, G)
$$

即,可以分为两步理解:

  1. 在$G$为常数的情况下,选择合适的$D$使得$V(D,G)$能够最大化。
  2. 在这之后,选取合适的$G$来最小化$V(D, G)$,这个$G$就是我们想要的生成模型。

在每一步的训练中:

  • 取$m$个真实数据:

    $$
    \{x^{(1)}, x^{(2)}, x^{(3)}, \cdots, x^{(m)}\}
    $$

    使用$G$和$m$组随机数(服从于噪音分布$P_G$,一般使用服从正态分布的随机数)

    $$
    \{\boldsymbol {z}^{(1)}, \boldsymbol {z}^{(2)}, \boldsymbol {z}^{(3)}, \cdots, \boldsymbol {z}^{(m)}\}
    $$

    生成$m$个假数据,其中

    $$
    \forall i \in [1, m], i \in \mathbb Z \Rightarrow x^{(i)} \sim P_{\text{data}}, \boldsymbol {z} ^ {(i)} \sim P_z
    $$

  • 根据$\max$部分的目标使用随机梯度上升(Stochastic Gradient Ascent)更新$D$的参数,提高$D$的分辨能力

    $$
    \theta_d := \theta_d + \alpha_d \nabla_{\theta_d} \frac 1 m \sum _ {i = 1} ^ m \Big [ \log D\big ( x^{(i)} \big) + \log \big ( 1 - D \big ( G(\boldsymbol{z}^{(i)}) \big ) \big ) \Big]
    $$

  • 根据$\min$部分的目标使用随机梯度下降(Stochastic Gradient Descent)更新$G$的参数,使$G$生成的数据更有迷惑性

    $$
    \theta_g := \theta_g - \alpha_g\nabla_{\theta_g} \frac 1 m \sum _{i = 1} ^ m\log \Big ( 1 - D\big ( G(\boldsymbol{z} ^ {(i)})\big ) \Big )
    $$

GAN的数学原理

Prerequisites

信息量(自信息)

信息量是指信息多少的量度,即,对于一条信息,传达这条信息所需的最少信息长度为自信息。

信息论创始人C.E.Shannon,1938年首次使用比特(bit)概念:1(bit)= $\log_2 2$。它相当于对二个可能结局所作的一次选择量。信息论采用对随机分布概率取对数的办法,解决了不定度的度量问题。

定义:符合分布$P$的某一事件$x$出现,传达出这条信息的信息量记为:

$$
I = \log \frac 1 {P(x)} = - \log P(x)
$$

香农熵

从离散分布$P$中随机抽选一个事件,传达这条信息所需的最优平均信息长度为香农熵,表达为:

$$
H(P) = \sum_x P(x) \log \frac 1 {P(x)} = - \sum_x P(x) \log P(x)
$$

若分布是连续的,则:

$$
H(P) = \int_x P(x) \log \frac 1 {P(x)} \mathrm dx = -\int_x P(x) \log P(x) \mathrm dx
$$

交叉熵

用分布$P$的最佳信息传递方式来传达分布$Q$中随机抽选的一个事件,所需的平均信息长度为交叉熵,表达为

$$
H_P(Q) = \sum_x Q(x) \log \frac 1 {P(x)} = - \sum_x Q(x) \log P(x)
$$

$$
H_P(Q) = \int_x Q(x) \log \frac 1 {P(x)} \mathrm dx = - \int_x Q(x) \log P(x) \mathrm dx
$$

$KL$ Divergence

$KL$散度:用分布$P$的最佳信息传递方式来传达分布$Q$,比用分布$Q$自己的最佳信息传递方式来传达分布$Q$,平均多耗费的信息长度为$KL$散度,表达为$D_P(Q)$或$D_{KL}(Q||P)$,$KL$散度衡量了两个分布之间的差异。

$$
\begin{aligned}
D_{KL}(Q||P) = D_P(Q) &= H_P(Q) - H(Q) \\
&= \sum_x Q(x) \log \frac 1 {P(x)} - \sum _x Q(x) \log \frac 1 {Q(x)} \\
&= \sum_x Q(x) \log \frac {Q(x)} {P(x)}
\end{aligned}
$$

对于连续分布:

$$
D_{KL}(Q||P) = D_P(Q) = \int_{-\infty} ^ {\infty} P(x)\log \frac {P(x)}{Q(x)} \mathrm dx
$$

KL Divergence越大,两个分布差异越大,反之差异越小。

数学原理

看完Prerequisites,我们回归正题讨论GAN的原理。我们现在想要做的事情,其实就是将一个服从$P_G$的随机噪声$\boldsymbol z$通过一个生成网络$G$得到一个和真实数据分布$P_{\text {data}}(x)$差不多的生成分布$P_G(x;\theta_g)$,其中$\theta_g$为生成网络$G$的参数。我们希望找到一个$\theta_g$使得两个分布$P_{\text {data}}(x)$与$P_G(x;\theta)$尽可能地相似(使得他们地KL散度尽可能得小)。

我们从真实数据分布$P_\text{data}(x)$中取$m$个样本,记作:

$$
\{x^{(1)}, x^{(2)}, x^{(3)}, \cdots, x^{(m)}\}
$$

根据生成网络的参数$\theta_g$,我们可以计算出这$m$个真实样本在生成网络中出现的概率$P_G(x^{(i)}; \theta_g)$,那么生成这样的$m$个样本数据的似然(likelihood)为:

$$
L = \prod_{i = 1} ^ m P_G(x^{(i)}; \theta_g)
$$

由于我们想要两个分布尽量相似,那么我们肯定希望这个似然$L$尽量大,即生成这样的真实数据的概率尽量大,遂我们最大化这个似然,找到$\theta_g^*$:

$$
\begin{aligned}
\theta_g^* &= \arg \max _ {\theta_g} \prod _ {i = 1} ^ m P_G(x^{(i)}; \theta_g) \\
&\Leftrightarrow \arg \max _ {\theta_g} \log \prod _ {i = 1} ^ m P_G(x^{(i)}; \theta_g) \\
&= \arg \max _ {\theta_g} \sum_{i = 1} ^ m \log P_G(x^{(i)}; \theta_g) \\
&\approx \arg \max _ {\theta_g} \mathbb E_{x \sim P_{\text {data}}} [\log P_G(x; \theta_g)] \\
&= \arg \max _ {\theta_g} \int _ x P_\text{data}(x) \log P_G(x; \theta_g)\mathrm dx \\
&= \arg \max _ {\theta_g} \int _ x P_\text{data}(x) \log P_G(x; \theta_g)\mathrm dx \\&~~~~- \int _ x P_\text{data}(x) \log P_\text {data}(x)\mathrm dx \\
&= \arg \max _ {\theta_g} \int _ x P_\text {data} (x) \log \frac {P_G(x; \theta_g)} {P_\text {data}(x)} \mathrm dx \\
&= \arg \max _ {\theta_g} - \int _ x P_\text {data} (x) \log \frac {P_\text {data}(x)} {P_G(x; \theta_g)} \mathrm dx \\
&= \arg \min _ {\theta_g} KL(P_\text {data} || P_G(x; \theta))
\end{aligned}
$$

所以可见,其实最大化这个似然,和最小化KL散度是基本相同的。

上述式子中,$P_G(x;\theta_g)$代表在生成分布中出现$x$的概率,也可以如下计算得到:

$$
P_G(x) = \int _ \boldsymbol{z} P_z(z)\cdot 1\{G(z) = x\} \mathrm dz
$$

注:$1\{\cdot\}$的含义是若打括号内的逻辑运算为真则取$1$,假则取$0$. 即

$$
1\{\text {True}\} = 1, 1 \{\text {False}\} = 0
$$

但是我们发现,上述的过程是难以进行计算的,甚至完全没办法求$P_G(x)$,这只是模型的想法而已。

现在我们看回之前我们提到的目标函数:

$$
V(D, G) = \mathbb E_{x \sim P_{\text{data}} } [\log D(x)] + \mathbb E_{\boldsymbol {z} \sim P_z}[\log(1 - D(G(\boldsymbol {z})))]
$$

与最优化生成模型:

$$
G^* = \arg \min _G \max _D V(D, G)
$$

我们接下来分步解释。

首先,我们不妨解释一下$\max_D V(G, D)$,这部分的含义之前也解释过,是在给定$G$的情况下,最大化$V(G, D)$。观察发现,其形式其实与交叉熵损失函数非常相似:

$$
\mathcal L(\hat y, y) = -(y\log \hat y + (1 - y) \log (1 - \hat y))
$$

其实他们表达的目的也差不多。我们先化简一下$V(G, D)$看看能得到什么结果:

$$
\begin{aligned}
V(G, D) &= \mathbb E_{x \sim P_{\text{data}} } [\log D(x)] + \mathbb E_{\boldsymbol {z} \sim P_z}[\log(1 - D(G(\boldsymbol {z})))] \\
&= \int_x P_\text {data}(x) \log D(x) \mathrm dx + \int_\boldsymbol{z} P_z(\boldsymbol{z}) \log (1 - D(G(\boldsymbol{z})))\mathrm dz \\
&= \int_x P_\text {data}(x) \log D(x) \mathrm dx + \int_x P_G(x)\log (1 - D(x)) \mathrm dx \\
&= \int_x [P_\text {data}(x) \log D(x) + P_G(x)\log (1 - D(x))] \mathrm dx
\end{aligned}
$$

让我们考察积分内部的项,我们可以对它做指数运算,即:

$$
e^{P_\text {data}(x) \log D(x) + P_G(x)\log (1 - D(x))} = D(x) ^ {P_\text {data}}\times(1- D(x))^{P_G(x)}
$$

其想表达什么便不言而喻了,它表达的就是判别器判别是真的的正确率和判别是假的的正确率,总体来说就是衡量$D$的能力,所以我们想要最大化$V$,提高$D$的判别能力。

$$
f(D) = P_\text {data}(x) \log D(x) \mathrm dx + P_G(x)\log (1 - D(x))
$$

因为这里$P_\text {data} (x)$和$P_G(x)$都可以看作常数,所以

$$
\arg \max _ D \int_x f(D) \mathrm dx = \arg \max _D f(D)
$$

最大化$f(D)$,即令其导数为$0$:

$$
\frac {\mathrm d f(D)} {\mathrm dD } = \frac {P_\text {data}(x)}{D} - \frac {P_G(x)} {1 - D} = 0
$$

则:

$$
\frac {P_\text {data}(x)}{D} = \frac {P_G(x)} {1 - D}
$$

$$
D^*(x) = \frac {P_\text {data}(x)}{P_\text {data}(x) + P_G(x)}
$$

这样,我们就得到了那个状态下最优的$D^*$的表达式。我们将这个能够最大化$V$的$D$代入回$V(G, D)$:

$$
\begin{aligned}
\max V(G, D) &= V(G, D^*) \\
&= \mathbb E_{x \sim P_\text {data}} \Big [ \log \frac {P_\text {data}(x)}{P_\text {data}(x) + P_G(x)} \Big] + \mathbb E_{x \sim P_G} \Big [ \frac {P_G(x)}{P_\text {data}(x) + P_G(x)}\Big ] \\
&= \int_x P_\text {data} (x) \log \frac {\frac 1 2 P_\text {data}}{\frac 1 2 (P_\text{data}(x) + P_G(x))} \mathrm dx + \int _ x P_G(x) \log \frac {\frac 1 2 P_G(x)} {\frac 1 2 (P_\text {data} (x) + P_G(x))} \mathrm dx \\
&= -2\log 2 + D_{KL} (P_\text {data} || \frac 1 2 [P_\text{data}(x) + P_G(x)]) + D_{KL} (P_G(x) || \frac 1 2 [P_\text {data}(x) + P_G(x)]) \\
&= -2\log 2 + JSD (P_\text{data}(x) || P_G(x))
\end{aligned}
$$

其中,我们引入了$JS$ Divergence,定义如下:

$$
JSD(P||Q) = \frac 1 2 D_{KL} (P||M) + \frac 1 2 D_{KL} (Q || M), M = \frac 1 2 (P + Q)
$$

容易得到,KL Divergence是不对称的,而JS Divergence是对称的。他们都可以衡量两组分布建的差异。这里我们想要两组分布差异最小,故取$\min$

所以,这也就解释了为什么:

$$
\arg \min _G \max _D V(G, D)
$$

是我们的目标过程。

Canny边缘检测算法(基于OpenCV的Java实现)

2020年3月5日 00:21

Canny边缘检测算法(基于OpenCV的Java实现)

绪论

最近在学习ORB的过程中又仔细学习了Canny,故写下此篇笔记,以作总结。

Canny边缘检测算法的发展历史

Canny边缘检测于1986年由JOHN CANNY首次在论文《A Computational Approach to Edge Detection》中提出,就此拉开了Canny边缘检测算法的序幕。

Canny边缘检测是从不同视觉对象中提取有用的结构信息并大大减少要处理的数据量的一种技术,目前已广泛应用于各种计算机视觉系统。Canny发现,在不同视觉系统上对边缘检测的要求较为类似,因此,可以实现一种具有广泛应用意义的边缘检测技术。边缘检测的一般标准包括:

  • 以低的错误率检测边缘,也即意味着需要尽可能准确的捕获图像中尽可能多的边缘。
  • 检测到的边缘应精确定位在真实边缘的中心。
  • 图像中给定的边缘应只被标记一次,并且在可能的情况下,图像的噪声不应产生假的边缘。

为了满足这些要求,Canny使用了变分法。Canny检测器中的最优函数使用四个指数项的和来描述,它可以由高斯函数的一阶导数来近似。

在目前常用的边缘检测方法中,Canny边缘检测算法是具有严格定义的,可以提供良好可靠检测的方法之一。由于它具有满足边缘检测的三个标准和实现过程简单的优势,成为边缘检测最流行的算法之一。

Canny边缘检测算法的处理流程

Canny边缘检测算法可以分为以下5个步骤:

  • 使用高斯滤波器,以平滑图像,滤除噪声。
  • 计算图像中每个像素点的梯度强度和方向。
  • 应用非极大值(Non-Maximum Suppression)抑制,以消除边缘检测带来的杂散响应。
  • 应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。
  • 通过抑制孤立的弱边缘最终完成边缘检测。

下面详细介绍每一步的实现思路。

用高斯滤波器平滑图像

高斯滤波是一种线性平滑滤波,适用于消除高斯噪声,特别是对抑制或消除服从正态分布的噪声非常有效。滤波可以消除或降低图像中噪声的影响,使用高斯滤波器主要是基于在滤波降噪的同时也可以最大限度保留边缘信息的考虑。

高斯滤波实现步骤:

彩色RGB图像转换为灰度图像

边缘检测是基于对图像灰度差异运算实现的,所以如果输入的是RGB彩色图像,需要先进行灰度图的转换。
RGB转换成灰度图像的一个常用公式是:

$$
Gray = R*0.299 + G*0.587 + B*0.114
$$

注意一般情况下图像处理中彩色图像各分量的排列顺序是B、G、R。

RGB原图像:

转换后的灰度图:

Java代码调用系统库实现:

public static Mat RGB2Gray(Mat image) {
    // Gray = R*0.299 + G*0.587 + B*0.114
    Mat gray = new Mat();
    Imgproc.cvtColor(image, gray, Imgproc.COLOR_BGR2GRAY);
    return gray;
}

一维,二维高斯函数及分布

一维高斯函数表述为:

$$
G(x) = \frac {1}{\sqrt {2\pi}\sigma}\exp(-\frac {(x-\mu_x)^2}{2\sigma^2})
$$

对应图形:

二维高斯函数表述为:

$$
G(x) = \frac {1}{ {2\pi}\sigma^2}\exp(-\frac {(x-\mu_x)^2+(y-\mu_y)^2}{2\sigma^2})
$$

对应图形:

一些重要特性说明:

  1. 一维二维高斯函数中$μ$是服从正态分布的随机变量的均值,称为期望或均值影响正态分布的位置,实际的图像处理应用中一般取$μ=0;~σ$是标准差,$σ^2$是随机变量的方差,$σ$定义了正态分布数据的离散程度,$σ$越大,数据分布越分散,σ越小,数据分布越集中。

在图形或滤波效果上表现为:$σ$越大,曲线越扁平,高斯滤波器的频带就越宽,平滑程度就越好,$σ$越小,曲线越瘦高,高斯滤波的频带就越窄,平滑程度也越弱;

  1. 二维高斯函数具有旋转对称性,即滤波器在各个方向上的平滑程度是相同的.一般来说,一幅图像的边缘方向是事先不知道的,因此,在滤波前是无法确定一个方向上比另一方向上需要更多的平滑.旋转对称性意味着高斯平滑滤波器在后续边缘检测中不会偏向任一方向;
  2. 高斯函数是单值函数。这表明,高斯滤波器用像素邻域的加权均值来代替该点的像素值,而每一邻域像素点权值是随该点与中心点的距离单调增减的。这一性质是很重要的,因为边缘是一种图像局部特征,如果平滑运算对离算子中心很远的像素点仍然有很大作用,则平滑运算会使图像失真;
  3. 相同条件下,高斯卷积核的尺寸越大,图像的平滑效果越好,表现为图像越模糊,同时图像细节丢失的越多;尺寸越小,平滑效果越弱,图像细节丢失越少;

以下对比一下不同大小标准差$σ$(Sigma)对图像平滑的影响:

原图:

卷积核尺寸5*5,σ=0.1:

卷积核尺寸5*5,σ=1:

对比可以看到,Sigma(σ)越大,平滑效果越明显。

生成高斯滤波卷积核

滤波的主要目的是降噪,一般的图像处理算法都需要先进行降噪。而高斯滤波主要使图像变得平滑(模糊),同时也有可能增大了边缘的宽度。

高斯函数是一个类似与正态分布的中间大两边小的函数。

对于一个位置(m,n)的像素点,其灰度值(这里只考虑二值图)为f(m,n)。

那么经过高斯滤波后的灰度值将变为:

$$
g_\sigma(m,n)=\frac 1 {2\pi\sigma^2} \exp(-\frac {m^2+n^2}{2\sigma^2})f(m,n)
$$

其中,

$$
\frac 1 {{2\pi\sigma^2}} \exp(-\frac {m^2+n^2}{2\sigma^2})
$$

是二元高斯函数。

为了尽可能减少噪声对边缘检测结果的影响,所以必须滤除噪声以防止由噪声引起的错误检测。为了平滑图像,使用高斯滤波器与图像进行卷积,该步骤将平滑图像,以减少边缘检测器上明显的噪声影响。大小为(2k+1)x(2k+1)的高斯滤波器核的生成方程式由下式给出:

下面是一个sigma = 1.4,尺寸为3x3的高斯卷积核的例子(需要注意归一化):

若图像中一个3x3的窗口为A,要滤波的像素点为e,则经过高斯滤波之后,像素点e的亮度值为:

其中*为卷积符号,sum表示矩阵中所有元素相加求和。

重要的是需要理解,高斯卷积核大小的选择将影响Canny检测器的性能。尺寸越大,检测器对噪声的敏感度越低,但是边缘检测的定位误差也将略有增加。一般5x5是一个比较不错的trade off。

public static double[][] getGaussianArray(int size, double sigma) {
    double[][] array = new double[size][size];
    double center_i, center_j;
    if (size % 2 == 1) { center_i = (double) (size / 2); center_j = (double) (size / 2); }
    else { center_i = (double) (size / 2) - 0.5; center_j = (double) (size / 2) - 0.5; }
    double sum = 0.0;
    for (int i = 0; i < size; i ++)
            for (int j = 0; j < size; j ++) {
                array[i][j] = Math.exp((-1.0) * ((i - center_i) * (i - center_i) + (j - center_j) * (j - center_j)) / 2.0 / sigma / sigma);
                sum += array[i][j];
            }
    for (int i = 0; i < size; i ++)
            for (int j = 0; j < size; j++)
                array[i][j] /= sum;
    return array;
}

Sigma为1时,求得的3*3大小的高斯卷积核参数为:

Sigma为1,5*5大小的高斯卷积核参数为:

在计算出高斯滤波卷积核之后就是用这个卷积核来扫过整张图像,对每个像素点进行加权平均。

单色高斯滤波与彩色高斯滤波

加入了多线程的优化,代码实现:

public static Mat greyGaussianFilter(Mat src, double[][] array, int size) throws InterruptedException {
    Mat temp = src.clone();
    final CountDownLatch latch = new CountDownLatch(src.rows());
    for (int i = 0; i < src.rows(); i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 0; j < src.cols(); j ++)
                        if (finalI > (size / 2) - 1 && j > (size / 2) - 1 &&
                            finalI < src.rows() - (size / 2) && j < src.cols() - (size / 2)) {
                            double sum = 0.0;
                            for (int k = 0; k < size; k++)
                                    for (int l = 0; l < size; l++)
                                        sum += src.get(finalI - k + size / 2, j - l + size / 2)[0] * array[k][l];
                            temp.put(finalI, j, sum);
                        }
                latch.countDown();
            }).start();
    }
    latch.await();
    return temp;
}

public static Mat colorGaussianFilter(Mat src, int size, double sigma) throws InterruptedException {
    // return variable
    Mat ret = new Mat();
    // list for merge and split
    List<Mat> channels = new ArrayList<>();
    List<Mat> new_channels = new ArrayList<>();
    Map<Integer, Mat> temp_channels = new TreeMap<>();
    // split into 3 channels (r, g, b)
    Core.split(src, channels);
    // get gaussion array
    double[][] array = SmartGaussian.getGaussianArray(size, sigma);
    // multi-thread
    final CountDownLatch latch = new CountDownLatch(channels.size());
    channels.forEach(mat -> {
            new Thread(() -> {
                Mat tmp = new Mat();
                try {
                        tmp = SmartGaussian.greyGaussianFilter(mat, array, size);
                } catch (InterruptedException e) {
                        e.printStackTrace();
                }
                temp_channels.put(channels.indexOf(mat), tmp);
                latch.countDown();
            }).start();
    });
    latch.await();
    for (int i = 0; i < channels.size(); i ++)
          new_channels.add(temp_channels.get(i));
    Core.merge(new_channels, ret);
    return ret;
}

效果图:(左图为原图,中间为上面的实现,右边是OpenCV实现)

用Sobel等梯度算子计算梯度幅值和方向

梯度

梯度的本意是一个向量(矢量),表示某一函数在该点处的方向导数沿着该方向取得最大值,即函数在该点处沿着该方向(此梯度的方向)变化最快,变化率最大(为该梯度的模)。

设二元函数:

$$
z=f(x,y)
$$

在平面区域D上具有一阶连续偏导数,则对于每一个点$P(x,y)$都可定出一个向量:

$$
\Big\{\frac {\partial f}{\partial x},\frac {\partial f} {\partial y} \Big\}
=f_x(x,y)\vec i + f_y(x,y)\vec j
$$

该函数就称为函数$z=f(x,y)$在点$P(x,y)$的梯度,记作$\text{grad}~f(x,y)$或$\nabla f(x,y)$,即有:

$$
\text{grad}~f(x,y)=\nabla f(x,y)=\Big\{\frac {\partial f}{\partial x},\frac {\partial f} {\partial y} \Big\}
=f_x(x,y)\vec i + f_y(x,y)\vec j
$$

其中$\nabla=\displaystyle\frac{\partial}{\partial x}\vec i+\displaystyle\frac{\partial}{\partial y}\vec j$称为(二维的)向量微分算子或Nabla算子。

设$e=\{\cos\alpha,\cos\beta\}$是方向$l$上的单位向量,则

$$
\begin{aligned}
&\displaystyle\frac{\partial f}{\partial l} \\ =& \displaystyle\frac{\partial f}{\partial x} \cos \alpha + \displaystyle\frac{\partial f}{\partial y} \cos \beta \\ =& \Big \{\displaystyle\frac{\partial f}{\partial x},\displaystyle\frac{\partial f}{\partial y} \Big\} \{\cos\alpha, \cos \beta\} \\ =& \operatorname{grad}f(x,y)e \\ =& |\operatorname{grad}f(x,y)||e|\cos[\operatorname{grad}f(x,y),e]
\end{aligned}
$$

由于当方向$l$与梯度方向一致时,有$\cos[\operatorname{grad}f(x,y),e]=1$。

所以当$l$与梯度方向一致时,方向导数$\displaystyle\frac{\partial f}{\partial x}$有最大值,且最大值为梯度的模,即$|\operatorname{grad}f(x,y)| = \sqrt{\Big(\displaystyle\frac{\partial f}{\partial x}\Big)^2+\Big(\displaystyle\frac{\partial f}{\partial y}\Big)^2}$。

因此说,函数在一点沿梯度方向的变化率最大,最大值为该梯度的

图像灰度值的梯度的简单求法

图像灰度值的梯度可以使用最简单的一阶有限差分来进行近似,使用以下图像在x和y方向上偏导数的两个矩阵:

计算公式为:

其中$f$为图像灰度值,$P[i,j]$代表$[i,j]$在$X$方向梯度幅值,$Q[i,j]$代表$[i,j]$在$Y$方向的梯度幅值,$M[i,j]$是该点幅值,$\theta[i,j]$是梯度方向,也就是角度。

使用Sobel算子来计算梯度的大小及方向:

图像中的边缘可以指向各个方向,因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。边缘检测的算子(如Roberts,Prewitt,Sobel等)返回水平Gx和垂直Gy方向的一阶导数值,由此便可以确定像素点的梯度G和方向theta 。

其中G为梯度强度, theta表示梯度方向,arctan为反正切函数。下面以Sobel算子为例讲述如何计算梯度强度和方向。

x和y方向的Sobel算子分别为:

其中Sx表示x方向的Sobel算子,用于检测y方向的边缘; Sy表示y方向的Sobel算子,用于检测x方向的边缘(边缘方向和梯度方向垂直)。在直角坐标系中,Sobel算子的方向如下图所示。

若图像中一个3x3的窗口为A,要计算梯度的像素点为e,则和Sobel算子进行卷积之后,像素点e在x和y方向的梯度值分别为:

其中*为卷积符号,sum表示矩阵中所有元素相加求和。根据公式(3-2)便可以计算出像素点e的梯度和方向。

下面是Sobel算子求梯度的java实现:

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.Core;
import org.opencv.core.CvType;
import org.opencv.core.Mat;

import java.util.concurrent.CountDownLatch;

public class SmartSobel {

    private Mat gradientX = new Mat(), gradientY = new Mat(), gradientXY = new Mat();
    private double[][] pointDirection;

    public SmartSobel() {}

    public void compute(Mat image) throws InterruptedException {
        pointDirection = new double[image.rows()][image.cols()];
        for (int i = 0; i < image.rows(); i ++)
            for (int j = 0; j < image.cols(); j ++)
                pointDirection[i][j] = 0;
        gradientX = Mat.zeros(image.size(), CvType.CV_32SC1);
        gradientY = Mat.zeros(image.size(), CvType.CV_32SC1);
        gradientXY = Mat.zeros(image.size(), CvType.CV_32SC1);
        final CountDownLatch latch = new CountDownLatch(image.rows() - 2);
        for (int i = 1; i < image.rows() - 1; i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 1; j < image.cols() - 1; j++) {
                    double gX = (-1) * image.get(finalI - 1, j - 1)[0] +
                            1 * image.get(finalI - 1, j + 1)[0] +
                            (-2) * image.get(finalI, j - 1)[0] +
                            2 * image.get(finalI, j + 1)[0] +
                            (-1) * image.get(finalI + 1, j - 1)[0] +
                            1 * image.get(finalI + 1, j + 1)[0];
                    double gY = 1 * image.get(finalI - 1, j - 1)[0] +
                            2 * image.get(finalI - 1, j)[0] +
                            1 * image.get(finalI - 1, j + 1)[0] +
                            (-1) * image.get(finalI + 1, j - 1)[0] +
                            (-2) * image.get(finalI + 1, j)[0] +
                            (-1) * image.get(finalI + 1, j + 1)[0];
                    gradientY.put(finalI, j, Math.abs(gY));
                    gradientX.put(finalI, j, Math.abs(gX));
                    gradientXY.put(finalI, j, Math.sqrt(gX * gX + gY * gY));
                    // 防止除以0的情况发生
                    if (gX == 0) gX = 0.00000000000000001;
                    pointDirection[finalI][j] = Math.atan(gY / gX);
                }
                latch.countDown();
            }).start();
        }
        latch.await();
    }

    public void convert() {
        Core.convertScaleAbs(gradientX, gradientX);
        Core.convertScaleAbs(gradientY, gradientY);
        Core.convertScaleAbs(gradientXY, gradientXY);
    }

    public Mat getGradientX() { return this.gradientX; }

    public Mat getGradientY() { return this.gradientY; }

    public Mat getGradientXY() { return this.gradientXY; }

    public double[][] getPointDirection() { return this.pointDirection; }

}

由于这里使用的是$Math.tan()$,所以最终的角度是映射在$[-\frac \pi 2, \frac \pi 2]$的范围之内的。如果使用$Math.tan2()$会映射到$[-\pi,\pi]$的范围内,并且无需考虑符号影响,更加精确。但是这里我们并不关心另外的一个$\pi$的情况,我们只关心其所在直线(这在后文中会提到,也就是非极大值抑制),所以无需多考虑。

X方向梯度图:

Y方向梯度图:

X、Y方向梯度融合效果:

Opencv Sobel函数效果:

对梯度幅值进行非极大值抑制

非极大值抑制是一种边缘稀疏技术,非极大值抑制的作用在于“瘦”边。对图像进行梯度计算后,仅仅基于梯度值提取的边缘仍然很模糊。对于标准3,对边缘有且应当只有一个准确的响应。而非极大值抑制则可以帮助将局部最大值之外的所有梯度值抑制为0,对梯度图像中每个像素进行非极大值抑制的算法是:

  1. 将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。
  2. 如果当前像素的梯度强度与另外两个像素相比最大,则该像素点保留为边缘点,否则该像素点将被抑制。

通常为了更加精确的计算,在跨越梯度方向的两个相邻像素之间使用线性插值来得到要比较的像素梯度,现举例如下:

          图3-2 梯度方向分割

如图3-2所示,将梯度分为8个方向,分别为E、NE、N、NW、W、SW、S、SE,其中0代表$0^\circ\sim45^\circ$,1代表$45^\circ\sim90^\circ$,2代表$-90^\circ\sim-45^\circ$,3代表$-45^\circ\sim0^\circ$。像素点P的梯度方向为$\theta$,则像素点P1和P2的梯度线性插值为:

$$
\begin{aligned}
\tan \theta = G_y ~/~G_x \\
G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW
\end{aligned}
$$

上面也只是图中的情况,具体情况如下:

$$
\begin{aligned}
\theta \in [0, \frac \pi 4]:
\begin {cases}
G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW
\end {cases}\\\theta \in [\frac \pi 4, \frac \pi 2]:
\begin {cases}
G_{p1} = (1-\tan\theta)\times N + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times S + \tan\theta \times SW
\end {cases}\\\theta \in [-\frac \pi 4, 0]:
\begin {cases}
G_{p1} = (1-\tan(-\theta))\times E + \tan(-\theta) \times SE \\
G_{p2} = (1-\tan(-\theta))\times W + \tan(-\theta) \times NW
\end {cases}\\\theta \in [-\frac \pi 2, -\frac \pi 4]:
\begin {cases}
G_{p1} = (1-\tan(-\theta))\times S + \tan(-\theta) \times SE \\
G_{p2} = (1-\tan(-\theta))\times N + \tan(-\theta) \times NW
\end {cases}
\end{aligned}
$$

因此非极大值抑制的伪代码描写如下:

需要注意的是,如何标志方向并不重要,重要的是梯度方向的计算要和梯度算子的选取保持一致。

Java实现:

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.CvType;
import org.opencv.core.Mat;

import java.util.concurrent.CountDownLatch;

public class SmartNMS {

    public static Mat NMS(Mat gradientImage, double[][] pointDirection) throws InterruptedException {
        Mat outputImage = gradientImage.clone();
        final CountDownLatch latch = new CountDownLatch(gradientImage.rows() - 2);
        for (int i = 1; i < gradientImage.rows() - 1; i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 1; j < gradientImage.cols() - 1; j ++) {
                    double GP = gradientImage.get(finalI, j)[0],
                            E = gradientImage.get(finalI, j + 1)[0],
                            NE = gradientImage.get(finalI - 1, j + 1)[0],
                            N = gradientImage.get(finalI - 1, j)[0],
                            NW = gradientImage.get(finalI - 1, j - 1)[0],
                            W = gradientImage.get(finalI, j - 1)[0],
                            SW = gradientImage.get(finalI + 1, j - 1)[0],
                            S = gradientImage.get(finalI + 1, j)[0],
                            SE = gradientImage.get(finalI + 1, j + 1)[0];
                    double GP1 = 0, GP2 = 0;
                    double theta = pointDirection[finalI][j];
                    if (theta >= 0 && theta <= Math.PI / 4) {
                        GP1 = E * (1 - Math.tan(theta)) + NE * Math.tan(theta);
                        GP2 = W * (1 - Math.tan(theta)) + SW * Math.tan(theta);
                    } else if (theta > Math.PI / 4) {
                        GP1 = N * (1 - 1 / Math.tan(theta)) + NE * 1 / Math.tan(theta);
                        GP2 = S * (1 - 1 / Math.tan(theta)) + SW * 1 / Math.tan(theta);
                    } else if (theta < 0 && theta >= -Math.PI / 4) {
                        GP1 = E * (1 - Math.tan(-theta)) + SE * Math.tan(-theta);
                        GP2 = W * (1 - Math.tan(-theta)) + NW * Math.tan(-theta);
                    } else {
                        GP1 = S * (1 - 1 / Math.tan(-theta)) + SE * 1 / Math.tan(-theta);
                        GP2 = N * (1 - 1 / Math.tan(-theta)) + NW * 1 / Math.tan(-theta);
                    }
                    if (GP < GP1 || GP < GP2) outputImage.put(finalI, j, 0);
                }
                latch.countDown();
            }).start();
        }
        latch.await();
        return outputImage;
    }

}

双阈值检测

在施加非极大值抑制之后,剩余的像素可以更准确地表示图像中的实际边缘。然而,仍然存在由于噪声和颜色变化引起的一些边缘像素。为了解决这些杂散响应,必须用弱梯度值过滤边缘像素,并保留具有高梯度值的边缘像素,可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值,则将其标记为强边缘像素;如果边缘像素的梯度值小于高阈值并且大于低阈值,则将其标记为弱边缘像素;如果边缘像素的梯度值小于低阈值,则会被抑制。阈值的选择取决于给定输入图像的内容。

双阈值检测的伪代码描写如下:

抑制孤立低阈值点

到目前为止,被划分为强边缘的像素点已经被确定为边缘,因为它们是从图像中的真实边缘中提取出来的。然而,对于弱边缘像素,将会有一些争论,因为这些像素可以从真实边缘提取也可以是因噪声或颜色变化引起的。为了获得准确的结果,应该抑制由后者引起的弱边缘。通常,由真实边缘引起的弱边缘像素将连接到强边缘像素,而噪声响应未连接。为了跟踪边缘连接,通过查看弱边缘像素及其8个邻域像素,只要其中一个为强边缘像素,则该弱边缘点就可以保留为真实的边缘。

抑制孤立边缘点的伪代码描述如下:

实现:(使用递归)

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.Mat;

import java.util.ArrayList;
import java.util.List;

public class SmartCanny {

    private List<Integer[]> highPoints = new ArrayList<Integer[]>();

    private void DoubleThreshold(Mat image, double lowThreshold, double highThreshold) {
        for (int i = 1; i < image.rows() - 1; i ++)
            for (int j = 1; j < image.cols() - 1; j ++)
                if (image.get(i, j)[0] >= highThreshold) {
                    image.put(i, j, 255);
                    Integer[] p = new Integer[2];
                    p[0] = i; p[1] = j;
                    highPoints.add(p);
                } else if (image.get(i, j)[0] < lowThreshold)
                    image.put(i, j, 0);
    }

    private void DoubleThresholdLink(Mat image, double lowThreshold) {
        for (Integer[] p : highPoints) {
            DoubleThresholdLinkRecurrent(image, lowThreshold, p[0], p[1]);
        }
        for (int i = 1; i < image.rows() - 1; i ++)
            for (int j = 1; j < image.cols() - 1; j ++)
                if (image.get(i, j)[0] < 255)
                    image.put(i, j, 0);
    }

    private void DoubleThresholdLinkRecurrent(Mat image, double lowThreshold, int i, int j) {
        if (i <= 0 || j <= 0 || i >= image.rows() - 1 || j >= image.cols() - 1) return;
        if (image.get(i - 1, j - 1)[0] >= lowThreshold && image.get(i - 1, j - 1)[0] < 255) {
            image.put(i - 1, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j - 1);
        }
        if (image.get(i - 1, j)[0] >= lowThreshold && image.get(i - 1, j)[0] < 255) {
            image.put(i - 1, j, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j);
        }
        if (image.get(i - 1, j + 1)[0] >= lowThreshold && image.get(i - 1, j + 1)[0] < 255) {
            image.put(i - 1, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j + 1);
        }
        if (image.get(i, j - 1)[0] >= lowThreshold && image.get(i, j - 1)[0] < 255) {
            image.put(i, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i, j - 1);
        }
        if (image.get(i, j + 1)[0] >= lowThreshold && image.get(i, j + 1)[0] < 255) {
            image.put(i, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i, j + 1);
        }
        if (image.get(i + 1, j - 1)[0] >= lowThreshold && image.get(i + 1, j - 1)[0] < 255) {
            image.put(i + 1, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j - 1);
        }
        if (image.get(i + 1, j)[0] >= lowThreshold && image.get(i + 1, j)[0] < 255) {
            image.put(i + 1, j, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j);
        }
        if (image.get(i + 1, j + 1)[0] >= lowThreshold && image.get(i + 1, j + 1)[0] < 255) {
            image.put(i + 1, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j + 1);
        }
    }

    public Mat Canny(Mat image, int size, double sigma, double lowThreshold, double highThreshold) throws InterruptedException {
        Mat tmp = SmartConverter.RGB2Gray((SmartGaussian.colorGaussianFilter(image, size, sigma)));
        SmartSobel ss = new SmartSobel();
        ss.compute(tmp);
        ss.convert();
        Mat ret = SmartNMS.NMS(ss.getGradientXY(), ss.getPointDirection());
        this.DoubleThreshold(ret, lowThreshold, highThreshold);
        this.DoubleThresholdLink(ret, lowThreshold);
        return ret;
    }

}

Reference

[1] 高斯滤波及高斯卷积核C++实现 https://blog.csdn.net/dcrmg/article/details/52304446

[2] 边缘检测之Canny https://www.cnblogs.com/techyan1990/p/7291771.html

[3] Canny边缘检测及C++实现 https://blog.csdn.net/dcrmg/article/details/52344902

[4] Canny边缘检测算法 https://zhuanlan.zhihu.com/p/42122107

[5] Sobel算子及C++实现 https://blog.csdn.net/dcrmg/article/details/52280768

深入理解 PBFT 协议

2023年2月12日 04:08

PBFT - Practical Byzantine Fault Tolerance

论文链接:https://pmg.csail.mit.edu/papers/osdi99.pdf

BFT

BFT (Byzantine Fault Tolerant) 拜占庭容错,指的是在一个分布式系统当中,在有节点可能出现任意行为(包括但不限于崩溃、发送错误数据、出现恶意行为等)的情况下,仍可以正常运行即达成共识。

Wikipedia Link.

PBFT 的重要概念

PBFT 可以保证在一个大小为 $3f+1$ 的分布式网络下,容忍至多 $f$ 个节点的拜占庭错误。

  • Primary (Leader): PBFT 是一个 leader-based 的协议,即在任意时刻,都有一个 leader 主导协议的走向。
  • 视图 View: 在 Primary 不发生改变的一段时间我们称之为 View。注意到:View 的编号是单调递增的。给定一个视图编号 $v$,我们就可以唯一确定这个视图的 Primary 编号 $p = v \bmod |\mathcal R|$,其中 $\mathcal R$ 是所有节点的集合,即所有 Replica 的集合。
  • Replica: 可以理解为服务器节点。
  • 目标: 达成 Consensus。即在每一个 Replica 上,所有的 transaction 内容和顺序都保证一致。
  • View Change: 视图改变。正因为协议由一个 Leader 主导,那在 leader 出现 malicious behaviour 的时候,协议可以保证触发 View Change,即更换 leader,来保证 liveness。具体内容后文会详细说明。

PBFT 的大致流程

Normal Case Operation 即 Primary 不出问题的情况下

  1. [REQUEST] 客户端向 Primary 发送 $\langle Request, o, t, c \rangle _{\sigma_c}$ 请求
  • $o$: Operation
  • $c$: Client
  • $t$: Timestamp - 这是用来保证执行的 exactly-once semantics
  1. [PREPREPARE] Primary 将请求 Multicast 给所有 Replica: $\langle \langle Preprepare, v, n, d \rangle _{\sigma_p}, m \rangle$
  • $v$: view
  • $n$: Primary 给这个请求分配的 Sequence Number 【关于 Sequence Number: 一个正常的 Replica 只会按顺序 (Sequence Number) 处理消息。即若其收到了 127 号消息但始终没有收到 125、126 号消息,127 号消息便一直不会处理。】
  • $d$: $d = D(m)$ 消息 $m$ 的 digest (摘要)
  1. [PREPARE] 每个 Replica $i$ 收到请求之后,如果选择接受请求(原文有更详细的描述),就会向所有 Replica 广播 Prepare: $\langle Prepare, v, n, d, i \rangle _ {\sigma_i}$
  2. [COMMIT] 接下来,如果每个 Replica $i$ 收到了其他 Replica (包括自己在内) 一共 $2f+1$ 个 Prepare 请求,这个时候认为其进入了 Prepared 状态,并广播 Commit 消息: $\langle Commit, v, n, d, i \rangle _ {\sigma_i}$

Checkpoint Protocol

见原文,不赘述。

View Change Protocol

  • 在怀疑 Primary 作恶的情况下,Replica 会广播 View Change 消息: $\langle ViewChange, v+1, n, \mathcal C, \mathcal P, i \rangle _{\sigma_i}$
    • $\mathcal C$ 是用来证明上一个 Checkpoint 的有效性的
    • $\mathcal P$ : 对于所有编号在 $n$ 之前,上一个 Checkpoint 之后的消息 $m$,$\mathcal P_m \in \mathcal P$,其中 $\mathcal P_m$ 包含了 $m$ 的 Preprepare 消息和 $2f+1$ Matching。
    • 注意到之所以包含了这么多消息,主要就是为了防止消息伪造,全部演算一遍数字签名,出现纰漏就直接忽略,拒不执行。
  • 如果一个 Replica 收到了 $f+1$ 个 View Change 请求,那他也会广播 View Change 请求,加速 View Change 进程。
  • New Primary 收到 $2f$ 个 View Change 请求之后(不包括自己),广播 $\langle NewView, v+1, \mathcal V, \mathcal O \rangle _{\sigma_i}$,具体的计算细节原文写的很详细了。

Fault 的分类

见:https://dinhtta.github.io/pbft/

Network Failure Model

Byzantine fault tolerant protocols tolerate at most $f$ Byzantine failures. In a distributed system, however, there are two types of failures: node and network. It is important to distinguish them, especially since they determine guarantees about the protocol’s safety and liveness.

Node failure means a node (or server, peer, entity, etc.) behaves arbitrarily. This captures the strongest adversary model. A network failure refers to network partitions which can last for an unbounded amount of time. It can be quantified as the number of nodes being isolated from the rest of the network, although they can still communicate with each other.

Given this distinction, what are being counted toward $f$? In PBFT

  • $f$ refers to node failures. The protocol guarantees safety upto $f$ node failures. Safety is independent of network failure. That is, even if the network is severely partition, namely more than $f$ nodes are isolated (but across all partitions there are still fewer than $f$ node failures), the protocol is still safe.
  • However, liveness is only guaranteed with fewer than $f$ total failures, i.e. counting both node and network failures. This means at least $2f+1$ nodes must be reachable.

Note that safety being independent of network failures is a strong guarantee, and it is common in most variants of PBFT. Not until Liu et al. recent work, namely XFT, is it relaxed in a way that separates out node and network failures.

为什么达成 Prepared 状态需要 $2f+1$ 个 Matching

Prepared 保证了什么?

要回答上面的大问题首先要理解我们想要 Preprepare + Prepare Phase 保证什么。论文 Section 4.2 已经回答了这个问题,即

$$
prepared(m, v, n, i) \Rightarrow \lnot prepared(m', v, n, j), \forall i,j,m,m', s.t. D(m') \neq D(m)
$$

即我们要保证在同一个 view 内,每一个 sequence number 只对应唯一一个 prepared message。

反例构造

首先我们可以构造一个如果只有 $2f$ 个 matching 就能打破上述原则的反例。

证明

接下来我们尝试用 Quorum 的思想来证明。使用反证法:如果每个节点都有 $2f+1$ 个 matching,但是最后仍然出现了两个 Replica $a \neq b$,使得

$$
prepared(m, v, n, a) \land prepared(m', v, n, b) \land D(m) \neq D(m')
$$

根据题设,有 $2f+1$ 个节点向 $a$ 发送了 $\langle Prepare, v, n, d, i \rangle _{\sigma_i}$,有 $2f+1$ 个节点向 $b$ 发送了 $\langle Prepare, v, n, d', i \rangle _{\sigma_i}$。根据容斥原理,至少有

$$
2 \times (2f + 1) - (3f+1) = f + 1
$$

个节点向 $a, b$ 分别发送了不同的消息。这样的行为毫无疑问是 malicious 的,所以至少存在 $f+1$ 个恶意节点,与至多有 $f$ 个恶意节点的协议假设矛盾。证毕。

Commit Phase 的必要性

Commit Phase 能够保证什么?

首先,Commit Phase 是 View Change 的副产物。后面的章节会详细解释 View Change,在这个阶段我们只需要知道 View Change 用来解决 Primary 作恶时的 Liveness 问题。而 Commit Phase 就是用来解决 View Change 过程中 Safety 的问题,即在 View 切换的前后:

  • 一个 Sequence Number 对应的 Message 必定在所有正常节点之间达成 Consensus
  • 不会出现一部分节点 Commit 了,但是另一部分节点完全不知道这件事,并继续执行下一个 Commit 或者重复使用 Sequence Number 的情况。
  • 更加明确一点就是原文中的表达:

$$
Committed(n, m, v, i) \Rightarrow \lnot Committed(n, m', v', j), \forall D(m) \neq D(m'), v' > v
$$

反例构造

现在我们假设 Commit Phase 不存在,即在节点达成 Prepared 的时候就直接给 Client 发送 Reply。考虑下图的情况:

可以发现反例非常容易构造,我们甚至构造出了 $f$ 个节点使用了 $n+1$ 的 Sequence Number,但是剩余节点以及 New Primary 没有的情况。

理解

从上面这个例子已经可以看出,Commit Phase 主要就是需要保障一个 Sequence Number 在跨 View 的时候不能被 Malicious Primary 或者 Network Partition 影响,使得被重复使用,破坏 Safety。

增加一个 Commit Phase 的主要功能,就是 只有在知道有一定数量的其他 Replica 已经处于 Prepared 状态之后再进行 Commit,这个数量一定要保证能够在 View Change 的时候至少有一个节点能够证明这个请求已经 Prepared。

因为注意到上图:在节点 Commit 的时候,很多节点没有达到 Prepared 状态,无法在 Network Partition 的时候提供证据向新的 Primary 证明这个请求已经被 Commit。综上所述,这就是 Commit Phase 存在的意义。

为什么达成 Committed 状态需要 $2f+1$ 个 Prepared

估算

上一节提到,在 Commit 前,我们至少要保证有足够多的 Replica 已经处于 Prepared 阶段,使得在 Network Partition 出现的时候仍然能够向 Primary 提供信息,告诉 Primary 一个请求已经被 Commit。那这个数量到底应该是多少呢?通过简单的估算其实就可以得到。

已知 New Primary 需要 $2f$ 个除了自己以外的 View Change 请求。最悲观地考虑,出现了最大限度的 Network Partition,即有 $f$ 个节点无法通信。在这个极限情况下的极限情况是这 $f$ 个 Replica 都处于 Prepared 状态却无法与 Primary 通信。而其余 $2f+1$ 个节点中又有 $f$ 个 malicious node,可以假装不知道自己 Prepared 了。

所以,在那 $2f+1$ 个节点中至少需要 $f+1$ 个节点处于 Prepared 状态才能保证 New Primary 收到 Commit 记录,再算上 Network Partition 的 $f$ 个 Replica,一共得有 $2f+1$ 个节点处于 Prepared 状态才能保证信息能够在 Network Partition 和 View Change 的情况下送达 New Primary。

反例构造

上述文字选取的是极端情况,有一条不成立即为反例。所以 $2f$ 个 Prepared 却仍然失败的反例和没有 Commit Phase 一样非常好构造,这里不再赘述。

为什么 New Primary 需要 $2f$ 个 View Change

理解了 Commit Phase 的 $2f+1$ 便不难发现,这里 $2f$ 这个数目的选取是和 Committed 数目的选取是相关的。

两个 $2f+1$ (实际上 View Change 算上 New Primary 也是 $2f+1$) 在总共 $3f+1$ 的池子里,Intersection (Quorum) 的大小总是不小于 $f+1$ 个 Replica 的。这样就保证了总有一个 Non-faulty Replica 位于最终的 Quorum 当中(faulty replica 可以隐瞒不报,装作不知道,所以一定需要 non-faulty replica),即既是 Prepared,又没有被 Network Partition 和 Malicious Replica 影响到。

为什么 Reply 需要 $f+1$ 个 Matching

这个就非常 Intuitive 了,因为这样必定能有一个 Non-faulty Replica 回复。

Checkpoint (Garbage Collection) 的必要性

注意到 View Change 的时候需要附上所有的记录来证明自己是正常节点,这所带来的开销随着时间的推移显然是不可接受的,所以我们需要想办法回收一部分 log。

Checkpoint 为什么在数字签名模式下只需要 $f+1$?

因为 $f+1$ 保证了至少有一个 Non-faulty Node 进行了 Checkpoint,也就是在没有 Network Partition 的情况下必定可以从网络中访问到这个 Checkpoint。

但是这里其实有一个问题:如果 Network Partition 了呢?注意到 PBFT 只在没有 Network Fault 的情况下 Guarantee Liveness,所以在那种情况下,需要请求 Checkpoint 的节点就 In-dark 了。如果这样的情况足够多,Liveness 也就无法保证了。

View Change 和 Liveness 之间的关系

View Change 主要是为了解决没有 Network Fault 的情况下的 Liveness 问题。当一个 Replica 的请求迟迟无法完成的时候就会 Trigger View Change。因为请求无法完成无非几种情况:

  • Primary 作恶,比方说只发给一部分 Replica Preprepare 消息
  • Network Partition

View Change 在第二种情况下其实并不能解决问题,但是可以解决第一种情况的问题。这就是 View Change 的 Intuition。

Primary 作恶带来的 In-Dark 问题

注意到 Primary 可以故意不给某些节点发请求,即 Primary 作恶。在这种情况下,这些节点既不会触发 View Change,也不会做任何事情,也就是字面意义上的 In-Dark。

由于前面的论证,PBFT 已经保证了 Safety,同时在没有 Network fault 的情况下也保证了 Liveness,所以这样的 In-Dark 并不会有什么特别大的影响。

只有当 In-Dark Replica 的数量多到使得 Consensus 无法达成了,才会 trigger View Change,使得协议能够推进下去。

Network Recovery 的流程

  • 首先向网络请求 Checkpoint
  • 如果能顺利 Catch Up,那就无缝衔接
  • 反之,迎来的结果必然是请求积压无法处理,timeout 之后便会 Request View Change,然后在 View Change 的过程中完成 Catch Up

Misc

  • Q: 如果 Primary 无中生有虚构信息或者篡改 Client 信息怎么办?
  • A: 注意到 $\langle Request, o, t, c \rangle _{\sigma_c}$ 的角标 $\sigma_c$,客户端签名保证了 Primary 不能篡改消息也不能伪造消息。
  • Q: 如果 Client 不知道当前的 Primary 怎么办?
  • A: 如果客户端一段时间内没有收到 Reply,触发 Timeout,重新 Multicast 请求。
  • Q: 如果因为 Network Partition 正常节点没有把回复送达到 Client,Client 认为请求失败,重传怎么办?这不就等于会执行两次了吗?
  • A: 注意到 $\langle Request, o, t, c \rangle _{\sigma_c}$ 中的 timestamp 可以解决这个问题

其他小的点:

  • View Change 后面 New View 这里应该也是有 Time out 机制的,否则 New Primary 迟迟不发 New View 的话 Liveness 就寄了
  • Client 发请求之后应该也有 Timeout 机制然后 Multicast,用来规避:(1)Primary 作恶,故意不发 Preprocess(2)发送的 Replica 不是 Primary

西瓜书 绪论 习题

2022年5月17日 19:44

注:本文为《机器学习》(周志华)第一章的作业笔记。

1.1

表 1.1 的编号为 1, 4 的两个样例:

色泽 根蒂 敲声 好瓜
青绿 蜷缩 浊响
乌黑 稍蜷 沉闷

因为只有两个样本,我们只要能够保证:

$$
s_好 \in S, s_坏 \notin S
$$

1.2

首先我们可以将决策看作一棵树,而每一个取式其实就是一棵子树或者叶子节点。会过来看合取式的析合范式,其实就是在一整棵树中取 $k$ 个互相不包含的子树或者叶子节点。

这里不能直接确定答案,因为题目并没有说决策树能不能给出形如下面的节点:

$$
(色泽\in\{青绿,乌黑\}; 根蒂=*; 敲声=*)
$$

我们先假设题目不允许上述的情况出现,那我们可以断言,包含 $n$ 个并列条件的决策,整棵树的深度为 $n + 1$,因为每一个父子关系其实就是析构 $*$ 和具体内容。

这个问题其实不能简单估算,但是我们可以用树形动态规划来解决。设状态 $f[u][p]$ 代表以节点 $u$ 为根的子决策树有 $p$ 个合取式的析合范式情况数。同时我们注意到书上的旁注:会有冗余的等价情况,那么我们还需要去除这样的等价情况,即在本题中我们需要去除

$$
\bigcup_{v \in u.sons} v = u
$$

的情况。

下面是递推方程:

$$
f[u][p] = \sum_{\sum_{1 \leq i \leq |u.sons|} p_i = p, p_i \geq 0} f[v_i][p_i]
$$

减去冗余情况:

$$
f[u][|u.sons|] := f[u][|u.sons|] - 1
$$

1.3

举例:线性回归。

1.4

$$
E_{ote} (\mathfrak L_a | X, f) = \sum_h \sum_{\boldsymbol x \in \mathcal X - X} P(\boldsymbol x) \ell (h(\boldsymbol x), f(\boldsymbol x)) P(h | X, \mathfrak L_a)
$$

$$
\begin{aligned}
\sum_f E_{ote} &= \sum_f \sum_h \sum_{\boldsymbol x \in \mathcal X - X} P(\boldsymbol x) \ell (h(\boldsymbol x), f(\boldsymbol x)) P(h | X, \mathfrak L_a) \\
&=\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x)\sum_h P(h | X, \mathfrak L_a) \sum_f \ell (h(\boldsymbol x), f(\boldsymbol x))
\end{aligned}
$$

注意到最后一部分为与 $\mathfrak L$ 无关的常数。不妨设

$$
L = \sum_f \ell (h(\boldsymbol x), f(\boldsymbol x))
$$

故原式等于

$$
L\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x)\sum_h P(h | X, \mathfrak L_a) = L\sum_{\boldsymbol x \in \mathcal X - X}P(\boldsymbol x) \cdot 1
$$

Q.E.D.

1.5

略。

P3572 [POI2014] PTA-Little Bird - DP 单调队列

2022年4月8日 23:50

题面

https://www.luogu.com.cn/problem/P3572

题解

首先,不难列出方程:$f[i]$ 代表跳到 $i$ 的最小疲劳值

$$
f[i] = \min_{i - k \leq j < i} \{f[j] + (a[i] \geq a[j])\}
$$

时间复杂度为 $O(qn^2)$,不理想。但是不难发现,这道题可以用单调队列进行优化。为什么?

如下状态一定是更优的:随着 index 递增

  • $f$ 不下降(第一维比较)
  • $f$ 若相同,同时有 $a$ 不上升(第二维比较)

所以我们就可以用单调队列来维护 index,时间复杂度为 $O(qn)$

注意:

  • 初值的设定,详见代码
  • 写的顺序和单调队列的模板稍微有一点不同,因为 $i - k \leq j < i$,右边的小于号不是小于等于号,即 $f[i]$ 不能被自己更新,所以调整了一下逻辑块的顺序。

代码

#include <iostream>

using namespace std;

const int maxn = 1e6+5;

int qq, n, k, head, tail, a[maxn], f[maxn], q[maxn];

int main() {
    cin >> n;
    for (int i = 1; i <= n; i ++) cin >> a[i];
    cin >> qq;
    while (qq --> 0) {
        cin >> k;
        head = 1; tail = 1; q[1] = 1;
        for (int i = 2; i <= n; i ++) {
            while (head <= tail && q[head] < i - k) head ++;
            f[i] = f[q[head]] + (a[i] >= a[q[head]]);
            while (head <= tail && (
                f[i] < f[q[tail]] || (
                    f[i] == f[q[tail]] && a[i] >= a[q[tail]]
                )
            )) tail --;
            q[++ tail] = i;
        }
        cout << f[n] << endl;
    }
    return 0;
}

POJ 2823 滑动窗口 单调队列 - 致逝去的青春

2022年4月8日 23:03

Preface

三年前的四月,第一次接触 OI。初三下,小教室里听着 fcx 讲着滑动窗口单调队列,放学后去本部蹭了第一节斜率优化的课,随即拉开了我长达两年的竞赛摆烂之旅。

记忆犹新。第一个注册的 OJ 账户就是 POJ,第一道 AC 的题便是 POJ 2823 滑动窗口。

多年过去,今天再看这个基础中的基础,茫然。只能说这是究极摆烂的下场罢。总之就好好写一下题解。

题面

https://www.luogu.com.cn/problem/P1886

题解

接下来只讨论求最小值的情况。设 $f[i]$ 代表终止于位置 $i$ 的滑动窗口的最小值:

$$
f[i] = \min_{i - k < j \leq i} a[j]
$$

如果暴力做的话时间复杂度为 $O(n^2)$。这里我们需要使用单调队列来优化。

对于最小值:单调队列内的元素值递增,同时 index 也是递增的。为什么这样做?思考一个情况:index 递增的情况下,如果一个 $a[i] > a[j], i < j$,那么在 $i,j$ 都有效的范围内,$a[i]$ 永远不可能作为答案。

这里有一句经典的话:如果一个选手年龄比你小,但却比你强,那你就可以退役了。我们把 index 比作出生年份,$a[]$ 比作能力,那么这句话就维护了一个求最大值的单调队列。

接下来是代码中需要注意的一些问题:

  • head, tail 不是左闭右开,而是 $[head, tail]$。这样做的原因是访问 headtail 的元素就可以简单的:q[head], q[tail]
  • 队列里存的是 index (pos)

注意:第二个 while 可以修改为 if,因为每次 i ++q 内的 pos 递增,写成 while 是习惯。

代码

#include <iostream>

using namespace std;

const int maxn = 1e6+5;
int n, k, a[maxn], q[maxn], head, tail;

int main() {
    cin >> n >> k;
    for (int i = 1; i <= n; i ++) cin >> a[i];
    // min
    head = 1, tail = 0;
    for (int i = 1; i <= n; i ++) {
        while (head <= tail && a[i] <= a[q[tail]]) tail --;
        q[++ tail] = i;
        while (head < tail && q[head] <= i - k) head ++;
        if (i >= k) cout << a[q[head]] << " ";
    } cout << endl;
    // max
    head = 1, tail = 0;
    for (int i = 1; i <= n; i ++) {
        while (head <= tail && a[i] >= a[q[tail]]) tail --;
        q[++ tail] = i;
        while (head < tail && q[head] <= i - k) head ++;
        if (i >= k) cout << a[q[head]] << " ";
    } cout << endl;
    return 0;
}

leetcode 1787 使所有区间的异或结果为零 - DP - 随机跳题计划

2022年4月7日 22:16

Preface

好耶,少见的难题。

题面

https://leetcode-cn.com/problems/make-the-xor-of-all-segments-equal-to-zero/

题解

大部分题解已经提到了:就是下面这个结论。

结论:最终构造出来的数列以 $k$ 为一个周期

证明:

由题设:

$$
a[1] \oplus a[2] \oplus \cdots \oplus a[k] = a[2] \oplus a[3] \oplus \cdots \oplus a[k + 1]
$$

所以我们有:

$$
(a[1] \oplus a[2] \oplus \cdots \oplus a[k]) \oplus (a[2] \oplus a[3] \oplus \cdots \oplus a[k + 1]) = 0
$$

根据异或的交换律与结合律:

$$
a[1] \oplus a[k + 1] = 0
$$

故 $a[1] = a[k + 1]$。易证,此结论可以将 $1$ 替换为任意定义域内的数。

接下来思考一下最终构造出来的数列满足什么特征:

  1. 以 $k$ 为一个周期
  2. 第一个周期异或和为 0

接下来看图说话:

在上图中 $n = 13, k = 3$,我们就是需要让

$$
a[1] = a[4] = a[7] = a[10] = a[13]
$$

$$
a[2] = a[5] = a[8] = a[11]
$$

$$
a[3] = a[6] = a[9] = a[12]
$$

$$
a[1] \oplus a[2] \oplus a[3] = 0
$$

那么我们不妨一组一组按列处理,接下来是状态的设法:$f[i][j]$ 代表处理完第 $i$ 列,使得前 $i$ 列的异或和为 $j$ 的最少修改次数。

思考转移方程:都是从上一列转移,但是有两种转移方式:

  1. 第 $i$ 列全部修改
  2. 第 $i$ 列部分修改,部分保留

设 $cnt$ 为第 $i$ 列的个数,$t[j]$ 代表第 $i$ 列 $j$ 出现的次数,那么:

$$
f[i][j] = \min (f[i][j], cnt + \min_{k} f[i - 1][k])
$$

$$
f[i][j] = \min_{k \in \text{column }i} \{f[i][j], f[i - 1][j \oplus k] + cnt - t[k]\}
$$

注意:全部修改仍然有可能比部分修改更优,这是因为本题不是简单的贪心。具体可以见下一小节,我们可以构造出反例,证明贪心是错误的。

此外,我们还需要注意初值的设定。

构造贪心的反例

贪心?怎么贪心?那当然是每一列选择相同数量最多的,再把这个数量排序,相同数量最少的一列全部改变,使得最后异或值为 0。

但是这个反例是非常容易构造的。

第一列:

$$
a, a, b, b, b, c, c
$$

第二列:

$$
d, d, d, d, e, e, e
$$

第三列:

$$
f,f,g,g,h,h,i
$$

按照之前说的贪心,我们会选保留 $b, d$,然后第三列适应答案使得异或为 $0$。但是如果很不巧,我们的构造在这样的情况下其他元素全部需要改变,那么此时的改动数量为 $4 + 5 + 7 = 16$。然而如果 $a \oplus e \oplus f = 0$,那么如果我们选择 $a, e, f$,改动的数量就是 $5 + 4 + 5 = 14 < 16$,更优。构造完毕。

代码

const int maxn = 2e3+5;
const int maxl = 10;

class Solution {

int f[maxn][1 << maxl];
int g[maxn], a[maxn], t[1 << maxl];

public:
    int minChanges(vector<int>& nums, int k) {
        int n = nums.size();
        for (int i = 0; i < n; i ++) a[i + 1] = nums[i];
        memset(f, 0x7f, sizeof(f));
        memset(g, 0x7f, sizeof(g));

        memset(t, 0, sizeof(t));
        int cnt = 0;
        for (int j = 1; j <= n; j += k)
            t[a[j]] ++, cnt ++;
        for (int j = 0; j < 1024; j ++) {
            f[1][j] = min(f[1][j], cnt - t[j]);
            g[1] = min(g[1], f[1][j]);
        }

        for (int i = 2; i <= k; i ++) {
            memset(t, 0, sizeof(t));
            cnt = 0;
            for (int j = i; j <= n; j += k)
                t[a[j]] ++, cnt ++;
            for (int j = 0; j < 1024; j ++) {
                f[i][j] = g[i - 1] + cnt;
                for (int m = 0; m < cnt; m ++) {
                    f[i][j] = min(f[i][j], f[i - 1][j ^ a[m * k + i]] + cnt - t[a[m * k + i]]);
                }
                g[i] = min(g[i], f[i][j]);
            }
        }
        return f[k][0];
    }
};

P2210 Haywire - 状压 DP

2022年4月7日 00:26

题面

https://www.luogu.com.cn/problem/P2210

题解

一开始没反应出来是状压,看了少数题解才知道可以状压。今天姑且 A 掉状压的写法,明天写写看模拟退火。

首先,我们来设状压的状态方程:

$$
f[i]
$$

  • $i$ 代表 $i$ 表示成二进制的时候有 $1$ 那些位所代表的牛组成的集合
  • 比方说 $i = 6 = (110)_2$ 代表有 $2, 3$ 号牛组成的集合
  • $f[i]$ 代表牛的集合为 $i$ 时对答案产生的贡献。

什么是贡献?这个说法很模棱两可。其实是这样的:两个朋友都在集合里面,那答案必然要包括这对朋友连线的贡献。但是如果只有一头牛在集合里呢?很简单,贡献就是:按照某种最佳排列时,这头牛所处的位置到集合队尾的距离。这么说很抽象,我们举个例子:

$$
~ * ~ * ~ * ~ \triangle ~ * ~ *
$$

上图中,$*$ 代表无所谓的某头牛,$\triangle$ 代表我们现在考察的单个朋友,在这个状态下,$\triangle$ 对答案的贡献就是 $2$,因为 $\triangle$ 之后还有两头牛.

而如果在这个队列之后紧跟就是 $\triangle$ 的朋友:另一个 $\triangle$,那么在增加前 $\triangle$ 给出的贡献 $2$ 就是这对朋友对答案的贡献了:

$$
~ * ~ * ~ * ~ \triangle ~ * ~ * ~ \triangle
$$

好。赘述完关于状态的定义,我们来思考状态怎么转移。对于集合 $i$,我们可以枚举位于队列尾端的到底是哪头牛。每次转移我们只需要加上原来的集合中孤立朋友的数量即可。最终给出的代码是做了一些小的优化,稍微难以理解,所以下面给出一个等价的形式,但是时间复杂度多一个 $n$,为 $O(n^22^n)$:

for (int i = 1; i < (1 << n); i ++) {
    for (int j = 0; j < n; j ++) {
        if ((1 << j) & i) {
            int pending_links = 0;
            int ii = i & ~(1 << j);
            for (int k = 0; k < n; k ++)
                if ((1 << k) & ii)
                    pending_links += 3 - 
                        (((ii >> nbr[k][0]) & 1) + 
                        ((ii >> nbr[k][1]) & 1) + 
                        ((ii >> nbr[k][2]) & 1));
            f[i] = min(f[i], f[ii] + pending_links);
        }
    }
}

代码

最终给出的代码时间复杂度为 $O(n2^n)$:

#include <iostream>
#include <cstring>

using namespace std;

const int maxn = 12;

int f[1 << maxn], nbr[maxn + 5][3];

int main() {
    int n; cin >> n;
    for (int i = 0; i < n; i ++) {
        cin >> nbr[i][0] >> nbr[i][1] >> nbr[i][2];
        nbr[i][0] --; nbr[i][1] --; nbr[i][2] --;
    }
    memset(f, 0x7f, sizeof(f));
    f[0] = 0;
    for (int i = 1; i < (1 << n); i ++) {
        int pending_links = 0;
        for (int j = 0; j < n; j ++)
            if ((1 << j) & i)
                pending_links += 3 - 
                    (((i >> nbr[j][0]) & 1) + 
                    ((i >> nbr[j][1]) & 1) + 
                    ((i >> nbr[j][2]) & 1));
        for (int j = 0; j < n; j ++) {
            if ((1 << j) & i) {
                f[i] = min(f[i], f[i & ~(1 << j)] + pending_links - (3 - 
                    (((i >> nbr[j][0]) & 1) + 
                    ((i >> nbr[j][1]) & 1) + 
                    ((i >> nbr[j][2]) & 1))) + 
                    ((i >> nbr[j][0]) & 1) + 
                    ((i >> nbr[j][1]) & 1) + 
                    ((i >> nbr[j][2]) & 1));
            }
        }
    }
    cout << f[(1 << n) - 1] << endl;
    return 0;
}

Canny边缘检测算法(基于OpenCV的Java实现)

2020年3月5日 00:21

Canny边缘检测算法(基于OpenCV的Java实现)

绪论

最近在学习ORB的过程中又仔细学习了Canny,故写下此篇笔记,以作总结。

Canny边缘检测算法的发展历史

Canny边缘检测于1986年由JOHN CANNY首次在论文《A Computational Approach to Edge Detection》中提出,就此拉开了Canny边缘检测算法的序幕。

Canny边缘检测是从不同视觉对象中提取有用的结构信息并大大减少要处理的数据量的一种技术,目前已广泛应用于各种计算机视觉系统。Canny发现,在不同视觉系统上对边缘检测的要求较为类似,因此,可以实现一种具有广泛应用意义的边缘检测技术。边缘检测的一般标准包括:

  • 以低的错误率检测边缘,也即意味着需要尽可能准确的捕获图像中尽可能多的边缘。
  • 检测到的边缘应精确定位在真实边缘的中心。
  • 图像中给定的边缘应只被标记一次,并且在可能的情况下,图像的噪声不应产生假的边缘。

为了满足这些要求,Canny使用了变分法。Canny检测器中的最优函数使用四个指数项的和来描述,它可以由高斯函数的一阶导数来近似。

在目前常用的边缘检测方法中,Canny边缘检测算法是具有严格定义的,可以提供良好可靠检测的方法之一。由于它具有满足边缘检测的三个标准和实现过程简单的优势,成为边缘检测最流行的算法之一。

Canny边缘检测算法的处理流程

Canny边缘检测算法可以分为以下5个步骤:

  • 使用高斯滤波器,以平滑图像,滤除噪声。
  • 计算图像中每个像素点的梯度强度和方向。
  • 应用非极大值(Non-Maximum Suppression)抑制,以消除边缘检测带来的杂散响应。
  • 应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。
  • 通过抑制孤立的弱边缘最终完成边缘检测。

下面详细介绍每一步的实现思路。

用高斯滤波器平滑图像

高斯滤波是一种线性平滑滤波,适用于消除高斯噪声,特别是对抑制或消除服从正态分布的噪声非常有效。滤波可以消除或降低图像中噪声的影响,使用高斯滤波器主要是基于在滤波降噪的同时也可以最大限度保留边缘信息的考虑。

高斯滤波实现步骤:

彩色RGB图像转换为灰度图像

边缘检测是基于对图像灰度差异运算实现的,所以如果输入的是RGB彩色图像,需要先进行灰度图的转换。
RGB转换成灰度图像的一个常用公式是:

$$
Gray = R*0.299 + G*0.587 + B*0.114
$$

注意一般情况下图像处理中彩色图像各分量的排列顺序是B、G、R。

RGB原图像:

转换后的灰度图:

Java代码调用系统库实现:

public static Mat RGB2Gray(Mat image) {
    // Gray = R*0.299 + G*0.587 + B*0.114
    Mat gray = new Mat();
    Imgproc.cvtColor(image, gray, Imgproc.COLOR_BGR2GRAY);
    return gray;
}

一维,二维高斯函数及分布

一维高斯函数表述为:

$$
G(x) = \frac {1}{\sqrt {2\pi}\sigma}\exp(-\frac {(x-\mu_x)^2}{2\sigma^2})
$$

对应图形:

二维高斯函数表述为:

$$
G(x) = \frac {1}{ {2\pi}\sigma^2}\exp(-\frac {(x-\mu_x)^2+(y-\mu_y)^2}{2\sigma^2})
$$

对应图形:

一些重要特性说明:

  1. 一维二维高斯函数中$μ$是服从正态分布的随机变量的均值,称为期望或均值影响正态分布的位置,实际的图像处理应用中一般取$μ=0;~σ$是标准差,$σ^2$是随机变量的方差,$σ$定义了正态分布数据的离散程度,$σ$越大,数据分布越分散,σ越小,数据分布越集中。

在图形或滤波效果上表现为:$σ$越大,曲线越扁平,高斯滤波器的频带就越宽,平滑程度就越好,$σ$越小,曲线越瘦高,高斯滤波的频带就越窄,平滑程度也越弱;

  1. 二维高斯函数具有旋转对称性,即滤波器在各个方向上的平滑程度是相同的.一般来说,一幅图像的边缘方向是事先不知道的,因此,在滤波前是无法确定一个方向上比另一方向上需要更多的平滑.旋转对称性意味着高斯平滑滤波器在后续边缘检测中不会偏向任一方向;
  2. 高斯函数是单值函数。这表明,高斯滤波器用像素邻域的加权均值来代替该点的像素值,而每一邻域像素点权值是随该点与中心点的距离单调增减的。这一性质是很重要的,因为边缘是一种图像局部特征,如果平滑运算对离算子中心很远的像素点仍然有很大作用,则平滑运算会使图像失真;
  3. 相同条件下,高斯卷积核的尺寸越大,图像的平滑效果越好,表现为图像越模糊,同时图像细节丢失的越多;尺寸越小,平滑效果越弱,图像细节丢失越少;

以下对比一下不同大小标准差$σ$(Sigma)对图像平滑的影响:

原图:

卷积核尺寸5*5,σ=0.1:

卷积核尺寸5*5,σ=1:

对比可以看到,Sigma(σ)越大,平滑效果越明显。

生成高斯滤波卷积核

滤波的主要目的是降噪,一般的图像处理算法都需要先进行降噪。而高斯滤波主要使图像变得平滑(模糊),同时也有可能增大了边缘的宽度。

高斯函数是一个类似与正态分布的中间大两边小的函数。

对于一个位置(m,n)的像素点,其灰度值(这里只考虑二值图)为f(m,n)。

那么经过高斯滤波后的灰度值将变为:

$$
g_\sigma(m,n)=\frac 1 {2\pi\sigma^2} \exp(-\frac {m^2+n^2}{2\sigma^2})f(m,n)
$$

其中,

$$
\frac 1 {{2\pi\sigma^2}} \exp(-\frac {m^2+n^2}{2\sigma^2})
$$

是二元高斯函数。

为了尽可能减少噪声对边缘检测结果的影响,所以必须滤除噪声以防止由噪声引起的错误检测。为了平滑图像,使用高斯滤波器与图像进行卷积,该步骤将平滑图像,以减少边缘检测器上明显的噪声影响。大小为(2k+1)x(2k+1)的高斯滤波器核的生成方程式由下式给出:

下面是一个sigma = 1.4,尺寸为3x3的高斯卷积核的例子(需要注意归一化):

若图像中一个3x3的窗口为A,要滤波的像素点为e,则经过高斯滤波之后,像素点e的亮度值为:

其中*为卷积符号,sum表示矩阵中所有元素相加求和。

重要的是需要理解,高斯卷积核大小的选择将影响Canny检测器的性能。尺寸越大,检测器对噪声的敏感度越低,但是边缘检测的定位误差也将略有增加。一般5x5是一个比较不错的trade off。

public static double[][] getGaussianArray(int size, double sigma) {
    double[][] array = new double[size][size];
    double center_i, center_j;
    if (size % 2 == 1) { center_i = (double) (size / 2); center_j = (double) (size / 2); }
    else { center_i = (double) (size / 2) - 0.5; center_j = (double) (size / 2) - 0.5; }
    double sum = 0.0;
    for (int i = 0; i < size; i ++)
            for (int j = 0; j < size; j ++) {
                array[i][j] = Math.exp((-1.0) * ((i - center_i) * (i - center_i) + (j - center_j) * (j - center_j)) / 2.0 / sigma / sigma);
                sum += array[i][j];
            }
    for (int i = 0; i < size; i ++)
            for (int j = 0; j < size; j++)
                array[i][j] /= sum;
    return array;
}

Sigma为1时,求得的3*3大小的高斯卷积核参数为:

Sigma为1,5*5大小的高斯卷积核参数为:

在计算出高斯滤波卷积核之后就是用这个卷积核来扫过整张图像,对每个像素点进行加权平均。

单色高斯滤波与彩色高斯滤波

加入了多线程的优化,代码实现:

public static Mat greyGaussianFilter(Mat src, double[][] array, int size) throws InterruptedException {
    Mat temp = src.clone();
    final CountDownLatch latch = new CountDownLatch(src.rows());
    for (int i = 0; i < src.rows(); i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 0; j < src.cols(); j ++)
                        if (finalI > (size / 2) - 1 && j > (size / 2) - 1 &&
                            finalI < src.rows() - (size / 2) && j < src.cols() - (size / 2)) {
                            double sum = 0.0;
                            for (int k = 0; k < size; k++)
                                    for (int l = 0; l < size; l++)
                                        sum += src.get(finalI - k + size / 2, j - l + size / 2)[0] * array[k][l];
                            temp.put(finalI, j, sum);
                        }
                latch.countDown();
            }).start();
    }
    latch.await();
    return temp;
}

public static Mat colorGaussianFilter(Mat src, int size, double sigma) throws InterruptedException {
    // return variable
    Mat ret = new Mat();
    // list for merge and split
    List<Mat> channels = new ArrayList<>();
    List<Mat> new_channels = new ArrayList<>();
    Map<Integer, Mat> temp_channels = new TreeMap<>();
    // split into 3 channels (r, g, b)
    Core.split(src, channels);
    // get gaussion array
    double[][] array = SmartGaussian.getGaussianArray(size, sigma);
    // multi-thread
    final CountDownLatch latch = new CountDownLatch(channels.size());
    channels.forEach(mat -> {
            new Thread(() -> {
                Mat tmp = new Mat();
                try {
                        tmp = SmartGaussian.greyGaussianFilter(mat, array, size);
                } catch (InterruptedException e) {
                        e.printStackTrace();
                }
                temp_channels.put(channels.indexOf(mat), tmp);
                latch.countDown();
            }).start();
    });
    latch.await();
    for (int i = 0; i < channels.size(); i ++)
          new_channels.add(temp_channels.get(i));
    Core.merge(new_channels, ret);
    return ret;
}

效果图:(左图为原图,中间为上面的实现,右边是OpenCV实现)

用Sobel等梯度算子计算梯度幅值和方向

梯度

梯度的本意是一个向量(矢量),表示某一函数在该点处的方向导数沿着该方向取得最大值,即函数在该点处沿着该方向(此梯度的方向)变化最快,变化率最大(为该梯度的模)。

设二元函数:

$$
z=f(x,y)
$$

在平面区域D上具有一阶连续偏导数,则对于每一个点$P(x,y)$都可定出一个向量:

$$
\Big\{\frac {\partial f}{\partial x},\frac {\partial f} {\partial y} \Big\}
=f_x(x,y)\vec i + f_y(x,y)\vec j
$$

该函数就称为函数$z=f(x,y)$在点$P(x,y)$的梯度,记作$\text{grad}~f(x,y)$或$\nabla f(x,y)$,即有:

$$
\text{grad}~f(x,y)=\nabla f(x,y)=\Big\{\frac {\partial f}{\partial x},\frac {\partial f} {\partial y} \Big\}
=f_x(x,y)\vec i + f_y(x,y)\vec j
$$

其中$\nabla=\displaystyle\frac{\partial}{\partial x}\vec i+\displaystyle\frac{\partial}{\partial y}\vec j$称为(二维的)向量微分算子或Nabla算子。

设$e=\{\cos\alpha,\cos\beta\}$是方向$l$上的单位向量,则

$$
\begin{aligned}
&\displaystyle\frac{\partial f}{\partial l} \\ =& \displaystyle\frac{\partial f}{\partial x} \cos \alpha + \displaystyle\frac{\partial f}{\partial y} \cos \beta \\ =& \Big \{\displaystyle\frac{\partial f}{\partial x},\displaystyle\frac{\partial f}{\partial y} \Big\} \{\cos\alpha, \cos \beta\} \\ =& \operatorname{grad}f(x,y)e \\ =& |\operatorname{grad}f(x,y)||e|\cos[\operatorname{grad}f(x,y),e]
\end{aligned}
$$

由于当方向$l$与梯度方向一致时,有$\cos[\operatorname{grad}f(x,y),e]=1$。

所以当$l$与梯度方向一致时,方向导数$\displaystyle\frac{\partial f}{\partial x}$有最大值,且最大值为梯度的模,即$|\operatorname{grad}f(x,y)| = \sqrt{\Big(\displaystyle\frac{\partial f}{\partial x}\Big)^2+\Big(\displaystyle\frac{\partial f}{\partial y}\Big)^2}$。

因此说,函数在一点沿梯度方向的变化率最大,最大值为该梯度的

图像灰度值的梯度的简单求法

图像灰度值的梯度可以使用最简单的一阶有限差分来进行近似,使用以下图像在x和y方向上偏导数的两个矩阵:

计算公式为:

其中$f$为图像灰度值,$P[i,j]$代表$[i,j]$在$X$方向梯度幅值,$Q[i,j]$代表$[i,j]$在$Y$方向的梯度幅值,$M[i,j]$是该点幅值,$\theta[i,j]$是梯度方向,也就是角度。

使用Sobel算子来计算梯度的大小及方向:

图像中的边缘可以指向各个方向,因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。边缘检测的算子(如Roberts,Prewitt,Sobel等)返回水平Gx和垂直Gy方向的一阶导数值,由此便可以确定像素点的梯度G和方向theta 。

其中G为梯度强度, theta表示梯度方向,arctan为反正切函数。下面以Sobel算子为例讲述如何计算梯度强度和方向。

x和y方向的Sobel算子分别为:

其中Sx表示x方向的Sobel算子,用于检测y方向的边缘; Sy表示y方向的Sobel算子,用于检测x方向的边缘(边缘方向和梯度方向垂直)。在直角坐标系中,Sobel算子的方向如下图所示。

若图像中一个3x3的窗口为A,要计算梯度的像素点为e,则和Sobel算子进行卷积之后,像素点e在x和y方向的梯度值分别为:

其中*为卷积符号,sum表示矩阵中所有元素相加求和。根据公式(3-2)便可以计算出像素点e的梯度和方向。

下面是Sobel算子求梯度的java实现:

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.Core;
import org.opencv.core.CvType;
import org.opencv.core.Mat;

import java.util.concurrent.CountDownLatch;

public class SmartSobel {

    private Mat gradientX = new Mat(), gradientY = new Mat(), gradientXY = new Mat();
    private double[][] pointDirection;

    public SmartSobel() {}

    public void compute(Mat image) throws InterruptedException {
        pointDirection = new double[image.rows()][image.cols()];
        for (int i = 0; i < image.rows(); i ++)
            for (int j = 0; j < image.cols(); j ++)
                pointDirection[i][j] = 0;
        gradientX = Mat.zeros(image.size(), CvType.CV_32SC1);
        gradientY = Mat.zeros(image.size(), CvType.CV_32SC1);
        gradientXY = Mat.zeros(image.size(), CvType.CV_32SC1);
        final CountDownLatch latch = new CountDownLatch(image.rows() - 2);
        for (int i = 1; i < image.rows() - 1; i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 1; j < image.cols() - 1; j++) {
                    double gX = (-1) * image.get(finalI - 1, j - 1)[0] +
                            1 * image.get(finalI - 1, j + 1)[0] +
                            (-2) * image.get(finalI, j - 1)[0] +
                            2 * image.get(finalI, j + 1)[0] +
                            (-1) * image.get(finalI + 1, j - 1)[0] +
                            1 * image.get(finalI + 1, j + 1)[0];
                    double gY = 1 * image.get(finalI - 1, j - 1)[0] +
                            2 * image.get(finalI - 1, j)[0] +
                            1 * image.get(finalI - 1, j + 1)[0] +
                            (-1) * image.get(finalI + 1, j - 1)[0] +
                            (-2) * image.get(finalI + 1, j)[0] +
                            (-1) * image.get(finalI + 1, j + 1)[0];
                    gradientY.put(finalI, j, Math.abs(gY));
                    gradientX.put(finalI, j, Math.abs(gX));
                    gradientXY.put(finalI, j, Math.sqrt(gX * gX + gY * gY));
                    // 防止除以0的情况发生
                    if (gX == 0) gX = 0.00000000000000001;
                    pointDirection[finalI][j] = Math.atan(gY / gX);
                }
                latch.countDown();
            }).start();
        }
        latch.await();
    }

    public void convert() {
        Core.convertScaleAbs(gradientX, gradientX);
        Core.convertScaleAbs(gradientY, gradientY);
        Core.convertScaleAbs(gradientXY, gradientXY);
    }

    public Mat getGradientX() { return this.gradientX; }

    public Mat getGradientY() { return this.gradientY; }

    public Mat getGradientXY() { return this.gradientXY; }

    public double[][] getPointDirection() { return this.pointDirection; }

}

由于这里使用的是$Math.tan()$,所以最终的角度是映射在$[-\frac \pi 2, \frac \pi 2]$的范围之内的。如果使用$Math.tan2()$会映射到$[-\pi,\pi]$的范围内,并且无需考虑符号影响,更加精确。但是这里我们并不关心另外的一个$\pi$的情况,我们只关心其所在直线(这在后文中会提到,也就是非极大值抑制),所以无需多考虑。

X方向梯度图:

Y方向梯度图:

X、Y方向梯度融合效果:

Opencv Sobel函数效果:

对梯度幅值进行非极大值抑制

非极大值抑制是一种边缘稀疏技术,非极大值抑制的作用在于“瘦”边。对图像进行梯度计算后,仅仅基于梯度值提取的边缘仍然很模糊。对于标准3,对边缘有且应当只有一个准确的响应。而非极大值抑制则可以帮助将局部最大值之外的所有梯度值抑制为0,对梯度图像中每个像素进行非极大值抑制的算法是:

  1. 将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。
  2. 如果当前像素的梯度强度与另外两个像素相比最大,则该像素点保留为边缘点,否则该像素点将被抑制。

通常为了更加精确的计算,在跨越梯度方向的两个相邻像素之间使用线性插值来得到要比较的像素梯度,现举例如下:

          图3-2 梯度方向分割

如图3-2所示,将梯度分为8个方向,分别为E、NE、N、NW、W、SW、S、SE,其中0代表$0^\circ\sim45^\circ$,1代表$45^\circ\sim90^\circ$,2代表$-90^\circ\sim-45^\circ$,3代表$-45^\circ\sim0^\circ$。像素点P的梯度方向为$\theta$,则像素点P1和P2的梯度线性插值为:

$$
\begin{aligned}
\tan \theta = G_y ~/~G_x \\
G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW
\end{aligned}
$$

上面也只是图中的情况,具体情况如下:

$$
\begin{aligned}
\theta \in [0, \frac \pi 4]:
\begin {cases}
G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW
\end {cases}\\\theta \in [\frac \pi 4, \frac \pi 2]:
\begin {cases}
G_{p1} = (1-\tan\theta)\times N + \tan\theta \times NE \\
G_{p2} = (1-\tan\theta)\times S + \tan\theta \times SW
\end {cases}\\\theta \in [-\frac \pi 4, 0]:
\begin {cases}
G_{p1} = (1-\tan(-\theta))\times E + \tan(-\theta) \times SE \\
G_{p2} = (1-\tan(-\theta))\times W + \tan(-\theta) \times NW
\end {cases}\\\theta \in [-\frac \pi 2, -\frac \pi 4]:
\begin {cases}
G_{p1} = (1-\tan(-\theta))\times S + \tan(-\theta) \times SE \\
G_{p2} = (1-\tan(-\theta))\times N + \tan(-\theta) \times NW
\end {cases}
\end{aligned}
$$

因此非极大值抑制的伪代码描写如下:

需要注意的是,如何标志方向并不重要,重要的是梯度方向的计算要和梯度算子的选取保持一致。

Java实现:

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.CvType;
import org.opencv.core.Mat;

import java.util.concurrent.CountDownLatch;

public class SmartNMS {

    public static Mat NMS(Mat gradientImage, double[][] pointDirection) throws InterruptedException {
        Mat outputImage = gradientImage.clone();
        final CountDownLatch latch = new CountDownLatch(gradientImage.rows() - 2);
        for (int i = 1; i < gradientImage.rows() - 1; i ++) {
            int finalI = i;
            new Thread(() -> {
                for (int j = 1; j < gradientImage.cols() - 1; j ++) {
                    double GP = gradientImage.get(finalI, j)[0],
                            E = gradientImage.get(finalI, j + 1)[0],
                            NE = gradientImage.get(finalI - 1, j + 1)[0],
                            N = gradientImage.get(finalI - 1, j)[0],
                            NW = gradientImage.get(finalI - 1, j - 1)[0],
                            W = gradientImage.get(finalI, j - 1)[0],
                            SW = gradientImage.get(finalI + 1, j - 1)[0],
                            S = gradientImage.get(finalI + 1, j)[0],
                            SE = gradientImage.get(finalI + 1, j + 1)[0];
                    double GP1 = 0, GP2 = 0;
                    double theta = pointDirection[finalI][j];
                    if (theta >= 0 && theta <= Math.PI / 4) {
                        GP1 = E * (1 - Math.tan(theta)) + NE * Math.tan(theta);
                        GP2 = W * (1 - Math.tan(theta)) + SW * Math.tan(theta);
                    } else if (theta > Math.PI / 4) {
                        GP1 = N * (1 - 1 / Math.tan(theta)) + NE * 1 / Math.tan(theta);
                        GP2 = S * (1 - 1 / Math.tan(theta)) + SW * 1 / Math.tan(theta);
                    } else if (theta < 0 && theta >= -Math.PI / 4) {
                        GP1 = E * (1 - Math.tan(-theta)) + SE * Math.tan(-theta);
                        GP2 = W * (1 - Math.tan(-theta)) + NW * Math.tan(-theta);
                    } else {
                        GP1 = S * (1 - 1 / Math.tan(-theta)) + SE * 1 / Math.tan(-theta);
                        GP2 = N * (1 - 1 / Math.tan(-theta)) + NW * 1 / Math.tan(-theta);
                    }
                    if (GP < GP1 || GP < GP2) outputImage.put(finalI, j, 0);
                }
                latch.countDown();
            }).start();
        }
        latch.await();
        return outputImage;
    }

}

双阈值检测

在施加非极大值抑制之后,剩余的像素可以更准确地表示图像中的实际边缘。然而,仍然存在由于噪声和颜色变化引起的一些边缘像素。为了解决这些杂散响应,必须用弱梯度值过滤边缘像素,并保留具有高梯度值的边缘像素,可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值,则将其标记为强边缘像素;如果边缘像素的梯度值小于高阈值并且大于低阈值,则将其标记为弱边缘像素;如果边缘像素的梯度值小于低阈值,则会被抑制。阈值的选择取决于给定输入图像的内容。

双阈值检测的伪代码描写如下:

抑制孤立低阈值点

到目前为止,被划分为强边缘的像素点已经被确定为边缘,因为它们是从图像中的真实边缘中提取出来的。然而,对于弱边缘像素,将会有一些争论,因为这些像素可以从真实边缘提取也可以是因噪声或颜色变化引起的。为了获得准确的结果,应该抑制由后者引起的弱边缘。通常,由真实边缘引起的弱边缘像素将连接到强边缘像素,而噪声响应未连接。为了跟踪边缘连接,通过查看弱边缘像素及其8个邻域像素,只要其中一个为强边缘像素,则该弱边缘点就可以保留为真实的边缘。

抑制孤立边缘点的伪代码描述如下:

实现:(使用递归)

package edu.sfls.Jeff.JavaDev.CVLib;

import org.opencv.core.Mat;

import java.util.ArrayList;
import java.util.List;

public class SmartCanny {

    private List<Integer[]> highPoints = new ArrayList<Integer[]>();

    private void DoubleThreshold(Mat image, double lowThreshold, double highThreshold) {
        for (int i = 1; i < image.rows() - 1; i ++)
            for (int j = 1; j < image.cols() - 1; j ++)
                if (image.get(i, j)[0] >= highThreshold) {
                    image.put(i, j, 255);
                    Integer[] p = new Integer[2];
                    p[0] = i; p[1] = j;
                    highPoints.add(p);
                } else if (image.get(i, j)[0] < lowThreshold)
                    image.put(i, j, 0);
    }

    private void DoubleThresholdLink(Mat image, double lowThreshold) {
        for (Integer[] p : highPoints) {
            DoubleThresholdLinkRecurrent(image, lowThreshold, p[0], p[1]);
        }
        for (int i = 1; i < image.rows() - 1; i ++)
            for (int j = 1; j < image.cols() - 1; j ++)
                if (image.get(i, j)[0] < 255)
                    image.put(i, j, 0);
    }

    private void DoubleThresholdLinkRecurrent(Mat image, double lowThreshold, int i, int j) {
        if (i <= 0 || j <= 0 || i >= image.rows() - 1 || j >= image.cols() - 1) return;
        if (image.get(i - 1, j - 1)[0] >= lowThreshold && image.get(i - 1, j - 1)[0] < 255) {
            image.put(i - 1, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j - 1);
        }
        if (image.get(i - 1, j)[0] >= lowThreshold && image.get(i - 1, j)[0] < 255) {
            image.put(i - 1, j, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j);
        }
        if (image.get(i - 1, j + 1)[0] >= lowThreshold && image.get(i - 1, j + 1)[0] < 255) {
            image.put(i - 1, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j + 1);
        }
        if (image.get(i, j - 1)[0] >= lowThreshold && image.get(i, j - 1)[0] < 255) {
            image.put(i, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i, j - 1);
        }
        if (image.get(i, j + 1)[0] >= lowThreshold && image.get(i, j + 1)[0] < 255) {
            image.put(i, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i, j + 1);
        }
        if (image.get(i + 1, j - 1)[0] >= lowThreshold && image.get(i + 1, j - 1)[0] < 255) {
            image.put(i + 1, j - 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j - 1);
        }
        if (image.get(i + 1, j)[0] >= lowThreshold && image.get(i + 1, j)[0] < 255) {
            image.put(i + 1, j, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j);
        }
        if (image.get(i + 1, j + 1)[0] >= lowThreshold && image.get(i + 1, j + 1)[0] < 255) {
            image.put(i + 1, j + 1, 255);
            DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j + 1);
        }
    }

    public Mat Canny(Mat image, int size, double sigma, double lowThreshold, double highThreshold) throws InterruptedException {
        Mat tmp = SmartConverter.RGB2Gray((SmartGaussian.colorGaussianFilter(image, size, sigma)));
        SmartSobel ss = new SmartSobel();
        ss.compute(tmp);
        ss.convert();
        Mat ret = SmartNMS.NMS(ss.getGradientXY(), ss.getPointDirection());
        this.DoubleThreshold(ret, lowThreshold, highThreshold);
        this.DoubleThresholdLink(ret, lowThreshold);
        return ret;
    }

}

Reference

[1] 高斯滤波及高斯卷积核C++实现 https://blog.csdn.net/dcrmg/article/details/52304446

[2] 边缘检测之Canny https://www.cnblogs.com/techyan1990/p/7291771.html

[3] Canny边缘检测及C++实现 https://blog.csdn.net/dcrmg/article/details/52344902

[4] Canny边缘检测算法 https://zhuanlan.zhihu.com/p/42122107

[5] Sobel算子及C++实现 https://blog.csdn.net/dcrmg/article/details/52280768

❌
❌