Преобразовать облако точек в воксели с помощью усреднения

Question

У меня есть следующие данные:

N = 10^3;
x = randn(N,1);
y = randn(N,1);
z = randn(N,1);
f = x.^2+y.^2+z.^2;

Теперь я хочу разбить это непрерывное трехмерное пространство на nB корзины.

nB = 20;
[~,~,x_bins] = histcounts(x,nB);
[~,~,y_bins] = histcounts(y,nB);
[~,~,z_bins] = histcounts(z,nB);

И положить в каждый куб среднее значение f или nan если в кубе не наблюдалось никаких наблюдений:

F = nan(50,50,50);

for iX = 1:20
    for iY = 1:20
        for iZ = 1:20
            idx = (x_bins==iX)&(y_bins==iY)&(z_bins==iZ);
            F(iX,iY,iZ) = mean(f(idx));
        end
    end
end
isosurface(F,0.5)

Этот код выполняет то, что я хочу.Моя проблема в скорости.Этот код очень медленный, когда N > 10^5 и nB = 100.

Как я могу ускорить этот код?

---

Я также попробовал функцию accumarray():

subs=([x_bins,y_bins,z_bins]);
F2 = accumarray(subs,f,[],@mean);
all(F(:) == F2(:)) % false

Однако этот код дает другой результат.

Answer 1

Проблема с кодом в OP заключается в том, что он проверяет все элементы данных для каждого элемента в выходном массиве. Выходной массив содержит nB^3 элементов, данные содержат N элементов, поэтому алгоритм имеет значение O (N*nB^3). Вместо этого можно зациклить элементы N на входе и установить соответствующий элемент в выходном массиве, который является операцией O (N) (2-й блок кода ниже).

Для решения accumarray в OP необходимо использовать параметр fillvals, установите для него значение NaN (третий кодовый блок ниже).

Чтобы сравнить результаты, нужно явно проверить, что оба массива имеют NaN в тех же местах и ​​имеют равные не-NaN значения в других местах:

all( ( isnan(F(:)) & isnan(F2(:)) ) | ( F(:) == F2(:) ) )
%    \-------same NaN values------/   \--same values--/

Вот код. Все три версии дают одинаковые результаты. Времена в Octave 4.4.1 (без JIT), в MATLAB код цикла должен быть быстрее. (Используя входные данные из OP, с N=10^3 и nB=20).

%% OP's code, O(N*nB^3)
tic
F = nan(nB,nB,nB);
for iX = 1:nB
    for iY = 1:nB
        for iZ = 1:nB
            idx = (x_bins==iX)&(y_bins==iY)&(z_bins==iZ);
            F(iX,iY,iZ) = mean(f(idx));
        end
    end
end
toc
% Elapsed time is 1.61736 seconds.

%% Looping over input, O(N)
tic
s = zeros(nB,nB,nB);
c = zeros(nB,nB,nB);
ind = sub2ind([nB,nB,nB],x_bins,y_bins,z_bins);
for ii=1:N
   s(ind(ii)) = s(ind(ii)) + f(ii);
   c(ind(ii)) = c(ind(ii)) + 1;
end
F2 = s ./ c;
toc
% Elapsed time is 0.0606539 seconds.

%% Other alternative, using accumarray
tic
ind = sub2ind([nB,nB,nB],x_bins,y_bins,z_bins);
F3 = accumarray(ind,f,[nB,nB,nB],@mean,NaN);
toc
% Elapsed time is 0.14113 seconds.