LCOV - code coverage report
Current view: top level - external/barnes_hut_sne - quadtree.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 120 120 100.0 %
Date: 2013-05-24 Functions: 10 10 100.0 %
Branches: 108 122 88.5 %

           Branch data     Line data    Source code
       1                 :            : /**
       2                 :            :  * Copyright (c) 2013, Laurens van der Maaten (Delft University of Technology)
       3                 :            :  * All rights reserved.
       4                 :            :  *
       5                 :            :  * Redistribution and use in source and binary forms, with or without
       6                 :            :  * modification, are permitted provided that the following conditions are met:
       7                 :            :  * 1. Redistributions of source code must retain the above copyright
       8                 :            :  *    notice, this list of conditions and the following disclaimer.
       9                 :            :  * 2. Redistributions in binary form must reproduce the above copyright
      10                 :            :  *    notice, this list of conditions and the following disclaimer in the
      11                 :            :  *    documentation and/or other materials provided with the distribution.
      12                 :            :  * 3. All advertising materials mentioning features or use of this software
      13                 :            :  *    must display the following acknowledgement:
      14                 :            :  *    This product includes software developed by the Delft University of Technology.
      15                 :            :  * 4. Neither the name of the Delft University of Technology nor the names of
      16                 :            :  *    its contributors may be used to endorse or promote products derived from
      17                 :            :  *    this software without specific prior written permission.
      18                 :            :  *
      19                 :            :  * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
      20                 :            :  * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
      21                 :            :  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
      22                 :            :  * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
      23                 :            :  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
      24                 :            :  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
      25                 :            :  * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
      26                 :            :  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
      27                 :            :  * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
      28                 :            :  * OF SUCH DAMAGE.
      29                 :            :  *
      30                 :            :  */
      31                 :            : 
      32                 :            : #include <math.h>
      33                 :            : #include <float.h>
      34                 :            : #include <stdlib.h>
      35                 :            : #include <stdio.h>
      36                 :            : #include <algorithm>
      37                 :            : 
      38                 :            : #ifndef QUADTREE_H
      39                 :            : #define QUADTREE_H
      40                 :            : 
      41                 :            : namespace tsne
      42                 :            : {
      43                 :            : 
      44                 :            : class Cell {
      45                 :            : 
      46                 :            : public:
      47                 :            : 
      48                 :            :         double x;
      49                 :            :         double y;
      50                 :            :         double hw;
      51                 :            :         double hh;
      52                 :            : 
      53                 :     621192 :         bool containsPoint(double point[])
      54                 :            :         {
      55         [ +  + ]:     621192 :                 if(x - hw > point[0]) return false;
      56         [ +  + ]:     573077 :                 if(x + hw < point[0]) return false;
      57         [ -  + ]:     391773 :                 if(y - hh > point[1]) return false;
      58         [ +  + ]:     391773 :                 if(y + hh < point[1]) return false;
      59                 :     621192 :                 return true;
      60                 :            :         }
      61                 :            : 
      62                 :            : };
      63                 :            : 
      64                 :            : 
      65                 :            : class QuadTree
      66                 :            : {
      67                 :            : 
      68                 :            :         // Fixed constants    
      69                 :            :         static const int QT_NO_DIMS = 2;
      70                 :            :         static const int QT_NODE_CAPACITY = 1;
      71                 :            : 
      72                 :            :         // A buffer we use when doing force computations
      73                 :            :         double buff[QT_NO_DIMS];
      74                 :            : 
      75                 :            :         // Properties of this node in the tree
      76                 :            :         QuadTree* parent;
      77                 :            :         bool is_leaf;
      78                 :            :         int size;
      79                 :            :         int cum_size;
      80                 :            : 
      81                 :            :         // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
      82                 :            :         Cell boundary;
      83                 :            : 
      84                 :            :         // Indices in this quad tree node, corresponding center-of-mass, and list of all children
      85                 :            :         double* data;
      86                 :            :         double center_of_mass[QT_NO_DIMS];
      87                 :            :         int index[QT_NODE_CAPACITY];
      88                 :            : 
      89                 :            :         // Children
      90                 :            :         QuadTree* northWest;
      91                 :            :         QuadTree* northEast;
      92                 :            :         QuadTree* southWest;
      93                 :            :         QuadTree* southEast;
      94                 :            : 
      95                 :            : public:
      96                 :            : 
      97                 :            :         // Default constructor for quadtree -- build tree, too!
      98                 :       1000 :         QuadTree(double* inp_data, int N) : 
      99                 :            :                 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
     100                 :       1000 :                 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
     101                 :            :         {
     102                 :            :                 // Compute mean, width, and height of current map (boundaries of quadtree)
     103         [ +  + ]:       3000 :                 double* mean_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] = .0;
     104         [ +  + ]:       3000 :                 double*  min_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++)  min_Y[d] =  DBL_MAX;
     105         [ +  + ]:       3000 :                 double*  max_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++)  max_Y[d] = -DBL_MAX;
     106         [ +  + ]:      51000 :                 for(int n = 0; n < N; n++) {
     107         [ +  + ]:     150000 :                         for(int d = 0; d < QT_NO_DIMS; d++) {
     108                 :     100000 :                                 mean_Y[d] += inp_data[n * QT_NO_DIMS + d];
     109         [ +  + ]:     100000 :                                 if(inp_data[n * QT_NO_DIMS + d] < min_Y[d]) min_Y[d] = inp_data[n * QT_NO_DIMS + d];
     110         [ +  + ]:     100000 :                                 if(inp_data[n * QT_NO_DIMS + d] > max_Y[d]) max_Y[d] = inp_data[n * QT_NO_DIMS + d];
     111                 :            :                         }
     112                 :            :                 }
     113         [ +  + ]:       3000 :                 for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] /= (double) N;
     114                 :            :                 
     115                 :            :                 // Construct quadtree
     116                 :       1000 :                 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,
     117                 :       2000 :                                                            std::max(max_Y[1] - mean_Y[1], mean_Y[1] - min_Y[1]) + 1e-5);
     118                 :       1000 :                 fill(N);
     119 [ +  - ][ +  - ]:       1000 :                 delete[] mean_Y; delete[] max_Y; delete[] min_Y;
                 [ +  - ]
     120                 :       1000 :         }
     121                 :            : 
     122                 :            :         // Constructor for quadtree with particular size and parent -- build the tree, too!
     123                 :            :         QuadTree(double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
     124                 :            :                 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
     125                 :            :                 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
     126                 :            :         {
     127                 :            :                 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
     128                 :            :         }
     129                 :            : 
     130                 :            :         // Constructor for quadtree with particular size and parent -- build the tree, too!
     131                 :            :         QuadTree(double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
     132                 :            :                 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
     133                 :            :                 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
     134                 :            :         {
     135                 :            :                 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
     136                 :            :                 fill(N);
     137                 :            :         }
     138                 :            : 
     139                 :            :         // Constructor for quadtree with particular size (do not fill the tree)
     140                 :            :         QuadTree(QuadTree* inp_parent, double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
     141                 :            :                 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
     142                 :            :                 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
     143                 :            :         {
     144                 :            :                 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
     145                 :            :                 fill(N);
     146                 :            :         }
     147                 :            : 
     148                 :            :         // Constructor for quadtree with particular size and parent (do not fill the tree)
     149                 :     161904 :         QuadTree(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
     150                 :            :                 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
     151                 :     161904 :                 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
     152                 :            :         {
     153                 :     161904 :                 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
     154                 :     161904 :         }
     155                 :            : 
     156                 :            :         // Destructor for quadtree
     157                 :     162904 :         ~QuadTree()
     158                 :            :         {
     159         [ +  + ]:     162904 :                 delete northWest;
     160         [ +  + ]:     162904 :                 delete northEast;
     161         [ +  + ]:     162904 :                 delete southWest;
     162         [ +  + ]:     162904 :                 delete southEast;
     163                 :     162904 :         }
     164                 :            : 
     165                 :            :         void setData(double* inp_data)
     166                 :            :         {
     167                 :            :                 data = inp_data;
     168                 :            :         }
     169                 :            : 
     170                 :            :         QuadTree* getParent()
     171                 :            :         {
     172                 :            :                 return parent;
     173                 :            :         }
     174                 :            : 
     175                 :            :         //void construct(Cell boundary);
     176                 :            : 
     177                 :            :         // Insert a point into the QuadTree
     178                 :     621192 :         bool insert(int new_index)
     179                 :            :         {
     180                 :            :                 // Ignore objects which do not belong in this quad tree
     181                 :     621192 :                 double* point = data + new_index * QT_NO_DIMS;
     182         [ +  + ]:     621192 :                 if(!boundary.containsPoint(point))
     183                 :     346179 :                         return false;
     184                 :            :                 
     185                 :            :                 // Online update of cumulative size and center-of-mass
     186                 :     275013 :                 cum_size++;
     187                 :     275013 :                 double mult1 = (double) (cum_size - 1) / (double) cum_size;
     188                 :     275013 :                 double mult2 = 1.0 / (double) cum_size;
     189         [ +  + ]:     825039 :                 for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] *= mult1;
     190         [ +  + ]:     825039 :                 for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] += mult2 * point[d];
     191                 :            :                 
     192                 :            :                 // If there is space in this quad tree and it is a leaf, add the object here
     193 [ +  + ][ +  + ]:     275013 :                 if(is_leaf && size < QT_NODE_CAPACITY) {
     194                 :      90476 :                         index[size] = new_index;
     195                 :      90476 :                         size++;
     196                 :      90476 :                         return true;
     197                 :            :                 }
     198                 :            :                 
     199                 :            :                 // Don't add duplicates for now (this is not very nice)
     200                 :     184537 :                 bool any_duplicate = false;
     201         [ +  + ]:     225013 :                 for(int n = 0; n < size; n++) {
     202                 :      40476 :                         bool duplicate = true;
     203         [ +  - ]:      40476 :                         for(int d = 0; d < QT_NO_DIMS; d++) {
     204         [ +  - ]:      40476 :                                 if(point[d] != data[index[n] * QT_NO_DIMS + d]) { duplicate = false; break; }
     205                 :            :                         }
     206                 :      40476 :                         any_duplicate = any_duplicate | duplicate;
     207                 :            :                 }
     208         [ -  + ]:     184537 :                 if(any_duplicate) return true;
     209                 :            :                 
     210                 :            :                 // Otherwise, we need to subdivide the current cell
     211         [ +  + ]:     184537 :                 if(is_leaf) subdivide();
     212                 :            :                 
     213                 :            :                 // Find out where the point can be inserted
     214         [ +  + ]:     184537 :                 if(northWest->insert(new_index)) return true;
     215         [ +  + ]:     130354 :                 if(northEast->insert(new_index)) return true;
     216         [ +  + ]:      95722 :                 if(southWest->insert(new_index)) return true;
     217         [ +  - ]:      58107 :                 if(southEast->insert(new_index)) return true;
     218                 :            :                 
     219                 :            :                 // Otherwise, the point cannot be inserted (this should never happen)
     220                 :     621192 :                 return false;
     221                 :            :         }
     222                 :            : 
     223                 :            :         // Create four children which fully divide this cell into four quads of equal area
     224                 :      40476 :         void subdivide()
     225                 :            :         { 
     226                 :            :                 // Create four children
     227         [ +  - ]:      40476 :                 northWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
     228         [ +  - ]:      40476 :                 northEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
     229         [ +  - ]:      40476 :                 southWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
     230         [ +  - ]:      40476 :                 southEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
     231                 :            :                 
     232                 :            :                 // Move existing points to correct children
     233         [ +  + ]:      80952 :                 for(int i = 0; i < size; i++) {
     234                 :      40476 :                         bool success = false;
     235         [ +  - ]:      40476 :                         if(!success) success = northWest->insert(index[i]);
     236         [ +  + ]:      40476 :                         if(!success) success = northEast->insert(index[i]);
     237         [ +  + ]:      40476 :                         if(!success) success = southWest->insert(index[i]);
     238         [ +  + ]:      40476 :                         if(!success) success = southEast->insert(index[i]);
     239                 :      40476 :                         index[i] = -1;
     240                 :            :                 }
     241                 :            :                 
     242                 :            :                 // Empty parent node
     243                 :      40476 :                 size = 0;
     244                 :      40476 :                 is_leaf = false;
     245                 :      40476 :         }
     246                 :            : 
     247                 :            :         // Checks whether the specified tree is correct
     248                 :            :         bool isCorrect()
     249                 :            :         {
     250                 :            :                 for(int n = 0; n < size; n++) {
     251                 :            :                         double* point = data + index[n] * QT_NO_DIMS;
     252                 :            :                         if(!boundary.containsPoint(point)) return false;
     253                 :            :                 }
     254                 :            :                 if(!is_leaf) return northWest->isCorrect() &&
     255                 :            :                                                         northEast->isCorrect() &&
     256                 :            :                                                         southWest->isCorrect() &&
     257                 :            :                                                         southEast->isCorrect();
     258                 :            :                 else return true;
     259                 :            :         }
     260                 :            : 
     261                 :            :         // Rebuilds a possibly incorrect tree (LAURENS: This function is not tested yet!)
     262                 :            :         void rebuildTree()
     263                 :            :         {
     264                 :            :                 for(int n = 0; n < size; n++) {
     265                 :            :                         // Check whether point is erroneous
     266                 :            :                         double* point = data + index[n] * QT_NO_DIMS;
     267                 :            :                         if(!boundary.containsPoint(point)) {
     268                 :            :                                 
     269                 :            :                                 // Remove erroneous point
     270                 :            :                                 int rem_index = index[n];
     271                 :            :                                 for(int m = n + 1; m < size; m++) index[m - 1] = index[m];
     272                 :            :                                 index[size - 1] = -1;
     273                 :            :                                 size--;
     274                 :            :                                 
     275                 :            :                                 // Update center-of-mass and counter in all parents
     276                 :            :                                 bool done = false;
     277                 :            :                                 QuadTree* node = this;
     278                 :            :                                 while(!done) {
     279                 :            :                                         for(int d = 0; d < QT_NO_DIMS; d++) {
     280                 :            :                                                 node->center_of_mass[d] = ((double) node->cum_size * node->center_of_mass[d] - point[d]) / (double) (node->cum_size - 1);
     281                 :            :                                         }
     282                 :            :                                         node->cum_size--;
     283                 :            :                                         if(node->getParent() == NULL) done = true;
     284                 :            :                                         else node = node->getParent();
     285                 :            :                                 }
     286                 :            :                                 
     287                 :            :                                 // Reinsert point in the root tree
     288                 :            :                                 node->insert(rem_index);
     289                 :            :                         }
     290                 :            :                 }    
     291                 :            :                 
     292                 :            :                 // Rebuild lower parts of the tree
     293                 :            :                 northWest->rebuildTree();
     294                 :            :                 northEast->rebuildTree();
     295                 :            :                 southWest->rebuildTree();
     296                 :            :                 southEast->rebuildTree();
     297                 :            :         }
     298                 :            : 
     299                 :            :         // Build a list of all indices in quadtree
     300                 :            :         void getAllIndices(int* indices)
     301                 :            :         {
     302                 :            :                 getAllIndices(indices, 0);
     303                 :            :         }
     304                 :            : 
     305                 :            :         int getDepth()
     306                 :            :         {
     307                 :            :                 if(is_leaf) return 1;
     308                 :            :                 return 1 + std::max(std::max(northWest->getDepth(),
     309                 :            :                                              northEast->getDepth()),
     310                 :            :                                     std::max(southWest->getDepth(),
     311                 :            :                                              southEast->getDepth()));                
     312                 :            :         }
     313                 :            : 
     314                 :            :         // Compute non-edge forces using Barnes-Hut algorithm
     315                 :    2417084 :         void computeNonEdgeForces(int point_index, double theta, double neg_f[], double* sum_Q)
     316                 :            :         {
     317                 :            :                 
     318                 :            :                 // Make sure that we spend no time on empty nodes or self-interactions
     319 [ +  + ][ +  + ]:    3884672 :                 if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index)) return;
         [ +  - ][ +  + ]
     320                 :            :                 
     321                 :            :                 // Compute distance between point and center-of-mass
     322                 :    1467588 :                 double D = .0;
     323                 :    1467588 :                 int ind = point_index * QT_NO_DIMS;
     324         [ +  + ]:    4402764 :                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = data[ind + d];
     325         [ +  + ]:    4402764 :                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= center_of_mass[d];
     326         [ +  + ]:    4402764 :                 for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
     327                 :            :                 
     328                 :            :                 // Check whether we can use this node as a "summary"
     329 [ +  + ][ +  + ]:    1467588 :                 if(is_leaf || std::max(boundary.hh, boundary.hw)/sqrt(D) < theta) {
                 [ +  + ]
     330                 :            :                 
     331                 :            :                         // Compute and add t-SNE force between point and current node
     332                 :     875817 :                         double Q = 1.0 / (1.0 + D);
     333                 :     875817 :                         *sum_Q += cum_size * Q;
     334                 :     875817 :                         double mult = cum_size * Q * Q;
     335         [ +  + ]:    2627451 :                         for(int d = 0; d < QT_NO_DIMS; d++) neg_f[d] += mult * buff[d];
     336                 :            :                 }
     337                 :            :                 else {
     338                 :            : 
     339                 :            :                         // Recursively apply Barnes-Hut to children
     340                 :     591771 :                         northWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
     341                 :     591771 :                         northEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
     342                 :     591771 :                         southWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
     343                 :     591771 :                         southEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
     344                 :            :                 }
     345                 :            :         }
     346                 :            : 
     347                 :            :         // Computes edge forces
     348                 :       1000 :         void computeEdgeForces(int* row_P, int* col_P, double* val_P, int N, double* pos_f)
     349                 :            :         {
     350                 :            :                 // Loop over all edges in the graph
     351                 :            :                 int ind1, ind2;
     352                 :            :                 double D;
     353         [ +  + ]:      51000 :                 for(int n = 0; n < N; n++) {
     354                 :      50000 :                         ind1 = n * QT_NO_DIMS;
     355         [ +  + ]:    1794000 :                         for(int i = row_P[n]; i < row_P[n + 1]; i++) {
     356                 :            :                         
     357                 :            :                                 // Compute pairwise distance and Q-value
     358                 :    1744000 :                                 D = .0;
     359                 :    1744000 :                                 ind2 = col_P[i] * QT_NO_DIMS;
     360         [ +  + ]:    5232000 :                                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = data[ind1 + d];
     361         [ +  + ]:    5232000 :                                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= data[ind2 + d];
     362         [ +  + ]:    5232000 :                                 for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
     363                 :    1744000 :                                 D = val_P[i] / (1.0 + D);
     364                 :            :                                 
     365                 :            :                                 // Sum positive force
     366         [ +  + ]:    5232000 :                                 for(int d = 0; d < QT_NO_DIMS; d++) pos_f[ind1 + d] += D * buff[d];
     367                 :            :                         }
     368                 :            :                 }
     369                 :       1000 :         }
     370                 :            : 
     371                 :            :         // Print out tree
     372                 :            :         void print()
     373                 :            :         {
     374                 :            :                 if(cum_size == 0) {
     375                 :            :                         printf("Empty node\n");
     376                 :            :                         return;
     377                 :            :                 }
     378                 :            : 
     379                 :            :                 if(is_leaf) {
     380                 :            :                         printf("Leaf node; data = [");
     381                 :            :                         for(int i = 0; i < size; i++) {
     382                 :            :                                 double* point = data + index[i] * QT_NO_DIMS;
     383                 :            :                                 for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", point[d]);
     384                 :            :                                 printf(" (index = %d)", index[i]);
     385                 :            :                                 if(i < size - 1) printf("\n");
     386                 :            :                                 else printf("]\n");
     387                 :            :                         }        
     388                 :            :                 }
     389                 :            :                 else {
     390                 :            :                         printf("Intersection node with center-of-mass = [");
     391                 :            :                         for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", center_of_mass[d]);
     392                 :            :                         printf("]; children are:\n");
     393                 :            :                         northEast->print();
     394                 :            :                         northWest->print();
     395                 :            :                         southEast->print();
     396                 :            :                         southWest->print();
     397                 :            :                 }
     398                 :            :         }
     399                 :            : 
     400                 :            : private:
     401                 :            : 
     402                 :            :         QuadTree(const QuadTree&);
     403                 :            :         QuadTree& operator=(const QuadTree&);
     404                 :            : 
     405                 :     162904 :         void init(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
     406                 :            :         {
     407                 :     162904 :                 parent = inp_parent;
     408                 :     162904 :                 data = inp_data;
     409                 :     162904 :                 is_leaf = true;
     410                 :     162904 :                 size = 0;
     411                 :     162904 :                 cum_size = 0;
     412                 :     162904 :                 boundary.x  = inp_x;
     413                 :     162904 :                 boundary.y  = inp_y;
     414                 :     162904 :                 boundary.hw = inp_hw;
     415                 :     162904 :                 boundary.hh = inp_hh;
     416                 :     162904 :                 northWest = NULL;
     417                 :     162904 :                 northEast = NULL;
     418                 :     162904 :                 southWest = NULL;
     419                 :     162904 :                 southEast = NULL;
     420         [ +  + ]:     488712 :                 for(int i = 0; i < QT_NO_DIMS; i++) center_of_mass[i] = .0;
     421                 :     162904 :         }
     422                 :            : 
     423                 :            :         // Build quadtree on dataset
     424                 :       1000 :         void fill(int N)
     425                 :            :         {
     426         [ +  + ]:      51000 :                 for(int i = 0; i < N; i++) insert(i);
     427                 :       1000 :         }
     428                 :            : 
     429                 :            :         // Build a list of all indices in quadtree
     430                 :            :         int getAllIndices(int* indices, int loc)
     431                 :            :         {
     432                 :            :                 
     433                 :            :                 // Gather indices in current quadrant
     434                 :            :                 for(int i = 0; i < size; i++) indices[loc + i] = index[i];
     435                 :            :                 loc += size;
     436                 :            :                 
     437                 :            :                 // Gather indices in children
     438                 :            :                 if(!is_leaf) {
     439                 :            :                         loc = northWest->getAllIndices(indices, loc);
     440                 :            :                         loc = northEast->getAllIndices(indices, loc);
     441                 :            :                         loc = southWest->getAllIndices(indices, loc);
     442                 :            :                         loc = southEast->getAllIndices(indices, loc);
     443                 :            :                 }
     444                 :            :                 return loc;
     445                 :            :         }
     446                 :            : 
     447                 :            :         //bool isChild(int test_index, int start, int end);
     448                 :            : };
     449                 :            : 
     450                 :            : }
     451                 :            : 
     452                 :            : #endif

Generated by: LCOV version 1.9