% 立方体范围内找点

X = [8028.14600000000        3378.20500000000        36.5120000000000

8028.22800000000        3378.18000000000        36.5260000000000

8028.80400000000        3377.69100000000        36.5000000000000

8028.10600000000        3378.32300000000        36.5150000000000

8028.81300000000        3377.81500000000        36.4950000000000

8028.95100000000        3377.74100000000        36.5120000000000

8029.00100000000        3377.74200000000        36.5140000000000

8028.99800000000        3377.79100000000        36.4990000000000

8028.93900000000        3377.88200000000        36.4940000000000

8029.02600000000        3377.85100000000        36.4920000000000

8029.27900000000        3378.02700000000        36.5010000000000

8029.38400000000        3378.06900000000        36.4930000000000

8029.41500000000        3378.13300000000        36.5040000000000

8029.45200000000        3378.14400000000        36.5040000000000

8029.17900000000        3378.36500000000        36.4820000000000

8029.46800000000        3378.17400000000        36.5220000000000

8029.30700000000        3378.30500000000        36.4950000000000

8029.45000000000        3378.23100000000        36.5110000000000

8029.30900000000        3378.34500000000        36.4840000000000

8029.48300000000        3378.24600000000        36.5160000000000

8029.49600000000        3378.27700000000        36.5060000000000

8029.58200000000        3378.29100000000        36.4970000000000

8026.79200000000        3379.27600000000        36.5170000000000

8026.19900000000        3379.86100000000        36.5100000000000

8026.33400000000        3379.78700000000        36.5200000000000

8026.12900000000        3379.96100000000        36.5000000000000

8026.27400000000        3379.87500000000        36.5080000000000];

X = repmat(X,100000,1); % 矩阵扩充100000倍

targetPt = X(3,:); % 选取第三个点为目标点

r = 0.5; % 半径为0.5

% 立方体

tic;

xyzMin = targetPt - r;

xyzMax = targetPt + r;

xInBox = (X(:,1) >= xyzMin(1) & X(:,1)<= xyzMax(1) );

yInBox = (X(:,2) >= xyzMin(2) & X(:,2)<= xyzMax(2) );

zInBox = (X(:,3) >= xyzMin(3) & X(:,3)<= xyzMax(3) );

ptsInBox = xInBox & yInBox & zInBox;

allPtsInBox = X(ptsInBox,:);

toc;

% 立方体(循环方法，优点：慢到没朋友！)

% tic;

% nRows = size(X,1);

% xyzMin = targetPt - r;

% xyzMax = targetPt + r;

%

% allPtsInBox = [];

% % ptsInBox = zeros(nRows,1);

% for i = 1:nRows

%     if (X(i,1) >= xyzMin(1) && X(i,1)<= xyzMax(1))

%         if (X(i,2) >= xyzMin(2) && X(i,2)<= xyzMax(2) )

%             if (X(i,3) >= xyzMin(3) && X(i,3)<= xyzMax(3) )

%                 allPtsInBox = [allPtsInBox; X(i,:)];

% %                 ptsInBox(i) = 1;

%             end

%         end

%     end

% end

% toc;

% 球体

tic;

tarPtMat = repmat(targetPt,size(X,1),1);

distP2PMat = sqrt(sum((X - tarPtMat).^2,2));

ptsInSphere = find(distP2PMat <= r);

allPtsInSphere = X(ptsInSphere,:);

toc;