橢圓體的相交量

Question

我繪製了橢球體的兩個橢球體，並使用旋轉手柄。爲了計算相交體積（即透鏡）和橢圓體之間沒有旋轉的相交體積的情況，我使用球體相交體積的分析代碼。 但是我被卡住了，如果橢圓體處於某種旋轉狀態，我將如何計算它們之間的相交體積（即透鏡）。另外我需要計算相交直徑。的鏡頭。 橢圓體相對於其沿z軸的最大軸具有不同的旋轉角度。這裏是一個代碼：

draw ellipsoid, [x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,n); 
(xc,yc,zc)=center; semi-axis lengths = (Fmax,Fmin,Fmin);n 
X=[0,2]; two ellipsoid x coordinates i.e 0 is first ellipsoid 
      and 2 is second ellipsoid respectively 
Y=[0,2]; two ellipsoid y coordinates 
Z=[0,2]; two ellipsoid z coordinates; 
ROTATIONANGLE=[90,30]; 
RMIN=[1,2]; two ellipsoid minor axis radius 
RMAX=[3,5]; two ellipsoid major axis radius. 
for a = 1:2 display ellipsoid 
     [x, y, z] = ellipsoid(X(a),Y(a),Z(a),RMIN(a),RMAX(a),3.25,30); 
S = surfl(x, y, z); 
    rotate(S,[0,0,1],ROTATIONANGLE(a)) 
colormap copper 
axis equal 
xlabel('X') 
ylabel('Y')`enter code here 
zlabel('Z') 
check the boolean condition 
hold on 
end 

Answer 1

據我所知在Matlab中沒有這樣的功能。我會使用Monte Carlo method以找到交點的近似體積。

您需要從某個框中生成隨機點並檢查它們是否屬於交叉點。你計算有多少點落入交點。知道盒子的體積，交點的點數和生成點的總數就可以計算出所需的體積。

圖片上可以看到路口內的紅點：

它找到一個迭代次數這是你的問題不夠好是很重要的。我通過計算一個小橢圓體的體積來驗證這個方法，這個橢圓體放在一個更大的橢圓體內。在這種情況下，交叉點就是小橢圓本身，其體積可以通過分析找到。

你可以看到相對誤差爲下次畫面的迭代次數的函數：

下面是代碼：

的代碼適用於旋轉僅繞Z軸

X=[0,2]; %two ellipsoid x coordinates 
Y=[0,2]; %two ellipsoid y coordinates 
Z=[0,2]; %two ellipsoid z coordinates 
ROTATIONANGLE=[90,30]; 
AXIS_A=[1,2]; %two ellipsoid minor axis radius 
AXIS_B=[3,5]; %two ellipsoid major axis radius 
AXIS_C=[3.25,3.25]; %two ellipsoid major axis radius 

ranges = zeros(3, 2); 

do_plot = 1; %either plot ellipsoids or not 

step_number = 100000; 

for i = 1:2 %display ellipsoid 

    if (do_plot == 1) 
     [x, y, z] = ellipsoid(X(i),Y(i),Z(i),AXIS_A(i),AXIS_B(i),AXIS_C(i),20); 
     S = surf(x, y, z); 
     alpha(.1); 
     rotate(S,[0,0,1], ROTATIONANGLE(i), [X(i),Y(i),Z(i)]); 
    end 

    %calculate the ranges for the simulation box 
    ranges(1, 1) = min(ranges(1, 1), X(i) - max(AXIS_A(i), AXIS_B(i)) ); 
    ranges(1, 2) = max(ranges(1, 2), X(i) + max(AXIS_A(i), AXIS_B(i)) ); 

    ranges(2, 1) = min(ranges(2, 1), Y(i) - max(AXIS_A(i), AXIS_B(i)) ); 
    ranges(2, 2) = max(ranges(2, 2), Y(i) + max(AXIS_A(i), AXIS_B(i)) ); 

    ranges(3, 1) = min(ranges(3, 1), Z(i) - AXIS_C(i)); 
    ranges(3, 2) = max(ranges(3, 2), Z(i) + AXIS_C(i)); 

    if (do_plot == 1) 
     hold on; 
    end 
end 

counter = 0; %how many points targeted the intersection 

for i = 1:step_number 

    R = rand(3, 1).*(ranges(:, 2) - ranges(:, 1)) + ranges(:, 1); %a random point 

    n = 1; 
    val = calc_ellipsoid(R(1), R(2), R(3), X(n),Y(n),Z(n),AXIS_A(n),AXIS_B(n),AXIS_C(n),ROTATIONANGLE(n)*pi/180); 

    if (val <= 1.0) 
     n = 2; 
     val = calc_ellipsoid(R(1), R(2), R(3), X(n),Y(n),Z(n),AXIS_A(n),AXIS_B(n),AXIS_C(n),ROTATIONANGLE(n)*pi/180); 

     if (val <= 1.0) 

      if (do_plot == 1) 
       plot3(R(1), R(2), R(3), 'or', 'MarkerSize', 1, 'MarkerFaceColor','r'); 
      end 

      counter = counter + 1; 
     end 
    end 

end 

cube_vol = 1; %the volume of the simulation box 
for i=1:3 
    cube_vol = cube_vol * (ranges(i, 2) - ranges(i, 1)); 
end 

%approximated volume of the intersection 
section_vol = cube_vol * counter/step_number; 

display(['Cube volume: ', num2str(cube_vol)]); 
display(['Targeted points: ', num2str(counter), ' from ', num2str(step_number)]); 
display(['Section volume: ', num2str(section_vol)]); 

if (do_plot == 1) 
    hold off; 
end 

函數來查找一個點是否屬於第橢圓體（它的確如果val是< = 1。0）：

function [ val ] = calc_ellipsoid(x, y, z, x0, y0, z0, a, b, c, theta) 

    x_cmp = ((x - x0)*cos(theta) + (y - y0)*sin(theta))^2 /(a^2); 
    y_cmp = ((x - x0)*sin(theta) - (y - y0)*cos(theta))^2 /(b^2); 
    z_cmp = (z - z0)^2/(c^2); 

    val = x_cmp + y_cmp + z_cmp; 
end 

UPDATE

爲了找到你可以做以下的共同體內的最大軸：

1. 保存所有點落入相交處一個數組（代碼中的R_arr）

2. 從結果數組構造凸包：

   [K, v] = convhull(R_arr(:, 1), R_arr(:, 2), R_arr(:, 3));

   它會給你一個近似體積v和點的指數，用於構建船體K

3. 現在你有你的觀點的一個子集，裏面躺着的表面上相交體。所有點都在陣列中存在好幾次。讓我們擺脫重複的：

   K = K(:); K = unique(K, 'stable');

4. 現在陣列是非常小（上面的例子中約300點），你可以通過它找到的最長距離。我爲它寫了功能findAxis。

   [A, B, d] = findAxis(R_arr(K, :));

   它返回的最遠點A和B，以及它們之間的距離d。

5. 現在，當您知道這兩點時，您可以定義屬於該軸的另一個點。使用功能calc_coord：

   C = calc_coord(A, B, ranges(3, 1));

   D = calc_coord(A, B, ranges(3, 2));

   予以指定軸線的座標 - z使用從陣列ranges值。

6. 繪製軸！

對於情節上面的例子如下所示：

的凸包的體積爲5.933。蒙特卡洛方法給了我們6.1889，所以結果非常接近。

最長軸的長度爲4.3264。

要獲得直覺它的體積近似是更好的，我計算出的相對誤差：

正如你所看到的蒙特卡羅方法即使經過幾次迭代提供更好的精度。

這裏是更新的代碼與功能軸計算：

X=[0,2]; %two ellipsoid x coordinates 
Y=[0,2]; %two ellipsoid y coordinates 
Z=[0,2]; %two ellipsoid z coordinates 
ROTATIONANGLE=[90,30]; 
AXIS_A=[1,2]; %two ellipsoid minor axis radius 
AXIS_B=[3,5]; %two ellipsoid major axis radius 
AXIS_C=[3.25,13.25]; %two ellipsoid major axis radius 

ranges = zeros(3, 2); 

do_plot = 1; %either plot ellipsoids or not 

step_number = 1000000; 

for i = 1:2 %display ellipsoid 

    if (do_plot == 1) 
     [x, y, z] = ellipsoid(X(i),Y(i),Z(i),AXIS_A(i),AXIS_B(i),AXIS_C(i),20); 
     S = surf(x, y, z); 
     alpha(.05); 
     rotate(S,[0,0,1], ROTATIONANGLE(i), [X(i),Y(i),Z(i)]); 
    end 

    %calculate the ranges for the simulation box 
    ranges(1, 1) = min(ranges(1, 1), X(i) - max(AXIS_A(i), AXIS_B(i)) ); 
    ranges(1, 2) = max(ranges(1, 2), X(i) + max(AXIS_A(i), AXIS_B(i)) ); 

    ranges(2, 1) = min(ranges(2, 1), Y(i) - max(AXIS_A(i), AXIS_B(i)) ); 
    ranges(2, 2) = max(ranges(2, 2), Y(i) + max(AXIS_A(i), AXIS_B(i)) ); 

    ranges(3, 1) = min(ranges(3, 1), Z(i) - AXIS_C(i)); 
    ranges(3, 2) = max(ranges(3, 2), Z(i) + AXIS_C(i)); 

    if (do_plot == 1) 
     hold on; 
    end 
end 

counter = 0; %how many points targeted the intersection 

R_arr = []; 

for i = 1:step_number 

    R = rand(3, 1).*(ranges(:, 2) - ranges(:, 1)) + ranges(:, 1); %a random point 

    n = 1; 
    val = calc_ellipsoid(R(1), R(2), R(3), X(n),Y(n),Z(n),AXIS_A(n),AXIS_B(n),AXIS_C(n),ROTATIONANGLE(n)*pi/180); 

    if (val <= 1.0) 
     n = 2; 
     val = calc_ellipsoid(R(1), R(2), R(3), X(n),Y(n),Z(n),AXIS_A(n),AXIS_B(n),AXIS_C(n),ROTATIONANGLE(n)*pi/180); 

     if (val <= 1.0) 

      if (do_plot == 1) 
       %plot3(R(1), R(2), R(3), 'or', 'MarkerSize', 1, 'MarkerFaceColor','r'); 
      end 

      counter = counter + 1; 

      R_arr = [R_arr; R']; 
     end 
    end 

end 

cube_vol = 1; %the volume of the simulation box 
for i=1:3 
    cube_vol = cube_vol * (ranges(i, 2) - ranges(i, 1)); 
end 

%approximated volume of the intersection 
section_vol = cube_vol * counter/step_number; 

display(['Cube volume: ', num2str(cube_vol)]); 
display(['Targeted points: ', num2str(counter), ' from ', num2str(step_number)]); 
display(['Section volume: ', num2str(section_vol)]); 

if (counter > 0) 

    %hull 
    [K, v] = convhull(R_arr(:, 1), R_arr(:, 2), R_arr(:, 3)); 


    display(['Approximated volume: ', num2str(v)]); 
    trisurf(K, R_arr(:, 1), R_arr(:, 2), R_arr(:, 3), 'FaceColor','cyan') 
    hold on; 
    axis equal 

    K = K(:); 
    K = unique(K, 'stable'); 

    [A, B, d] = findAxis(R_arr(K, :)); 

    plot3(A(1, 1), A(1, 2), A(1, 3), 'or', 'MarkerSize', 10, 'MarkerFaceColor','r'); 
    plot3(B(1, 1), B(1, 2), B(1, 3), 'or', 'MarkerSize', 10, 'MarkerFaceColor','r'); 

    %axis 
    C = calc_coord(A, B, ranges(3, 1)); 
    D = calc_coord(A, B, ranges(3, 2)); 

    p = [C; D]; 

    plot3(p(:, 1), p(:, 2), p(:, 3), '--g', 'LineWidth', 3); 

else 
    display('There is no intersection!'); 
end 

hold off; 

功能：

function [ C] = calc_coord(A, B, z) 
    x1 = A(1, 1); 
    x2 = B(1, 1); 

    y1 = A(1, 2); 
    y2 = B(1, 2); 

    z1 = A(1, 3); 
    z2 = B(1, 3); 

    y = (y2 - y1)*(z - z1)/(z2 - z1) + y1; 
    x = (x2 - x1)*(y - y1)/(y2 - y1) + x1; 

    C = [x, y, z]; 
end 

function [ A, B, d ] = findAxis(point_arr) 

    A_ind = 0; B_ind = 0; d = -1; 

    n = size(point_arr, 1); 
    i = 1; 

    while (i < n) 

     for j=(i+1):n 

      cur_d = norm(point_arr(i, :) - point_arr(j, :)); 

      if (cur_d > d) 
       d = cur_d; 
       A_ind = i; 
       B_ind = j; 
      end 
     end 

     i = i + 1; 
    end 

    A = point_arr(A_ind, :); 
    B = point_arr(B_ind, :); 
end 

解釋

功能calc_ellipsoid

該函數是基於所述橢球方程（橢球由角度theta繞z軸）：

的函數定義，如果一個給定的點的橢圓體的內部規定或不

%the function calculates a value for some given point, which shows if the 
%point is inside of the ellipsoid or not. 
%for a point to be inside, the value has to be <=1 

function [ val ] = calc_ellipsoid(x, y, z, x0, y0, z0, a, b, c, theta) 

    %x, y, z - coordinates of the point to be checked 
    %x0, y0, z0 - center of the ellipsoid 
    %a, b, c - axes of the ellipsoid 
    %theta - angle of the rotation about the Z-axis 

    x_cmp = ((x - x0)*cos(theta) + (y - y0)*sin(theta))^2 /(a^2); 
    y_cmp = ((x - x0)*sin(theta) - (y - y0)*cos(theta))^2 /(b^2); 
    z_cmp = (z - z0)^2/(c^2); 

    val = x_cmp + y_cmp + z_cmp; 
end 

如果您需要使用其他形狀，您需要找到這樣的方程，以檢查您的點是否屬於該形狀。橢球方程是最簡單的情況之一。

功能findAxis

兩個點的函數搜索從分，這讓他們之間的距離最長的雲。在這裏，我們假設交點形狀的軸線貫穿這兩點（這是一個非常粗略的假設，但它的工作原理）。

%the function is given an array of all points which lay on the constructed 
%body. we assume, that the axis of the body is the longest line, connecting 
%some two points of this array. the function goes through all the points 
%and calculates the euclidian distance between them. 
%the function returns points A and B, and the distance between them 

function [ A, B, d ] = findAxis(point_arr) 

    %point_arr - points from the bodies surface 
    %A, B - end points of the found exis 
    %d - distance between A and B 

    A_ind = 0; B_ind = 0; d = -1; 

    n = size(point_arr, 1); 
    i = 1; 

    while (i < n) 
     for j=(i+1):n 

      cur_d = norm(point_arr(i, :) - point_arr(j, :)); 

      if (cur_d > d) 
       d = cur_d; 
       A_ind = i; 
       B_ind = j; 
      end 
     end 

     i = i + 1; 
    end 

    A = point_arr(A_ind, :); 
    B = point_arr(B_ind, :); 
end 

功能calc_coord

的函數基於在三維空間中的直線方程：

方程式定義了一個線，其通過2點座標(x1, y1, z1)並進入(x2, y2, z2)。我們使用這個函數找出兩個交點形狀外的其他點，以便繪製出更加明顯的軸。如果你不打算拍出漂亮的照片，你並不需要這個功能。

%the function is given two points A and B, which were found in the 
%'findAxis' function. we want to find a new point C, which lay on the axis, 
%and have z-coordinate equal to 'z'. it is done merely to visualize the 
%axis in a better way. 

function [C] = calc_coord(A, B, z) 
    x1 = A(1, 1); 
    x2 = B(1, 1); 

    y1 = A(1, 2); 
    y2 = B(1, 2); 

    z1 = A(1, 3); 
    z2 = B(1, 3); 

    y = (y2 - y1)*(z - z1)/(z2 - z1) + y1; 
    x = (x2 - x1)*(y - y1)/(y2 - y1) + x1; 

    C = [x, y, z]; 
end