Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Determine which points of a grid belong to a certain volume

Subject: Determine which points of a grid belong to a certain volume

From: Andreas Frölich

Date: 14 Oct, 2010 14:29:04

Message: 1 of 5

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 z-direction and its short axis are of equal length rx and lie in the x-y-plane.

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

Subject: Determine which points of a grid belong to a certain volume

From: Sean

Date: 14 Oct, 2010 14:53:04

Message: 2 of 5

"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 z-direction and its short axis are of equal length rx and lie in the x-y-plane.
>
> 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

Well you know the formula for a 3D ellipsoid so this is a basic logical scheme:

[xx yy zz] = meshgrid(1:Nx,1:Ny,1:Nz);

These are three matrices which contain the coordinates to every point in your block.
Now write the formula for the ellipsoid in terms of the center and these three matrices with some maximum radii.

See this thread for the equivalent with just a 2D circle:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/291115#777679

Subject: Determine which points of a grid belong to a certain volume

From: ImageAnalyst

Date: 14 Oct, 2010 16:36:23

Message: 3 of 5

Andreas:
That is done by a morphological dilation.
http://en.wikipedia.org/wiki/Dilation_%28morphology%29
You just need to set up the kernel that defines the ellipsoid, and
draw your path into your volume (set those pixels to 1), and then call
imdilate().
-ImageAnalyst

Subject: Determine which points of a grid belong to a certain volume

From: Roger Stafford

Date: 14 Oct, 2010 19:10:05

Message: 4 of 5

"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 z-direction and its short axis are of equal length rx and lie in the x-y-plane.
>
> 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 (infinitely-long) 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

 (x-sx)^2/a^2 + (y-sy)^2/b^2 + (z-sz)^2/c^2 <= 1 (1)

 (x-ex)^2/a^2 + (y-ey)^2/b^2 + (z-ez)^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*(x-sx)/a^2, 2*(y-sy)/b^2, 2*(z-sz)/c^2, so for such points we must have

 (x-sx)/a^2*(ex-sx) + (y-sy)/b^2*(ey-sy) + (z-sz)/c^2*(ez-sz) = 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

 (x-sx)*(ex-sx)/a^2 + (y-sy)*(ey-sy)/b^2 + (z-sz)*(ez-sz)/c^2 >= 0 (3)

 (x-ex)*(ex-sx)/a^2 + (y-ey)*(ey-sy)/b^2 + (z-ez)*(ez-sz)/c^2 <= 0 (4)

  If we define dx = ex-sx, dy = ey-sy, and dz = ez-sz, and vary the parameter t as we move the ellipsoid, it satisfies

 (x-sx-dx*t)^2/a^2 + (y-sy-dy*t)^2/b^2 + (z-sz-dz*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 above-mentioned 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 non-negative:

 B^2-4*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

Subject: Determine which points of a grid belong to a certain volume

From: Roger Stafford

Date: 15 Oct, 2010 00:41:03

Message: 5 of 5

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i97kic$8p4$1@fred.mathworks.com>...
> 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 (infinitely-long) 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
> ........
- - - - - - - - - - -
  It has occurred to me that you could simplify the analysis of your problem by first transforming all the z coordinates of your mesh points by

 Z = a/c*z

Remembering that your a and b are equal, this transforms your ellipsoids into spheres and you could then test for the transformed points lying in a circular cylinder with spherical hemispheres at each end, an easier analysis to perform.

Roger Stafford

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us