创建滤波器的拉普拉斯矩阵并求解图像滤波的线性方程

Question

为了过滤图像，我有一个优化问题要解决。
我创建了一个处理稀疏矩阵的问题的线性方程。

首先我会展示问题。
一、问题的拉普拉斯（邻接）矩阵：

矩阵 Dx / Dy 是前向差分算子 -> 因此它的转置是后向差分算子。
矩阵 Ax / Ay 是对角矩阵，其权重是图像梯度的函数（逐点，即该值仅取决于该像素本身的梯度）。
权重为：

其中 Ix(i) 是输入图像在第 i 个像素处的水平梯度（当对输入图像进行矢量化时）。

假设输入图像 G -> g = vec(G) = G(:)。
我想找到并成像 U -> u = vec(U) = U(:) s.t.:

我的问题是：

1. 如何有效地构建矩阵 Dx / Dy / Ax / Ay（它们都是稀疏的）？
2. 通过设置M = (I + \lambda * {L}_{g})，有没有优化的方式直接创建M？
3. 在 MATLAB 中解决此线性问题的最佳方法是什么？有没有办法绕过内存限制（即处理大图像并且仍然能够解决它）？
4. 有没有开源库可以在内存资源有限的情况下解决？任何带有 MATLAB API 的库？

谢谢。

Answer 1

鉴于您的 cmets，让我们在概要中回答每个问题，然后从那里开始：

1. 我将使用sparse 和其他相关函数回答下面的问题
2. 使用 (1)，我们绝对可以以优化的方式构建 M。
3. 简单地说，\ 运算符是求解逆时的最佳选择。 MathWorks 花了很多时间来优化它，而且它在后台几乎使用了 LAPACK 和 BLAS，如果你不使用它会发疯的。在 (4) 中回答了您唯一无法使用它的情况。
4. 有一些 MATLAB 脚本可以迭代地处理矩阵求解，例如连续过度松弛技术，但您应该只在内存不足时使用这些脚本（即如果 \ 没有给您答案）。由于矩阵的稀疏表示，这不应该（希望）发生，所以让我们暂时避免使用这些函数。

回到您的问题，我们可以很好地生成L_g 的稀疏表示。给定Dx 和Dy 的定义，我们可以使用eye 命令的sparse 版本，称为speye。因此，Dx 和 Dy 可以通过 Dx = diff(speye(size(inputImage))); 计算得到。例如，如果您尝试在 7 x 5 图像上执行此操作，将会产生这样的结果。

>> diff(speye(7,5))

ans =

(1,1)       -1
(1,2)        1
(2,2)       -1
(2,3)        1
(3,3)       -1
(3,4)        1
(4,4)       -1
(4,5)        1
(5,5)       -1

如您所见，我们仅引用非零条目。第 1 行第 1 列的系数为 -1，第 1 行第 2 列的系数为 1，依此类推。至于你的Ax 和Ay，这也很容易做到。我们有一个对角矩阵，我们可以手动设置每个条目。我们要做的就是指定一组行索引、列索引以及每个点的值。因此，我们可以这样做：

inputImage = im2double(inputImage); %//Important
rows = 1 : numel(inputImage); %// Assuming a 2D matrix
cols = rows; % // Row and column indices are the same
valuesDx = exp(-(gradX(rows).^2 / 2*sigma*sigma ));
valuesDy = exp(-(gradY(rows).^2 / 2*sigma*sigma ));

第一次调用的原因是因为我们要确保像素是双精度的，因为在 MATLAB 中求逆需要您这样做。它还确保我们不会溢出类型，因为我们正在规范 0 和 1 之间的强度。您可能需要调整标准偏差以反映这一点。现在我们只需要构造我们的Ax 和Ay 矩阵，我们把它和Dx 和Dy 放在一起：

numberElements = numel(inputImage);
Ax = sparse(rows, cols, valuesDx, numberElements, numberElements);
Ay = sparse(rows, cols, valuesDy, numberElements, numberElements);
identity = speye(numberElements, numberElements); 
Dx = diff(identity);
Dy = Dx.'; %// Transpose

我之所以将Dx 转置为Dy，是因为垂直方向上的差异运算符应该只是转置（对我来说很有意义）。这些都应该是您想要的每个矩阵的稀疏表示。矩阵运算也可以在稀疏矩阵上执行，包括乘法和逆。因此：

Lg = Dx.' * Ax * Dx + Dy.' * Ay * Dy;

您现在可以通过以下方式求解u：

u = (identity + lambda*Lg) \ g;

这假定g 的结构是您的图像中的像素采用column-major 格式。我对像素进行采样以构建 Ax 和 Ay 的方式自然遵循这一点。因此，请执行g = inputImage(:);，假设我们已转换为双精度并在 0 和 1 之间进行归一化。

当您最终解决 u 时，您可以通过以下操作将其重新整形为图像：

u = reshape(u, size(inputImage, 1), size(inputImage, 2));

u 也可能是稀疏的，所以如果你想要返回原始图像，请使用full() 投射它：

u = full(u);

希望这会有所帮助！