Есть много вещей, которые вы могли бы сделать, чтобы улучшить свою программу - как алгоритм, так и код.
На стороне кода, одна из вещей, которые действительно замедлят является тот факт, что не только вы используете for цикл (который медленно), но в линии
P = [P;P1];
вы добавляете элементы к массиву. Каждый раз, когда это происходит, Matlab нужно найти новое место для размещения данных, копируя все точки в процессе. Это быстро становится очень медленным. Предварительное выделение массива с
P = zeros(1000000, 3);
отслеживающих число N точек вы нашли до сих пор, и изменение расчета расстояния до
D = pdist2(P1, P(1:N, :), 'euclidean');
бы по крайней мере, адрес ...
Другая проблема заключается в том, что вы проверяете новые очки против всех найденных ранее очков - поэтому, когда у вас есть 100 очков, вы проверяете около 100x100, за 1000 - 1000x1000. Вы можете видеть, что этот алгоритм O (N^3), по крайней мере ... не считая того, что вы получите больше «промахов» по мере роста плотности. Алгоритм O (N^3) с N = 10E6 занимает не менее 10E18 циклов; если бы у вас была машина 4 ГГц с одним тактовым циклом для сравнения, вам понадобилось бы 2.5E9 секунд = примерно 80 лет. Вы можете попробовать параллельную обработку, но это просто грубая сила - кто этого хочет?
Я рекомендую вам подумать о том, чтобы разбить проблему на более мелкие кусочки (довольно буквально): например, если вы разделите сферу на маленькие квадратики размером примерно с вашего максимального расстояния, и для каждой коробки вы отслеживаете какие точки в нем, тогда вам нужно только проверить очки в «этом» поле и его ближайших соседях - всего 27 ящиков. Если ваши коробки имеют ширину 2,5 мм, у вас будет 100x100x100 = 1M боксов. Это похоже на много, но теперь ваше время вычисления будет резко сокращено, так как у вас будет (по окончанию алгоритма) только 1 балл в среднем за поле ... Конечно, с критерием расстояния, который вы используете, вы будете есть больше точек рядом с центром, но это детали.
Структура данных, которая вам нужна, будет массивом ячеек 100x100x100, и каждая ячейка содержит индекс хороших точек, найденных до сих пор «в этой ячейке».Проблема с массивом ячеек заключается в том, что он не поддается векторизации. Если вместо этого у вас есть память, вы можете назначить ее как 4D-массив размером 10x100x100x100, если у вас будет не более 10 очков на ячейку (если вы это сделаете, вам придется обрабатывать это отдельно, работайте со мной здесь ...) , Используйте индекс -1 точек пока не найдено
Тогда ваш чек будет что-то вроде этого:
% initializing:
bigList = zeros(10,102,102,102)-1; % avoid hitting the edge...
NPlist = zeros(102, 102, 102); % track # valid points in each box
bottomcorner = [-25.5, -25.5, -25.5]; % boxes span from -25.5 to +25.5
cellSize = 0.5;
.
% in your loop:
P1= [x, y, z];
cellCoords = ceil(P1/cellSize);
goodFlag = true;
pointsSoFar = bigList(:, cellCoords(1)+(-1:1), cellCoords(2)+(-1:1), cellCoords(3)+(-1:1));
pointsToCheck = find(pointsSoFar>0); % this is where the big gains come...
r=sum(P1.^2);
D = pdist2(P1,P(pointsToCheck, :),'euclidean'); % euclidean distance
if D>0.146*r^(2/3)
P(k,:) = P1;
% check the syntax of this line...
cci = ind2sub([102 102 102], cellCoords(1), cellCoords(2), cellCoords(3));
NP(cci)=NP(cci)+1; % increasing number of points in this box
% you want to handle the case where this > 10!!!
bigList(NP(cci), cci) = k;
k=k+1;
end
....
Я не знаю, можете ли вы взять его отсюда; если вы не можете сказать это в заметках, и у меня может быть какое-то время в эти выходные, чтобы более подробно описать это. Есть способы ускорить его с помощью некоторой векторизации, но с этим быстро становится сложно справиться.
Я думаю, что случайное размещение в пространстве большего количества точек, а затем использование выше для гигантского векторизованного отбраковки, может быть, путь. Но сначала я рекомендую сделать небольшие шаги ... если вы можете заставить все работать хорошо, вы можете затем оптимизировать его (размер массива и т. Д.).
Я использую эту программу: –
Nmax = 10000; R = 25; P = rand (1,3); k = 1; , в то время как k 0,146 * r^(2/3) P = [P; P1]; k = k + 1; конец i = i + 1; конец x = P (:, 1); y = P (:, 2); z = P (:, 3); plot3 (x, y, z, '.'); –