Implementing the KD Tree

[ LiB ]

In order to implement the KD Tree algorithm, you need a structure that will be used for the nodes of the tree. In this section, you'll write a few methods that will be located in the PhotonMap class for subdivision and searching of the KD Tree. The KD Tree algorithm doesn't require any methods or variables from the cScene class, so it will be defined only in the PhotonMap class. The only element that the KD Tree requires is the list of photons in the photon map array.

The KD Tree is composed of a set of smaller elements called KD nodes . A KD node is a component of the KD Tree that dictates the type of splitting that occurs. You will have a KD node for each photon in the photon map array to give each photon a splitting dimension and sorting value. The KD nodes are used to construct the hierarchy of point sets by partitioning photons with each other and flagging each KD node of the type of dimensional splitting that occurs. This in essence preserves the splitting information or the hierarchy of the point list so that you can search through the KD nodes based on their assigned sorting values.

Because you already have a method for locating the nearest photons found using the photon map array, you need a way to toggle between using the KD Tree and the linear search method of retrieving the closest photons. This will be a simple pre-compiled macro located in the PhotonMap class. If this macro is defined, your code can use the KD Tree when building and searching for photons.

Adding the KD Node

The KD node structure is the most important element of the KD Treeit contains a photon and two scalar values. The iSort value is used to build the photon map as well as to determine which kind of splitting should be conducted . The dValue holds the dimensional value for the partitioning plane.

For example, a photon located in space (0,0,10) that is used as the splitting node would be listed in the structure as follows :

You have a splitting node at location 0,0, and 10 in space for the pho-ton position.
If you decipher the point, you know that a partition should be conducted positive 10 units on the Z axis, indicating that the dValue is set to 10 units into the Z axis.
The iSort or splitting dimension is set to the Z plane in world space.

The node structure is defined as follows:

 typedef struct {    photon_t tPhoton;    int iSort;    float dValue; } KDNode;

Adding the KD Tree

Typically, you add the KD Tree algorithm to photon mapping by using a macro to switch the logic of the code to the KD Tree's methods and variables. The KD node array is listed in the PhotonMap class. The pre-compiled macro is listed as follows:

 #define   USE_KD_TREE

The KD Tree is listed in the PhotonMap class as follows:

 class  PhotonMap { // rest of code assumed KDNode   * pKDNodes; }

Building the KD Tree

You'll now write a method that will take the photon map array and create a local photon array to mangle photons. A new KD node is created for each photon listed in the photon map array. A variable, called sortKey , is defined for partitioning photons on a specific dimension. If the sort key is 1, the split was conducted on the X dimension. If it's 2, it's for the Y dimension, whereas 3 is for the Z dimension. Each node of the KD Tree is unsorted in the beginning. The function ends by calling MakeSubTree() , which begins the process of building and spacing the partitioning photons.

 photon_t *tLocalPhotons; long lLocalPhotonCount; int sortKey; void PhotonMap::MakeKDTree() {    long i;   // allocate memory to hold photons in the photon map array    lLocalPhotonCount = lNumGlobalPhotons;    tLocalPhotons = (photon_t *) malloc(sizeof(photon_t)           * lNumGlobalPhotons);    // copy elements over to local array    for (i = 0; i < lNumGlobalPhotons; i++)       tLocalPhotons[i] = pGlobalPhoton[i];    // allocate memory for each KD node for    // each photon stored in photon map array    pKDNodes = (KDNode *) malloc(sizeof(KDNode) * lNumGlobalPhotons);    // set each KD node to unsorted    for (i = 0; i < lNumGlobalPhotons;  i++)    {       pKDNodes[i].iSort = -1;    }   // build tree from available photons    MakeSubTree(0,lNumGlobalPhotons);

Using the SortProc Method

If you recall, when building the KD Tree, points need to be relative to one another. Sometimes during the build, the elements of the photon-mapping array need to be placed in an orderly fashion. In order to split a set of points on a desired plane, you need to sort the points based on their positions in the array. For example, if the sort key was the Y dimension, all points in the local array must be in Y distance order, from closest to farthest. This example uses the quick sort method for reordering photons in the local array. Using quick sort is an easy way of sorting. All you need to do is define a recursive function for the algorithm to call. The recursive function will be called constantly during the re-ordering process. After quick sort completes, the list of photons are in their appropriate order based on a specific dimension.

 int SortProc( const void *a, const void *b ) {    float r1, r2;    photon_t pa = *(photon_t *)a;    photon_t pb = *(photon_t *)b;    // if the sort key was on the X Axis    if (sortKey == 0)    {       r1 = pa.tPosition.x;       r2 = pb.tPosition.x;    }    // if the sort key was on the Y Axis    else if (sortKey == 1)    {       r1 = pa.tPosition.y;       r2 = pb.tPosition.y;    }    // if the sort key was on the Z Axis    else if (sortKey == 2)    {       r1 = pa.tPosition.z;       r2 = pb.tPosition.z;    }    else       return 0;    if( r1 < r2 )       return -1;    else if( r1 > r2 )        return 1;    else       return 0; }

Using the MakeSubTree() Method

The MakeSubTree() method's main purpose is to recursively subdivide the input point list. Because this is a KD Tree in a 3D environment, you use cubes (or bounding boxes) when partitioning photons. Each bounding box is sliced through the center divided by X, Y, and Z planes. If the photon map array has only one photon, one KD node is returned unsorted and unchanged. The first time the algorithm begins, a few bogus values are used to set up the maximum distance for the large bounding box encompassing the scene. This is conducted so as to accommodate the whole scene without leaving out any photons. The new correct bounding box is then computed by the list of photons in the photon map array. The farthest photons on each side of the box are used to create the boundaries of the box.

The process then cycles through each photon in the array and decides how to conduct the specific type of split. The job of the algorithm is to yield a dimension to split on for every photon. Subtracting the maximum (in terms of X, Y, and Z) from the minimum (in terms of X, Y, and Z) determines a range; the largest range of the three dictates the dimension to subdivide the remaining photons to split on. The current point is used as the median split where two new child lists will be spawned. The algorithm calls the recursive function again for each new child list. Every time a splitting occurs, two lists (photons to the left and photons to the right) are generated. The whole process begins once again with the child as the parent node.

 void PhotonMap::MakeSubTree(long start, long n) { long left, nleft, right, nright; long maxRangeIndex; double maxRange; long j; photon_t p; double rmin_x, rmin_y, rmin_z; double rmax_x, rmax_y, rmax_z;    // only one photon if (n == 1) {    pKDNodes[start].tPhoton = tLocalPhotons[start];    //pKDNodes[start].dValue = 0;    //pKDNodes[start].iSort = 0;    return; } // Initialize tree stuff - left   = start; nleft  = n / 2; right  = left + nleft; nright = n - nleft; maxRangeIndex = -1; maxRange      = -100000000000000; // set some bogus ranges so that we don't leave out any photons // PS: The objects in the scene cannot // possibly fall outside this big bounding box rmin_x  = 100000000000000; rmax_x  = -100000000000000; rmin_y  = 100000000000000; rmax_y  = -100000000000000; rmin_z  = 100000000000000; rmax_z  = -100000000000000; // get bounding box sides from available list of photons for(j = start; j < start+n; j++ ) {    p = tLocalPhotons[j];    if( p.tPosition.x < rmin_x )       rmin_x = p.tPosition.x;    else if( p.tPosition.x > rmax_x )       rmax_x = p.tPosition.x;    if( p.tPosition.y < rmin_y )       rmin_y = p.tPosition.y;    else if( p.tPosition.y > rmax_y )       rmax_y = p.tPosition.y;    if( p.tPosition.z < rmin_z )        rmin_z = p.tPosition.z;    else if( p.tPosition.z > rmax_z )       rmax_z = p.tPosition.z;  } // select the axis that has the maximum range // of the farthest distance of available photons  if (rmax_x-rmin_x > maxRange)  {     maxRangeIndex = 0;     maxRange = rmax_x-rmin_x;  }  if (rmax_y-rmin_y > maxRange)  {      maxRangeIndex = 1;      maxRange = rmax_y-rmin_y;  }  if (rmax_z-rmin_z > maxRange)  {     maxRangeIndex = 2;     maxRange = rmax_z-rmin_z;  } // the sort key is set to the dimensional split  sortKey   = maxRangeIndex; // sort available photons in list to selected dimension  qsort( (void *)&(tLocalPhotons[start]), n, sizeof(photon_t), SortProc );   // Median split   pKDNodes[right].tPhoton = tLocalPhotons[right];  pKDNodes[right].iSort   = sortKey;    // what kind of dimensional split? X, Y, or Z?    if (sortKey == 0)       pKDNodes[right].dValue = tLocalPhotons[right].tPosition.x;    else if (sortKey == 1)       pKDNodes[right].dValue = tLocalPhotons[right].tPosition.y;    else       pKDNodes[right].dValue = tLocalPhotons[right].tPosition.z;    // the children to the left of me    MakeSubTree( left, nleft );    // the children to the right of me    MakeSubTree( right, nright ); }

Checking the Photons

As you run through different nodes in the tree and create a photon pool, you need a method to track how many photons are accumulated in the temporary distance array. This requires a few global variables for retrieving the relevant photons because the distance array is used to hold relevant photons found by comparing the distance between found photons and the point of intersection. A counter is used to tally the number of photons found. By testing new photons with the point of intersection, the program can then decide whether to add the photon to the distance array or not.

The function searches through the photon map, finds the closest number of photons (based on their distance to the intersection point), and saves them in the distance array. The function finally returns with the number of pho-tons found as well as the distance array filled with the closest photons. It also lists the photons in distance order. For example, the nearest photons to the point of intersection are listed first in the distance array and farther photons are listed towards the end of the array. The farthest photon is listed as the last element in the array; it creates the expanding sphere. The radius of the expanding sphere is the distance from the sphere's center to the last photon in the array. This photon is used to calculate the area of the bounding sphere.

 photon_t    *tLocalPhotons; long        lLocalPhotonCount; cVector3    tRefPoint; float       maxDistance; long        lFoundCount; photonlist_t  *tLocalList; long        tLocalListSize;     void PhotonMap::CheckPhoton(photon_t p)    {       float dCurrDist;       long j;       cVector3 Dummy;       Dummy     = p.tPosition - tRefPoint;       dCurrDist = Dummy.Mag();       if (dCurrDist < maxDistance)       {          for (j = ((lFoundCount < tLocalListSize-1) ?             lFoundCount : tLocalListSize-1); (j > 0)          && (tLocalList[ j - 1].dDistance > dCurrDist); j)             tLocalList[j] = tLocalList[j-1];          tLocalList[j].dDistance = dCurrDist;          tLocalList[j].tPhoton = p;      if (lFoundCount != tLocalListSize)         lFoundCount ++;      if (lFoundCount == tLocalListSize)         maxDistance = tLocalList[tLocalListSize-1].dDistance;    } }

Using the SearchTree Function

The SearchTree function is recursively called until n photons are found in the KD Tree for the point of intersection within a given volume. You begin searching from the root of the tree and then through each node respectively. As the algorithm runs through a series of nodes in a specific area, a local photon pool is allocated with photons that are within a certain distance from the hit point. When a photon is found, it is tested to determine what kind of split was conducted in the build process. The intersection point is constantly tested with each splitting node and is either branched down into the left or right subtrees. The global distance array for the closest neighbors is constantly re-ordered with photons in distance order. This is done so that, as the algorithm runs through a series of nodes in the KD Tree, it would be easy to replace the farthest photons with newly found closer photons.

 void PhotonMap::SearchTree( int start, int n ) {    int left, nleft, right, nright;    int sk;    photon_t ph;    float d, d2;    // Check terminator     if( n <= 2)    {      CheckPhoton(pKDNodes[start].tPhoton);      if (n == 2)         CheckPhoton(pKDNodes[start+1].tPhoton);      return;    }    // Initialize tree stuff -    left = start;    nleft = n / 2;    right = left + nleft;    nright = n - nleft;    //- Do search     sk = pKDNodes[right].iSort;    ph = pKDNodes[right].tPhoton;    if (sk == 0)       d = tRefPoint.x - pKDNodes[right].dValue;    else if (sk == 1)       d = tRefPoint.y - pKDNodes[right].dValue;    else if (sk == 2)       d = tRefPoint.z - pKDNodes[right].dValue;    else    {       CheckPhoton(pKDNodes[right].tPhoton);       return;    }    d2 = fabsf(d);    if( d <= 0.0 )    {       SearchTree( left, nleft);       if( d2 <= maxDistance )          SearchTree( right, nright );     }    else    {       SearchTree( right, nright );       if( d2 <= maxDistance )          SearchTree( left, nleft );    } }

Using the GetNearestPhotonsKD Method

This method replaces the previous method but also utilizes the KD Tree. The searching variables are initialized for optimal performance. The reference point is set and the found count is set to zero. The local distance array is set to the local return array passed as a parameter. The maximum photons to return are set and the maximum search distance is set to a very large number so that it includes all photons in the scene. This value will, of course, be updated when the algorithm begins. The SearchTree() method is called and begins searching through the relevant nodes of the KD Tree. When the function is done calling SearchTree() , the nearest count is returned and the passed array is updated with the closest photons found.

 long PhotonMap::GetNearestPhotonsKD(photonlist_t tList[],       long lLength, cVector3 tInter) {    tRefPoint        = tInter;    lFoundCount      = 0;    tLocalList       = tList;    tLocalListSize   = lLength;    maxDistance      = 100000000;    SearchTree(0,lLocalPhotonCount);    return lFoundCount; }

Updating the Photon- Contribution Method

You can now update the photon-contribution method so that it encapsulates the KD Tree algorithm. Based on the preprocessor macro, the KD Tree is utilized or discarded. The new code is listed in bold and is fairly straightforward. The change to the method is minimal with only a few lines of new code.

 color3 cScene::PhotonContribution(        int iObjectIndex,cVector3 tDirection, cVector3 tInter) {     // allocate memory     photonlist_t *pList = new photonlist_t[lMaxPhotonNeighbors];     long lNumPhotons;     long i;     color3 tColor, tTmp;     float dDot;     float dPiR2;     tColor = color3::Black; #ifdef USE__KD__TREE     lNumPhotons = GetNearestPhotonsKD((pList,lMaxPhotonNeighbors,tInter); #else     lNumPhotons = GetNearestPhotons(pList,lMaxPhotonNeighbors,tInter); #endif     // no photons found, omg very bad     if (lNumPhotons == 0)        return tColor;     // for the amount of friends found     for (i = 0; i < lNumPhotons; i++)     {         // is this friend facing the point of intersection         dDot = pObjectList[ iObjectIndex ].Normal(tInter) *                  pList[i].tPhoton.tDirection;         // if not, let's go to next array node         if (dDot <= 0.0)         continue;         // add the power of the photon and the cosine shading method         tTmp = pList[i].tPhoton.tPower *          (float) dDot * pObjectList[iObjectIndex ].fDiffuseFactor;         // sum up everything         tColor.r = tColor.r + tTmp.r;         tColor.g = tColor.g + tTmp.g;         tColor.b = tColor.b + tTmp.b;      }      // divide the summed photons by the area of the circle found     dPiR2 = 3.1415926f *            pList[lNumPhotons-1].dDistance *            pList[lNumPhotons-1].dDistance;     tColor.r = tColor.r / (float) dPiR2;     tColor.g = tColor.g / (float) dPiR2;     tColor.b = tColor.b / (float) dPiR2;     // used to color the photon with the object     tColor.r *= pObjectList[ iObjectIndex ].Color.r;     tColor.g *= pObjectList[ iObjectIndex ].Color.g;     tColor.b *= pObjectList[ iObjectIndex ].Color.b;     // saturate colors out of unit range     tColor.Sat();     // free memory if(pList) delete [] pList;    // return final solution     return tColor; }

[ LiB ]