1 #ifndef LIBGEODECOMP_GEOMETRY_REGION_H
2 #define LIBGEODECOMP_GEOMETRY_REGION_H
3 
4 #include <libgeodecomp/geometry/coordbox.h>
5 #include <libgeodecomp/geometry/regionstreakiterator.h>
6 #include <libgeodecomp/geometry/streak.h>
7 #include <libgeodecomp/misc/stdcontaineroverloads.h>
8 
9 namespace LibGeoDecomp {
10 
11 class RegionTest;
12 
13 namespace RegionHelpers {
14 
15 class RegionCommonHelper
16 {
17 public:
pairCompareFirst(const std::pair<int,int> & a,const std::pair<int,int> & b)18     static inline bool pairCompareFirst(const std::pair<int, int>& a, const std::pair<int, int>& b)
19     {
20         return a.first < b.first;
21     }
22 
pairCompareSecond(const std::pair<int,int> & a,const std::pair<int,int> & b)23     static inline bool pairCompareSecond(const std::pair<int, int>& a, const std::pair<int, int>& b)
24     {
25         return a.second < b.second;
26     }
27 
28 protected:
29     typedef std::pair<int, int> IntPair;
30     typedef std::vector<IntPair> VecType;
31 
incRemainder(const VecType::iterator & start,const VecType::iterator & end,const int & inserts)32     inline void incRemainder(const VecType::iterator& start, const VecType::iterator& end, const int& inserts)
33     {
34         if (inserts == 0) {
35             return;
36         }
37 
38         for (VecType::iterator incrementer = start;
39              incrementer != end; ++incrementer) {
40             incrementer->second += inserts;
41         }
42     }
43 };
44 
45 template<int DIM>
46 class ConstructStreakFromIterators
47 {
48 public:
49     typedef std::pair<int, int> IntPair;
50     typedef std::vector<IntPair> VecType;
51 
52     template<int STREAK_DIM>
operator()53     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators)
54     {
55         ConstructStreakFromIterators<DIM - 1>()(streak, iterators);
56         streak->origin[DIM] = iterators[DIM]->first;
57     }
58 };
59 
60 template<>
61 class ConstructStreakFromIterators<0>
62 {
63 public:
64     typedef std::pair<int, int> IntPair;
65     typedef std::vector<IntPair> VecType;
66 
67     template<int STREAK_DIM>
operator()68     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators)
69     {
70         streak->origin[0] = iterators[0]->first;
71         streak->endX      = iterators[0]->second;
72     }
73 };
74 
75 template<int DIM>
76 class StreakIteratorInitSingleOffset
77 {
78 public:
79     typedef std::pair<int, int> IntPair;
80     typedef std::vector<IntPair> VecType;
81 
StreakIteratorInitSingleOffset(const std::size_t & offset)82     explicit StreakIteratorInitSingleOffset(const std::size_t& offset) :
83         offset(offset)
84     {}
85 
86     template<int STREAK_DIM, typename REGION>
operator()87     inline std::size_t operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
88     {
89         StreakIteratorInitSingleOffset<DIM - 1> delegate(offset);
90         std::size_t newOffset = delegate(streak, iterators, region);
91 
92         VecType::const_iterator upperBound =
93             std::upper_bound(region.indicesBegin(DIM),
94                              region.indicesEnd(DIM),
95                              IntPair(0, newOffset),
96                              RegionHelpers::RegionCommonHelper::pairCompareSecond);
97         iterators[DIM] = upperBound - 1;
98         newOffset =  iterators[DIM] - region.indicesBegin(DIM);
99 
100         return newOffset;
101     }
102 
103 private:
104     const std::size_t& offset;
105 };
106 
107 template<>
108 class StreakIteratorInitSingleOffset<0>
109 {
110 public:
111     typedef std::pair<int, int> IntPair;
112     typedef std::vector<IntPair> VecType;
113 
StreakIteratorInitSingleOffset(const std::size_t & offset)114     explicit StreakIteratorInitSingleOffset(const std::size_t& offset) :
115         offset(offset)
116     {}
117 
118     template<int STREAK_DIM, typename REGION>
operator()119     inline std::size_t operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
120     {
121         iterators[0] = region.indicesBegin(0) + offset;
122         return offset;
123     }
124 
125 private:
126     const std::size_t& offset;
127 };
128 
129 template<int DIM>
130 class StreakIteratorInitSingleOffsetWrapper
131 {
132 public:
133     typedef std::pair<int, int> IntPair;
134     typedef std::vector<IntPair> VecType;
135 
StreakIteratorInitSingleOffsetWrapper(const std::size_t & offset)136     explicit StreakIteratorInitSingleOffsetWrapper(const std::size_t& offset) :
137         offset(offset)
138     {}
139 
140     template<int STREAK_DIM, typename REGION>
operator()141     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
142     {
143         StreakIteratorInitSingleOffset<DIM> delegate(offset);
144         delegate(streak, iterators, region);
145         ConstructStreakFromIterators<DIM>()(streak, iterators);
146     }
147 
148 private:
149     const std::size_t& offset;
150 };
151 
152 template<int DIM, int COORD_DIM>
153 class StreakIteratorInitOffsets
154 {
155 public:
156     typedef std::pair<int, int> IntPair;
157     typedef std::vector<IntPair> VecType;
158 
StreakIteratorInitOffsets(const Coord<COORD_DIM> & offsets)159     explicit StreakIteratorInitOffsets(const Coord<COORD_DIM>& offsets) :
160         offsets(offsets)
161     {}
162 
163     template<int STREAK_DIM, typename REGION>
operator()164     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
165     {
166         iterators[DIM] = region.indicesBegin(DIM) + offsets[DIM];
167 
168         StreakIteratorInitOffsets<DIM - 1, COORD_DIM> delegate(offsets);
169         delegate(streak, iterators, region);
170     }
171 
172 private:
173     const Coord<COORD_DIM>& offsets;
174 };
175 
176 template<int COORD_DIM>
177 class StreakIteratorInitOffsets<0, COORD_DIM>
178 {
179 public:
180     typedef std::pair<int, int> IntPair;
181     typedef std::vector<IntPair> VecType;
182 
StreakIteratorInitOffsets(const Coord<COORD_DIM> & offsets)183     explicit StreakIteratorInitOffsets(const Coord<COORD_DIM>& offsets) :
184         offsets(offsets)
185     {}
186 
187     template<int STREAK_DIM, typename REGION>
operator()188     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
189     {
190         iterators[0] = region.indicesBegin(0) + offsets[0];
191 
192         if (int(region.indicesSize(0)) > offsets[0]) {
193             ConstructStreakFromIterators<STREAK_DIM - 1>()(streak, iterators);
194         }
195     }
196 
197 private:
198     const Coord<COORD_DIM>& offsets;
199 };
200 
201 template<int DIM>
202 class StreakIteratorInitBegin
203 {
204 public:
205     typedef std::pair<int, int> IntPair;
206     typedef std::vector<IntPair> VecType;
207 
208     template<int STREAK_DIM, typename REGION>
operator()209     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
210     {
211         iterators[DIM] = region.indicesBegin(DIM);
212         StreakIteratorInitBegin<DIM - 1>()(streak, iterators, region);
213     }
214 };
215 
216 template<>
217 class StreakIteratorInitBegin<0>
218 {
219 public:
220     typedef std::pair<int, int> IntPair;
221     typedef std::vector<IntPair> VecType;
222 
223     template<int STREAK_DIM, typename REGION>
operator()224     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
225     {
226         iterators[0] = region.indicesBegin(0);
227 
228         if (region.indicesSize(0) > 0) {
229             ConstructStreakFromIterators<STREAK_DIM - 1>()(streak, iterators);
230         }
231     }
232 };
233 
234 template<int DIM>
235 class StreakIteratorInitEnd
236 {
237 public:
238     typedef std::pair<int, int> IntPair;
239     typedef std::vector<IntPair> VecType;
240 
241     template<int STREAK_DIM, typename REGION>
operator()242     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
243     {
244         StreakIteratorInitEnd<DIM - 1>()(streak, iterators, region);
245         iterators[DIM] = region.indicesEnd(DIM);
246     }
247 };
248 
249 template<>
250 class StreakIteratorInitEnd<0>
251 {
252 public:
253     typedef std::pair<int, int> IntPair;
254     typedef std::vector<IntPair> VecType;
255 
256     template<int STREAK_DIM, typename REGION>
operator()257     inline void operator()(Streak<STREAK_DIM> *streak, VecType::const_iterator *iterators, const REGION& region) const
258     {
259         iterators[0] = region.indicesEnd(0);
260     }
261 };
262 
263 template<int DIM>
264 class RegionIntersectHelper
265 {
266 public:
267     template<int MY_DIM>
intersects(const Streak<MY_DIM> & s1,const Streak<MY_DIM> & s2)268     static inline bool intersects(const Streak<MY_DIM>& s1, const Streak<MY_DIM>& s2)
269     {
270         if (s1.origin[DIM] != s2.origin[DIM]) {
271             return false;
272         }
273 
274         return RegionIntersectHelper<DIM - 1>::intersects(s1, s2);
275     }
276 
277     template<int MY_DIM>
lessThan(const Streak<MY_DIM> & s1,const Streak<MY_DIM> & s2)278     static inline bool lessThan(const Streak<MY_DIM>& s1, const Streak<MY_DIM>& s2)
279     {
280         if (s1.origin[DIM] < s2.origin[DIM]) {
281             return true;
282         }
283         if (s1.origin[DIM] > s2.origin[DIM]) {
284             return false;
285         }
286 
287         return RegionIntersectHelper<DIM - 1>::lessThan(s1, s2);
288     }
289 };
290 
291 template<>
292 class RegionIntersectHelper<0>
293 {
294 public:
295     template<int MY_DIM>
intersects(const Streak<MY_DIM> & s1,const Streak<MY_DIM> & s2)296     static inline bool intersects(const Streak<MY_DIM>& s1, const Streak<MY_DIM>& s2)
297     {
298         return ((s1.origin.x() <= s2.origin.x() && s2.origin.x() < s1.endX) ||
299                 (s2.origin.x() <= s1.origin.x() && s1.origin.x() < s2.endX));
300 
301     }
302 
303     template<int MY_DIM>
lessThan(const Streak<MY_DIM> & s1,const Streak<MY_DIM> & s2)304     static inline bool lessThan(const Streak<MY_DIM>& s1, const Streak<MY_DIM>& s2)
305     {
306         return s1.endX < s2.endX;
307     }
308 };
309 
310 template<int DIM>
311 class RegionLookupHelper;
312 
313 template<int DIM>
314 class RegionInsertHelper;
315 
316 template<int DIM>
317 class RegionRemoveHelper;
318 
319 }
320 
321 /**
322  * Region stores a set of coordinates. It performs a run-length
323  * coding. Instead of storing complete Streak objects, these objects
324  * get split up and are stored implicitly in the hierarchical indices
325  * vectors.
326  */
327 template<int DIM>
328 class Region
329 {
330 public:
331     friend class Serialization;
332 
333     template<int MY_DIM> friend class RegionHelpers::RegionLookupHelper;
334     template<int MY_DIM> friend class RegionHelpers::RegionInsertHelper;
335     template<int MY_DIM> friend class RegionHelpers::RegionRemoveHelper;
336     friend class LibGeoDecomp::RegionTest;
337 
338     typedef std::pair<int, int> IntPair;
339     typedef std::vector<IntPair> VecType;
340     typedef RegionStreakIterator<DIM, Region<DIM> > StreakIterator;
341 
342     class Iterator : public std::iterator<std::forward_iterator_tag,
343                                           const Coord<DIM> >
344     {
345     public:
Iterator(const StreakIterator & streakIterator)346         inline explicit Iterator(const StreakIterator& streakIterator) :
347             streakIterator(streakIterator),
348             cursor(streakIterator->origin)
349         {}
350 
351         inline void operator++()
352         {
353             cursor.x()++;
354             if (cursor.x() >= streakIterator->endX) {
355                 ++streakIterator;
356                 cursor = streakIterator->origin;
357             }
358         }
359 
360         inline bool operator==(const Iterator& other) const
361         {
362             return streakIterator == other.streakIterator &&
363                 (streakIterator.endReached() || cursor == other.cursor);
364         }
365 
366         inline bool operator!=(const Iterator& other) const
367         {
368             return !(*this == other);
369         }
370 
371         inline const Coord<DIM>& operator*() const
372         {
373             return cursor;
374         }
375 
376         inline const Coord<DIM> *operator->() const
377         {
378             return &cursor;
379         }
380 
381     private:
382         StreakIterator streakIterator;
383         Coord<DIM> cursor;
384     };
385 
Region()386     inline Region() :
387         mySize(0),
388         geometryCacheTainted(false)
389     {}
390 
391     template<class ITERATOR1, class ITERATOR2>
Region(const ITERATOR1 & start,const ITERATOR2 & end)392     inline Region(const ITERATOR1& start, const ITERATOR2& end) :
393         mySize(0),
394         geometryCacheTainted(false)
395     {
396         load(start, end);
397     }
398 
399     template<class ITERATOR1, class ITERATOR2>
load(const ITERATOR1 & start,const ITERATOR2 & end)400     inline void load(const ITERATOR1& start, const ITERATOR2& end)
401     {
402         for (ITERATOR1 i = start; i != end; ++i) {
403             *this << *i;
404         }
405     }
406 
clear()407     inline void clear()
408     {
409         for (int i = 0; i < DIM; ++i) {
410             indices[i].clear();
411         }
412         mySize = 0;
413         myBoundingBox = CoordBox<DIM>();
414         geometryCacheTainted = false;
415     }
416 
boundingBox()417     inline const CoordBox<DIM>& boundingBox() const
418     {
419         if (geometryCacheTainted) {
420             resetGeometryCache();
421         }
422         return myBoundingBox;
423     }
424 
size()425     inline const std::size_t& size() const
426     {
427         if (geometryCacheTainted) {
428             resetGeometryCache();
429         }
430         return mySize;
431     }
432 
numStreaks()433     inline std::size_t numStreaks() const
434     {
435         return indices[0].size();
436     }
437 
438     inline Region expand(const unsigned& width=1) const
439     {
440         Region ret;
441         Coord<DIM> dia = Coord<DIM>::diagonal(width);
442 
443         for (StreakIterator i = beginStreak(); i != endStreak(); ++i) {
444             Streak<DIM> streak = *i;
445 
446             Coord<DIM> boxOrigin = streak.origin - dia;
447             Coord<DIM> boxDim = Coord<DIM>::diagonal(2 * width + 1);
448             boxDim.x() += streak.length() - 1;
449 
450             CoordBox<DIM> box(boxOrigin, boxDim);
451             ret << box;
452         }
453 
454         return ret;
455     }
456 
457     /**
458      * does the same as expand, but will wrap overlap at edges
459      * correctly. The instance of the TOPOLOGY is actually unused, but
460      * without it g++ would complain...
461      */
462     template<typename TOPOLOGY>
expandWithTopology(const unsigned & width,const Coord<DIM> & dimensions,TOPOLOGY)463     inline Region expandWithTopology(
464         const unsigned& width,
465         const Coord<DIM>& dimensions,
466         TOPOLOGY /* unused */) const
467     {
468         Region ret;
469         Coord<DIM> dia = Coord<DIM>::diagonal(width);
470 
471         for (StreakIterator i = beginStreak(); i != endStreak(); ++i) {
472             Streak<DIM> streak = *i;
473 
474             Coord<DIM> boxOrigin = streak.origin - dia;
475             Coord<DIM> boxDim = Coord<DIM>::diagonal(2 * width + 1);
476             boxDim.x() += streak.length() - 1;
477 
478             CoordBox<DIM> box(boxOrigin, boxDim);
479 
480             for (typename CoordBox<DIM>::StreakIterator i = box.beginStreak(); i != box.endStreak(); ++i) {
481                 Streak<DIM> newStreak(*i);
482                 if (TOPOLOGY::template WrapsAxis<0>::VALUE) {
483                     splitStreak<TOPOLOGY>(newStreak, &ret, dimensions);
484                 } else {
485                     normalizeStreak<TOPOLOGY>(
486                         trimStreak(newStreak, dimensions), &ret, dimensions);
487                 }
488             }
489         }
490 
491         return ret;
492     }
493 
494     inline bool operator==(const Region<DIM>& other) const
495     {
496         for (int i = 0; i < DIM; ++i) {
497             if (indices[i] != other.indices[i]) {
498                 return false;
499             }
500         }
501 
502         return true;
503     }
504 
505     /**
506      * Checks whether the Region includes the given Streak.
507      */
count(const Streak<DIM> & s)508     bool count(const Streak<DIM>& s) const
509     {
510         return RegionHelpers::RegionLookupHelper<DIM - 1>()(*this, s);
511     }
512 
513     /**
514      * Is the Coord contained withing the Region?
515      */
count(const Coord<DIM> & c)516     bool count(const Coord<DIM>& c) const
517     {
518         return RegionHelpers::RegionLookupHelper<DIM - 1>()(*this, Streak<DIM>(c, c[0] + 1));
519     }
520 
521     /**
522      * Alias for operator<<
523      */
524     template<typename ADDEND>
insert(const ADDEND & a)525     inline void insert(const ADDEND& a)
526     {
527         *this << a;
528     }
529 
530     /**
531      * Add all coordinates of the Streak to this Region
532      */
533     inline Region& operator<<(const Streak<DIM>& s)
534     {
535         //ignore 0 length streaks
536         if (s.endX <= s.origin.x()) {
537             return *this;
538         }
539 
540         geometryCacheTainted = true;
541         RegionHelpers::RegionInsertHelper<DIM - 1>()(this, s);
542         return *this;
543     }
544 
545     inline Region& operator<<(const Coord<DIM>& c)
546     {
547         *this << Streak<DIM>(c, c.x() + 1);
548         return *this;
549     }
550 
551     inline Region& operator<<(const CoordBox<DIM>& box)
552     {
553         for (typename CoordBox<DIM>::StreakIterator i = box.beginStreak(); i != box.endStreak(); ++i) {
554             *this << *i;
555         }
556 
557         return *this;
558     }
559 
560     /**
561      * Remove the given Streak (or all of its coordinates) from the
562      * Region. Takes into account that not all (or none) of the
563      * coordinates may have been contained in the Region before.
564      */
565     inline Region& operator>>(const Streak<DIM>& s)
566     {
567         //ignore 0 length streaks and empty selves
568         if (s.endX <= s.origin.x() || empty()) {
569             return *this;
570         }
571 
572         geometryCacheTainted = true;
573         RegionHelpers::RegionRemoveHelper<DIM - 1>()(this, s);
574         return *this;
575     }
576 
577     inline Region& operator>>(const Coord<DIM>& c)
578     {
579         *this >> Streak<DIM>(c, c.x() + 1);
580         return *this;
581     }
582 
583     inline void operator-=(const Region& other)
584     {
585         if(other.empty()) return;
586         for (StreakIterator i = other.beginStreak(); i != other.endStreak(); ++i) {
587             *this >> *i;
588         }
589     }
590 
591     /**
592      * Equvalent to (A and (not B)) in sets, where other corresponds
593      * to B and *this corresponds to A.
594      */
595     inline Region operator-(const Region& other) const
596     {
597         Region ret(*this);
598         ret -= other;
599         return ret;
600     }
601 
602     inline void operator&=(const Region& other)
603     {
604         Region intersection = other & *this;
605         *this = intersection;
606     }
607 
608     /**
609      * Computes the intersection of both regions.
610      */
611     inline Region operator&(const Region& other) const
612     {
613         Region ret;
614         StreakIterator myIter = beginStreak();
615         StreakIterator otherIter = other.beginStreak();
616 
617         StreakIterator myEnd = endStreak();
618         StreakIterator otherEnd = other.endStreak();
619 
620         for (;;) {
621             if ((myIter == myEnd) ||
622                 (otherIter == otherEnd)) {
623                 break;
624             }
625 
626             if (RegionHelpers::RegionIntersectHelper<DIM - 1>::intersects(*myIter, *otherIter)) {
627                 Streak<DIM> intersection = *myIter;
628                 intersection.origin.x() = (std::max)(myIter->origin.x(), otherIter->origin.x());
629                 intersection.endX = (std::min)(myIter->endX, otherIter->endX);
630                 ret << intersection;
631             }
632 
633             if (RegionHelpers::RegionIntersectHelper<DIM - 1>::lessThan(*myIter, *otherIter)) {
634                 ++myIter;
635             } else {
636                 ++otherIter;
637             }
638         }
639 
640         return ret;
641     }
642 
643     inline void operator+=(const Region& other)
644     {
645         if(other.empty()) return;
646         for (StreakIterator i = other.beginStreak(); i != other.endStreak(); ++i) {
647             *this << *i;
648         }
649     }
650 
651     inline Region operator+(const Region& other) const
652     {
653         Region ret(*this);
654         ret += other;
655         return ret;
656     }
657 
toVector()658     inline std::vector<Streak<DIM> > toVector() const
659     {
660         std::vector<Streak<DIM> > ret(numStreaks());
661         std::copy(beginStreak(), endStreak(), ret.begin());
662         return ret;
663     }
664 
toString()665     inline std::string toString() const
666     {
667         std::ostringstream buf;
668         buf << "Region<" << DIM << ">(\n";
669         for (int dim = 0; dim < DIM; ++dim) {
670             buf << "  indices[" << dim << "] = "
671                 << indices[dim] << "\n";
672         }
673         buf << ")\n";
674 
675         return buf.str();
676 
677     }
678 
prettyPrint()679     inline std::string prettyPrint() const
680     {
681         std::ostringstream buf;
682         buf << "Region<" << DIM << ">(\n";
683 
684         for (StreakIterator i = beginStreak(); i != endStreak(); ++i) {
685             buf << "  " << *i << "\n";
686         }
687 
688         buf << ")\n";
689 
690         return buf.str();
691     }
692 
empty()693     inline bool empty() const
694     {
695         return (indices[0].size() == 0);
696     }
697 
beginStreak()698     inline StreakIterator beginStreak() const
699     {
700         return StreakIterator(this, RegionHelpers::StreakIteratorInitBegin<DIM - 1>());
701     }
702 
endStreak()703     inline StreakIterator endStreak() const
704     {
705         return StreakIterator(this, RegionHelpers::StreakIteratorInitEnd<DIM - 1>());
706     }
707 
708     /**
709      * Returns an iterator whose internal iterators are set to the
710      * given offsets from the corresponding array starts. Runs in O(1)
711      * time.
712      */
713     inline StreakIterator operator[](const Coord<DIM>& offsets) const
714     {
715         return StreakIterator(this, RegionHelpers::StreakIteratorInitOffsets<DIM - 1, DIM>(offsets));
716     }
717 
718     /**
719      * Yields an iterator to the offset'th Streak in the Region. Runs
720      * in O(log n) time.
721      */
722     inline StreakIterator operator[](std::size_t offset) const
723     {
724         if (offset == 0) {
725             return beginStreak();
726         }
727         if (offset >= numStreaks()) {
728             return endStreak();
729         }
730         return StreakIterator(this, RegionHelpers::StreakIteratorInitSingleOffsetWrapper<DIM - 1>(offset));
731     }
732 
begin()733     inline Iterator begin() const
734     {
735         return Iterator(beginStreak());
736     }
737 
end()738     inline Iterator end() const
739     {
740         return Iterator(endStreak());
741     }
742 
indicesSize(std::size_t dim)743     inline std::size_t indicesSize(std::size_t dim) const
744     {
745         return indices[dim].size();
746     }
747 
indicesAt(std::size_t dim,std::size_t offset)748     inline VecType::const_iterator indicesAt(std::size_t dim, std::size_t offset) const
749     {
750         return indices[dim].begin() + offset;
751     }
752 
indicesBegin(std::size_t dim)753     inline VecType::const_iterator indicesBegin(std::size_t dim) const
754     {
755         return indices[dim].begin();
756     }
757 
indicesEnd(std::size_t dim)758     inline VecType::const_iterator indicesEnd(std::size_t dim) const
759     {
760         return indices[dim].end();
761     }
762 
763 private:
764     VecType indices[DIM];
765     mutable CoordBox<DIM> myBoundingBox;
766     mutable std::size_t mySize;
767     mutable bool geometryCacheTainted;
768 
determineGeometry()769     inline void determineGeometry() const
770     {
771         if (empty()) {
772             myBoundingBox = CoordBox<DIM>();
773         } else {
774             Streak<DIM> someStreak = *beginStreak();
775             Coord<DIM> min = someStreak.origin;
776             Coord<DIM> max = someStreak.origin;
777 
778             mySize = 0;
779             for (StreakIterator i = beginStreak();
780                  i != endStreak(); ++i) {
781                 Coord<DIM> left = i->origin;
782                 Coord<DIM> right = i->origin;
783                 right.x() = i->endX - 1;
784 
785                 min = (min.min)(left);
786                 max = (max.max)(right);
787                 mySize += i->endX - i->origin.x();
788             }
789 
790             myBoundingBox =
791                 CoordBox<DIM>(min, max - min + Coord<DIM>::diagonal(1));
792         }
793     }
794 
resetGeometryCache()795     inline void resetGeometryCache() const
796     {
797         determineGeometry();
798         geometryCacheTainted = false;
799     }
800 
trimStreak(const Streak<DIM> & s,const Coord<DIM> & dimensions)801     inline Streak<DIM> trimStreak(
802         const Streak<DIM>& s,
803         const Coord<DIM>& dimensions) const
804     {
805         int width = dimensions.x();
806         Streak<DIM> buf = s;
807         buf.origin.x() = (std::max)(buf.origin.x(), 0);
808         buf.endX = (std::min)(width, buf.endX);
809         return buf;
810     }
811 
812     template<typename TOPOLOGY>
splitStreak(const Streak<DIM> & streak,Region * target,const Coord<DIM> & dimensions)813     void splitStreak(
814         const Streak<DIM>& streak,
815         Region *target,
816         const Coord<DIM>& dimensions) const
817     {
818         int width = dimensions.x();
819 
820         int currentX = streak.origin.x();
821         if (currentX < 0) {
822             Streak<DIM> section = streak;
823             section.endX = (std::min)(streak.endX, 0);
824             currentX = section.endX;
825 
826             // normalize left overhang
827             section.origin.x() += width;
828             section.endX += width;
829             normalizeStreak<TOPOLOGY>(section, target, dimensions);
830         }
831 
832         if (currentX < streak.endX) {
833             Streak<DIM> section = streak;
834             section.origin.x() = currentX;
835             section.endX = (std::min)(streak.endX, width);
836             currentX = section.endX;
837 
838             normalizeStreak<TOPOLOGY>(section, target, dimensions);
839         }
840 
841         if (currentX < streak.endX) {
842             Streak<DIM> section = streak;
843             section.origin.x() = currentX;
844 
845             // normalize right overhang
846             section.origin.x() -= width;
847             section.endX -= width;
848             normalizeStreak<TOPOLOGY>(section, target, dimensions);
849         }
850     }
851 
852     template<typename TOPOLOGY>
normalizeStreak(const Streak<DIM> & streak,Region * target,const Coord<DIM> & dimensions)853     void normalizeStreak(
854         const Streak<DIM>& streak,
855         Region *target,
856         const Coord<DIM>& dimensions) const
857     {
858         Streak<DIM> ret;
859         ret.origin = TOPOLOGY::normalize(streak.origin, dimensions);
860         ret.endX = ret.origin.x() + streak.length();
861 
862         // it's bad to use a magic value to check for out of bounds
863         // accesses, but throwing exceptions would be slower
864         if (ret.origin != Coord<DIM>::diagonal(-1)) {
865             (*target) << ret;
866         }
867     }
868 };
869 
870 namespace RegionHelpers {
871 
872 template<int DIM>
873 class RegionLookupHelper : public RegionCommonHelper
874 {
875 public:
876     typedef Region<1>::IntPair IntPair;
877     typedef Region<1>::VecType VecType;
878 
879     template<int MY_DIM>
operator()880     inline bool operator()(const Region<MY_DIM>& region, const Streak<MY_DIM>& s)
881     {
882         const VecType& indices = region.indices[DIM];
883         return (*this)(region, s, 0, indices.size());
884     }
885 
886     template<int MY_DIM>
operator()887     inline bool operator()(const Region<MY_DIM>& region, const Streak<MY_DIM>& s, const int& start, const int& end)
888     {
889         int c = s.origin[DIM];
890         const VecType& indices = region.indices[DIM];
891 
892         VecType::const_iterator i =
893             std::upper_bound(
894                 indices.begin() + start,
895                 indices.begin() + end,
896                 IntPair(c, 0),
897                 RegionCommonHelper::pairCompareFirst);
898 
899         int nextLevelStart = 0;
900         int nextLevelEnd = 0;
901 
902         if (i != (indices.begin() + start)) {
903             VecType::const_iterator entry = i;
904             --entry;
905 
906             // recurse if found
907             if (entry->first == c) {
908                 nextLevelStart = entry->second;
909                 nextLevelEnd = region.indices[DIM - 1].size();
910                 if (i != indices.end()) {
911                     nextLevelEnd = i->second;
912                 }
913 
914                 return RegionLookupHelper<DIM-1>()(
915                     region,
916                     s,
917                     nextLevelStart,
918                     nextLevelEnd);
919             }
920         }
921 
922         return false;
923     }
924 };
925 
926 template<>
927 class RegionLookupHelper<0> : public RegionCommonHelper
928 {
929 public:
930     typedef Region<1>::IntPair IntPair;
931     typedef Region<1>::VecType VecType;
932 
933     template<int MY_DIM>
operator()934     inline bool operator()(const Region<MY_DIM>& region, const Streak<MY_DIM>& s, const int& start, int end)
935     {
936         IntPair curStreak(s.origin.x(), s.endX);
937         const VecType& indices = region.indices[0];
938 
939         VecType::const_iterator cursor =
940             std::upper_bound(indices.begin() + start, indices.begin() + end,
941                              curStreak, RegionCommonHelper::pairCompareFirst);
942         // This will yield the streak AFTER the current origin
943         // c. We can't really use lower_bound() as this doesn't
944         // replace the < operator by >= but rather by <=, which is
945         // IMO really sick...
946         if (cursor != (indices.begin() + start)) {
947             // ...so we revert to landing one past the streak we're
948             // searching and moving back afterwards:
949             cursor--;
950         }
951 
952         return (cursor->first <= s.origin[0]) && (cursor->second >= s.endX);
953     }
954 
955 };
956 
957 template<int DIM>
958 class RegionInsertHelper : public RegionCommonHelper
959 {
960 public:
961     typedef Region<1>::IntPair IntPair;
962     typedef Region<1>::VecType VecType;
963 
964     template<int MY_DIM>
operator()965     inline void operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s)
966     {
967         VecType& indices = region->indices[DIM];
968         (*this)(region, s, 0, indices.size());
969     }
970 
971     template<int MY_DIM>
operator()972     int operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s, const int& start, const int& end)
973     {
974         int c = s.origin[DIM];
975         VecType& indices = region->indices[DIM];
976 
977         VecType::iterator i =
978             std::upper_bound(
979                 indices.begin() + start,
980                 indices.begin() + end,
981                 IntPair(c, 0),
982                 RegionCommonHelper::pairCompareFirst);
983 
984         int nextLevelStart = 0;
985         int nextLevelEnd = 0;
986 
987         if (i != (indices.begin() + start)) {
988             VecType::iterator entry = i;
989             --entry;
990 
991             // short-cut: no need to insert if index already present
992             if (entry->first == c) {
993                 nextLevelStart = entry->second;
994                 nextLevelEnd = region->indices[DIM - 1].size();
995                 if (i != indices.end()) {
996                     nextLevelEnd = i->second;
997                 }
998 
999                 int inserts = RegionInsertHelper<DIM - 1>()(
1000                     region,
1001                     s,
1002                     nextLevelStart,
1003                     nextLevelEnd);
1004                 incRemainder(i, indices.end(), inserts);
1005                 return 0;
1006             }
1007         }
1008 
1009         if (i != indices.end()) {
1010             nextLevelStart = i->second;
1011         } else {
1012             nextLevelStart = region->indices[DIM - 1].size();
1013         }
1014 
1015         nextLevelEnd = nextLevelStart;
1016 
1017         VecType::iterator followingEntries;
1018 
1019         if (i == indices.end()) {
1020             indices << IntPair(c, nextLevelStart);
1021             followingEntries = indices.end();
1022         } else {
1023             followingEntries = indices.insert(i, IntPair(c, nextLevelStart));
1024             ++followingEntries;
1025         }
1026 
1027         int inserts = RegionInsertHelper<DIM - 1>()(region, s, nextLevelStart, nextLevelEnd);
1028         incRemainder(followingEntries, indices.end(), inserts);
1029 
1030         return 1;
1031     }
1032 };
1033 
1034 template<>
1035 class RegionInsertHelper<0>
1036 {
1037 public:
1038     friend class LibGeoDecomp::RegionTest;
1039     typedef Region<1>::IntPair IntPair;
1040     typedef Region<1>::VecType VecType;
1041 
1042     template<int MY_DIM>
operator()1043     inline int operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s, int start, int end)
1044     {
1045         IntPair curStreak(s.origin.x(), s.endX);
1046         VecType& indices = region->indices[0];
1047 
1048         VecType::iterator cursor =
1049             std::upper_bound(indices.begin() + start, indices.begin() + end,
1050                              curStreak, RegionCommonHelper::pairCompareFirst);
1051         // This will yield the streak AFTER the current origin
1052         // c. We can't really use lower_bound() as this doesn't
1053         // replace the < operator by >= but rather by <=, which is
1054         // IMO really sick...
1055         if (cursor != (indices.begin() + start)) {
1056             // ...so we revert to landing one past the streak we're
1057             // searching and moving back afterwards:
1058             cursor--;
1059         }
1060 
1061         int inserts = 1;
1062 
1063         while ((cursor != (indices.begin() + end)) &&
1064                (curStreak.second >= cursor->first)) {
1065             if (intersectOrTouch(*cursor, curStreak)) {
1066                 curStreak = fuse(*cursor, curStreak);
1067                 cursor = indices.erase(cursor);
1068                 --end;
1069                 --inserts;
1070             } else {
1071                 cursor++;
1072             }
1073 
1074             if ((cursor == (indices.begin() + end)) ||
1075                 (!intersectOrTouch(*cursor, curStreak))) {
1076                 break;
1077             }
1078         }
1079 
1080         indices.insert(cursor, curStreak);
1081         return inserts;
1082     }
1083 
1084 private:
intersectOrTouch(const IntPair & a,const IntPair & b)1085     inline bool intersectOrTouch(const IntPair& a, const IntPair& b) const
1086     {
1087         return
1088             ((a.first <= b.first && b.first <= a.second) ||
1089              (b.first <= a.first && a.first <= b.second));
1090     }
1091 
fuse(const IntPair & a,const IntPair & b)1092     inline IntPair fuse(const IntPair& a, const IntPair& b) const
1093     {
1094         return IntPair((std::min)(a.first, b.first),
1095                        (std::max)(a.second, b.second));
1096     }
1097 };
1098 
1099 template<int DIM>
1100 class RegionRemoveHelper : public RegionCommonHelper
1101 {
1102 public:
1103     typedef Region<1>::IntPair IntPair;
1104     typedef Region<1>::VecType VecType;
1105 
1106     template<int MY_DIM>
operator()1107     inline void operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s)
1108     {
1109         VecType& indices = region->indices[DIM];
1110         (*this)(region, s, 0, indices.size());
1111     }
1112 
1113     /**
1114      * tries to remove a streak from the set. Returns the number of
1115      * inserted streaks (may be negative).
1116      */
1117     template<int MY_DIM>
operator()1118     int operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s, const int& start, const int& end)
1119     {
1120         int c = s.origin[DIM];
1121         VecType& indices = region->indices[DIM];
1122 
1123         VecType::iterator i =
1124             std::upper_bound(
1125                 indices.begin() + start,
1126                 indices.begin() + end,
1127                 IntPair(c, 0),
1128                 RegionCommonHelper::pairCompareFirst);
1129 
1130         // key is not present, so no need to remove it
1131         if (i == (indices.begin() + start)) {
1132             return 0;
1133         }
1134 
1135         VecType::iterator entry = i;
1136         --entry;
1137 
1138         // ditto
1139         if (entry->first != c) {
1140             return 0;
1141         }
1142 
1143         int nextLevelStart = entry->second;
1144         int nextLevelEnd = region->indices[DIM - 1].size();
1145         if (i != indices.end()) {
1146             nextLevelEnd = i->second;
1147         }
1148 
1149         int inserts = RegionRemoveHelper<DIM - 1>()(
1150             region,
1151             s,
1152             nextLevelStart,
1153             nextLevelEnd);
1154 
1155         int myInserts = 0;
1156 
1157         // current entry needs to be removed if no childs are left
1158         if ((nextLevelStart - nextLevelEnd) == inserts) {
1159             entry = indices.erase(entry);
1160             myInserts = -1;
1161         } else {
1162             ++entry;
1163         }
1164 
1165         incRemainder(entry, indices.end(), inserts);
1166         return myInserts;
1167     }
1168 };
1169 
1170 template<>
1171 class RegionRemoveHelper<0>
1172 {
1173 public:
1174     friend class LibGeoDecomp::RegionTest;
1175     typedef Region<1>::IntPair IntPair;
1176     typedef Region<1>::VecType VecType;
1177 
1178     template<int MY_DIM>
operator()1179     int operator()(Region<MY_DIM> *region, const Streak<MY_DIM>& s, const int& start, int end)
1180     {
1181         int c = s.origin[0];
1182         VecType& indices = region->indices[0];
1183         int inserts = 0;
1184 
1185         // This will yield the streak AFTER the current origin
1186         // c. We can't really use lower_bound() as this doesn't
1187         // replace the < operator by >= but rather by <=, which is
1188         // IMO really sick...
1189         VecType::iterator cursor =
1190             std::upper_bound(
1191                 indices.begin() + start,
1192                 indices.begin() + end,
1193                 IntPair(c, 0),
1194                 RegionCommonHelper::pairCompareFirst);
1195         if (cursor != (indices.begin() + start)) {
1196             // ...so we resort to landing one past the streak we're
1197             // searching and moving back afterwards:
1198             --cursor;
1199         }
1200 
1201         IntPair curStreak(s.origin.x(), s.endX);
1202 
1203         while (cursor != (indices.begin() + end)) {
1204             if (intersect(curStreak, *cursor)) {
1205                 VecType newStreaks(substract(*cursor, curStreak));
1206                 cursor = indices.erase(cursor);
1207                 int delta = newStreaks.size() - 1;
1208                 end += delta;
1209                 inserts += delta;
1210 
1211                 for (VecType::iterator i = newStreaks.begin(); i != newStreaks.end(); ++i) {
1212                     cursor = indices.insert(cursor, *i);
1213                     ++cursor;
1214                 }
1215             } else {
1216                 ++cursor;
1217             }
1218 
1219             if (cursor == (indices.begin() + end) || !intersect(*cursor, curStreak)) {
1220                 break;
1221             }
1222         }
1223 
1224         return inserts;
1225     }
1226 
1227 private:
intersect(const IntPair & a,const IntPair & b)1228     inline bool intersect(const IntPair& a, const IntPair& b) const
1229     {
1230         return
1231             ((a.first <= b.first && b.first < a.second) ||
1232              (b.first <= a.first && a.first < b.second));
1233     }
1234 
substract(const IntPair & base,const IntPair & minuend)1235     inline VecType substract(const IntPair& base, const IntPair& minuend) const
1236     {
1237         if (!intersect(base, minuend)) {
1238             return std::vector<IntPair>(1, base);
1239         }
1240 
1241         std::vector<IntPair> ret;
1242         IntPair s1(base.first, minuend.first);
1243         IntPair s2(minuend.second, base.second);
1244 
1245         if (s1.second > s1.first) {
1246             ret.push_back(s1);
1247         }
1248         if (s2.second > s2.first) {
1249             ret.push_back(s2);
1250         }
1251         return ret;
1252     }
1253 };
1254 
1255 }
1256 
1257 template<typename _CharT, typename _Traits, int _Dim>
1258 std::basic_ostream<_CharT, _Traits>&
1259 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1260            const LibGeoDecomp::Region<_Dim>& region)
1261 {
1262     __os << region.toString();
1263     return __os;
1264 }
1265 
1266 }
1267 
1268 #endif
1269