CODEWITHSUNDEEP
matlabrotation
Jawad Faisal
Jawad Faisal

Reputation: 43

Intersection volume of ellipsoids

I plot two ellipsoids with ellipsoid function and also use rotation handle. In order to compute the intersection volume (i.e. lens) and for the case of intersection volume between ellipsoids without rotation , i used analytical code for spheres intersection volume. But I am stuck ,if ellipsoids are at some rotations ,how i will compute the intersection volume(i.e.lens) between them. Also i need to compute intersection dia. of lens. ellipsoids have different rotational angles with respect to their max axis along z axis. here is a code:

 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

Upvotes: 3

Views: 3588

Answers (1)

Anton
Anton

Reputation: 4684

As far as I know there is no such function in Matlab. I would use the Monte Carlo method in order to find an approximated volume of the intersection.

You need to generate random points from some box and check if they belong to the intersection or not. You count how many of the points fell into intersection. Knowing the volume of the box, the number of the points fallen into intersection and the overall number of generated points you can then calculate the desired volume.

On the picture you can see the red points inside the intersection:

Monte Carlo method to find volume of the intersection

It is important to find a number of iterations which is good enough for your problem. I validated the method by calculating the volume of a small ellipsoid, which was put inside of a bigger one. In this case the intersection is the small ellipsoid itself, whose volume can be found analytically.

You can see the relative error as function of the iteration number on the next picture:

Relative error of the Monte Carlo method

Here is the code:

The code works for the rotation only around the Z-axis

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

Function to find if a point belongs to the ellipsoid (it does if val is <= 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

In order to find the maximum axis of the common body you can do the following:

  1. Save all points which fall into the intersection in an array (R_arr in the code)

  2. Construct a convex hull from the resulting array:

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

    It will give you an approximate volume v and indices of the points, used to build the hull K

  3. Now you have a subset of your points, which lay on the surface of the intersection body. All points exist in the array several times. Let's get rid of duplicates:

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

  4. Now the array is pretty small (about 300 points for the example above) and you can go through it and find the longest distance. I wrote function findAxis for it.

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

    It returns the furthest points A and B, and distance d between them.

  5. Now when you know these two points, you can define another points belonging to the axis. Use function calc_coord:

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

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

    I used values from the array ranges in order to specify the z-coordinates of the axis.

  6. Plot the axis!

For the example above the plot looks like this:

Axis of the intersection body

The volume of the convex hull is 5.933. The Monte Carlo method gave us 6.1889, so the results are pretty close to each other.

The length of the longest axis is 4.3264.

To get the intuition which volume approximation is better, I calculated the relative error:

Relative error for Monte Carlo method and the convex hull

As you can see the Monte Carlo method provides much better accuracy even after few iterations.

Here is the updated code with function for axis calculation:

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;

Functions:

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

EXPLANATIONS

function calc_ellipsoid

The function is based on the ellipsoid equation (ellipsoid rotated about the z-axis by angle theta):

ellipsoid equation

The function defines if a given point lays inside the ellipsoid or not

%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

If you need to work with other shapes, you need to find such equations, to check if your points belong to the shape or not. The ellipsoid equation is one of the easiest cases.

function findAxis

The function searches for two points from a cloud of points, which have the longest distance between them. Here we assume that the axis of the intersection shape goes through such two points (it is a very rough assumption, but it works).

%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

function calc_coord

The function is based on the line equation in 3D space:

line equation

The equation defines a line which goes through 2 points with coordinates (x1, y1, z1) and (x2, y2, z2). We use this function two find two other points outside the intersection shape, in order to draw the axis more demonstrative. You don't really need this function, if you are not going to produce nice pictures.

%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

Upvotes: 2

Related Questions