Векторизация минимального расстояния в ядре

У меня есть массив Nx2 K1 с расположением N ключевых точек и трехмерный массив WxHx3 Kart1(width,height,coordinates), который отображает координаты для каждого пикселя изображения. Для каждой ключевой точки в K1 я хочу прочитать местоположение пикселя в Kart1 и оценить координаты (поиск минимума/максимума или вычислить среднее значение) в ядре 3x3 вокруг него и присвоить значение текущему пикселю в KPCoor1.

Мой текущий подход выглядит следующим образом:

for ii=1:length(K1(:,1)) %for every keypoint in K1

    MinDist=sqrt(sum(Kart1(K1(ii,2)-1,K1(ii,1)-1,:).^2)); %Calculate distance
    xShift=0;
    yShift=0;
    for kk=-1:1 %for every pixel in a 3x3 kernel...
        for ll=-1:1

            Distance=sqrt(sum(Kart1(K1(ii,2)+kk,K1(ii,1)+ll,:).^2));

            if Distance<MinDist   %... if the current distance is smaller than MinDist
                MinDist=Distance; %... update MinDist...
                xShift=kk; %... and take the kernel coordinate of the pixel
                yShift=ll;
            end

        end
    end

    KP1Coor(ii,:)=Kart1(K1(ii,2)+xShift,K1(ii,1)+yShift,:); %assign the coordinates of the pixel with the minimal distance in kernel.

end

и он работает, но уродлив, и я сомневаюсь, что он делает то, что я хочу. Я немного смущен "многомерностью" вопроса, не знаю многих функций для оценки ядер и не могу придумать, как использовать функции векторизации типа bsxfun() или логические операции (значит, я застрял, и мой мозг сухо :/)

Любые предложения о том, как устранить эти циклы/исправить код?


matlab lookup multidimensional-array vectorization bsxfun
person McMa    schedule 12.03.2014    source источник
comment
Я думаю, вам нужны xShift=-1 и yShift=-1 внутри первого вложенного цикла, не так ли? Кроме того, вам нужно MinDist в качестве вывода?   -  person Divakar    schedule 12.08.2014
comment
В этом случае меня больше интересует, «где» минимальное расстояние, и я читаю это по сдвигам x и y. Вот для чего я их использую. PS: я решил проблему более элегантным способом, все равно спасибо за комментарий к старому вопросу!   -  person McMa    schedule 13.08.2014
comment
Ну, я все же продолжил свое собственное решение, основанное на bsxfun! Вы сказали, что уже реализовали его эффективную версию, и я подумал, что не буду возражать против немного showdown of benchmarks по сравнению с представленной здесь, если вы тоже хотите :) Если вы это сделаете, возможно, вы можете отредактировать свой вопрос с ним или сделайте свое собственное решение здесь, хотя я думаю, что последнее было бы неплохо.   -  person Divakar    schedule 13.08.2014
comment
Я имел в виду, что мне не нужно больше искать это, но будьте уверены, что ваш ответ выиграл бы вскрытие с большим отрывом! ;)   -  person McMa    schedule 18.08.2014

Ответы (1)


arrow_upward
2
arrow_downward

Векторный подход: после подробного изучения векторной реализации это выглядело очень интересно look-up problem, поэтому, если вы все еще заинтересованы в векторных методах выполнения работы, вот подход с bsxfun -

%// Scalars to replace their repeated usages in the code
[W,H,NC]= size(Kart1);
d3a = 1:NC;

%// Indices for Kart1 at the first level of the nested loops
BI = bsxfun(@plus,K1(:,2)+(K1(:,1)-1)*W,(d3a-1)*W*H);

%// Indices for Kart1 in 3x3 kernel around the indices obtained at first level
BIW3 = bsxfun(@plus,[-1 0 1]',[-W 0 W]); %//'
%// Linear indices for the minimum value of Kart1 in the 3x3 neighborhood 
[~,MI3] = min(sqrt(sum(Kart1(bsxfun(@plus,...
    BI,permute(BIW3(:),[3 2 1]))).^2,2)),[],3);
%// X-Y indices
[xShift1,yShift1] = ind2sub([3 3],MI3);

%// Get Kart1 values corresponding to the indices for the minimum values
KP1Coor = Kart1(bsxfun(@plus,...
    K1(:,2)+xShift1-2 +(K1(:,1)+yShift1-2-1)*W,(d3a-1)*W*H));

Сравнительный анализ

Я также смог протестировать это с графическим процессором, используя gpuArray из Parallel Computing Toolbox, а затем провел несколько тестов с W = 1000, H = 1000 и использовал N в качестве размера данных, варьируя его с [1000 2000 5000 10000 20000]. Результаты кажутся достаточно сумасшедшими, хотя и не ненадежными, поскольку использовались утвержденные методы сравнительного анализа из Measure and Improve GPU Performance. Вот эталонный график для исходных и векторизованных кодов ЦП и ГП -

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

Тогда казалось уместным сравнить только векторизованные коды и даже большие объемы данных, график которых показан далее -

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

Выводы. Проблема в основном выглядит как задача поиска, где Kart1 — это данные, а K1 — массив индексов для поиска. Представленное здесь векторизованное решение в основном основано на грубой силе, и результаты тестов явно говорят в пользу этого с точки зрения производительности. Но было бы интересно посмотреть, может ли какой-либо подход без грубой силы, который может быть даже основан на цикле, но эффективно использует этот поиск, мог бы работать лучше.

person Divakar    schedule 13.08.2014
comment
Спасибо за такой замечательный ответ и за возвращение к старому вопросу! Несмотря на то, что вне актуальной проблемы, я многому учусь и черпаю новые идеи из ваших работ! - person McMa; 18.08.2014