2017-03-06 108 views
0

我實際上向量化了我的一個代碼,並且遇到了一些問題。Matlab:向量化三維矩陣的工藝

這是我最初的代碼:

CoordVorBd = random(N+1,3) 
CoordCP = random(N,3) 
v = random(1,3) 
for i = 1 : N 
    for j = 1 : N 
      ri1j = (-CoordVorBd (i,:) + CoordCP(j,:)); 

      vij(i,j,:) = cross(v,ri1j))/(norm(ri1j) 
    end  

end 

我已經開始向量化,創建一個包含3個* 1向量一些矩陣。我的矩陣大小是N * N * 3。

CoordVorBd1(1:N,:) = CoordVorBd(2:N+1,:); 
CoordCP_x= CoordCP(:,1); 
CoordCP_y= CoordCP(:,2); 
CoordCP_z= CoordCP(:,3); 
CoordVorBd_x = CoordVorBd([1:N],1); 
CoordVorBd_y = CoordVorBd([1:N],2); 
CoordVorBd_z = CoordVorBd([1:N],3); 
CoordVorBd1_x = CoordVorBd1(:,1); 
CoordVorBd1_y = CoordVorBd1(:,2); 
CoordVorBd1_z = CoordVorBd1(:,3); 
[X,Y] = meshgrid (1:N); 

ri1j_x = (-CoordVorBd_x(X) + CoordCP_x(Y)); 
ri1j_y = (-CoordVorBd_y(X) + CoordCP_y(Y)); 
ri1j_z = (-CoordVorBd_z(X) + CoordCP_z(Y)); 

ri1jmat(:,:,1) = ri1j_x(:,:); 
ri1jmat(:,:,2) = ri1j_y(:,:); 
ri1jmat(:,:,3) = ri1j_z(:,:); 
vmat(:,:,1) = ones(N)*v(1); 
vmat(:,:,2) = ones(N)*v(2); 
vmat(:,:,3) = ones(N)*v(3); 

此代碼的工作原理,但在變量創建方面很重。我沒有實現將矢量化應用於所有矩陣。

ri1jmat(X,Y,1:3) = (-CoordVorBd (X,:) + CoordCP(Y,:)); 

的一級方程式不工作... 如果有人有一些想法,有一些清潔劑。

在這一點上,我有一個N * N * 3矩陣ri1jmat與我所有的向量。

欲計算N * N矩陣rij1norm即矢量

rij1norm(i,j) = norm(ri1jmat(i,j,1:3)) 

能夠矢量化維吉矩陣的範數。

vij(:,:,1:3) = (cross(vmat(:,:,1:3),ri1jmat(:,:,1:3))/(ri1jmatnorm(:,:)); 

該跨產品的作品。

我嘗試了一些沒有實現這個rij1norm矩陣而沒有做雙重循環的方法。

如果有人有一些技巧,提前致謝。

回答

0

這是矢量化版本。請注意,您的原始循環不包括CoordVorBd的最後一列,所以如果這是故意的,您還需要從下面的代碼中刪除它。我認爲這是一個錯誤。

CoordVorBd = rand(N+1,3); 
CoordCP = rand(N,3); 
v = rand(1,3); 

repCoordVor=kron(CoordVorBd', ones(1,size(CoordCP,1)))'; %based on http://stackoverflow.com/questions/16266804/matlab-repeat-every-column-sequentially-n-times 
repCoordCP=repmat(CoordCP, size(CoordVorBd,1),1); %repeat matrix 
V2=-repCoordVor + repCoordCP; %your ri1j 
nrm123=sqrt(sum(V2.^2,2)); %vectorized norm for each row 
vij_unformatted=cat(3,(v(:,2).*V2(:,3) - V2(:,2).*v(:,3))./nrm123,(v(:,3).*V2(:,1) - V2(:,3).*v(:,1))./nrm123,(v(:,1).*V2(:,2) - V2(:,1).*v(:,2))./nrm123); % cross product, expanded, and each term divided by norm, could use bsxfun(@rdivide,cr123,nrm123) instead, if cr123 is same without divisions 
vij=permute(reshape(vij_unformatted,N,N+1,3),[2,1,3]); %reformat to match your vij 
0

這裏是另一種方式使用arrayfun

% Define a meshgrid of indices to run over 
[I, J] = meshgrid(1:N, 1:(N+1)); 
% Calculate ril for each index 
rilj = arrayfun(@(x, y) -CoordVorBd (y,:) + CoordCP(x,:), I, J, 'UniformOutput', false); 
%Calculate vij for each point 
temp_vij1 = arrayfun(@(x, y) cross(v, rilj{x, y})/norm(rilj{x, y}), J, I, 'UniformOutput', false); 
%Reshape the matrix into desired format 
temp_vij2 = cell2mat(temp_vij1); 
vij = cat(3, temp_vij2(:, 1:3:end), temp_vij2(:, 2:3:end), temp_vij2(:, 3:3:end));