2

I want to calculated the mass of a sphere based on a threedimensional discret inhomogenous density distribution. Lets say a set of 3x3x3 cubes of different densities is inscribed by a sphere. What is the fastest way to sum up the partitioned masses using Python?

I tried to calculate the volume under the mathematical equation for a sphere: x^2 +y^2 +z^2 = R^2 for the range of one of the cubes using scipy.integrate.dblquad. However, the result is only valid if the boundaries are smaller than the radius of the sphere and repetitive calculation for lets say 50'000 spheres with 27 cubes each would be quite slow.

On the other hand, the usual equation for CoM caluations could't be used in my opinion, due to the rather coarse and discrete mass distribution.

ProgramFOX
  • 6,131
  • 11
  • 45
  • 51
surcouf
  • 21
  • 2
  • 1
    "smaller then the radius" should be "smaller than the radius". Use "than" for comparisons, "then" for time sequence. – aaaa bbbb Dec 10 '10 at 19:03
  • If the cubes are homegeneous ...just add up the cubes mass. If the cubes are not homogeneous, let's see the density formula. You don't need the CoM for calculating the mass. – Dr. belisarius Dec 10 '10 at 19:08
  • I think this is not a real question unless the density of the cubes is specified. – Dr. belisarius Dec 10 '10 at 19:17
  • I think you'd have better luck asking this on http://math.stackexchange.com/ , first clarifying whether the cubelets are of each of uniform density, and also what you mean by 'inscribed by' - I think you mean the sphere inscribes the cube made up by the 27 cubelets, but that's just a guess. If both those are true you don't need numerical integration, just calculus. – AakashM Dec 10 '10 at 20:02
  • @AakashM:Your guess is correct. Each cubelet is described by a uniform density, but a different from the the adjecent cublets. With "inscribed by" I tried to describe the position of the sphere in the set of 3x3x3 cubes. The sphere could be somewhere in the space of the set. Thats the reason why I tried numerical integration, because the equation of the sphere is known. All I need - in my opinion - is a method to calculate the volume below the sphere, which belongs to each cublet. – surcouf Dec 12 '10 at 21:09
  • A second idea I had recently was to sum up subvolumes of the sphere and to use a lookup-function for the local density of each subvolume. This should be rather easy, but accuracy and speed are influenced by the choosen subvolume size. Any experience or hints on this approach? – surcouf Dec 12 '10 at 21:14

3 Answers3

1

Timing Experiment

You didn't specify your timing constraints, so I've done a little experiment with a nice integration package.

Without optimization, each integral in spherical coordinates can be evaluated in 0.005 secs in a standard laptop if the cubes densities are straightforward functions.

Just as a reference, this is the program in Mathematica: 

Clear@f;
(* Define a cuboid as density function *)
iP = IntegerPart;
f[{x_, y_, z_}, {lx_, ly_, lz_}] :=   iP[x - lx] + iP[y - ly] + iP[z - lz] /; 
   lx <= x <= lx + 3 && ly <= y <= ly + 3 && lz <= z <= lz + 3;

f[{x_, y_, z_}, {lx_, ly_, lz_}] := Break[] /; True;

Timing[Table[s = RandomReal[{0, 3}, 3]; (*sphere center random*)
   sphereRadius = Min[Union[s, 3 - s]]; (*max radius inside cuboid *)
   NIntegrate[(f[{x, y, z} - s, -s] /.  (*integrate in spherical coords *)
       {x -> r Cos@th Sin@phi, 
        y -> r Sin@th Sin@phi, 
        z -> r Cos@phi}) r^2 Sin@phi,
       {r, 0, sphereRadius}, {th, 0, 2 Pi}, {phi, 0, Pi}], 
         {10000}]][[1]]

The result is 52 secs for 10^4 iterations.

So perhaps you don't need to optimize a lot ...

Dr. belisarius
  • 60,527
  • 15
  • 115
  • 190
0

I cannot get your exact meaning of inscribed by a sphere. Also I havent tried the scipy.integrate. However, here are some though:

Set a 3x3x3 cube to unit density. Then take the integration for each cube respectively, so you should have the volume cube V_ijk here. Now for each of sphere, you can get the mass of each sphere by summing V_ijk*D_ijk, where the D_ijk is the density of the sphere.

It should be much faster because you do not need to do integration now.

unsym
  • 2,080
  • 1
  • 20
  • 19
0

You can obtain an analytic formula for the intersecting volume between a cube (or rectangular prism) and a sphere. It won't be easy, but it should be possible. I have done it for an arbitrary triangle and circle in 2D. The basic idea is to decompose the intersection into simpler pieces, like tetrahedra and volumetric spherical triangle sectors, for which relatively simple volume formulas are known. The main difficult is in considering all the possible cases of intersections. Luckily both objects are convex, so you are guaranteed a single convex intersection volume.

An approximate method might be to simply subdivide the cubes until your approximate numerical integration algorithm does work; this should still be relatively fast. Do you know about Pick's Theorem? That only works in 2D, but there are, I believe, 3D generalizations.

Victor Liu
  • 3,545
  • 2
  • 24
  • 37