"Andreas Frölich" <ThisPartBlocksSpam_andreas.froelich@kit.edu> wrote in message <i9743g$jip$1@fred.mathworks.com>...
> Hi everyone,
>
> I am looking for a good way to solve the following task:
>
> Given a 3d grid g=zeros(Nx,Ny,Nz)
> corresponding to a discretization of a block in physical space with the dimensions ax,ay and az I would like to determine, which points of g belong to the volume swept out by moving the center of an ellipsoid along a path inside the volume that g is a discretization of. The path is specified by a starting point [sx sy sz] and an endpoint [ex ey ez]. The ellipsoid has its long axis rz along the zdirection and its short axis are of equal length rx and lie in the xyplane.
>
> I hope this is a clearly put question, if not I will be happy to provide additional information.
>
> Thanks in advance for your help,
>
> Andreas
           
I will approach this problem from a mathematical point of view. If you had a sphere rather than an ellipsoid, the volume in question would consist of the inside of a circular cylinder with a hemispherical cap at either end. Testing that an arbitrary point (x,y,z) lies within that volume would be accomplished by determining if the point was a) inside either end sphere or b) simultaneously within the (infinitelylong) cylinder and between the two planes where the cylinder contacts each sphere.
With an ellipsoid the method is similar though more complicated. The points inside the "start" and "end" ellipsoids satisfy
(xsx)^2/a^2 + (ysy)^2/b^2 + (zsz)^2/c^2 <= 1 (1)
(xex)^2/a^2 + (yey)^2/b^2 + (zez)^2/c^2 <= 1 (2)
respectively. (I generalized to a, b, and c being all different.)
As the ellipsoid is moved from "start" to "end", an ellipsoidal cylinder is swept out by the ellipsoid. The points on either ellipsoid that make contact with this cylinder must have their surface normal orthogonal to a vector in the direction of motion from (sx,sy,sz) to (ex,ey,ez). For the start ellipsoid the surface normal points in the direction of the three partial derivatives 2*(xsx)/a^2, 2*(ysy)/b^2, 2*(zsz)/c^2, so for such points we must have
(xsx)/a^2*(exsx) + (ysy)/b^2*(eysy) + (zsz)/c^2*(ezsz) = 0
which defines a plane, just as with the above spheres. There is a similar parallel plane for the end ellipsoid, and the resulting inequalities for a point to lie between the planes would be
(xsx)*(exsx)/a^2 + (ysy)*(eysy)/b^2 + (zsz)*(ezsz)/c^2 >= 0 (3)
(xex)*(exsx)/a^2 + (yey)*(eysy)/b^2 + (zez)*(ezsz)/c^2 <= 0 (4)
If we define dx = exsx, dy = eysy, and dz = ezsz, and vary the parameter t as we move the ellipsoid, it satisfies
(xsxdx*t)^2/a^2 + (ysydy*t)^2/b^2 + (zszdz*t)^2/c^2 = 1
This equation can (after a little sweat or expeditious use of matlab's symbolic toolbox) be put in the form
A*t^2 + B*t + C = 0 (5)
where for example A = dx^2/a^2+dy^2/b^2+dz^2/c^2. (I'll let you work out B and C.) If a point (x,y,z) is to lie inside the abovementioned cylinder, then there must exist a value for t which satisfies the quadratic (5), and this is possible if and only if its discriminant is nonnegative:
B^24*A*C >= 0 (6)
Therefore the equation for the inside of the cylinder is (6) with A, B, and C determined as described above.
Thus we now have the criterion for an arbitrary point (x,y,z) to lie within the volume in question. It must either satisfy inequality (1), or satisfy inquality (2), or else it must satisfy inequalities (3), (4), and (6) simultaneously.
I realize this still requires some work on your part working out the details of inequality (6), but what we have here is a balance between the effort of doing this and the inefficiency of a more crude kind of matlab code that would place the ellipsoid in tiny successive steps along the path, testing all points in the grid each time. This latter could constitute a rather wasteful order n^3 type algorithm.
Roger Stafford
