矩阵"Zigzag"重新排序
fbr*_*eto 35 matlab loops matrix
我在MATLAB中有一个NxM矩阵,我想以类似于JPEG重新排序其子块像素的方式重新排序:
我希望算法是通用的,这样我就可以传入任何维度的二维矩阵.我是一名C++程序员,我非常想写一个旧的学校循环来实现这个目标,但我怀疑在MATLAB中有更好的方法.
我宁愿想要一个在NxN矩阵上运行的算法,然后从那里开始.
例:
1 2 3
4 5 6 --> 1 2 4 7 5 3 6 8 9
7 8 9
Amr*_*mro 25
考虑一下代码:
M = randi(100, [3 4]); %# input matrix
ind = reshape(1:numel(M), size(M)); %# indices of elements
ind = fliplr( spdiags( fliplr(ind) ) ); %# get the anti-diagonals
ind(:,1:2:end) = flipud( ind(:,1:2:end) ); %# reverse order of odd columns
ind(ind==0) = []; %# keep non-zero indices
M(ind) %# get elements in zigzag order
4x4矩阵的示例:
» M
M =
17 35 26 96
12 59 51 55
50 23 70 14
96 76 90 15
» M(ind)
ans =
17 35 12 50 59 26 96 51 23 96 76 70 55 14 90 15
和一个非方矩阵的例子:
M =
69 9 16 100
75 23 83 8
46 92 54 45
ans =
69 9 75 46 23 16 100 83 92 54 8 45
这种方法非常快:
X = randn(500,2000); %// example input matrix
[r, c] = size(X);
M = bsxfun(@plus, (1:r).', 0:c-1);
M = M + bsxfun(@times, (1:r).'/(r+c), (-1).^M);
[~, ind] = sort(M(:));
y = X(ind).'; %'// output row vector
标杆
下面的代码进行比较的运行时间与的荷银的出色答卷,使用timeit.它测试矩阵大小(条目数)和矩阵形状(行数与列数比率)的不同组合.
%// Amro's approach
function y = zigzag_Amro(M)
ind = reshape(1:numel(M), size(M));
ind = fliplr( spdiags( fliplr(ind) ) );
ind(:,1:2:end) = flipud( ind(:,1:2:end) );
ind(ind==0) = [];
y = M(ind);
%// Luis' approach
function y = zigzag_Luis(X)
[r, c] = size(X);
M = bsxfun(@plus, (1:r).', 0:c-1);
M = M + bsxfun(@times, (1:r).'/(r+c), (-1).^M);
[~, ind] = sort(M(:));
y = X(ind).';
%// Benchmarking code:
S = [10 30 100 300 1000 3000]; %// reference to generate matrix size
f = [1 1]; %// number of cols is S*f(1); number of rows is S*f(2)
%// f = [0.5 2]; %// plotted with '--'
%// f = [2 0.5]; %// plotted with ':'
t_Amro = NaN(size(S));
t_Luis = NaN(size(S));
for n = 1:numel(S)
X = rand(f(1)*S(n), f(2)*S(n));
f_Amro = @() zigzag_Amro(X);
f_Luis = @() zigzag_Luis(X);
t_Amro(n) = timeit(f_Amro);
t_Luis(n) = timeit(f_Luis);
end
loglog(S.^2*prod(f), t_Amro, '.b-');
hold on
loglog(S.^2*prod(f), t_Luis, '.r-');
xlabel('number of matrix entries')
ylabel('time')
下图是使用Matlab R2014b在Windows 7 64位上获得的.R2010b的结果非常相似.可以看出,新方法将运行时间缩短了2.5(对于小矩阵)和1.4(对于大矩阵).考虑到总条目数,结果被认为对矩阵形状几乎不敏感.
这是一个非循环解决方案zig_zag.m.它看起来很丑,但它有效!:
function [M,index] = zig_zag(M)
[r,c] = size(M);
checker = rem(hankel(1:r,r-1+(1:c)),2);
[rEven,cEven] = find(checker);
[cOdd,rOdd] = find(~checker.'); %'#
rTotal = [rEven; rOdd];
cTotal = [cEven; cOdd];
[junk,sortIndex] = sort(rTotal+cTotal);
rSort = rTotal(sortIndex);
cSort = cTotal(sortIndex);
index = sub2ind([r c],rSort,cSort);
M = M(index);
end
和一个测试矩阵:
>> M = [magic(4) zeros(4,1)];
M =
16 2 3 13 0
5 11 10 8 0
9 7 6 12 0
4 14 15 1 0
>> newM = zig_zag(M) %# Zig-zag sampled elements
newM =
16
2
5
9
11
3
13
10
7
4
14
6
8
0
0
12
15
1
0
0
这是一种如何做到这一点的方法.基本上,你的数组是一个hankel矩阵加上1:m的向量,其中m是每个对角线中元素的数量.也许其他人对于如何创建必须在没有循环的情况下添加到翻转的hankel阵列的对角线阵列有一个巧妙的想法.
我认为这应该可以推广到非正方形数组.
% for a 3x3 array
n=3;
numElementsPerDiagonal = [1:n,n-1:-1:1];
hadaRC = cumsum([0,numElementsPerDiagonal(1:end-1)]);
array2add = fliplr(hankel(hadaRC(1:n),hadaRC(end-n+1:n)));
% loop through the hankel array and add numbers counting either up or down
% if they are even or odd
for d = 1:(2*n-1)
if floor(d/2)==d/2
% even, count down
array2add = array2add + diag(1:numElementsPerDiagonal(d),d-n);
else
% odd, count up
array2add = array2add + diag(numElementsPerDiagonal(d):-1:1,d-n);
end
end
% now flip to get the result
indexMatrix = fliplr(array2add)
result =
1 2 6
3 5 7
4 8 9
之后,您只需调用reshape(image(indexMatrix),[],1)以获取重新排序元素的向量.
编辑
好的,从您的评论中看起来您需要sort像Marc建议的那样使用.
indexMatrixT = indexMatrix'; % ' SO formatting
[dummy,sortedIdx] = sort(indexMatrixT(:));
sortedIdx =
1 2 4 7 5 3 6 8 9
请注意,在索引之前,您需要首先转置输入矩阵,因为Matlab首先计数,然后计算正确.