Векторизованные суммы разных диагоналей в матрице - PullRequest
Новая
       21

Векторизованные суммы разных диагоналей в матрице

9 голосов
/ 19 мая 2010

Я хочу векторизовать следующий код MATLAB. Я думаю, это должно быть просто, но, тем не менее, меня это смущает. 

r = some constant less than m or n
[m,n] = size(C);
S = zeros(m-r,n-r);
for i=1:m-r+1
    for j=1:n-r+1
        S(i,j) = sum(diag(C(i:i+r-1,j:j+r-1)));
    end
end

Код вычисляет таблицу оценок, S , для алгоритма динамического программирования, из другой таблицы оценок, C .
Диагональное суммирование - это генерация баллов для отдельных фрагментов данных, использованных для генерации C , для всех возможных фрагментов (размера r).

Заранее спасибо за любые ответы! Извините, если это должно быть очевидно ...

Примечание
Встроенный conv2 оказался быстрее, чем convnfft, потому что мой глаз (r) довольно маленький (5 <= r <= 20). convnfft.m заявляет, что r должно быть> 20, чтобы проявить любую выгоду.

Ответы [ 4 ]

10 голосов
/ 19 мая 2010

Если я правильно понимаю, вы пытаетесь вычислить диагональную сумму каждого подмассива C, где вы удалили последнюю строку и столбец C (если вы не должны удалять строку / столбец, вам нужно выполнить цикл до m-r + 1, и вам нужно передать весь массив C функции в моем решении ниже).

Вы можете выполнить эту операцию с помощью свертки , например: 

S = conv2(C(1:end-1,1:end-1),eye(r),'valid');

Если C и r большие, вы можете взглянуть на CONVNFFT из Matlab File Exchange, чтобы ускорить вычисления.

3 голосов
/ 22 мая 2010

Исходя из идеи JS и, как указано в комментариях Джонас , это можно сделать в две строки, используя IM2COL с некоторыми манипуляциями с массивами: 

B = im2col(C, [r r], 'sliding');
S = reshape( sum(B(1:r+1:end,:)), size(C)-r+1 );

В основном B содержит элементы всех скользящих блоков размера r-на-r по матрице C. Затем мы берем элементы по диагонали каждого из этих блоков B(1:r+1:end,:), вычисляем их сумму и изменяем результат до ожидаемого размера.

Сравнивая это с решением на основе свертки Йонаса, он не выполняет никакого умножения матриц, а только индексирует ...

2 голосов
/ 19 мая 2010

Я думаю, вам может понадобиться переставить C в трехмерную матрицу, прежде чем суммировать ее по одному из измерений. Я отправлю с ответом в ближайшее время.

EDIT

Мне не удалось найти способ аккуратно векторизовать его, но я нашел функцию accumarray, которая могла бы помочь. Я посмотрю на это более подробно, когда буду дома.

EDIT # 2

Нашли более простое решение с помощью линейного индексирования, но это может потребовать много памяти.

При C (1,1) индексы, которые мы хотим суммировать, составляют 1+ [0, m + 1, 2 * m + 2, 3 * m + 3, 4 * m + 4, ...], или (0: r-1) + (0: m: (r-1) * m) 

sum_ind = (0:r-1)+(0:m:(r-1)*m);

создать S_offset, (mr) по (nr) по r матрице, так что S_offset (:,:, 1) = 0, S_offset (:,:, 2) = m + 1, S_offset (:, :, 3) = 2 * m + 2 и т. Д. 

S_offset = permute(repmat( sum_ind, [m-r, 1, n-r] ), [1, 3, 2]);

create S_base, матрица адресов базового массива, из которой будет рассчитываться смещение. 

S_base = reshape(1:m*n,[m n]);
S_base = repmat(S_base(1:m-r,1:n-r), [1, 1, r]);

Наконец, используйте S_base+S_offset для адресации значений C. 

S = sum(C(S_base+S_offset), 3);

Конечно, вы можете использовать bsxfun и другие методы, чтобы сделать его более эффективным; здесь я решил выложить это для ясности. Мне еще предстоит сравнить это с тем, чтобы сравнить его с методом двойной петли; Мне нужно сначала пойти домой на ужин!

1 голос
/ 17 октября 2013

Это то, что вы ищете? Эта функция добавляет диагонали и помещает их в вектор, аналогично тому, как функция sum суммирует все столбцы в матрице и помещает их в вектор. 

function [diagSum] = diagSumCalc(squareMatrix, LLUR0_ULLR1)
% 
% Input: squareMatrix: A square matrix.
%        LLUR0_ULLR1:  LowerLeft to UpperRight addition = 0     
%                      UpperLeft to LowerRight addition = 1
% 
% Output: diagSum: A vector of the sum of the diagnols of the matrix.
% 
% Example: 
% 
% >> squareMatrix = [1 2 3; 
%                    4 5 6;
%                    7 8 9];
% 
% >> diagSum = diagSumCalc(squareMatrix, 0);
% 
% diagSum = 
% 
%       1 6 15 14 9
% 
% >> diagSum = diagSumCalc(squareMatrix, 1);
% 
% diagSum = 
% 
%       7 12 15 8 3
% 
% Written by M. Phillips
% Oct. 16th, 2013
% MIT Open Source Copywrite
% Contact mphillips@hmc.edu fmi.
% 

if (nargin < 2)
    disp('Error on input. Needs two inputs.');
    return;
end

if (LLUR0_ULLR1 ~= 0 && LLUR0_ULLR1~= 1)
    disp('Error on input. Only accepts 0 or 1 as input for second condition.');
    return;
end

[M, N] = size(squareMatrix);

if (M ~= N)
    disp('Error on input. Only accepts a square matrix as input.');
    return;
end

diagSum = zeros(1, M+N-1);

if LLUR0_ULLR1 == 1
    squareMatrix = rot90(squareMatrix, -1);
end

for i = 1:length(diagSum)
    if i <= M
        countUp = 1;
        countDown = i;
        while countDown ~= 0
            diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i);
            countUp = countUp+1;
            countDown = countDown-1;
        end
    end
    if i > M
        countUp = i-M+1;
        countDown = M;
        while countUp ~= M+1
            diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i);
            countUp = countUp+1;
            countDown = countDown-1;
        end
    end
end

Приветствия

...