1 // An example of this class is at the bottom
2 // This uses two other snippets, PtrList and ValList . You need
3 // them present for this to compile. The PtrList and ValList
4 // snippets _must_ be split into multiple files for this to compile.
5 // However, _this_ snippet need not be split to compile the example
6 // at the bottom of this code. If you do split this snippet up,
7 // take care to uncomment the #include lines indicated.
8
9 // ------------------ Start of Interval.H
10 // Derive your own classes from the Interval class.
11
12 #ifndef Interval_H
13 #define Interval_H
14
15 class Interval {
16 protected:
17 virtual void set( double a, double b ) = 0 ;
18 public:
19 virtual double min() const = 0 ;
20 virtual double max() const = 0 ;
21 } ;
22
23
24 #endif // Interval_H
25
26 // ------------------- Start of Interval.C
27
28 bool operator== ( const Interval& a, const Interval& b )
29 {
30 if ( a.min() == b.min() &&
31 a.max() == b.max() )
32 return true ;
33 return false ;
34 }
35
36 // ------------------- Start of IntervalTree.H
37 #ifndef IntervalTree_H
38 #define IntervalTree_H
39
40 #include "PtrList.H"
41 #include "ValList.H"
42
43 /**@name An Interval Tree.
44 * This class stores lists of items that span a range of values (stored
45 * as doubles). The items stored must have member functions min() and max()
46 * defined as doubles so that they can be inserted into the tree. The items
47 * are stored by reference, not by value.
48 *
49 * Pseudocode for this algorithm is from Bajaj, Pascucci, and Schikore,
50 * "Accelerated Isocontouring of Scalar Fields", \emph{Data Visualization
51 * Techniques}, Wiley & Sons, 1998 (draft pp. 30-31)
52 */
53 //@{
54
55 /// The nodes of the IntervalTree
56 template< class T >
57 class IntervalTreeNode {
58 public:
59 /**@name Data */
60 //@{
61 /// Split value
62 double m_split ;
63 /// List of intervals spanning m_split sorted by increasing minimum value
64 PtrList<T> m_minlist ;
65 /// List of intervals spanning m_split sorted by decreasing maximum value
66 PtrList<T> m_maxlist ;
67 //@}
68
69 /**@name Data */
70 //@{
71 /// add an interval to the node
72 void insert( const T* item ) ;
73 //@}
74 } ;
75
76
77 typedef ValList< double > EndPtList ;
78
79 /// The interval tree interface
80 template< class T >
81 class IntervalTree {
82 protected:
83 /**@name Data*/
84 //@{
85 /// Endpoint vector of the interval tree
86 IntervalTreeNode<T>* m_unique ;
87 /// Number of items in the vector
88 long m_nu ;
89 /// Number of intervals stored in the whole tree
90 long m_ni ;
91 /// The isovalue for the current retrieve request
92 double m_value ;
93 /// The left boundary of the tree we are considering for the current traverse
94 long m_lf ;
95 /// The right boundary of the tree we are considering for the current traverse
96 long m_rt ;
97 /// The current traverse node for a retrieve request
98 IntervalTreeNode<T>* m_trav ;
99 /// Which interval list to report for m_trav
100 bool m_by_max ;
101 //@}
102
103 /**@name Members*/
104 //@{
105 /// Private recursive routine to insert an interval at the correct node
106 void p_insert( T* interval, long lf, long rt ) ;
107 //@}
108
109 public:
110 /**@name Members */
111 //@{
112 /// Create a tree for representing intervals with the given, sorted endpoints.
113 IntervalTree( EndPtList& unique ) ;
114 /// Destructor
115 ~IntervalTree() ;
116 /// Number of unique values
117 long numUnique() const { return m_nu ; } ;
118 /// Number of intervals stored
119 long numIntervals() const { return m_ni ; } ;
120 /// Insert an interval
121 void insert( T* interval ) ;
122 /// Prepare to traverse the list of nodes covering a given isovalue
123 void init( double value ) ;
124 /// Get the next node covering the isovalue from init(). NULL terminated.
125 const T* next() ;
126 ///
127 void dump() ;
128 //@}
129 } ;
130
131 //@}
132 #endif // IntervalTree_H
133
134 // ------------------ Start of IntervalTree.C
135 #include <stdlib.h>
136 // Uncomment this line if you split the snippet into multiple files
137 //#include "IntervalTree.H"
138
139
140 template< class T >
141 IntervalTree<T>::IntervalTree( EndPtList& unique )
142 {
143 if (unique.size() < 2) // If there's not enough points, barf and die
144 #ifdef THROW
145 throw new int(-1) ;
146 #else
147 exit(-1) ;
148 #endif
149
150 // sort the list
151 // Assume presorted or uncomment this: unique.sort() ;
152
153 m_nu = unique.size() + 1 ;
154 m_unique = new IntervalTreeNode<T>[ m_nu ] ;
155
156 // Since endpoints go to -\infty to \infty, we must take care
157 // to assign meaningful split values. In this case, we double
158 // the difference between the last two points and move that far away
159 // from the last one.
160 m_unique[ 0 ].m_split = 3.0 * unique[0] - 2.0 * unique[1] ;
161 m_unique[ m_nu - 1 ].m_split = 3.0 * unique[m_nu-2] - 2.0 * unique[m_nu-3] ;
162
163 // loop through list setting up IntervalTreeNodes with list values
164 unique.init() ;
165 double* current ;
166 double previous = *unique.next() ;
167 int i = 1 ;
168 while ( current = unique.next() ) {
169 m_unique[i].m_split = (*current + previous)/2.0 ;
170 previous = *current ;
171 i++ ;
172 }
173
174 m_lf = 0 ;
175 m_rt = m_nu - 1 ;
176 m_trav = 0 ;
177 m_value = 0.0 ;
178 m_ni = 0 ;
179
180 }
181
182 template< class T >
183 IntervalTree<T>::~IntervalTree()
184 {
185 delete [] m_unique ;
186 }
187
188 template< class T >
189 void IntervalTree<T>::init( double value )
190 {
191 m_value = value ;
192 m_lf = 0 ;
193 m_rt = m_nu - 1;
194 m_trav = 0 ;
195 }
196
197 template< class T >
198 const T* IntervalTree<T>::next()
199 {
200 T* next ;
201
202 if (!m_trav) {
203 if ( ! (m_lf <= m_rt) )
204 return 0 ; // done with the whole traverse, return null to signal
205
206 long mid = (m_lf + m_rt)/2 ;
207 if ( m_value > m_unique[ mid ].m_split ) {
208 m_trav = &m_unique[mid] ;
209 m_lf = mid + 1 ;
210 m_by_max = true ;
211 m_trav->m_maxlist.init() ;
212 } else { // ( m_value < m_unique[ mid ].m_split ) {
213 m_trav = &m_unique[mid] ;
214 m_rt = mid - 1 ;
215 m_by_max = false ;
216 m_trav->m_minlist.init() ;
217 }
218 }
219 if ( m_by_max ) {
220 next = m_trav->m_maxlist.next() ;
221 if ( next && next->max() < m_value ) // go no further
222 next = 0 ;
223 } else {
224 next = m_trav->m_minlist.next() ;
225 if ( next && next->min() > m_value ) // go no further
226 next = 0 ;
227 }
228 if ( !next ) {
229 m_trav = 0 ; // reset so we take the first case on the switch next time
230 return IntervalTree<T>::next() ; // do just that to get a return value
231 }
232 return next ;
233 }
234
235 template< class T >
236 void IntervalTree<T>::p_insert( T* interval, long lf, long rt )
237 {
238 long mid = (lf + rt)/2 ;
239 if ( m_unique[ mid ].m_split > interval->max() ) {
240 p_insert( interval, lf, mid - 1 ) ;
241 } else if ( m_unique[ mid ].m_split < interval->min() ) {
242 p_insert( interval, mid + 1, rt ) ;
243 } else {
244 m_unique[ mid ].insert( interval ) ;
245 }
246 }
247
248 template< class T >
249 void IntervalTree<T>::insert( T* interval )
250 {
251 m_ni++ ;
252 p_insert( interval, 0, m_nu - 1 ) ;
253 }
254
255 template< class T >
256 void IntervalTree<T>::dump()
257 {
258 for (int i=0; i<m_nu; i++) {
259 fprintf( stderr, "\nSplit value %f\n", m_unique[i].m_split ) ;
260 T* tmp ;
261 fprintf( stderr, "Min List:\n\t" ) ;
262 if (!m_unique[i].m_minlist.size() ) {
263 fprintf( stderr, "empty\n" ) ;
264 } else {
265 m_unique[i].m_minlist.init() ;
266 while (tmp = m_unique[i].m_minlist.next() )
267 fprintf( stderr, "[ %1f %1f ] ", tmp->min(), tmp->max() ) ;
268 fprintf( stderr, "\n" ) ;
269 }
270 fprintf( stderr, "Max List:\n\t" ) ;
271 if (!m_unique[i].m_minlist.size() ) {
272 fprintf( stderr, "empty\n" ) ;
273 } else {
274 m_unique[i].m_maxlist.init() ;
275 while (tmp = m_unique[i].m_maxlist.next() )
276 fprintf( stderr, "[ %1f %1f ] ", tmp->min(), tmp->max() ) ;
277 fprintf( stderr, "\n" ) ;
278 }
279 }
280 }
281
282 template< class T >
283 void IntervalTreeNode<T>::insert( const T* interval )
284 {
285 // Aiieeee. This is an ugly kludge. We should use an AVL tree or some such!!!
286
287 // Jimminy!!! another kludge already:
288 // stage has two significant bits.
289 // bit 0 == interval has not been placed in maxlist
290 // bit 1 == interval has not been placed in minlist
291 int stage = 3 ;
292 m_minlist.init() ;
293 m_maxlist.init() ;
294 T* tmin ;
295 T* tmax ;
296 bool stopme = false ;
297
298 while ( stage && !stopme ) {
299 if (tmin = (T*) m_minlist.next() ) {
300 if ( (stage & 2) && interval->min() < tmin->min() ) {
301 stage -= 2 ;
302 m_minlist.insert_before( interval ) ;
303 }
304 }
305
306 if (tmax = (T*) m_maxlist.next() ) {
307 if ( (stage & 1) && interval->max() > tmax->max() ) {
308 stage -= 1 ;
309 m_maxlist.insert_before( interval ) ;
310 }
311 }
312
313 if (!tmin && !tmax)
314 stopme = true ;
315 }
316
317 if (stage & 2) {
318 m_minlist.append( (T*) interval ) ;
319 }
320 if (stage & 1) {
321 m_maxlist.append( (T*) interval ) ;
322 }
323 }
324
325 // --------------- Start of test_IntervalTree.C
326 /* This is the sample program that shows how IntervalTree works.
327 * Just compile and link it. When you run it, pass a list of numbers and
328 * it will report a list of intervals that span each number you pass it.
329 * For instance, running it with an argument of "0 2 5" should yield:
330 Query at w = 0.000000
331 [ -inf 4.000000 ]
332 [ -inf 2.000000 ]
333 ---
334
335 Query at w = 2.000000
336 [ -inf 4.000000 ]
337 [ 1.000000 3.000000 ]
338 [ 1.000000 4.000000 ]
339 [ 1.000000 inf ]
340 [ 2.000000 4.000000 ]
341 [ 2.000000 3.000000 ]
342 [ -inf 2.000000 ]
343 [ 1.000000 2.000000 ]
344 ---
345
346 Query at w = 5.000000
347 [ 1.000000 inf ]
348 ---
349 */
350
351 #include <stdlib.h>
352
353 #include "PtrList.H"
354 #include "ValList.H"
355 // Uncomment this if you split the snippet into multiple files
356 //#include "Interval.H"
357 //#include "IntervalTree.H"
358
359 // These values cause compilation warnings, but that's OK.
360 double dtNaN = 0.0/0.0 ;
361 double dtPInf = 1.0/0.0 ;
362 double dtNInf = -1.0/0.0 ;
363
364 class MyInterval : public Interval {
365 public:
366 double data[2] ;
367 double min() const { return data[0] ; } ;
368 double max() const { return data[1] ; } ;
369 void set( double a, double b ) { data[0] = a; data[1] = b ; } ;
370 } ;
371
372 // Template instantiations
373 #ifdef __GNUC__
374 # ifndef NO_EQUALITY_OP
375 # define NO_EQUALITY_OP
376 # endif
377 #include "PtrList.C"
378 #include "ValList.C"
379 // Uncomment this if you split the snippet into multiple files
380 //#include "IntervalTree.C"
381 template class ValList< double > ;
382 template class PtrList< Interval > ;
383 template class IntervalTree< Interval > ;
384 template class IntervalTreeNode< Interval > ;
385 #undef NO_EQUALITY_OP
386 #endif // __GNUC__
387
388 void report_query( IntervalTree<MyInterval>& itree, double val )
389 {
390 Interval* j ;
391
392 fprintf( stdout, "Query at w = %f\n", val) ;
393 itree.init( val ) ;
394 while ( j = (Interval*) itree.next() ) {
395 fprintf( stdout, "\t[ %f %f ]\n", j->min(), j->max() ) ;
396 }
397 fprintf( stdout, "---\n\n", val) ;
398 }
399
400 int main(int argc, char**argv) {
401 ValList<double> endpoints ;
402 endpoints.append( 1.0 ) ;
403 endpoints.append( 2.0 ) ;
404 endpoints.append( 3.0 ) ;
405 endpoints.append( 4.0 ) ;
406
407 MyInterval* i = new MyInterval[9] ;
408 i[0].set( 1.0, 3.0 ) ;
409 i[1].set( 2.0, 4.0 ) ;
410 i[2].set( 1.0, 4.0 ) ;
411 i[3].set( 2.0, 3.0 ) ;
412 i[4].set( 1.0, 2.0 ) ;
413 i[5].set( 3.0, 4.0 ) ;
414 i[6].set( dtNInf, 4.0 ) ;
415 i[7].set( dtNInf, 2.0 ) ;
416 i[8].set( 1.0, dtPInf ) ;
417
418 IntervalTree<MyInterval> itree( endpoints ) ;
419 itree.insert( &i[0] ) ;
420 itree.insert( &i[1] ) ;
421 itree.insert( &i[2] ) ;
422 itree.insert( &i[3] ) ;
423 itree.insert( &i[4] ) ;
424 itree.insert( &i[5] ) ;
425 itree.insert( &i[6] ) ;
426 itree.insert( &i[7] ) ;
427 itree.insert( &i[8] ) ;
428
429 // itree.dump() ;
430
431 for( int counter = 1; counter < argc; counter++) {
432 double request = atof( argv[ counter ] ) ;
433 report_query( itree, request ) ;
434 }
435
436 if ( argc < 2 ) {
437 fprintf( stderr, "Hey, you should give a list of numbers to query!\n" ) ;
438 }
439
440 return 0 ;
441 }
442