I posted an article in povray.binaries.images about random placement
of non-intersecting spheres.
  I got 5000 spheres of diameter 0.2 placed randomly in a planar area
between <-10, 0, -10> and <10, 0, 10> in 15 seconds with my 3.4GHz P4.
When I increased the diameter of the spheres to 0.23 (a much tougher
job because there's less space to place all the 5000 spheres) the time
was 1:50.
  There's no trick in calculating the random location of the spheres,
ie. no speedups using repeating patterns or anything, but every single
sphere is placed with a honest
MinCoords+(MaxCoords-MinCoords)*<rand(Seed), 0, rand(Seed)>

  Basically I use the naive method of putting the center coordinate of
each new sphere in an array and each time I try to create a new sphere,
I first go through the array and check if it would be too close to any
of the other spheres (IOW. if the distance between the centers would
be smaller than the diameter of the spheres).

  Of course a linear search of the array is very slow and results in
an O(n^2) algorithm which gets increasingly slow as the number of
spheres grows. I thought how hard would it be to develop an algorithm
which could make this check faster than in linear time.
  Of course this is already difficult for spheres located in a plane,
not to talk about spheres freely located in space.
  Then I thought that it would require some kind of binary space
partitioning or any similar speeding-up algorithm.
  But then it just hit me: POV-Ray *already* has speeding algorithms
for checking points against objects much faster than in linear time.
Thus the only trick is how to get to use these already-existing
algorithms.

  There's a simple (at least in theory) way for this: Make a union of
spheres of double diameter from all the existing coordinates and then
check if the new center coordinate you calculated for the new sphere
is inside() this union. That's it. POV-Ray uses internally quite advanced
algorithms for speeding up this inside() test and it does it much faster
than in linear time.

  Of course there's a problem here: You just can't create the union of
spheres each time you want to create a new sphere. That would actually
be slower than the naive approach.

  One possible solution would be to build up a nested union of spheres.
In other words, each time you create a new sphere you would create a
union of this sphere and the existing union of spheres.
  However, the problem with this solution is, I believe, that it's really
resource-heavy and consumes lots of memory because each such union would
exist all at the same time in memory. It might also be somewhat slow to
create such a nested union each time you create a new sphere. (To be
honest, however, I didn't even try it to see what happens.)

  So what I ended up doing was a compromise between the naive solution
and the inside()-of-union solution:
  For the first 200 spheres (the number was searched by trial-and-error)
I create them by the naive method. Then I create a union of these 200
spheres and put them in an array.
  For the next 200 spheres I first test them against the union(s) in
that array (using inside()) and only if that test is negative I then
check it with the naive way with the newly-created coordinates. Then
I create another union of these new 200 spheres and put them in that
array.
  And so on for the next 200 spheres etc.

  The code below is not very optimized as I wrote it more or less hastily,
but it still resulted in the times I mentioned at the beginning.


#declare SpheresAmount = 5000;
#declare MinCoords = <-10, 0, -10>;
#declare MaxCoords = <10, 0, 10>;
#declare SphereDiameter = .2;
#declare Seed = seed(0);

#declare ChunkSize = 200;

#declare ChunksAmount = ceil(SpheresAmount/ChunkSize);
#declare Chunks = array[ChunksAmount];
#declare ChunkIndex = 0;

#declare CheckArray = array[ChunkSize];

#while(ChunkIndex < ChunksAmount)
  #declare CAInd = 0;
  #while(CAInd < ChunkSize)
    #declare CoordOk = false;
    #while(!CoordOk)
      #declare Coord =
      MinCoords+(MaxCoords-MinCoords)*<rand(Seed), 0, rand(Seed)>;
      #declare CoordOk = true;
      #declare Ind = 0;
      #while(Ind < ChunkIndex)
        #if(inside(Chunks[Ind], Coord))
          #declare CoordOk = false;
          #declare Ind = ChunkIndex;
        #end
        #declare Ind = Ind+1;
      #end
      #if(CoordOk)
        #declare Ind = CAInd-1;
        #while(Ind >= 0)
          #if(vlength(CheckArray[Ind]-Coord) < SphereDiameter)
            #declare CoordOk = false;
            #declare Ind = -1;
          #end
          #declare Ind = Ind-1;
        #end // #while(Ind >= 0)
      #end // #if
    #end // #while(!CoordOk)
    #declare CheckArray[CAInd] = Coord;
    #declare CAInd = CAInd+1;

    sphere
    { Coord, SphereDiameter/2
      pigment { rgb <rand(Seed), rand(Seed), rand(Seed)> }
      finish { specular .7 }
    }
  #end

  #declare ChunkIndex = ChunkIndex+1;
  #if(ChunkIndex < ChunksAmount)
    #declare Chunks[ChunkIndex-1] =
      union
      { #declare Ind = 0;
        #while(Ind < ChunkSize)
          sphere { CheckArray[Ind], SphereDiameter }
          #declare Ind = Ind+1;
        #end
      };
  #end

  #debug concat("Chunks: ", str(ChunkIndex,0,0), "/", str(ChunksAmount,0,0),
                "\n")
#end

camera { location <-.5, 25, -5>*1.5 look_at 0 angle 35 }
light_source { <100, 200, -50>, 1 }
plane
{ y, -SphereDiameter/2
  pigment { rgb <.8,.9,1> }
}

Post a reply to this message

Warp <war### [at] tagpovrayorg> wrote:
> I posted an article in povray.binaries.images about random placement
> of non-intersecting spheres.
> ... and so on..

Actually, your method is quite different from mine.

Again, the spheres are placed in the exact same way, random generation.
The methoid of checking intersections is also based on the naive approach,
but uses a different methodology.  Rather than checking each and every
sphere, only spheres that are within a certain distance are checked.
This is done by partitioning the space into square regions sized such that
only one sphere center may be present without intersection (D/sqrt(2)) and
build an array (what I like to call mailboxes).  Then as each position is
generated, first the mailbox that contains that position is checked, if it
is full, discard and generate a new point. If it is empty, check the
surrounding mailboxes that could contain a sphere that could intersect
(i.e, any mailbox within D of the position).  If any mailboxes are full,
check the distance between the two positions, if they intersect, discard
and move on. If no intersection is found, then the position is good and the
sphere is placed, the position value of this sphere is placed in the
corresponding mailbox and the mailbox is marked full.

Using this method avoids having to check each sphere.  Rather, each sphere
records its presence and position in the mailbox and each subsequent one
only needs to poll the necessary mailboxes for intersectiong spheres.
As more spheres get placed the method slows down as it gets harder to find a
random position that is good.  Eventually, when no positions exist, an
endless loop can be reached, however I haven't included any termination
criteria for this.
Partitions of any size can be used as long as they can only contain one
sphere (ie D/sqrt(2) or smaller) but I found that partition sizes that
equal to D/sqrt(2) work the fastest.  Larger partition sizes can be used,
but requires modification to the code to allow for multiple spheres within
the mailboxes, but I found it slows as the partition size gets bigger as
well.

Hope this is informative, if not boring.

code:
//start
camera{
  up y
  right x*image_width/image_height
  angle 60
  location <-8,30,1>
  look_at  0
}
light_source{<0,200,100> rgb 1 fade_power 2 fade_distance 200}

box{<-10,-1,-10> <10,0,10>
 texture{pigment{rgb 1} finish{ambient 0 diffuse 0.75}}
}



#local minPos=-10; #local maxPos= 10; //bounds
#local Rng=(maxPos-minPos);
#local D=0.23; //sphere radius
#local D2=D/sqrt(2); //partition size
#local nS=5000; //number of spheres
#local mn=ceil(Rng/D2); //number of partions needed
#local MBox=array[mn][mn]; //initialize mailboxes
#local R=seed(0);
#local xztomn=function(POS){ceil((POS-minPos)/D2)-1} //convert position to
mailbox number
#local D3=function(J,POSX)
{sqrt(pow(D,2)-pow(POSX-(J+select(J-xztomn(POSX),1,0))*D2-minPos,2))}
#local OBJ=sphere{0 D/2 texture{pigment{rgb <1,0,0>} finish{ambient 0
diffuse 0.75}}}

#local i=0; #while (i<nS)
  #local OK=false;
  #while (!OK)
    #local Pos1=<rand(R)*Rng+minPos,0,rand(R)*Rng+minPos>; //generate
position
    #ifndef(MBox[xztomn(Pos1.x)][xztomn(Pos1.z)]) #local OK=true; #end
//check mailbox. if defined then it is full, discard and generate new
position.
  #end
  #local j=max(xztomn(Pos1.x-D),0); #while (j<=min(xztomn(Pos1.x+D),mn-1) &
OK)//check surrounding mailboxes (only those within D of the position)
    #local DD3=D3(j,Pos1.x);
    #local k=max(xztomn(Pos1.z-DD3),0); #while
(k<=min(xztomn(Pos1.z+DD3),mn-1) & OK)
      #ifdef(MBox[j][k]) //if mailbox is defined, then it is full and
intersection needs to be checked
        #if (vlength(Pos1-MBox[j][k])<D) #local OK=false; #end //if
intersection found, break loop and discard
      #end
    #local k=k+1;#end
  #local j=j+1;#end
  #if (OK) //all mailboxes checked and no intersetion found
    #local MBox[xztomn(Pos1.x)][xztomn(Pos1.z)]=Pos1;//write position to
corresponding mailbox, this mailbox is now full
    object{OBJ translate Pos1} //place sphere
    #local i=i+1;
  #end
#end

//end

-tgq

Post a reply to this message

"Trevor G Quayle" <Tin### [at] hotmailcom> wrote:
> Warp <war### [at] tagpovrayorg> wrote:
> > I posted an article in povray.binaries.images about random placement
> > of non-intersecting spheres.
> > ... and so on..
>
> Actually, your method is quite different from mine.
>

Warp and Trevor:

Thank you both so much for posting these.  I learn quite a bit by following
your logic!

Dave Matthews

Post a reply to this message

Warp <war### [at] tagpovrayorg> wrote:
> I posted an article in povray.binaries.images about random placement
> of non-intersecting spheres.

Yes, it's an interesting little problem, Warp. I have C code (on my Amiga)
to place circles of semi-random radius in a plane. It generates a random
maximum radius for each circle, looks for an empty spot for the circle's
centre (by directly reading the screen buffer) and then determines the
actual radius for the circle to fit the space, which is obviously limited
by the nearest neighbour. The results are quite aesthetically pleasing,
IMNSHO.

>   There's a simple (at least in theory) way for this: Make a union of
> spheres of double diameter from all the existing coordinates and then
> check if the new center coordinate you calculated for the new sphere
> is inside() this union. That's it. POV-Ray uses internally quite advanced
> algorithms for speeding up this inside() test and it does it much faster
> than in linear time.
>
>   Of course there's a problem here: You just can't create the union of
> spheres each time you want to create a new sphere. That would actually
> be slower than the naive approach.

But how inefficient is it if you incrementally build your union, like this
snippet from one of my Penrose tiles scenes:

#debug "\nIgnore following warning, this union intentionally empty\n"
#declare Null = union{}
#debug "Ignore previous warning, this union intentionally empty\n"

#declare SetT = Null;
#declare SetP = array[2]{Null, Null};

#macro AddTri(A,B,C) #declare SetT = union{object{SetT} Tri(A,B,C)} #end
#macro AddTri0(A,B,C,P)
  #declare SetP[P] = union{object{SetP[P]}Tri0(A,B,C)}
#end

Post a reply to this message

This is quite fascinating! A "proximity detector" for placing random
objects. Thanks, guys, for going into such detail, both here and in
binaries:images. I see a lot of uses for the fundamental idea. The
challenge now (for me) is to adapt the concept to irregular surfaces like
HFs and isosurfaces. (Rather than use JC's  Repartition macro referenced
over in binaries:images, I think I'll try coming up with something of my
own--which is always more fun!) Right now, the optimization ideas aren't as
important to me as the basic concept.

Thanks again!

Ken W.

Post a reply to this message