如何在matlab中进行快速卷积小补丁
how to perform fast convolution small patches in matlab
我有大量的二维图像块和一组二维滤镜。它们的大小相同。当我应用大小为 dxd 的 2D 滤镜时，我想获得图像块 (dxd) 的滤镜响应 (1x1)。
最快的方法是什么？
图像块存储为 3D 矩阵。补丁和过滤器的大小为 15x15。我在 matlab 中试过 conv 但似乎很慢。 （5000 个补丁大约需要 20-30 秒）。
我不能通过一些矩阵乘法来实现吗？我不能将每个补丁和每个过滤器转换为一个 D 向量并进行乘法运算吗？
这取决于你有什么工具。据我所知，没有批处理函数，所以你必须在 for 循环中进行。但是，如果您安装了 Parallel Computing Toolbox，则可以使用 parfor 而不是 for，这应该会提高性能。
不要使用 conv2。这是为一般二维信号设计的。但是，根据您的描述，您在 3D 矩阵中存储了 15 x 15 图像块，以及一个 15 x 15 过滤器。 imfilter 专门用于过滤图像，但这对您也没有帮助，因为仅支持 full 或 same 输出。我的猜测是，在使用 conv2 时，您使用的是 valid 标志，因此这将减少到 1 x 1 输出。
我建议你做的是展开你的图像补丁，这样你就可以创建一个二维矩阵。每个补丁适合一列，所有图像补丁按列连接以创建 2D 矩阵。您还可以展开过滤器，使其适合单个列，然后使用 bsxfun 实现对每个补丁的过滤。
因此，假设您的图像补丁存储在 A 中并且您的滤镜存储在 K 中：
B = reshape(A, 15*15, []);
C = sum(bsxfun(@times, B, K(:)), 1);
B 将您的 3D 矩阵重塑为 2D 矩阵，其中每一列代表一个图像块。接下来，C 是每个带有过滤器 K 的图像块的输出响应。 bsxfun 允许您使用展开到列向量中的过滤器，并且它会为与 B 中一样多的列复制过滤器。这会产生另一个与 B 大小相同的二维矩阵，但每一列都是 相同的 （即 K）。使用此矩阵，我们进行逐元素乘法，然后沿行求和以获得滤波器的输出响应。因此，C 的每个元素都会为您提供图像块对感兴趣的特定滤镜的响应。
请注意，此代码假设已执行关联。如果 2D 滤波器是对称的，则卷积与相关性相同。但是，如果要进行卷积并且滤波器不是对称的，则需要在应用卷积之前翻转每个维度。因此，您必须为每个过滤器执行此操作：
K = K(end:-1:1, end:-1:1);
然后您将应用以上两行代码来实现过滤。
因为您有多个过滤器，所以您可以使用 for 循环为您拥有的每个过滤器执行上述操作。您当然可以对其进行矢量化并使用 bsxfun 方法，但我的直觉是这会更慢。我认为最好单独处理每个滤镜，而不是试图为每个图像补丁找到所有滤镜的响应。因此，如果您的过滤器存储在 3D 矩阵中...我们将再次称它为 K，您可以这样做：
B = reshape(A, 15*15, []);
num_patches = size(B, 2);
num_filters = size(K, 3);
C = zeros(num_filters, num_patches);
for idx = 1 : num_filters
filt = K(:,:,idx);
%filt = filt(end:-1:1, end:-1:1); %// Do this if the filter is not symmetric
C(idx,:) = sum(bsxfun(@times, B, filt(:)), 1);
end
此处，C 将是一个二维矩阵，其中每个 行 表示每个图像块对特定过滤器的响应。
另一种可能的方法是使用 blockproc 函数。
由于图像块和过滤器的大小相同，并且您似乎想要一个 'valid' 卷积，这减少到简单地将图像块和过滤器的相应元素 (.*) 相乘，然后求和结果。这可以通过简单的 sum(sum(image_patch.*filter_patch,1),2).
来完成
设置一些数据：
numimg = 1000;
numfilt = 100;
img = rand(15,15,numimg); % 1000 15x15 image patches
filt = rand(15,15,numfilt); % 100 15x15 filters
为了使用blockproc，我们首先将 3-D 图像展平
img_flat = reshape(img,15,[]);
为blockproc定义一个处理函数。
filter_func = @(block_struct) squeeze( sum(sum(repmat(block_struct.data,1,1,numfilt) .* filt,1),2));
处理使用 blockproc.
out = blockproc(img_flat,[15 15],filter_func);
结果是一个 numfilt x numimg 矩阵，每个图像块上的每个过滤器都有过滤器输出。
在我的机器上，运行 用了不到一秒的时间。如果 numfilt 或 numimg 大得多，您可以考虑打开 blockproc 的 'useParallel' 标志。
这是对@rayryeng 回答的跟进。我（大部分）采用相同的表示法。
可以应用单个矩阵乘法来解决这个问题。这是可能的原因是，要获得图像 I 和滤波器 F 的滤波器响应，我们只需取 I 和 Fflipped 的标量积，我们在这里定义 Fflipped = F(end:-1:1,end:-1:1);。该标量积可以通过多种方式计算，例如使用：sum(Fflipped(:).*I(:))，但使用 transpose(Fflipped(:))*I(:) 的矩阵乘法更简单。这样做的最大好处是，您可以通过以正确的方式堆叠向量来一次计算多个标量积：
M = [V(:,1).'; V(:,2).'; V(:,3).'; V(:,4).'] * [W(:,1) W(:,2) W(:,3)]
将为您提供 V(:,i) 和 W(:,j).
组合的所有标量积的矩阵 M(i,j) == dot(V(:,i), W(:,j))
所以，首先我们通过以下方式翻转所有过滤器：
K = K(end:-1:1,end:-1:1,:);
然后我们重塑图像的 3D 数组并将过滤器翻转为将按列存储数据的 2D 数组，正如我们已经看到的那样。
psize = size(A,1); % Patch size
A = reshape(A, psize^2, []); % Images
K = reshape(K, psize^2, []); % Flipped filters
现在我们转置过滤器并通过使用单个矩阵乘法来执行我们所有的标量乘积以获得我们的响应矩阵：
C = K.'*A;
这将实现以下功能：
function Responses = getFilterResponses(Images, Filters)
Filters = Filters(end:-1:1,end:-1:1,:); %// Flip filters
psize = size(Images,1); %// Patch size
Images = reshape(Images, psize^2, []); %// Images
Filters = reshape(Filters, psize^2, []); %// Flipped filters
%%// Compute all scalar products via nifty matrix multiplication trick
Responses = Filters.'*Images;
对于 5000 15-by-15 个过滤器，这将比循环过滤器快大约 25 倍。
我有大量的二维图像块和一组二维滤镜。它们的大小相同。当我应用大小为 dxd 的 2D 滤镜时，我想获得图像块 (dxd) 的滤镜响应 (1x1)。
最快的方法是什么？
图像块存储为 3D 矩阵。补丁和过滤器的大小为 15x15。我在 matlab 中试过 conv 但似乎很慢。 （5000 个补丁大约需要 20-30 秒）。
我不能通过一些矩阵乘法来实现吗？我不能将每个补丁和每个过滤器转换为一个 D 向量并进行乘法运算吗？
这取决于你有什么工具。据我所知，没有批处理函数，所以你必须在 for 循环中进行。但是，如果您安装了 Parallel Computing Toolbox，则可以使用 parfor 而不是 for，这应该会提高性能。
不要使用 conv2。这是为一般二维信号设计的。但是，根据您的描述，您在 3D 矩阵中存储了 15 x 15 图像块，以及一个 15 x 15 过滤器。 imfilter 专门用于过滤图像，但这对您也没有帮助，因为仅支持 full 或 same 输出。我的猜测是，在使用 conv2 时，您使用的是 valid 标志，因此这将减少到 1 x 1 输出。
我建议你做的是展开你的图像补丁，这样你就可以创建一个二维矩阵。每个补丁适合一列，所有图像补丁按列连接以创建 2D 矩阵。您还可以展开过滤器，使其适合单个列，然后使用 bsxfun 实现对每个补丁的过滤。
因此，假设您的图像补丁存储在 A 中并且您的滤镜存储在 K 中：
B = reshape(A, 15*15, []);
C = sum(bsxfun(@times, B, K(:)), 1);
B 将您的 3D 矩阵重塑为 2D 矩阵，其中每一列代表一个图像块。接下来，C 是每个带有过滤器 K 的图像块的输出响应。 bsxfun 允许您使用展开到列向量中的过滤器，并且它会为与 B 中一样多的列复制过滤器。这会产生另一个与 B 大小相同的二维矩阵，但每一列都是 相同的 （即 K）。使用此矩阵，我们进行逐元素乘法，然后沿行求和以获得滤波器的输出响应。因此，C 的每个元素都会为您提供图像块对感兴趣的特定滤镜的响应。
请注意，此代码假设已执行关联。如果 2D 滤波器是对称的，则卷积与相关性相同。但是，如果要进行卷积并且滤波器不是对称的，则需要在应用卷积之前翻转每个维度。因此，您必须为每个过滤器执行此操作：
K = K(end:-1:1, end:-1:1);
然后您将应用以上两行代码来实现过滤。
因为您有多个过滤器，所以您可以使用 for 循环为您拥有的每个过滤器执行上述操作。您当然可以对其进行矢量化并使用 bsxfun 方法，但我的直觉是这会更慢。我认为最好单独处理每个滤镜，而不是试图为每个图像补丁找到所有滤镜的响应。因此，如果您的过滤器存储在 3D 矩阵中...我们将再次称它为 K，您可以这样做：
B = reshape(A, 15*15, []);
num_patches = size(B, 2);
num_filters = size(K, 3);
C = zeros(num_filters, num_patches);
for idx = 1 : num_filters
filt = K(:,:,idx);
%filt = filt(end:-1:1, end:-1:1); %// Do this if the filter is not symmetric
C(idx,:) = sum(bsxfun(@times, B, filt(:)), 1);
end
此处，C 将是一个二维矩阵，其中每个 行 表示每个图像块对特定过滤器的响应。
另一种可能的方法是使用 blockproc 函数。
由于图像块和过滤器的大小相同，并且您似乎想要一个 'valid' 卷积，这减少到简单地将图像块和过滤器的相应元素 (.*) 相乘，然后求和结果。这可以通过简单的 sum(sum(image_patch.*filter_patch,1),2).
设置一些数据：
numimg = 1000;
numfilt = 100;
img = rand(15,15,numimg); % 1000 15x15 image patches
filt = rand(15,15,numfilt); % 100 15x15 filters
为了使用blockproc，我们首先将 3-D 图像展平
img_flat = reshape(img,15,[]);
为blockproc定义一个处理函数。
filter_func = @(block_struct) squeeze( sum(sum(repmat(block_struct.data,1,1,numfilt) .* filt,1),2));
处理使用 blockproc.
out = blockproc(img_flat,[15 15],filter_func);
结果是一个 numfilt x numimg 矩阵，每个图像块上的每个过滤器都有过滤器输出。
在我的机器上，运行 用了不到一秒的时间。如果 numfilt 或 numimg 大得多，您可以考虑打开 blockproc 的 'useParallel' 标志。
这是对@rayryeng 回答的跟进。我（大部分）采用相同的表示法。
可以应用单个矩阵乘法来解决这个问题。这是可能的原因是，要获得图像 I 和滤波器 F 的滤波器响应，我们只需取 I 和 Fflipped 的标量积，我们在这里定义 Fflipped = F(end:-1:1,end:-1:1);。该标量积可以通过多种方式计算，例如使用：sum(Fflipped(:).*I(:))，但使用 transpose(Fflipped(:))*I(:) 的矩阵乘法更简单。这样做的最大好处是，您可以通过以正确的方式堆叠向量来一次计算多个标量积：
M = [V(:,1).'; V(:,2).'; V(:,3).'; V(:,4).'] * [W(:,1) W(:,2) W(:,3)]
将为您提供 V(:,i) 和 W(:,j).
M(i,j) == dot(V(:,i), W(:,j))
所以，首先我们通过以下方式翻转所有过滤器：
K = K(end:-1:1,end:-1:1,:);
然后我们重塑图像的 3D 数组并将过滤器翻转为将按列存储数据的 2D 数组，正如我们已经看到的那样。
psize = size(A,1); % Patch size
A = reshape(A, psize^2, []); % Images
K = reshape(K, psize^2, []); % Flipped filters
现在我们转置过滤器并通过使用单个矩阵乘法来执行我们所有的标量乘积以获得我们的响应矩阵：
C = K.'*A;
这将实现以下功能：
function Responses = getFilterResponses(Images, Filters)
Filters = Filters(end:-1:1,end:-1:1,:); %// Flip filters
psize = size(Images,1); %// Patch size
Images = reshape(Images, psize^2, []); %// Images
Filters = reshape(Filters, psize^2, []); %// Flipped filters
%%// Compute all scalar products via nifty matrix multiplication trick
Responses = Filters.'*Images;
对于 5000 15-by-15 个过滤器，这将比循环过滤器快大约 25 倍。