Эффективный расчет смежности 3D-доменов в Matlab? - PullRequest
Новая
       23

Эффективный расчет смежности 3D-доменов в Matlab?

0 голосов
/ 01 мая 2019

У меня есть скрипт Matlab, который генерирует аддитивные домены Вороного в топологии трехмерного тора. Поскольку в Matlab нет нативной взвешенной диаграммы Вороного, я просто вычисляю ее как матрицу, в которой воксели принадлежат разным ячейкам, индексированным целыми числами. A slice through the partition Тем не менее, я хотел бы получить данные о том, какие ячейки находятся рядом друг с другом, чтобы я мог проанализировать детали топологии.

Наиболее очевидный способ сделать это - выполнить итерацию проверки вокселей, принадлежат ли 6 или 26 соседних вокселей к другому домену. Это будет выполняться за O (N ^ 3) времени и не будет выглядеть слишком параллелизуемым.

Можно также использовать круговое смещение, чтобы найти места, где есть граница (например, d=abs(im-circshift(im,[1,0]))>0; [i,j]=find(d);), перебрать граничные точки и заполнить массив отношений смежности (один проход для каждого направления в массиве вокселей). Поскольку find не очень хорошо обрабатывает трехмерные массивы, мне нужен специальный случай для направления z (пример кода см. Ниже). Это лучше и не слишком медленно, но грязно и не распараллеливается.

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

Есть идеи, как эффективно рассчитать смежность? 

% Input 3D NxNxN array gg of domain membership
w=zeros(Nciv+1); % Adjacency matrix; index 1 corresponds to empty domains gg(i,j,k)=0
for z=1:N
    im=gg(:,:,z); % local slize
    imagesc(im);
    drawnow

    m1=abs(im-circshift(im,[1,0]))>0; % Borders in x-direction
    [i,j]=find(m1);
    for k=1:length(i)
        if (i(k)<N) % if not spanning edge of array
            w(1+im(i(k),j(k)),1+im(i(k)+1,j(k)))=1;
        else
            w(1+im(i(k),j(k)),1+im(1,j(k)))=1;
        end
    end

    m2=abs(im-circshift(im,[0,1]))>0; % Borders in y-direction
    [i,j]=find(m2);
    for k=1:length(i)
        if (j(k)<N) % if not spanning edge of array
            w(1+im(i(k),j(k)),1+im(i(k),j(k)+1))=1;
        else
            w(1+im(i(k),j(k)),1+im(i(k),1))=1;
        end
    end

    % z-direction case
    if (z<N)
        im2=gg(:,:,z+1);
    else
        im2=gg(:,:,1);
    end
    [i,j]=find(abs(im-im2)>0);
    for k=1:length(i)
        w(1+im(i(k),j(k)),1+im2(i(k),j(k)))=1;
    end
end
% Make adjacency matrix nice and display it 
w=w+w';
w=w-diag(diag(w));
w=(w>0);
spy(w)
...