Векторизовать циклы for, которые вызывают другие функции

У меня есть следующий фрагмент кода в Matlab с двумя циклами for: «I» — это двоичное изображение, которое было предварительно выделено.

введите здесь описание изображения

 ...
    [x,y] = find(bwmorph(I,'endpoints'));
    n=numel(x);
    m=numel(x)-1;
    n=m+1;
    r=i+1;
    for i= 1:m
        for j = r:n
            I=linept(I, x(i), y(i), x(j), y(j));
        end;
    end;
    ...

Функция linept приведена ниже. Она взята из файлового обмена Matlab:

function result=linept(matrix, X1, Y1, X2, Y2)
result = matrix;
a=max(1, X1);b=sign(X2 - X1);c=max(1, X2);
for x=a:b:c
    y = round(f(x, X1, Y1, X2, Y2));
    if y > 0
        result(x, y) = 1;
    end
end
d=max(1, Y1);e=sign(Y2 - Y1);g=max(1, Y2);
for y=d:e:g
    x = round(f2(y, X1, Y1, X2, Y2));
    if x > 0
        result(x, y) = 1;
    end
end

function y=f(x, X1, Y1, X2, Y2)
a = (Y2 - Y1)/(X2 - X1);
b = Y1 - X1 * a;
y = a * x + b;

function x=f2(y, X1, Y1, X2, Y2)
if X1==X2
    x = X1;
else
    a = (Y2 - Y1)/(X2 - X1);
    b = Y1 - X1 * a;
    x = (y - b)/a;
end

Из-за большого количества циклов for и вызовов функций этот код выполняется очень медленно. Он выполняется быстро для простого изображения с несколькими конечными точками, но занимает много времени, когда количество ребер больше. Это немного быстрее. если размер изображения уменьшен. Я попытался векторизовать его и предварительно распределил некоторые переменные, но особых улучшений нет. Может ли кто-нибудь помочь мне в том, как векторизовать коды, которые вызывают функции в циклах. Спасибо


person Matte    schedule 05.02.2015    source источник
comment
Не могли бы вы также поделиться кодами функций func2 и func3? Кроме того, что такое m в n=m+1; и i в r=i+1; в начале? Точно так же эти a,b,c.. не могли бы вы добавить, что они означают?   -  person Divakar    schedule 05.02.2015
comment
Да, сейчас я отредактирую вопрос   -  person Matte    schedule 05.02.2015
comment
Прошу прощения, редактирование занимает время. На самом деле это связано с моими предыдущими вопросами. Я скоро закончу редактирование.   -  person Matte    schedule 05.02.2015
comment
Вы можете попробовать использовать циклы parfor для разных изображений, чтобы улучшить время работы. Он будет распределять работу по разным ядрам на вашем компьютере.   -  person C.Colden    schedule 05.02.2015
comment
Благодарю вас! я посмотрю   -  person Matte    schedule 05.02.2015


Ответы (1)


Это была одна большая проблема!!

Краткое обсуждение и коды решений

Далее в списке указан векторизованный подход, в котором активно используется bsxfun в разных местах, чтобы позаботиться о соединении всех точек со всеми другими точками с помощью expansions, что в основном и является тем, как работает bsxfun. Код в нем использует пример входных данных для демонстрационных целей. Взглянем -

%// Create demo input data
img = false(20);
img(2,5:15) = 1;
img(12,5:15) = 1;
figure,imagesc(img), colormap gray, title('Starting Binary image')

%// Find endpoints
[X,Y] = find(bwmorph(img,'endpoints'));

%// Make a new binary image with only endpoints in it
I = false(size(img));
I(sub2ind(size(I),X,Y)) = 1;
figure,imagesc(I), colormap gray, title('Endpoints detected')

%-------- Vectorized code starts here ...
[nrows,ncols] = size(I);  %// Parameters
npts = numel(X);

twopts = nchoosek(1:npts,2);
slopes = (Y(twopts(:,1)) - Y(twopts(:,2)))./(X(twopts(:,1)) - X(twopts(:,2)));

%// Find the connecting indices with X-Y as they are, to work with 
%// slopes within [-1 1]
stage1 = abs(slopes)<=1;
[colID,rowID] = connecting_idx(X,Y,twopts(stage1,:),slopes(stage1));
valid = colID>0 & rowID>0 & colID<=ncols & rowID<=nrows;
I((colID(valid)-1)*nrows + rowID(valid))=1;

%// Find the connecting indices with X-Y interchanged, to work with 
%// slopes outside [-1 1]
[rowID,colID] = connecting_idx(Y,X,twopts(~stage1,:),1./slopes(~stage1));
valid = colID>0 & rowID>0 & colID<=ncols & rowID<=nrows;
I((colID(valid)-1)*nrows + rowID(valid))=1;
figure,imagesc(I),colormap gray,title('Final output : Endpoints connected')

Связанный код функции (наиболее важная часть кодовой базы) -

function [y_all,x_all] = connecting_idx(X,Y,twopts,slopes)

%// Find XY indices that connects the X, Y anchor points given the two points
%// combinations and the corresponding slopes

X_twopts = reshape(X(twopts),size(twopts));
Y_twopts = reshape(Y(twopts),size(twopts));

[sortedX_pairs1,sorted_idx1] = sort(X_twopts,2);
X_starts1 = sortedX_pairs1(:,1);
Y_pairs = Y_twopts;
Y_starts = Y_pairs(sorted_idx1==1);
offsets = Y_starts - slopes.*X_starts1;
max_X_len = max(diff(sortedX_pairs1,[],2));

x_all = bsxfun(@plus,X_starts1,[0:max_X_len]);
x_lens = diff(sortedX_pairs1,[],2);
mask = bsxfun(@gt,x_lens+1,0:max_X_len);
y_all = round(bsxfun(@plus,bsxfun(@times,x_all,slopes),offsets));

y_all = y_all(mask);
x_all = x_all(mask);

return;

Отладка изображений после запуска кода —

введите здесь описание изображения

введите здесь описание изображения

введите здесь описание изображения

person Divakar    schedule 05.02.2015
comment
Не могли бы вы подробнее объяснить, что вы сделали в этой части:stage1 = abs(slopes)<=1; [colID,rowID] = connecting_idx(X,Y,twopts(stage1,:),slopes(stage1)); valid = colID>0 & rowID>0 & colID<=ncols & rowID<=nrows; I((colID(valid)-1)*nrows + rowID(valid))=1; - person Matte; 05.02.2015
comment
@Matt Well slopes - это массив, в котором хранятся наклоны для всех точек попарно. С помощью stage1 я выбираю пары, которые имеют наклоны в пределах [-1 1], а затем нахожу связь между этими парами. На следующем шаге, если вы видите, я использую ~stage1, чтобы выбрать остальные пары и найти их соединительные индексы. Таким образом, после получения связующих индексов с обоих шагов мы устанавливаем их во входной матрице как единицы. - person Divakar; 05.02.2015
comment
@Matt Да, я только что заметил эту проблему и с некоторыми другими изображениями, думаю, мне нужно еще раз отладить то, что происходит внутри кодов! - person Divakar; 05.02.2015
comment
Я пройдусь по коду еще раз и постараюсь понять его в совершенстве! - person Matte; 05.02.2015
comment
Не могли бы вы указать мне ошибку в коде, если вы его отладили? Я пытался найти его, но мне потребовалось много времени, чтобы понять векторизованный код, и я не мог найти ошибку. - person Matte; 06.02.2015