Tapkee
quadtree.hpp
Go to the documentation of this file.
00001 
00032 #include <math.h>
00033 #include <float.h>
00034 #include <stdlib.h>
00035 #include <stdio.h>
00036 #include <algorithm>
00037 
00038 #ifndef QUADTREE_H
00039 #define QUADTREE_H
00040 
00041 namespace tsne
00042 {
00043 
00044 class Cell {
00045 
00046 public:
00047 
00048     double x;
00049     double y;
00050     double hw;
00051     double hh;
00052 
00053     bool containsPoint(double point[])
00054     {
00055         if(x - hw > point[0]) return false;
00056         if(x + hw < point[0]) return false;
00057         if(y - hh > point[1]) return false;
00058         if(y + hh < point[1]) return false;
00059         return true;
00060     }
00061 
00062 };
00063 
00064 
00065 class QuadTree
00066 {
00067 
00068     // Fixed constants    
00069     static const int QT_NO_DIMS = 2;
00070     static const int QT_NODE_CAPACITY = 1;
00071 
00072     // A buffer we use when doing force computations
00073     double buff[QT_NO_DIMS];
00074 
00075     // Properties of this node in the tree
00076     QuadTree* parent;
00077     bool is_leaf;
00078     int size;
00079     int cum_size;
00080 
00081     // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
00082     Cell boundary;
00083 
00084     // Indices in this quad tree node, corresponding center-of-mass, and list of all children
00085     double* data;
00086     double center_of_mass[QT_NO_DIMS];
00087     int index[QT_NODE_CAPACITY];
00088 
00089     // Children
00090     QuadTree* northWest;
00091     QuadTree* northEast;
00092     QuadTree* southWest;
00093     QuadTree* southEast;
00094 
00095 public:
00096 
00097     // Default constructor for quadtree -- build tree, too!
00098     QuadTree(double* inp_data, int N) : 
00099         parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
00100         northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
00101     {
00102         // Compute mean, width, and height of current map (boundaries of quadtree)
00103         double* mean_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] = .0;
00104         double*  min_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++)  min_Y[d] =  DBL_MAX;
00105         double*  max_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++)  max_Y[d] = -DBL_MAX;
00106         for(int n = 0; n < N; n++) {
00107             for(int d = 0; d < QT_NO_DIMS; d++) {
00108                 mean_Y[d] += inp_data[n * QT_NO_DIMS + d];
00109                 if(inp_data[n * QT_NO_DIMS + d] < min_Y[d]) min_Y[d] = inp_data[n * QT_NO_DIMS + d];
00110                 if(inp_data[n * QT_NO_DIMS + d] > max_Y[d]) max_Y[d] = inp_data[n * QT_NO_DIMS + d];
00111             }
00112         }
00113         for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] /= (double) N;
00114         
00115         // Construct quadtree
00116         init(NULL, inp_data, mean_Y[0], mean_Y[1], std::max(max_Y[0] - mean_Y[0], mean_Y[0] - min_Y[0]) + 1e-5,
00117                                                    std::max(max_Y[1] - mean_Y[1], mean_Y[1] - min_Y[1]) + 1e-5);
00118         fill(N);
00119         delete[] mean_Y; delete[] max_Y; delete[] min_Y;
00120     }
00121 
00122     // Constructor for quadtree with particular size and parent -- build the tree, too!
00123     QuadTree(double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
00124         parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
00125         northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
00126     {
00127         init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
00128     }
00129 
00130     // Constructor for quadtree with particular size and parent -- build the tree, too!
00131     QuadTree(double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
00132         parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
00133         northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
00134     {
00135         init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
00136         fill(N);
00137     }
00138 
00139     // Constructor for quadtree with particular size (do not fill the tree)
00140     QuadTree(QuadTree* inp_parent, double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
00141         parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
00142         northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
00143     {
00144         init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
00145         fill(N);
00146     }
00147 
00148     // Constructor for quadtree with particular size and parent (do not fill the tree)
00149     QuadTree(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
00150         parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
00151         northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
00152     {
00153         init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
00154     }
00155 
00156     // Destructor for quadtree
00157     ~QuadTree()
00158     {
00159         delete northWest;
00160         delete northEast;
00161         delete southWest;
00162         delete southEast;
00163     }
00164 
00165     void setData(double* inp_data)
00166     {
00167         data = inp_data;
00168     }
00169 
00170     QuadTree* getParent()
00171     {
00172         return parent;
00173     }
00174 
00175     //void construct(Cell boundary);
00176 
00177     // Insert a point into the QuadTree
00178     bool insert(int new_index)
00179     {
00180         // Ignore objects which do not belong in this quad tree
00181         double* point = data + new_index * QT_NO_DIMS;
00182         if(!boundary.containsPoint(point))
00183             return false;
00184         
00185         // Online update of cumulative size and center-of-mass
00186         cum_size++;
00187         double mult1 = (double) (cum_size - 1) / (double) cum_size;
00188         double mult2 = 1.0 / (double) cum_size;
00189         for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] *= mult1;
00190         for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] += mult2 * point[d];
00191         
00192         // If there is space in this quad tree and it is a leaf, add the object here
00193         if(is_leaf && size < QT_NODE_CAPACITY) {
00194             index[size] = new_index;
00195             size++;
00196             return true;
00197         }
00198         
00199         // Don't add duplicates for now (this is not very nice)
00200         bool any_duplicate = false;
00201         for(int n = 0; n < size; n++) {
00202             bool duplicate = true;
00203             for(int d = 0; d < QT_NO_DIMS; d++) {
00204                 if(point[d] != data[index[n] * QT_NO_DIMS + d]) { duplicate = false; break; }
00205             }
00206             any_duplicate = any_duplicate | duplicate;
00207         }
00208         if(any_duplicate) return true;
00209         
00210         // Otherwise, we need to subdivide the current cell
00211         if(is_leaf) subdivide();
00212         
00213         // Find out where the point can be inserted
00214         if(northWest->insert(new_index)) return true;
00215         if(northEast->insert(new_index)) return true;
00216         if(southWest->insert(new_index)) return true;
00217         if(southEast->insert(new_index)) return true;
00218         
00219         // Otherwise, the point cannot be inserted (this should never happen)
00220         return false;
00221     }
00222 
00223     // Create four children which fully divide this cell into four quads of equal area
00224     void subdivide()
00225     { 
00226         // Create four children
00227         northWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
00228         northEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
00229         southWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
00230         southEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
00231         
00232         // Move existing points to correct children
00233         for(int i = 0; i < size; i++) {
00234             bool success = false;
00235             if(!success) success = northWest->insert(index[i]);
00236             if(!success) success = northEast->insert(index[i]);
00237             if(!success) success = southWest->insert(index[i]);
00238             if(!success) success = southEast->insert(index[i]);
00239             index[i] = -1;
00240         }
00241         
00242         // Empty parent node
00243         size = 0;
00244         is_leaf = false;
00245     }
00246 
00247     // Checks whether the specified tree is correct
00248     bool isCorrect()
00249     {
00250         for(int n = 0; n < size; n++) {
00251             double* point = data + index[n] * QT_NO_DIMS;
00252             if(!boundary.containsPoint(point)) return false;
00253         }
00254         if(!is_leaf) return northWest->isCorrect() &&
00255                             northEast->isCorrect() &&
00256                             southWest->isCorrect() &&
00257                             southEast->isCorrect();
00258         else return true;
00259     }
00260 
00261     // Rebuilds a possibly incorrect tree (LAURENS: This function is not tested yet!)
00262     void rebuildTree()
00263     {
00264         for(int n = 0; n < size; n++) {
00265             // Check whether point is erroneous
00266             double* point = data + index[n] * QT_NO_DIMS;
00267             if(!boundary.containsPoint(point)) {
00268                 
00269                 // Remove erroneous point
00270                 int rem_index = index[n];
00271                 for(int m = n + 1; m < size; m++) index[m - 1] = index[m];
00272                 index[size - 1] = -1;
00273                 size--;
00274                 
00275                 // Update center-of-mass and counter in all parents
00276                 bool done = false;
00277                 QuadTree* node = this;
00278                 while(!done) {
00279                     for(int d = 0; d < QT_NO_DIMS; d++) {
00280                         node->center_of_mass[d] = ((double) node->cum_size * node->center_of_mass[d] - point[d]) / (double) (node->cum_size - 1);
00281                     }
00282                     node->cum_size--;
00283                     if(node->getParent() == NULL) done = true;
00284                     else node = node->getParent();
00285                 }
00286                 
00287                 // Reinsert point in the root tree
00288                 node->insert(rem_index);
00289             }
00290         }    
00291         
00292         // Rebuild lower parts of the tree
00293         northWest->rebuildTree();
00294         northEast->rebuildTree();
00295         southWest->rebuildTree();
00296         southEast->rebuildTree();
00297     }
00298 
00299     // Build a list of all indices in quadtree
00300     void getAllIndices(int* indices)
00301     {
00302         getAllIndices(indices, 0);
00303     }
00304 
00305     int getDepth()
00306     {
00307         if(is_leaf) return 1;
00308         return 1 + std::max(std::max(northWest->getDepth(),
00309                                      northEast->getDepth()),
00310                             std::max(southWest->getDepth(),
00311                                      southEast->getDepth()));                
00312     }
00313 
00314     // Compute non-edge forces using Barnes-Hut algorithm
00315     void computeNonEdgeForces(int point_index, double theta, double neg_f[], double* sum_Q)
00316     {
00317         
00318         // Make sure that we spend no time on empty nodes or self-interactions
00319         if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index)) return;
00320         
00321         // Compute distance between point and center-of-mass
00322         double D = .0;
00323         int ind = point_index * QT_NO_DIMS;
00324         for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = data[ind + d];
00325         for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= center_of_mass[d];
00326         for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
00327         
00328         // Check whether we can use this node as a "summary"
00329         if(is_leaf || std::max(boundary.hh, boundary.hw)/sqrt(D) < theta) {
00330         
00331             // Compute and add t-SNE force between point and current node
00332             double Q = 1.0 / (1.0 + D);
00333             *sum_Q += cum_size * Q;
00334             double mult = cum_size * Q * Q;
00335             for(int d = 0; d < QT_NO_DIMS; d++) neg_f[d] += mult * buff[d];
00336         }
00337         else {
00338 
00339             // Recursively apply Barnes-Hut to children
00340             northWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
00341             northEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
00342             southWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
00343             southEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
00344         }
00345     }
00346 
00347     // Computes edge forces
00348     void computeEdgeForces(int* row_P, int* col_P, double* val_P, int N, double* pos_f)
00349     {
00350         // Loop over all edges in the graph
00351         int ind1, ind2;
00352         double D;
00353         for(int n = 0; n < N; n++) {
00354             ind1 = n * QT_NO_DIMS;
00355             for(int i = row_P[n]; i < row_P[n + 1]; i++) {
00356             
00357                 // Compute pairwise distance and Q-value
00358                 D = .0;
00359                 ind2 = col_P[i] * QT_NO_DIMS;
00360                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = data[ind1 + d];
00361                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= data[ind2 + d];
00362                 for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
00363                 D = val_P[i] / (1.0 + D);
00364                 
00365                 // Sum positive force
00366                 for(int d = 0; d < QT_NO_DIMS; d++) pos_f[ind1 + d] += D * buff[d];
00367             }
00368         }
00369     }
00370 
00371     // Print out tree
00372     void print()
00373     {
00374         if(cum_size == 0) {
00375             printf("Empty node\n");
00376             return;
00377         }
00378 
00379         if(is_leaf) {
00380             printf("Leaf node; data = [");
00381             for(int i = 0; i < size; i++) {
00382                 double* point = data + index[i] * QT_NO_DIMS;
00383                 for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", point[d]);
00384                 printf(" (index = %d)", index[i]);
00385                 if(i < size - 1) printf("\n");
00386                 else printf("]\n");
00387             }        
00388         }
00389         else {
00390             printf("Intersection node with center-of-mass = [");
00391             for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", center_of_mass[d]);
00392             printf("]; children are:\n");
00393             northEast->print();
00394             northWest->print();
00395             southEast->print();
00396             southWest->print();
00397         }
00398     }
00399 
00400 private:
00401 
00402     QuadTree(const QuadTree&);
00403     QuadTree& operator=(const QuadTree&);
00404 
00405     void init(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
00406     {
00407         parent = inp_parent;
00408         data = inp_data;
00409         is_leaf = true;
00410         size = 0;
00411         cum_size = 0;
00412         boundary.x  = inp_x;
00413         boundary.y  = inp_y;
00414         boundary.hw = inp_hw;
00415         boundary.hh = inp_hh;
00416         northWest = NULL;
00417         northEast = NULL;
00418         southWest = NULL;
00419         southEast = NULL;
00420         for(int i = 0; i < QT_NO_DIMS; i++) center_of_mass[i] = .0;
00421     }
00422 
00423     // Build quadtree on dataset
00424     void fill(int N)
00425     {
00426         for(int i = 0; i < N; i++) insert(i);
00427     }
00428 
00429     // Build a list of all indices in quadtree
00430     int getAllIndices(int* indices, int loc)
00431     {
00432         
00433         // Gather indices in current quadrant
00434         for(int i = 0; i < size; i++) indices[loc + i] = index[i];
00435         loc += size;
00436         
00437         // Gather indices in children
00438         if(!is_leaf) {
00439             loc = northWest->getAllIndices(indices, loc);
00440             loc = northEast->getAllIndices(indices, loc);
00441             loc = southWest->getAllIndices(indices, loc);
00442             loc = southEast->getAllIndices(indices, loc);
00443         }
00444         return loc;
00445     }
00446 
00447     //bool isChild(int test_index, int start, int end);
00448 };
00449 
00450 }
00451 
00452 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines