《最优化理论》期末作业2：像素聚类

maltadirk · 发表于 2022-1-11 20:39

问题描述

给定一幅高宽为 $w \times h$ 的图像，在其上随机撒上个不同的种子点 $\bm{p}_i$ 。将图像的每个像素分配给离各自最近的种子点，则这组种子点将图像分割成个区域。在每个区域采用最小二乘法求得最佳的逼近多项式 $g_i(x, y, a_i, b_i)$ ，则每个区域的逼近误差为：

$\sum_{x \in \Omega_i}\sum_{y \in \Omega_i} (f(x, y) - g_i(x, y, a_i, b_i))^2\tag{1}$

个区域（所有像素值）的总误差为：

$\sum_{i = 1}^m\sum_{x \in \Omega_i}\sum_{y \in \Omega_i} (f(x, y) - g_i(x, y, a_i, b_i))^2\tag{2}$
求：优化种子点位置以获得最优的像素聚类结果 (此时的区域边界和图像特征比较吻合)，等价于平均像素值的误差最小。

$e = \min\frac{1}{wh}\sum_{i = 1}^m\sum_{x \in \Omega_i}\sum_{y \in \Omega_i} (f(x, y) - g_i(x, y, a_i, b_i))^2\tag{3}$
要求：至少采用两种优化算法，每种优化算法测试至少两幅图像，即总共四个结果，自定。
问题分析

为什么要像素聚类？

反证推理：如果不像素聚类，那么仅仅通过一个最小二乘法拟合，几乎很难完成。每个区域都拟合一个最小二乘法，那么难度大大降低。通常，单个区域中像素函数关系也更好表达。这体现了一种分治的思想：分而治之，大问题变小问题，从而可以做得了。
像素的聚类的本质就是划为区域，从而检测出有效的图像边缘。

要点分析：

图像分为灰度图像和伪彩色图像，灰度图像是关于位置坐标的实函数，伪彩色图像是关于位置坐标的向量函数。为了简化求解过程，假定图像为灰度图像，即 $f: \mathrm{R}^2 \to \mathrm{R}$ 。
图像的像素值通常是离散的无符号整数，且在0~255范围内，通常需要标准化到稠密的实数。又因为使用均方误差，所以标准化到-1~1之间内。
位置坐标是离散的无符号整数，所以在判断属于哪个使用距离时，要标准化到稠密的0~1实数范围内。
$p_{i, x} \in [0, w]$ ， $p_{i, y} \in [0, h]$ 。在要点3的基础上，也标准化到0~1之间，即 $p_{i, x} \gets \frac{p_{i, x}}{w}, p_{i, y} \gets \frac{p_{i, y}}{h}$ 。
采用非线性最小二乘法。
当越大，所有像素值拟合的效果就越好，但是太大，区域又划分太多，不能分割出有效的图像边缘。通常，一幅图像就几条边缘。因此，的设置也很关键，不能太大，也不能太小。

为了求解方便，我们改写公式(3)，则有：

$e = \mathop{\arg \min}\limits_{\bm{P}, \bm{A}}\frac{1}{wh}\sum_{k = 1}^{wh}(f(\bm{x}^{(k)}) - h(\bm{x}^{(k)}, \bm{P}, \bm{A}))^2 \tag{4}$
其中， $\bm{x}^{(k)} \in \mathrm{R}^2$ ，  为  的矩阵， $\bm{A}$ 为  的矩阵。约束条件（已标准化）有：

$p_{1, i} \in [0, 1], p_{2, i} \in [0, 1]\tag{5}$
对公式(5)使用拉格朗日乘子法有：

$\sum_{i = 1}^m (s_i - r_i)p_{1, i} + \sum_{i = 1}^m(u_i - t_i)p_{2, i} - \sum_{i = 1}^m(s_i + u_i) \tag{6}$
构造  维向量 $\bm{r}, \bm{s}, \bm{t}, \bm{u}$ ，从而有：

$\left<\bm{s} - \bm{r}, \bm{p}_1^\mathrm{T}\right> + \left<\bm{u} - \bm{t}, \bm{p}_2^\mathrm{T}\right> -\sum(\bm{s} + \bm{u}) \tag{7}$
联加公式(4)和(7)，从而构造无约束优化：

$l = \mathop{\arg \min}\limits_{\bm{P}, \bm{A}, \bm{r}, \bm{s}, \bm{t}, \bm{u}}\frac{1}{wh}\sum_{k = 1}^{wh}(f(\bm{x}^{(k)}) - h(\bm{x}^{(k)}, \bm{P}, \bm{A}))^2 + \left<\bm{s} - \bm{r}, \bm{p}_1^\mathrm{T}\right> + \left<\bm{u} - \bm{t}, \bm{p}_2^\mathrm{T}\right> -\sum(\bm{s} + \bm{u}) \tag{8}$
为了减少优化符号，构造一个 $4 \times m$ 的矩阵 $\bm{M}$ ：

$\bm{m}_1^\mathrm{T} = \bm{s}, \bm{m}_2^\mathrm{T} = \bm{r}, \bm{m}_3^\mathrm{T} = \bm{u}, \bm{m}_4^\mathrm{T} = \bm{t}\tag{9}$
问题解决

接下来，只需优化无约束公式(9)了，可以采用两种方法：同时优化和交替优化（先固定  ，后优化  ;再固定，后优化；重复）。本文采用同时优化。
实验

代码实现要点：

公式(9)中的函数 $h$ 是很难写成解析表达式的，所以基于中心差分求梯度的数值解。基于中心差分求梯度对预设置的微小变量很敏感。简接验证了误差方向传播算法的优越性。
可视化，绘制出任意相邻区域 $\Omega_i \cap \Omega_j$ 中像素点，即图像边缘。像素聚类的本质就是所有像素都拟合得很好，所以拟合的函数 $h(\bm{x}, \bm{P}, \bm{A})$ 的输出分布相似于 $f(\bm{x})$ 的分布，可以近似说明像素聚类的效果好，因此采用比较两者的分布图像即可，简化可视化的难度。题目要求的可视化就留给读者自己了。
在要点2的基础上，就可以随机初始化一幅较小的图像，减少训练时间。
拟合函数需要对输出的像素值进行裁剪修正。
只能通过优化算法来改变待训练参数的值，其他地方只能使用待训练参数。

代码：
#! -*- coding:utf-8 -*-

&#39;&#39;&#39;
@Author:       ZM
@Date and Time: 2021/12/17 18:07
@File:       pixel_clustering.py
&#39;&#39;&#39;

import numpy as np
import matplotlib.pyplot as plt

def loss_fun(arg_X, arg_P, arg_A, arg_M):
&#39;&#39;&#39;
:param arg_X: 像素矩阵
:return: 无约束优化
&#39;&#39;&#39;

&#39;&#39;&#39;平均均方误差
&#39;&#39;&#39;
h, w = arg_X.shape
avg_errors = 0.
for i in range(w):
      for j in range(h):
         ds = np.sum(np.square(arg_P - np.array([i / w, j / h], dtype=arg_P.dtype)[:, None]), axis=0) # 求距离，不开平方
         min_indexs = np.where(ds == np.min(ds))[0]    # 所有最短距离的索引
         ls_pixel = 0
         for min_index in min_indexs:
            ls_pixel += (arg_A[0, min_index] - i / w) + (arg_A[1, min_index] - j / h) ** 2 # 非线性最小二乘法拟合的像素值
         ls_pixel /= min_indexs.size
         avg_errors += (arg_X[j, i] - ls_pixel) ** 2    # 均方误差
avg_errors /= arg_X.size # 平均均方误差

&#39;&#39;&#39; 拉格朗日乘子项
&#39;&#39;&#39;
Lagrange = (np.einsum(&#39;i, i&#39;, arg_M[0] - arg_M[1], arg_P[0]) +
            np.einsum(&#39;i, i&#39;, arg_M[2] - arg_M[3], arg_P[1]) -
            np.sum(arg_M[0, 2]))

return 0.999 * avg_errors + 0.001 * Lagrange    # 无约束，系数分配让所有像素值拟合的效果更好一点

def grad_fun(arg_X, arg_P, arg_A, arg_M, f, epsilon=1e-14): # 基于中心差分求梯度对epsilon敏感
&#39;&#39;&#39;
:param arg_X: 像素矩阵
基于中心差分求梯度
&#39;&#39;&#39;

grad_P = np.zeros_like(arg_P)
grad_A = np.zeros_like(arg_A)
grad_M = np.zeros_like(arg_M)

for i in range(arg_P.shape[0]):
      for j in range(arg_P.shape[1]):
         t = arg_P[i, j]
         arg_P[i, j] = t + epsilon
         e1 = f(arg_X, arg_P, arg_A, arg_M)
         arg_P[i, j] = t - epsilon
         e2 = f(arg_X, arg_P, arg_A, arg_M)
         grad_P[i, j] = (e1 - e2) / (2 * epsilon)
         arg_P[i, j] = t

for i in range(arg_A.shape[0]):
      for j in range(arg_A.shape[1]):
         t = arg_A[i, j]
         arg_A[i, j] = t + epsilon
         e1 = f(arg_X, arg_P, arg_A, arg_M)
         arg_A[i, j] = t - epsilon
         e2 = f(arg_X, arg_P, arg_A, arg_M)
         grad_A[i, j] = (e1 - e2) / (2 * epsilon)
         arg_A[i, j] = t

for i in range(arg_M.shape[0]):
      for j in range(arg_M.shape[1]):
         t = arg_M[i, j]
         arg_M[i, j] = t + epsilon
         e1 = f(arg_X, arg_P, arg_M, arg_M)
         arg_M[i, j] = t - epsilon
         e2 = f(arg_X, arg_P, arg_M, arg_M)
         grad_M[i, j] = (e1 - e2) / (2 * epsilon)
         arg_M[i, j] = t

return grad_P, grad_A, grad_M

def fitting_fun(i, j, arg_P, arg_A):
&#39;&#39;&#39;
:param i: 横坐标
:param j: 纵坐标
:return: 像素值
&#39;&#39;&#39;
ds = np.sum(np.square(arg_P - np.array([i / w, j / h], dtype=arg_P.dtype)[:, None]), axis=0)  # 求距离，不开平方
min_indexs = np.where(ds == np.min(ds))[0]  # 所有最短距离的索引
ls_pixel = 0
for min_index in min_indexs:
      ls_pixel += (arg_A[0, min_index] - i / w) + (arg_A[1, min_index] - j / h) ** 2  # 非线性最小二乘法拟合的像素值
ls_pixel /= min_indexs.size
return max(min(ls_pixel, 1), 0)    # 裁剪修正

h, w = 10, 10
X = np.random.uniform(size=(h, w)) + np.random.random(size=(h, w)) + np.random.normal(size=(h, w))  # 随机初始化一幅图像
X = 2 * (np.max(X) - X) / (np.max(X) - np.min(X)) - 1 # 将像素值标准化到-1~1范围内
m = int(np.sqrt(h * w))    # 种子点的个数,既不能太大，也不能太小

&#39;&#39;&#39;待优化的参数
&#39;&#39;&#39;
Px = np.random.uniform(size=(1, m))    # 标准化到0~1范围内
Py = np.random.uniform(size=(1, m))    # 标准化到0~1范围内
P = np.concatenate([Px, Py], axis=0)
A = np.random.uniform(low=-np.sqrt(6 / m), high=np.sqrt(6 / m), size=(2, m))
M = np.random.uniform(low=-np.sqrt(6 / m), high=np.sqrt(6 / m), size=(4, m))

&#39;&#39;&#39; P的初始输入
&#39;&#39;&#39;
s = [f&#39;{P[0, i] * w:.2f}, {P[1, i] * h:.2f}\n&#39; for i in range(m)]
s[-1] = s[-1].strip()
with open(&#39;logs_init.txt&#39;, &#39;w&#39;, encoding=&#39;utf-8&#39;) as f:
f.writelines(s)

total_steps = 5000    # 总迭代次数
lr = init_lr = 0.01 # 学习步长
min_loss = 1e32

&#39;&#39;&#39;采用牛顿加速梯度下降法
&#39;&#39;&#39;
vP, vA, vM = np.zeros_like(P), np.zeros_like(A), np.zeros_like(M)
momentum = 0.9
for step in range(total_steps):
&#39;&#39;&#39;只能通过优化算法来改变带训练参数的值，其他地方只能使用待训练参数
&#39;&#39;&#39;
grad_P, grad_A, grad_M = grad_fun(X, P, A, M, loss_fun)
vP = momentum * vP + grad_P
vA = momentum * vA + grad_A
vM = momentum * vM + grad_M
P -= lr * vP
A -= lr * vA
M -= lr * vM

tmp_val = loss_fun(X, P, A, M)
print(f&#39;{lr=}, {tmp_val=}&#39;) # 输出每一步的训练日志

if min_loss > tmp_val:
      min_loss = tmp_val

if step < 2000:
      lr = init_lr
elif step < 3000:
      lr = init_lr * 0.1
elif step < 4000:
      lr = init_lr * 0.01
else:
      lr = init_lr * 0.001

&#39;&#39;&#39; P的输出
&#39;&#39;&#39;
P = np.clip(P, 0, 1) # 裁剪修正
s = [f&#39;{P[0, i] * w:.2f}, {P[1, i] * h:.2f}\n&#39; for i in range(m)]
s[-1] = s[-1].strip()
with open(&#39;logs_clip.txt&#39;, &#39;w&#39;, encoding=&#39;utf-8&#39;) as f:
f.writelines(s)

&#39;&#39;&#39;可视化
&#39;&#39;&#39;
fitting_X = np.empty_like(X)
for i in range(w):
for j in range(h):
      fitting_X[j, i] = fitting_fun(i, j, P, A)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(X.flatten(), &#39;o&#39;, label=&#39;Real fun&#39;)
ax.plot(fitting_X.flatten(), &#39;d&#39;, label=&#39;Fitting fun&#39;)
for i in range(m):
ax.plot(np.array([P[0, i] * P[1, i]] * 1000), np.random.uniform(size=1000))
ax.legend()
fig.savefig(&#39;rst.png&#39;)
plt.close()图1表明像素值拟合的效果还是挺好的。

图1 使用拉格朗日约束：像素值分布比较。

使用拉格朗日约束：约束条件的训练结果并没有出现过多  ，表明拉格朗日乘子项起到了约束作用。

的初始输入：
3.46, 4.82
2.04, 7.73
2.25, 4.06
6.91, 9.51
4.87, 5.49
6.86, 7.10
0.74, 9.39
3.15, 0.59
2.14, 7.19
4.83, 3.18
的输出：
4.18, 7.01
0.00, 2.76
6.92, 2.00
3.66, 3.72
8.52, 6.82
3.85, 2.71
2.23, 4.81
5.88, 2.92
8.57, 0.00
2.55, 10.00

由于拉格朗日乘子法引入过多参数，本文还尝试了只用梯度裁剪来防止优化的超出的可行解。而且，公式(9)是公式(5)的必要非充分条件，所以为了防止优化的 超出 的可行解，还可以组合拉格朗日乘子法和梯度裁剪。
图2表明梯度裁剪的拟合效果更好一点。

图2 使用梯度裁剪：像素值分布比较。

使用梯度裁剪：约束条件的训练结果也没有出现过多，表明梯度裁剪也起到了约束作用。

的初始输入：
4.63, 9.90
6.64, 9.96
9.56, 0.89
2.52, 2.07
2.11, 7.17
4.05, 7.90
0.34, 1.28
4.03, 5.24
7.34, 8.28
4.85, 7.07
的输出：
4.63, 9.90
6.64, 9.96
9.56, 0.89
2.52, 2.07
2.11, 7.17
4.05, 7.90
0.34, 1.28
4.03, 5.24
7.34, 8.28
4.85, 7.07

组合拉格朗日乘子法和梯度裁剪。

的输入：
4.49, 1.98
5.70, 8.08
9.49, 3.07
3.08, 1.39
7.77, 4.90
8.09, 2.79
8.50, 0.83
7.24, 8.99
7.86, 3.21
4.20, 1.82
的输出：
4.41, 0.00
6.72, 9.50
9.63, 3.13
3.74, 1.19
7.37, 7.12
10.00, 4.56
6.76, 0.00
8.05, 8.97
8.16, 0.35
5.20, 1.55

		自动登录	找回密码
密码			立即注册

《最优化理论》期末作业2：像素聚类

本帖子中包含更多资源