2014-09-28 54 views
0

我有一些形式的圓柱線源,我用它來通用3D核心;我通過首先生成一個簡單的二維核心來做到這一點,如下所示:在MATLAB多維數組中旋轉元素?

A = 625; 
B = 25; 

%2D kernal 
grid2d = zeros(101,101); 
c = 51; %centre 

for m = 1:101 
    for n = 1:101 

    r(m,n) = sqrt((m - c).^2 + (n - c).^2); %distance of point on grid to centre 

    %Populating the grid as a kernal 

    if r(m,n) <= 5 

     grid2d(m,n) = 100; 

    elseif r(m,n) >= 25 

     grid2d(m,n) = 0; 

    else 

     grid2d(m,n) = A./r(m,n) - B; 

    end 
    end 
end 

這給了我一個2D核心。現在,如果我將3D版本定義爲沿着更大網格內z軸的9個元素,我可以通過以下創建3D內核;

gz = 147:155; %9 elements in the z axis 
H = length(gz); 
kernel3D = zeros(301,301,301); % 

for n = 1:H 

    kernel3D(151-50:1:151+50,151-50:1:151+50,gz(n)) = grid2d; 
end 

這工作完全,如果我用垂直線源,但我很好奇它是否可以旋轉這個數組的元素在任何需要的方向,所以我可能會爲斜源內核;例如,假設我想將這個數組相對於XY平面旋轉45度,並且在(151,151,151)處通過線源中心將XZ平面旋轉60度?

有沒有一種優雅的方式來做到這一點,也許使用旋轉矩陣?

回答

0

我想我找到了一個工作但不雅的解決方案 - 在這個版本中,我生成了2D內核,並將它嵌套在一個3D數組中。然後我使用旋轉操作器將適當的值映射到旋轉變換的位置。這是醜陋的,但似乎工作。

%generate a test 2d kernel as before... 

A = 625; B = 25; c = 51; cen = 151; grid2d = zeros(101,101); 

for m = 1:101 
for n = 1:101 


     r(m,n) = sqrt((m - c).^2 + (n - c).^2); %distance of point on grid to centre 

%Populating the grid as a kernal 

if r(m,n) <= 5 

    grid2d(m,n) = 100; 

elseif r(m,n) >= 25 

    grid2d(m,n) = 0; 

else 

    grid2d(m,n) = A./r(m,n) - B; 

    end 
    end 
    end 

%make a bigger array for this to nest in 
flat3d = zeros(301,301,301); 
range = 101:201; %centre range of new grid 

flat3d(range,range,cen) = grid2d; % inserts grid2d at centre of new flat3d grid! 

%now the incline angles 
theta = input('Enter incline in XZ (vertical) plane in degrees : '); 
phi = input('Enter incline in XY (horizontal) plane in degrees : '); 

%Now a loop which inclines each point relative to theta around (x,y,z) = (cen,cen,cen); 

g = range(1); 

inc3d = zeros(301,301,301); %new inclined grid initialised 

for m = 1:101 
for n = 1:101 


    %the g - cen value is added to the index to map arrays sync up correctly and 
    %keep loop reasonable/avoid overlap outside index ! 

    xpart(m,n) = cen + cosd(theta).*(g + m -cen) - sind(theta).*(g + n - cen); 
    xn(m,n) = cen + cosd(phi).*(xpart(m,n) - cen) - sind(phi).*(g + n - cen); %this maps a value of x to new position theta and phi 
    yn(m,n) = cen + sind(phi).*(xn(m,n) - cen) + cosd(phi).*(g + n - cen); %maps a value of y to new position relative to phi 
    zn(m,n) = cen + sind(theta).*(g + m - cen) + cosd(theta).*(g + n - cen); %maps z to new position relative to theta (doesnt change relative to phi) 

    xn = round(xn); 
    yn = round(yn); %round to avoid integer problems with indexing 
    zn = round(zn); 


    inc3d(xn(m,n),yn(m,n),zn(m,n)) = flat3d(g + m,g + n,cen); 

end 
end 

這似乎工作;例如,如果我生成theta = 30和phi = 15的plot3圖片,我會得到類似的東西;

enter image description here

從這個二維平面圖像,與直線的方程,我覺得一個可以簡單地總結和生成所需要的這曾經地區內核。這很混亂,但似乎工作。如果有人知道更清潔的方式,通過一切手段添加它:)