SISCone  2.0.6
quadtree.cpp
1 
2 // File: quadtree.cpp //
3 // Description: source file for quadtree management (Cquadtree class) //
4 // This file is part of the SISCone project. //
5 // For more details, see http://projects.hepforge.org/siscone //
6 // //
7 // Copyright (c) 2006 Gavin Salam and Gregory Soyez //
8 // //
9 // This program is free software; you can redistribute it and/or modify //
10 // it under the terms of the GNU General Public License as published by //
11 // the Free Software Foundation; either version 2 of the License, or //
12 // (at your option) any later version. //
13 // //
14 // This program is distributed in the hope that it will be useful, //
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
17 // GNU General Public License for more details. //
18 // //
19 // You should have received a copy of the GNU General Public License //
20 // along with this program; if not, write to the Free Software //
21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
22 // //
23 // $Revision:: 320 $//
24 // $Date:: 2011-11-15 09:54:50 +0100 (Tue, 15 Nov 2011) $//
26 
27 #include "quadtree.h"
28 #include <math.h>
29 #include <stdio.h>
30 #include <iostream>
31 
32 namespace siscone{
33 
34 using namespace std;
35 
36 /*******************************************************************
37  * Cquadtree implementation *
38  * Implementation of a 2D quadtree. *
39  * This class implements the traditional two-dimensional quadtree. *
40  * The elements at each node are of 'Cmomentum' type. *
41  *******************************************************************/
42 
43 // default ctor
44 //--------------
46  v = NULL;
47 
48  children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL;
49  has_child = false;
50 }
51 
52 
53 // ctor with initialisation (see init for details)
54 //--------------------------
55 Cquadtree::Cquadtree(double _x, double _y, double _half_size_x, double _half_size_y){
56  v = NULL;
57 
58  children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL;
59  has_child = false;
60 
61  init(_x, _y, _half_size_x, _half_size_y);
62 }
63 
64 
65 // default destructor
66 // at destruction, everything is destroyed except
67 // physical values at the leaves
68 //------------------------------------------------
70  if (has_child){
71  if (v!=NULL) delete v;
72  delete children[0][0];
73  delete children[0][1];
74  delete children[1][0];
75  delete children[1][1];
76  }
77 }
78 
79 
80 /*
81  * init the tree.
82  * By initializing the tree, we mean setting the cell parameters
83  * and preparing the object to act as a seed for a new tree.
84  * - _x x-position of the center
85  * - _y y-position of the center
86  * - half_size_x half x-size of the cell
87  * - half_size_y half y-size of the cell
88  * return 0 on success, 1 on error. Note that if the cell
89  * is already filled, we return an error.
90  ******************************************************************/
91 int Cquadtree::init(double _x, double _y, double _half_size_x, double _half_size_y){
92  if (v!=NULL)
93  return 1;
94 
95  centre_x = _x;
96  centre_y = _y;
97  half_size_x = _half_size_x;
98  half_size_y = _half_size_y;
99 
100  return 0;
101 }
102 
103 
104 /*
105  * adding a particle to the tree.
106  * This method adds one vector to the quadtree structure which
107  * is updated consequently.
108  * - v vector to add
109  * return 0 on success 1 on error
110  ******************************************************************/
112  // Description of the method:
113  // --------------------------
114  // the addition process goes as follows:
115  // 1. check if the cell is empty, in which case, add the particle
116  // here and leave.
117  // 2. If there is a unique particle already inside,
118  // (a) create children
119  // (b) forward the existing particle to the appropriate child
120  // 3. Add current particle to this cell and forward to the
121  // adequate child
122  // NOTE: we assume in the whole procedure that the particle is
123  // indeed inside the cell !
124 
125  // step 1: the case of empty cells
126  if (v==NULL){
127  v = v_add;
128  return 0;
129  }
130 
131  // step 2: additional work if 1! particle already present
132  // we use the fact that only 1-particle systems have no child
133  if (!has_child){
134  double new_half_size_x = 0.5*half_size_x;
135  double new_half_size_y = 0.5*half_size_y;
136  // create children
137  children[0][0] = new Cquadtree(centre_x-new_half_size_x, centre_y-new_half_size_y,
138  new_half_size_x, new_half_size_y);
139  children[0][1] = new Cquadtree(centre_x-new_half_size_x, centre_y+new_half_size_y,
140  new_half_size_x, new_half_size_y);
141  children[1][0] = new Cquadtree(centre_x+new_half_size_x, centre_y-new_half_size_y,
142  new_half_size_x, new_half_size_y);
143  children[1][1] = new Cquadtree(centre_x+new_half_size_x, centre_y+new_half_size_y,
144  new_half_size_x, new_half_size_y);
145 
146  has_child = true;
147 
148  // forward to child
149  //? The following line assumes 'true'==1 and 'false'==0
150  // Note: v being a single particle, eta and phi are correct
151  children[v->eta>centre_x][v->phi>centre_y]->add(v);
152 
153  // copy physical params
154  v = new Cmomentum(*v);
155  }
156 
157  // step 3: add new particle
158  // Note: v_add being a single particle, eta and phi are correct
159  children[v_add->eta>centre_x][v_add->phi>centre_y]->add(v_add);
160  *v+=*v_add;
161 
162  return 0;
163 }
164 
165 
166 /*
167  * circle intersection.
168  * computes the intersection with a circle of given centre and radius.
169  * The output takes the form of a quadtree with all squares included
170  * in the circle.
171  * - cx circle centre x coordinate
172  * - cy circle centre y coordinate
173  * - cR2 circle radius SQUARED
174  * return the checksum for the intersection
175  ******************************************************************/
176 Creference Cquadtree::circle_intersect(double cx, double cy, double cR2){
177  // Description of the method:
178  // --------------------------
179  // 1. check if cell is empty => no intersection
180  // 2. if cell has 1! particle, check if it is inside the circle.
181  // If yes, add it and return, if not simply return.
182  // 3. check if the circle intersects the square. If not, return.
183  // 4. check if the square is inside the circle.
184  // If yes, add it to qt and return.
185  // 5. check intersections with children.
186 
187  // step 1: if there is no particle inside te square, no reason to go further
188  if (v==NULL)
189  return Creference();
190 
191  double dx, dy;
192 
193  // step 2: if there is only one particle inside the square, test if it is in
194  // the circle, in which case return associated reference
195  if (!has_child){
196  // compute the distance
197  // Note: v has only one particle => eta and phi are defined
198  dx = cx - v->eta;
199  dy = fabs(cy - v->phi);
200  if (dy>M_PI)
201  dy -= 2.0*M_PI;
202 
203  // test distance
204  if (dx*dx+dy*dy<cR2){
205  return v->ref;
206  }
207 
208  return Creference();
209  }
210 
211  // step 3: check if there is an intersection
212  //double ryp, rym;
213  double dx_c, dy_c;
214 
215  // store distance with the centre of the square
216  dx_c = fabs(cx-centre_x);
217  dy_c = fabs(cy-centre_y);
218  if (dy_c>M_PI) dy_c = 2.0*M_PI-dy_c;
219 
220  // compute (minimal) the distance (pay attention to the periodicity in phi).
221  dx = dx_c-half_size_x;
222  if (dx<0) dx=0;
223  dy = dy_c-half_size_y;
224  if (dy<0) dy=0;
225 
226  // check the distance
227  if (dx*dx+dy*dy>=cR2){
228  // no intersection
229  return Creference();
230  }
231 
232  // step 4: check if included
233 
234  // compute the (maximal) distance
235  dx = dx_c+half_size_x;
236  dy = dy_c+half_size_y;
237  if (dy>M_PI) dy = M_PI;
238 
239  // compute the distance
240  if (dx*dx+dy*dy<cR2){
241  return v->ref;
242  }
243 
244  // step 5: the square is not fully in. Recurse to children
245  return children[0][0]->circle_intersect(cx, cy, cR2)
246  + children[0][1]->circle_intersect(cx, cy, cR2)
247  + children[1][0]->circle_intersect(cx, cy, cR2)
248  + children[1][1]->circle_intersect(cx, cy, cR2);
249 }
250 
251 
252 /*
253  * output a data file for drawing the grid.
254  * This can be used to output a data file containing all the
255  * grid subdivisions. The file contents is as follows:
256  * first and second columns give center of the cell, the third
257  * gives the size.
258  * - flux opened stream to write to
259  * return 0 on success, 1 on error
260  ******************************************************************/
261 int Cquadtree::save(FILE *flux){
262 
263  if (flux==NULL)
264  return 1;
265 
266  if (has_child){
267  fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y);
268  children[0][0]->save(flux);
269  children[0][1]->save(flux);
270  children[1][0]->save(flux);
271  children[1][1]->save(flux);
272  }
273 
274  return 0;
275 }
276 
277 
278 /*
279  * output a data file for drawing the tree leaves.
280  * This can be used to output a data file containing all the
281  * tree leaves. The file contents is as follows:
282  * first and second columns give center of the cell, the third
283  * gives the size.
284  * - flux opened stream to write to
285  * return 0 on success, 1 on error
286  ******************************************************************/
287 int Cquadtree::save_leaves(FILE *flux){
288 
289  if (flux==NULL)
290  return 1;
291 
292  if (has_child){
293  if (children[0][0]!=NULL) children[0][0]->save_leaves(flux);
294  if (children[0][1]!=NULL) children[0][1]->save_leaves(flux);
295  if (children[1][0]!=NULL) children[1][0]->save_leaves(flux);
296  if (children[1][1]!=NULL) children[1][1]->save_leaves(flux);
297  } else {
298  fprintf(flux, "%e\t%e\t%e\t%e\n", centre_x, centre_y, half_size_x, half_size_y);
299  }
300 
301  return 0;
302 }
303 
304 }
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Mon May 6 2013 11:30:35 for SISCone by  Doxygen 1.8.3.1