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

Surface area of a plane within a cube

Asked by Cole on 21 Aug 2014
Latest activity Commented on by Cole on 28 Aug 2014

Hi,

I'm trying to solve a geometric problem. I have grid blocks which are represented as a cube. I also have rectangles on a plane that are orientated in and around the cube. I need to solve how much of the plane is within the cube.

If these planes were infinite, it would be easy. Since they are not infinite, it's not so trivial. Sometimes the corner points lie within the cube, sometimes outside the cube.

I've posted an example of one possible scenario.

This is a representation of a quartz vein in a block of rock. I'm simulating lots of these so the scenario can be varied.

So far, I check if all 4 points are within the bounds of the cube. If they are, the area is just calculated as is. If one point falls outside the cube (usually this is the case) then I project planes normal to x, normal to y and normal to z and this gives me the line of intersection. From here, I don't know how to solve for area outside the cube.

I'm really stuck on how to solve this issue. Some help would be great.

0 Comments

Cole

Tags

Products

No products are associated with this question.

2 Answers

Answer by Mike Garrity on 22 Aug 2014
Accepted answer

I don't have a code fragment handy, but a simple way to do this is with the Sutherland-Hodgman clipping algorithm.

Wikipedia article on Sutherland-Hodgman

Basically, you treat your plane as a ordered list of points and your cube as a set of 6 planes. For each plane, you loop through the list of points keeping track of when they cross the plane, and make a new list of points. You then use that new list for the next plane. When you're done, you've got the polygon which is the intersection of the cube and the plane. In this case, that polygon could be anywhere from empty to 6-sided.

If you google, you'll probably find MATLAB implementations of Sutherland-Hodgman, but it's the simplest of the polygon clipping algorithms to implement if you'd like to give it a go yourself. Note that it's limited to the case where the clipping volume is convex, but that's true of your cube in this case.

1 Comment

Cole on 22 Aug 2014

Hi Mike,

Thanks for the answer.

I found a code snippet from Rosetta Code:

Rosetta Code

It appears to be in 2D only. Do you think I could adapt this code to work in 3D?

As of now, I've solved the problem by brute force. I discretized each plane into a number of points and I count the number of points inside each block. It's kind of a Monte Carlo integration.

It works well but it's computationally expensive. With enough discretization points, the error is negligible. This is fine on a small synthetic example but for anything real world, the number of points would be in the trillions and it would just take forever to solve it all.

Cole

Mike Garrity
Answer by Mike Garrity on 28 Aug 2014

> Do you think I could adapt this code to work in 3D?

Sure, basically you just change the function Intersection so that it takes a plane equation and 3D point as input instead of a line equation and 2D point. Everything else stays exactly the same.

It would look something like this:

function t = Intersection(pt1,pt2,plane)
  t = nan;
  u = dot(pt1,plane(1:3)) - plane(4);    
  v = dot(pt2,plane(1:3)) - plane(4);    
  if (sign(u) ~= sign(v))
      t = (0-u) / (v-u);
  end
end

and here's an example of using it:

% define the 6 planes of a box
xrange = [-2 3];
yrange = [-2 4];
zrange = [-2 5];
planes = [1  0  0  xrange(2);
         -1  0  0 -xrange(1);
          0  1  0  yrange(2);
          0 -1  0 -yrange(1);
          0  0  1  zrange(2);
          0  0 -1 -zrange(1)];
% intesect a bunch of random lines against the box
for i=1:1000
  pt1 = rand(1,3);
  pt2 = 10*randn(1,3);
  t = zeros(1,6);
  for i=1:6
    t(i) = Intersection(pt1,pt2,planes(i,:));
  end
  startpt = pt1;
  endpt = pt1+min(t)*(pt2-pt1);
  line([startpt(1) endpt(1)],[startpt(2) endpt(2)],[startpt(3) endpt(3)])
end
view(3)
axis square

1 Comment

Cole on 28 Aug 2014

Thank you! Huge help!

Mike Garrity

Contact us