问如何加快在“移动”范围内的多个数组之间搜索数据的代码运行时

EN

% x = -5:2:5;       % x - position data
x = linspace(-5,5,120);
% y = -2:5;    % y - position data
y = linspace(-2,5,24);
% z = -5:2:5;       % z - position data
z = linspace(-5,5,120);
% my grid is much bigger actually

[X,Y,Z] = meshgrid(x,y,z);    % mesh for plotting data

dX = diff(x)';   dX(end+1) = dX(end);    % x grid intervals
dY = diff(y)';   dY(end+1) = dY(end);    % y grid intervals
dZ = diff(z)';   dZ(end+1) = dZ(end);    % z grid intervals

TestPoints = combvec(x,y,z)';    % you need the Matlab Neural Network Toolbox to run this
dXYZ = combvec(dX',dY',dZ')';

% TestPoints is the unrotated grid

M = length(x);   % size of grid x - direction
N = length(y);   % size of grid y - direction
P = length(z);   % size of grid z - direction

D = randi([-10,10],N,M,P,3);    % placeholder for data for 3 time steps (I have more than 3, this is a subset)
D2{3,M*N*P} = [];
PosAll{3,2} = [];

[xSph,ySph,zSph] = sphere(50);

c = 0.01;   % 1 cm
nu = 8e-6;  % 8 cSt
s = 3*c;    % span for Aspect Ratio 3
r_g = s/sqrt(3);
U_g = 110*nu/c; % velocity for Reynolds number 110
Omega = U_g/r_g;    % angular velocity
T = (2*pi)/Omega;   % period
dt = 40*T/1920;     % time interval
DeltaRotAngle = ((2*pi)/T)*dt;  % angle interval

timesteps = 121:123;   % time steps 121, 122, and 123
for ti=timesteps
    tj = find(ti==timesteps);
    Theta = ti*DeltaRotAngle;
    Rotate = [cos(Theta),0,sin(Theta);...
        0,1,0;...
        -sin(Theta),0,cos(Theta)];
    PosAll{tj,1} = (Rotate*TestPoints')';
end
tic
for i=1:M*N*P

    aa = TestPoints(i,1);
    bb = TestPoints(i,2);
    cc = TestPoints(i,3);
    rs = 0.8*sqrt(dXYZ(i,1)^2 + dXYZ(i,2)^2 + dXYZ(i,3)^2);

%     handles.H = figure;
%     hs = surf(xSph*rs+aa,ySph*rs+bb,zSph*rs+cc);
%     [Fs,Vs,~] = surf2patch(hs,'triangle');
%     close(handles.H)

    for ti=timesteps
        tj = find(timesteps==ti);

%         f = inpolyhedron(Fs,Vs,PosAll{tj,1},'FlipNormals',false);
        f = sqrt(sum((PosAll{tj,1}-[aa,bb,cc]).^2,2))<rs;
        TestPointsR_ti = PosAll{tj,1};
        PointsInSphere = TestPointsR_ti(f,:);

        p1 = [aa,bb,cc];
        p2 = [PointsInSphere(:,1),...
            PointsInSphere(:,2),...
            PointsInSphere(:,3)];

        w = 1./sqrt(sum(...
            (p2-repmat(p1,size(PointsInSphere,1),1))...
            .^2,2));

        D_ti = reshape(D(:,:,:,tj),M*N*P,1);
        D2{tj,i} = [D_ti(f),w];
    end
    if ~mod(i,10)
        toc
    end
end

D3{M*N*P,1} = [];
for i=1:M*N*P
    D3{i} = vertcat(D2{:,i});
end

D4 = zeros(M*N*P,1);
for i=1:M*N*P
    D4(i) = sum(D3{i}(:,1).*D3{i}(:,2))/...
                sum(D3{i}(:,2));
end
D_ta = reshape(D4,N,M,P);