Link to home
Start Free TrialLog in
Avatar of krusho
krushoFlag for United States of America

asked on

Closest Pair: Divide & Conquer

Ok, what I'm trying to do is take my function which is basically closest pair, and change it from the brute force method it currently is and turn it into a divide and conquer approach.  I'm having some issues understanding the descriptions of the problem I've found so far, so I was hoping to get some pseudo-code for c/c++ or possibly some code.  In the end this should get the minimum distance between any 2 points in an array of 3D points.  What I currently have works, but it's basically O( (N*(N+1))/2 ) and I've read with
divide and conquer it can be made into O( N*log(N) ).

I'm trying to write this in only C, so try to keep C++ stuff out for now.

typedef struct Vect // Vector struct for positions
{
      float x;
      float y;
      float z;
} Vect_t;

float findMinDistance_old(const int nPlayers, const Vect_t *playerArray)
{      // O( (N *(N + 1)) / 2) )
      Vect_t tmpVect;
      float tmpDist, minDist = 0.0f;
      int i, j;

      if(nPlayers <= 1) return 0; // error
      else if(nPlayers == 2)
      {
            // need this check because my loops won't work for 2 players
            tmpVect.x = playerArray[0].x - playerArray[1].x;
            tmpVect.y = playerArray[0].y - playerArray[1].y;
            tmpVect.z = playerArray[0].z - playerArray[1].z;

            tmpDist = (tmpVect.x * tmpVect.x)
                              + (tmpVect.y * tmpVect.y)
                              + (tmpVect.z * tmpVect.z);

            return (float)sqrt(tmpDist);
      }

      for(i = nPlayers; --i; )
            for(j = i; --j; )
            {
                  tmpVect.x = playerArray[i].x - playerArray[j].x;
                  tmpVect.y = playerArray[i].y - playerArray[j].y;
                  tmpVect.z = playerArray[i].z - playerArray[j].z;

                  tmpDist = (tmpVect.x * tmpVect.x)
                                    + (tmpVect.y * tmpVect.y)
                                    + (tmpVect.z * tmpVect.z);

                  if(tmpDist < minDist || minDist == 0) minDist = tmpDist;
            }

      return (float)sqrt(minDist);
}
Avatar of ozo
ozo
Flag of United States of America image

Avatar of krusho

ASKER

Yeah, I've seen that link as well as a couple others.  I guess the part that confuses me most about the description on that page is, how do I determine where l should be?  and how do I adapt that to a 3D point?
ASKER CERTIFIED SOLUTION
Avatar of Kashra
Kashra

Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of krusho

ASKER

Ok, so I've been working on this I got it somewhat working for 2D, but it seems to crash when I give it any more than 4 players.  It's getting a stack issue.  Do you have any tips on fixing this and making my code match a little more of what you are saying?

float findDistance(const Vect_t *v1, const Vect_t *v2)
{
      Vect_t tmpVect;
      float tmpDist;

      tmpVect.x = v1->x - v2->x;
      tmpVect.y = v1->y - v2->y;
      tmpVect.z = v1->z - v2->z;

      tmpDist = (tmpVect.x * tmpVect.x)
                        + (tmpVect.y * tmpVect.y)
                        + (tmpVect.z * tmpVect.z);

      return tmpDist;
}

float findMinDistance(const int nPlayers, const Vect_t *playerArray)
{
      Vect_t tmpVect;
      int i, j;
      float tmpDist;
      float left, right, middle;
      //      Draw dividing plane (using longest dimension or whatever approach you wish)
      float divide = (playerArray[(nPlayers>>1)+1].x - playerArray[(nPlayers>>1)].x);

      switch(nPlayers)
      {
            case 0:
            case 1:
                  return 0;

            case 2:
                  tmpVect.x = playerArray[0].x - playerArray[1].x;
                  tmpVect.y = playerArray[0].y - playerArray[1].y;
                  tmpVect.z = playerArray[0].z - playerArray[1].z;

                  tmpDist = (tmpVect.x * tmpVect.x)
                                    + (tmpVect.y * tmpVect.y)
                                    + (tmpVect.z * tmpVect.z);

                  return (float)sqrt(tmpDist);

            default:                        
                  //      Call DivideAndConquer twice to retrieve shortest distance between points in each bin
                  left = DivideAndConquer(0, (nPlayers >> 1), playerArray);
                  right = DivideAndConquer((nPlayers >> 1), nPlayers, playerArray);

                  //      Determine which of the two distances is shorter
                  if(left < right)
                        tmpDist = left;
                  else
                        tmpDist = right;

                  //      Search for any points that are within that distance of the dividng plane and compare distances between points that are from separate bins.
                  for(i = 0; i < nPlayers; i++)
                  {
                        if(playerArray[i].x < divide)
                        {
                              for(j = 0; j < (nPlayers>>1); j++)
                              {
                                    middle = findDistance(&playerArray[i], &playerArray[j]);
                                    if(tmpDist > middle && middle != 0)
                                          tmpDist = middle;
                              }
                        }
                        else
                        {
                              for(j = (nPlayers>>1); j < nPlayers; j++)
                              {
                                    middle = findDistance(&playerArray[i], &playerArray[j]);
                                    if(tmpDist > middle && middle != 0)
                                          tmpDist = middle;
                              }
                        }
                  }

                  return (float)sqrt(tmpDist);
                  break;
      }
}

float DivideAndConquer(const int start, const int end, const Vect_t *playerArray)
{
      Vect_t tmpVect;
      float tmpDist, minDist = 0.0f;
      float cut, left, right, dist, dist2;
      int i, j, nPlayers = (end - start);

      switch(nPlayers)
      {
            case 0:
            case 1:
                  //IF there is only one point, return 0
                  return 0;

            case 2:
                  //IF there are only two points, measure distance between them
                  tmpVect.x = playerArray[start].x - playerArray[end].x;
                  tmpVect.y = playerArray[start].y - playerArray[end].y;
                  tmpVect.z = playerArray[start].z - playerArray[end].z;

                  tmpDist = (tmpVect.x * tmpVect.x)
                                    + (tmpVect.y * tmpVect.y)
                                    + (tmpVect.z * tmpVect.z);

                  return tmpDist;

            default:
                  //IF there are more than two points:

                  //      Draw dividing plane (using longest dimension or whatever approach you wish)
                  cut = (playerArray[start + (nPlayers>>1) + 1].x - playerArray[start + (nPlayers>>1)].x);

                  //      Call DivideAndConquer twice to retrieve shortest distance between points in each bin
                  left = DivideAndConquer(start, (nPlayers>>1), playerArray);
                  right = DivideAndConquer((nPlayers>>1), end, playerArray);

                  //      Determine which of the two distances is shorter
                  if(left < right)
                        dist = left;
                  else
                        dist = right;

                  //      Search for any points that are within that distance of the dividng plane and compare distances between points that are from separate bins.
                  for(i = start; i < end; i++)
                  {
                        if(playerArray[i].x < cut)
                        {
                              for(j = start; j < (nPlayers>>1); j++)
                              {
                                    dist2 = findDistance(&playerArray[i], &playerArray[j]);
                                    if(dist > dist2 && dist2 != 0)
                                          dist = dist2;
                              }
                        }
                        else
                        {
                              for(j = (nPlayers>>1); j < end; j++)
                              {
                                    dist2 = findDistance(&playerArray[i], &playerArray[j]);
                                    if(dist > dist2 && dist2 != 0)
                                          dist = dist2;
                              }
                        }
                  }

                  return dist;
                  break;
      }
}
SOLUTION
Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of krusho

ASKER

Ok, I handled 1 & 2 with this code.  3 is fine because as you can see I sort it in the x myself.  I still am not sure if I am doing your 4th point the best way, if you could help with that.  And also, it seems when I run my function with more than 4 players it causes it to crash because somehow it keeps getting a negative number of players.  it has something to do with the (nPlayers>>1) that I use for end in my function, any ideas?

float findDistance(const Vect_t *v1, const Vect_t *v2)
{
      return (((v1->x - v2->x) * (v1->x - v2->x)) + ((v1->y - v2->y) * (v1->y - v2->y)) + ((v1->z - v2->z) * (v1->z - v2->z)));
}

float findMinDistance(const int nPlayers, const Vect_t *playerArray)
{
      Vect_t* playerArray_sorted = malloc(sizeof(Vect_t) * nPlayers);

      quickSort_x(playerArray, nPlayers, playerArray_sorted);
      
      return (float)sqrt(DivideAndConquer(playerArray_sorted, 0, nPlayers));
}

float DivideAndConquer(const Vect_t *playerArray, const int start, const int end)
{
      Vect_t cut;
      float left, right, dist;
      int i, j;
      int nPlayers = (end - start);

      switch(nPlayers)
      {
            case 0:
            case 1:
                  //IF there is only one point, return 0
                  return 0;

            case 2:
                  //IF there are only two points, measure distance between them
                  return findDistance(&playerArray[start], &playerArray[(end - 1)]);

            default:
                  //IF there are more than two points:

                  //      Draw dividing plane (using longest dimension or whatever approach you wish)
                  cut.x = playerArray[(nPlayers>>1)].x;
                  cut.y = 0;
                  cut.z = 0;

                  //      Call DivideAndConquer twice to retrieve shortest distance between points in each bin
                  left = DivideAndConquer(playerArray, start, (nPlayers>>1));
                  right = DivideAndConquer(playerArray, (nPlayers>>1), end);

                  //      Determine which of the two distances is shorter
                  if(left < right && left != 0)
                        dist = left;
                  else
                        dist = right;

                  //      Search for any points that are within that distance of the dividng plane and compare distances between points that are from separate bins.
                  for(i = start; i < end; i++)
                  {
                        cut.y = playerArray[i].y;
                        cut.z = 0; // temp

                        if( (playerArray[i].x <= cut.x) && (findDistance(&playerArray[i], &cut) < dist) )
                        {
                              for(j = start; j < end; j++)
                              {
                                    if( (playerArray[j].x >= cut.x) && (findDistance(&playerArray[j], &cut) < dist) )
                                    {
                                          left = findDistance(&playerArray[i], &playerArray[j]);
                                          if(left < dist && left != 0)
                                                dist = left;
                                    }
                              }
                        }
                        else// if( (playerArray[i].x >= cut.x) && (findDistance(&playerArray[i], &cut) < dist) )
                        {
                              for(j = start; j < end; j++)
                              {
                                    if( (playerArray[j].x <= cut.x) && (findDistance(&playerArray[j], &cut) < dist) )
                                    {
                                          right = findDistance(&playerArray[i], &playerArray[j]);
                                          if(right < dist && right != 0)
                                                dist = right;
                                    }
                              }
                        }
                  }

                  return dist;
                  break;
      }
}
SOLUTION
Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of Kashra
Kashra

Ugh, I was trying to be quick and dirty with the math, but as a side note, comparing the squared values doesn't actually work unless the actual values are all positive. Its still possible to do this "correctly", but you might have to take a square root of the dist once per function instead of squaring the x values.
Avatar of krusho

ASKER

I'm still working on this I just got busy with a couple other things.  I won't forget to give out my points when it's working.  I'll post some new code later this week.
Avatar of krusho

ASKER

Ok, I fixed the crash and it correctly finds the minimum distance.  The question is this the O( N log(N) ) ?
Avatar of krusho

ASKER

forgot to post the code

float findDistance(const Vect_t *v1, const Vect_t *v2)
{
      return (((v1->x - v2->x) * (v1->x - v2->x)) + ((v1->y - v2->y) * (v1->y - v2->y)) + ((v1->z - v2->z) * (v1->z - v2->z)));
}

float findMinDistance(const int nPlayers, const Vect_t *playerArray)
{
      Vect_t* playerArray_sorted = malloc(sizeof(Vect_t) * nPlayers);

      quickSort_x(playerArray, nPlayers, playerArray_sorted);
      
      return (float)sqrt(DivideAndConquer(playerArray_sorted, 0, nPlayers));
}

float DivideAndConquer(const Vect_t *playerArray, const int start, const int end)
{
      Vect_t cut;
      float left = 0, right = 0, dist = 0;
      int i, j;
      int nPlayers;

      if(start > end)
            nPlayers = 0;
      else
            nPlayers = (end - start);

      switch(nPlayers)
      {
            case 0:
            case 1:
                  //IF there is only one point, return 0
                  return 0;

            case 2:
                  //IF there are only two points, measure distance between them
                  return findDistance(&playerArray[start], &playerArray[(end - 1)]);

            default:
                  //IF there are more than two points:

                  //      Draw dividing plane (using longest dimension or whatever approach you wish)
                  cut.x = playerArray[(start + (nPlayers / 2))].x;
                  cut.y = 0;
                  cut.z = 0;

                  //      Call DivideAndConquer twice to retrieve shortest distance between points in each bin
                  left = DivideAndConquer(playerArray, start, (start + (nPlayers / 2)));
                  right = DivideAndConquer(playerArray, (start + (nPlayers / 2)), end);

                  //      Determine which of the two distances is shorter
                  if(left < right && left != 0)
                        dist = left;
                  else
                        dist = right;

                  //      Search for any points that are within that distance of the dividng plane and compare distances between points that are from separate bins.
                  for(i = start; i < end; i++)
                  {
                        cut.y = playerArray[i].y;
                        cut.z = 0; // temp

                        if( (playerArray[i].x <= cut.x) && (findDistance(&playerArray[i], &cut) < dist) )
                        {
                              for(j = start; j < end; j++)
                              {
                                    if( (playerArray[j].x >= cut.x) && (findDistance(&playerArray[j], &cut) < dist) )
                                    {
                                          left = findDistance(&playerArray[i], &playerArray[j]);
                                          if(left < dist && left != 0)
                                                dist = left;
                                    }
                              }
                        }
                        else if( (playerArray[i].x >= cut.x) && (findDistance(&playerArray[i], &cut) < dist) )
                        {
                              for(j = start; j < end; j++)
                              {
                                    if( (playerArray[j].x <= cut.x) && (findDistance(&playerArray[j], &cut) < dist) )
                                    {
                                          right = findDistance(&playerArray[i], &playerArray[j]);
                                          if(right < dist && right != 0)
                                                dist = right;
                                    }
                              }
                        }
                  }

                  return dist;
                  break;
      }
}
SOLUTION
Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of krusho

ASKER

Appreciate all the help, I'm guessing now it should be in O( N log( N ) ).  What would be the steps to convert this to handle 3D points, I know you were talking about Axis-Aligned Bounding Boxes (AABB), but not exactly sure the steps to getting this working.

                  for(i = start; i < end; i++)
                  {
                        cut.y = playerArray[i].y;
                        cut.z = 0; // temp

                        idist = (cut.x - playerArray[i].x);
                        if( (playerArray[i].x <= cut.x) && ((idist * idist) < dist) )
                        {
                              for(j = (i + 1); j < end; j++)
                              {
                                    idist = (playerArray[j].x - cut.x);
                                    if( (playerArray[j].x >= cut.x) && ((idist * idist) < dist) )
                                    {
                                          left = findDistance(&playerArray[i], &playerArray[j]);
                                          if(left < dist && left != 0)
                                                dist = left;
                                    }
                              }
                        }
                  }
By my understanding, these are 3D points. You have x, y, and z. Right now, it should work for 3D points.

The only thing to do now is to get "smarter" about where you decide to draw your dividing plane. As it stands, you're just arbitrarily choosing the X axis and drawing a plane through the center point. You could be a little smarter by just finding which axis is the longest (X, Y or Z) and using that one, instead.

Or you could be even more advanced and start drawing AABB around your points. That's a level of complication I'd rather not get into, here. Its been explained elsewhere and its not necessarily going to help your speed. Its a question of trading one set of calculations for another. If your points are strangely distributed, drawing AABB might help reduce the number of recursions you have to take, at the expense of doing all the math to figure out how to draw the AABB. If your points aren't strangely shaped, then drawing AABB's probably won't make much of a difference.
Avatar of krusho

ASKER

Oh you're right, I was just testing it in 2D because I wasn't sure if it would work in 3D yet.  Appreciate the help.