「移動する」球内の複数の配列間でデータを検索するコードのランタイムを高速化する方法

時間ステップのサブセットにわたってCFDデータ（スカラーN x M x P配列の形式です。NはY、Mからx、およびPからzに対応します）を平均化しようとしています。以下の希望する平均化プロセスの説明を簡略化しようとしました。

1. 各タイムステップでグリッドを指定された角度で回転させます（これは、フローが各タイムステップで回転して形状/サイズを変更するコヒーレント構造を持っているためです。それらを重ね合わせて、次のような構造の時間平均形式を見つけたいためです。時間の経過に伴う形状/サイズの変化を説明します）
2. 元の回転していないグリッドを中心とする球を描画する
3. 球内にあるすべての回転グリッドからグリッドポイントを特定する
4. 回転した各グリッドのグリッドポイントのインデックスを特定します
5. インデックスを使用して、球内の回転したグリッドポイントでスカラーデータを見つけます
6. 球内の値の平均を取る
7. その新しい平均値を回転していないグリッド上の場所に配置します

自分がやりたいことを正しく実行しているように見えるコードがありますが、計算を完了するのに時間がかかりすぎます。実行速度を上げたいので、必要に応じてコードを変更できます。以下は、小さいバージョンのデータで機能する私のコードのバージョンです。

x = -5:5;       % x - position data
y = -2:.5:5;    % y - position data
z = -5:5;       % z - position data
% 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

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);

        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
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);

各インデックスが、回転していないグリッド内の特定の位置にあるすべてのタイムステップをカバーするすべてのポイントの加重平均であるN x M xP配列を取得することを期待しています。ご覧のとおり、これはまさに私が得たものです。ただし、主な問題は、「実際の」データのより大きなセットを使用する場合に、そうするのにかかる時間の長さです。上記のコードの実行には数分しかかかりませんが、M = 120、N = 24、およびP = 120で、タイムステップ数が24の場合、これにはさらに時間がかかる可能性があります。私の見積もりによると、計算全体が完了するまでに約25日以上かかります。