1 /*
2     Copyright (c) 2005-2021 Intel Corporation
3 
4     Licensed under the Apache License, Version 2.0 (the "License");
5     you may not use this file except in compliance with the License.
6     You may obtain a copy of the License at
7 
8         http://www.apache.org/licenses/LICENSE-2.0
9 
10     Unless required by applicable law or agreed to in writing, software
11     distributed under the License is distributed on an "AS IS" BASIS,
12     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13     See the License for the specific language governing permissions and
14     limitations under the License.
15 */
16 
17 #ifndef NOMINMAX
18 #define NOMINMAX
19 #endif // NOMINMAX
20 
21 // Polygon overlay
22 //
23 #include <cstdlib>
24 #include <cstring>
25 #include <cassert>
26 
27 #include <iostream>
28 #include <algorithm>
29 
30 #include "oneapi/tbb/tick_count.h"
31 #include "oneapi/tbb/blocked_range.h"
32 #include "oneapi/tbb/parallel_for.h"
33 #include "oneapi/tbb/spin_mutex.h"
34 #include "oneapi/tbb/global_control.h"
35 
36 #include "common/utility/get_default_num_threads.hpp"
37 
38 #include "polyover.hpp"
39 #include "polymain.hpp"
40 #include "pover_video.hpp"
41 
42 /*!
43 * @brief intersects a polygon with a map, adding any results to output map
44 *
45 * @param[out] resultMap output map (must be allocated)
46 * @param[in] polygon to be intersected
47 * @param[in] map intersected against
48 * @param[in] lock to use when adding output polygons to result map
49 *
50 */
51 void OverlayOnePolygonWithMap(Polygon_map_t *resultMap,
52                               RPolygon *myPoly,
53                               Polygon_map_t *map2,
54                               oneapi::tbb::spin_mutex *rMutex) {
55     int r1, g1, b1, r2, g2, b2;
56     int myr = 0;
57     int myg = 0;
58     int myb = 0;
59     int p1Area = myPoly->area();
60     for (unsigned int j = 1; (j < map2->size()) && (p1Area > 0); j++) {
61         RPolygon *p2 = &((*map2)[j]);
62         RPolygon *pnew;
63         int newxMin, newxMax, newyMin, newyMax;
64         myPoly->getColor(&r1, &g1, &b1);
65         if (PolygonsOverlap(myPoly, p2, newxMin, newyMin, newxMax, newyMax)) {
66             p2->getColor(&r2, &g2, &b2);
67             myr = r1 + r2;
68             myg = g1 + g2;
69             myb = b1 + b2;
70             p1Area -= (newxMax - newxMin + 1) * (newyMax - newyMin + 1);
71             if (rMutex) {
72                 oneapi::tbb::spin_mutex::scoped_lock lock(*rMutex);
73                 resultMap->push_back(RPolygon(newxMin, newyMin, newxMax, newyMax, myr, myg, myb));
74             }
75             else {
76                 resultMap->push_back(RPolygon(newxMin, newyMin, newxMax, newyMax, myr, myg, myb));
77             }
78         }
79     }
80 }
81 
82 /*!
83 * @brief Serial version of polygon overlay
84 * @param[out] output map
85 * @param[in]  first map (map that individual polygons are taken from)
86 * @param[in]  second map (map passed to OverlayOnePolygonWithMap)
87 */
88 void SerialOverlayMaps(Polygon_map_t **resultMap, Polygon_map_t *map1, Polygon_map_t *map2) {
89     std::cout << "SerialOverlayMaps called"
90               << "\n";
91     *resultMap = new Polygon_map_t;
92 
93     RPolygon *p0 = &((*map1)[0]);
94     int mapxSize, mapySize, ignore1, ignore2;
95     p0->get(&ignore1, &ignore2, &mapxSize, &mapySize);
96     (*resultMap)->reserve(mapxSize * mapySize); // can't be any bigger than this
97     // push the map size as the first polygon,
98     (*resultMap)->push_back(RPolygon(0, 0, mapxSize, mapySize));
99     for (unsigned int i = 1; i < map1->size(); i++) {
100         RPolygon *p1 = &((*map1)[i]);
101         OverlayOnePolygonWithMap(*resultMap, p1, map2, nullptr);
102     }
103 }
104 
105 /*!
106 * @class ApplyOverlay
107 * @brief Simple version of parallel overlay (make parallel on polygons in map1)
108 */
109 class ApplyOverlay {
110     Polygon_map_t *m_map1, *m_map2, *m_resultMap;
111     oneapi::tbb::spin_mutex *m_rMutex;
112 
113 public:
114     /*!
115     * @brief functor to apply
116     * @param[in] r range of polygons to intersect from map1
117     */
118     void operator()(const oneapi::tbb::blocked_range<int> &r) const {
119         PRINT_DEBUG("From " << r.begin() << " to " << r.end());
120         for (int i = r.begin(); i != r.end(); i++) {
121             RPolygon *myPoly = &((*m_map1)[i]);
122             OverlayOnePolygonWithMap(m_resultMap, myPoly, m_map2, m_rMutex);
123         }
124     }
125     ApplyOverlay(Polygon_map_t *resultMap,
126                  Polygon_map_t *map1,
127                  Polygon_map_t *map2,
128                  oneapi::tbb::spin_mutex *rmutex)
129             : m_resultMap(resultMap),
130               m_map1(map1),
131               m_map2(map2),
132               m_rMutex(rmutex) {}
133 };
134 
135 /*!
136 * @brief apply the parallel algorithm
137 * @param[out] result_map generated map
138 * @param[in] polymap1 first map to be applied (algorithm is parallel on this map)
139 * @param[in] polymap2 second map.
140 */
141 void NaiveParallelOverlay(Polygon_map_t *&result_map,
142                           Polygon_map_t &polymap1,
143                           Polygon_map_t &polymap2) {
144     // -----------------------------------
145     bool automatic_threadcount = false;
146 
147     if (gThreadsLow == THREADS_UNSET || gThreadsLow == utility::get_default_num_threads()) {
148         gThreadsLow = gThreadsHigh = utility::get_default_num_threads();
149         automatic_threadcount = true;
150     }
151     result_map = new Polygon_map_t;
152 
153     RPolygon *p0 = &(polymap1[0]);
154     int mapxSize, mapySize, ignore1, ignore2;
155     p0->get(&ignore1, &ignore2, &mapxSize, &mapySize);
156     result_map->reserve(mapxSize * mapySize); // can't be any bigger than this
157     // push the map size as the first polygon,
158     oneapi::tbb::spin_mutex *resultMutex = new oneapi::tbb::spin_mutex();
159     int grain_size = gGrainSize;
160 
161     for (int nthreads = gThreadsLow; nthreads <= gThreadsHigh; nthreads++) {
162         oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
163                                       nthreads);
164         if (gIsGraphicalVersion) {
165             RPolygon *xp =
166                 new RPolygon(0, 0, gMapXSize - 1, gMapYSize - 1, 0, 0, 0); // Clear the output space
167             delete xp;
168         }
169         // put size polygon in result map
170         result_map->push_back(RPolygon(0, 0, mapxSize, mapySize));
171 
172         oneapi::tbb::tick_count t0 = oneapi::tbb::tick_count::now();
173         oneapi::tbb::parallel_for(
174             oneapi::tbb::blocked_range<int>(1, (int)(polymap1.size()), grain_size),
175             ApplyOverlay(result_map, &polymap1, &polymap2, resultMutex));
176         oneapi::tbb::tick_count t1 = oneapi::tbb::tick_count::now();
177 
178         double naiveParallelTime = (t1 - t0).seconds() * 1000;
179         std::cout << "Naive parallel with spin lock and ";
180         if (automatic_threadcount)
181             std::cout << "automatic";
182         else
183             std::cout << nthreads;
184         std::cout << ((nthreads == 1) ? " thread" : " threads");
185         std::cout << " took " << naiveParallelTime << " msec : speedup over serial "
186                   << (gSerialTime / naiveParallelTime) << "\n";
187         if (gCsvFile.is_open()) {
188             gCsvFile << "," << naiveParallelTime;
189         }
190 #if _DEBUG
191         CheckPolygonMap(result_map);
192         ComparePolygonMaps(result_map, gResultMap);
193 #endif
194         result_map->clear();
195     }
196     delete resultMutex;
197     if (gCsvFile.is_open()) {
198         gCsvFile << "\n";
199     }
200     // -----------------------------------
201 }
202 
203 template <typename T>
204 void split_at(Flagged_map_t &in_map,
205               Flagged_map_t &left_out,
206               Flagged_map_t &right_out,
207               const T median) {
208     left_out.reserve(in_map.size());
209     right_out.reserve(in_map.size());
210     for (Flagged_map_t::iterator i = in_map.begin(); i != in_map.end(); ++i) {
211         RPolygon *p = i->p();
212         if (p->xmax() < median) {
213             // in left map
214             left_out.push_back(*i);
215         }
216         else if (p->xmin() >= median) {
217             right_out.push_back(*i);
218             // in right map
219         }
220         else {
221             // in both maps.
222             left_out.push_back(*i);
223             right_out.push_back(RPolygon_flagged(p, true));
224         }
225     }
226 }
227 
228 // range that splits the maps as well as the range.  the flagged_map_t are
229 // vectors of pointers, and each range owns its maps (has to free them on destruction.)
230 template <typename T>
231 class blocked_range_with_maps {
232     typedef oneapi::tbb::blocked_range<T> my_range_type;
233 
234 private:
235     my_range_type my_range;
236     Flagged_map_t my_map1;
237     Flagged_map_t my_map2;
238 
239 public:
240     blocked_range_with_maps(T begin,
241                             T end,
242                             typename my_range_type::size_type my_grainsize,
243                             Polygon_map_t *p1,
244                             Polygon_map_t *p2)
245             : my_range(begin, end, my_grainsize) {
246         my_map1.reserve(p1->size());
247         my_map2.reserve(p2->size());
248         for (int i = 1; i < p1->size(); ++i) {
249             my_map1.push_back(RPolygon_flagged(&((*p1)[i]), false));
250         }
251         for (int i = 1; i < p2->size(); ++i) {
252             my_map2.push_back(RPolygon_flagged(&(p2->at(i)), false));
253         }
254     }
255 
256     // copy-constructor required for deep copy of flagged maps.  One copy is done at the start of the
257     // parallel for.
258     blocked_range_with_maps(const blocked_range_with_maps &other)
259             : my_range(other.my_range),
260               my_map1(other.my_map1),
261               my_map2(other.my_map2) {}
262     bool empty() const {
263         return my_range.empty();
264     }
265     bool is_divisible() const {
266         return my_range.is_divisible();
267     }
268 
269 #if _DEBUG
270     void check_my_map() {
271         assert(my_range.begin() <= my_range.end());
272         for (Flagged_map_t::iterator i = my_map1.begin(); i != my_map1.end(); ++i) {
273             RPolygon *rp = i->p();
274             assert(rp->xmax() >= my_range.begin());
275             assert(rp->xmin() < my_range.end());
276         }
277         for (Flagged_map_t::iterator i = my_map2.begin(); i != my_map2.end(); ++i) {
278             RPolygon *rp = i->p();
279             assert(rp->xmax() >= my_range.begin());
280             assert(rp->xmin() < my_range.end());
281         }
282     }
283 
284     void dump_map(Flagged_map_t &mapx) {
285         std::cout << " ** MAP **\n";
286         for (Flagged_map_t::iterator i = mapx.begin(); i != mapx.end(); ++i) {
287             std::cout << *(i->p());
288             if (i->isDuplicate()) {
289                 std::cout << " -- is_duplicate";
290             }
291             std::cout << "\n";
292         }
293         std::cout << "\n";
294     }
295 #endif
296 
297     blocked_range_with_maps(blocked_range_with_maps &lhs_r, oneapi::tbb::split)
298             : my_range(my_range_type(lhs_r.my_range, oneapi::tbb::split())) {
299         // lhs_r.my_range makes my_range from [median, high) and rhs_r.my_range from [low, median)
300         Flagged_map_t original_map1 = lhs_r.my_map1;
301         Flagged_map_t original_map2 = lhs_r.my_map2;
302         lhs_r.my_map1.clear();
303         lhs_r.my_map2.clear();
304         split_at(original_map1, lhs_r.my_map1, my_map1, my_range.begin());
305         split_at(original_map2, lhs_r.my_map2, my_map2, my_range.begin());
306 #if _DEBUG
307         this->check_my_map();
308         lhs_r.check_my_map();
309 #endif
310     }
311 
312     const my_range_type &range() const {
313         return my_range;
314     }
315     Flagged_map_t &map1() {
316         return my_map1;
317     }
318     Flagged_map_t &map2() {
319         return my_map2;
320     }
321 };
322 
323 /*!
324 * @class ApplySplitOverlay
325 * @brief parallel by columnar strip
326 */
327 class ApplySplitOverlay {
328     Polygon_map_t *m_map1, *m_map2, *m_resultMap;
329     oneapi::tbb::spin_mutex *m_rMutex;
330 
331 public:
332     /*!
333     * @brief functor for columnar parallel version
334     * @param[in] r range of map to be operated on
335     */
336     void operator()(/*const*/ blocked_range_with_maps<int> &r) const {
337 #ifdef _DEBUG
338         // if we are debugging, serialize the method.  That way we can
339         // see what is happening in each strip without the interleaving
340         // confusing things.
341         oneapi::tbb::spin_mutex::scoped_lock lock(*m_rMutex);
342         std::cout << std::unitbuf << "From " << r.range().begin() << " to " << r.range().end() - 1
343                   << "\n";
344 #endif
345         // get yMapSize
346         int r1, g1, b1, r2, g2, b2;
347         int myr = -1;
348         int myg = -1;
349         int myb = -1;
350         int i1, i2, i3, yMapSize;
351         (*m_map1)[0].get(&i1, &i2, &i3, &yMapSize);
352 
353         Flagged_map_t &fmap1 = r.map1();
354         Flagged_map_t &fmap2 = r.map2();
355 
356         // When intersecting polygons from fmap1 and fmap2, if BOTH are flagged
357         // as duplicate, don't add the result to the output map.  We can still
358         // intersect them, because we are keeping track of how much of the polygon
359         // is left over from intersecting, and quitting when the polygon is
360         // used up.
361 
362         for (unsigned int i = 0; i < fmap1.size(); i++) {
363             RPolygon *p1 = fmap1[i].p();
364             bool is_dup = fmap1[i].isDuplicate();
365             int parea = p1->area();
366             p1->getColor(&r1, &g1, &b1);
367             for (unsigned int j = 0; (j < fmap2.size()) && (parea > 0); j++) {
368                 int xl, yl, xh, yh;
369                 RPolygon *p2 = fmap2[j].p();
370                 if (PolygonsOverlap(p1, p2, xl, yl, xh, yh)) {
371                     if (!(is_dup && fmap2[j].isDuplicate())) {
372                         p2->getColor(&r2, &g2, &b2);
373                         myr = r1 + r2;
374                         myg = g1 + g2;
375                         myb = b1 + b2;
376 #ifdef _DEBUG
377 #else
378                         oneapi::tbb::spin_mutex::scoped_lock lock(*m_rMutex);
379 #endif
380                         (*m_resultMap).push_back(RPolygon(xl, yl, xh, yh, myr, myg, myb));
381                     }
382                     parea -= (xh - xl + 1) * (yh - yl + 1);
383                 }
384             }
385         }
386     }
387 
388     ApplySplitOverlay(Polygon_map_t *resultMap,
389                       Polygon_map_t *map1,
390                       Polygon_map_t *map2,
391                       oneapi::tbb::spin_mutex *rmutex)
392             : m_resultMap(resultMap),
393               m_map1(map1),
394               m_map2(map2),
395               m_rMutex(rmutex) {}
396 };
397 
398 /*!
399 * @brief intersects two maps strip-wise
400 *
401 * @param[out] resultMap output map (must be allocated)
402 * @param[in] polymap1 map to be intersected
403 * @param[in] polymap2 map to be intersected
404 */
405 void SplitParallelOverlay(Polygon_map_t **result_map,
406                           Polygon_map_t *polymap1,
407                           Polygon_map_t *polymap2) {
408     int nthreads;
409     bool automatic_threadcount = false;
410     double domainSplitParallelTime;
411     oneapi::tbb::tick_count t0, t1;
412     oneapi::tbb::spin_mutex *resultMutex;
413     if (gThreadsLow == THREADS_UNSET || gThreadsLow == utility::get_default_num_threads()) {
414         gThreadsLow = gThreadsHigh = utility::get_default_num_threads();
415         automatic_threadcount = true;
416     }
417     *result_map = new Polygon_map_t;
418 
419     RPolygon *p0 = &((*polymap1)[0]);
420     int mapxSize, mapySize, ignore1, ignore2;
421     p0->get(&ignore1, &ignore2, &mapxSize, &mapySize);
422     (*result_map)->reserve(mapxSize * mapySize); // can't be any bigger than this
423     resultMutex = new oneapi::tbb::spin_mutex();
424 
425     int grain_size;
426 #ifdef _DEBUG
427     grain_size = gMapXSize / 4;
428 #else
429     grain_size = gGrainSize;
430 #endif
431     for (nthreads = gThreadsLow; nthreads <= gThreadsHigh; nthreads++) {
432         oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
433                                       nthreads);
434         if (gIsGraphicalVersion) {
435             RPolygon *xp =
436                 new RPolygon(0, 0, gMapXSize - 1, gMapYSize - 1, 0, 0, 0); // Clear the output space
437             delete xp;
438         }
439         // push the map size as the first polygon,
440         (*result_map)->push_back(RPolygon(0, 0, mapxSize, mapySize));
441         t0 = oneapi::tbb::tick_count::now();
442         oneapi::tbb::parallel_for(
443             blocked_range_with_maps<int>(0, (int)(mapxSize + 1), grain_size, polymap1, polymap2),
444             ApplySplitOverlay((*result_map), polymap1, polymap2, resultMutex));
445         t1 = oneapi::tbb::tick_count::now();
446         domainSplitParallelTime = (t1 - t0).seconds() * 1000;
447         std::cout << "Splitting parallel with spin lock and ";
448         if (automatic_threadcount)
449             std::cout << "automatic";
450         else
451             std::cout << nthreads;
452         std::cout << ((nthreads == 1) ? " thread" : " threads");
453         std::cout << " took " << domainSplitParallelTime << " msec : speedup over serial "
454                   << (gSerialTime / domainSplitParallelTime) << "\n";
455         if (gCsvFile.is_open()) {
456             gCsvFile << "," << domainSplitParallelTime;
457         }
458 #if _DEBUG
459         CheckPolygonMap(*result_map);
460         ComparePolygonMaps(*result_map, gResultMap);
461 #endif
462         (*result_map)->clear();
463     }
464     delete resultMutex;
465     if (gCsvFile.is_open()) {
466         gCsvFile << "\n";
467     }
468 }
469 
470 class ApplySplitOverlayCV {
471     Polygon_map_t *m_map1, *m_map2;
472     concurrent_Polygon_map_t *m_resultMap;
473 
474 public:
475     /*!
476     * @brief functor for columnar parallel version
477     * @param[in] r range of map to be operated on
478     */
479     void operator()(blocked_range_with_maps<int> &r) const {
480         // get yMapSize
481         int r1, g1, b1, r2, g2, b2;
482         int myr = -1;
483         int myg = -1;
484         int myb = -1;
485         int i1, i2, i3, yMapSize;
486         (*m_map1)[0].get(&i1, &i2, &i3, &yMapSize);
487 
488         Flagged_map_t &fmap1 = r.map1();
489         Flagged_map_t &fmap2 = r.map2();
490 
491         // When intersecting polygons from fmap1 and fmap2, if BOTH are flagged
492         // as duplicate, don't add the result to the output map.  We can still
493         // intersect them, because we are keeping track of how much of the polygon
494         // is left over from intersecting, and quitting when the polygon is
495         // used up.
496 
497         for (unsigned int i = 0; i < fmap1.size(); i++) {
498             RPolygon *p1 = fmap1[i].p();
499             bool is_dup = fmap1[i].isDuplicate();
500             int parea = p1->area();
501             p1->getColor(&r1, &g1, &b1);
502             for (unsigned int j = 0; (j < fmap2.size()) && (parea > 0); j++) {
503                 int xl, yl, xh, yh;
504                 RPolygon *p2 = fmap2[j].p();
505                 if (PolygonsOverlap(p1, p2, xl, yl, xh, yh)) {
506                     if (!(is_dup && fmap2[j].isDuplicate())) {
507                         p2->getColor(&r2, &g2, &b2);
508                         myr = r1 + r2;
509                         myg = g1 + g2;
510                         myb = b1 + b2;
511                         (*m_resultMap).push_back(RPolygon(xl, yl, xh, yh, myr, myg, myb));
512                     }
513                     parea -= (xh - xl + 1) * (yh - yl + 1);
514                 }
515             }
516         }
517     }
518 
519     ApplySplitOverlayCV(concurrent_Polygon_map_t *resultMap,
520                         Polygon_map_t *map1,
521                         Polygon_map_t *map2)
522             : m_resultMap(resultMap),
523               m_map1(map1),
524               m_map2(map2) {}
525 };
526 
527 /*!
528 * @brief intersects two maps strip-wise, accumulating into a concurrent_vector
529 *
530 * @param[out] resultMap output map (must be allocated)
531 * @param[in] polymap1 map to be intersected
532 * @param[in] polymap2 map to be intersected
533 */
534 void SplitParallelOverlayCV(concurrent_Polygon_map_t **result_map,
535                             Polygon_map_t *polymap1,
536                             Polygon_map_t *polymap2) {
537     int nthreads;
538     bool automatic_threadcount = false;
539     double domainSplitParallelTime;
540     oneapi::tbb::tick_count t0, t1;
541     if (gThreadsLow == THREADS_UNSET || gThreadsLow == utility::get_default_num_threads()) {
542         gThreadsLow = gThreadsHigh = utility::get_default_num_threads();
543         automatic_threadcount = true;
544     }
545     *result_map = new concurrent_Polygon_map_t;
546 
547     RPolygon *p0 = &((*polymap1)[0]);
548     int mapxSize, mapySize, ignore1, ignore2;
549     p0->get(&ignore1, &ignore2, &mapxSize, &mapySize);
550     // (*result_map)->reserve(mapxSize*mapySize); // can't be any bigger than this
551 
552     int grain_size;
553 #ifdef _DEBUG
554     grain_size = gMapXSize / 4;
555 #else
556     grain_size = gGrainSize;
557 #endif
558     for (nthreads = gThreadsLow; nthreads <= gThreadsHigh; nthreads++) {
559         oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
560                                       nthreads);
561         if (gIsGraphicalVersion) {
562             RPolygon *xp =
563                 new RPolygon(0, 0, gMapXSize - 1, gMapYSize - 1, 0, 0, 0); // Clear the output space
564             delete xp;
565         }
566         // push the map size as the first polygon,
567         (*result_map)->push_back(RPolygon(0, 0, mapxSize, mapySize));
568         t0 = oneapi::tbb::tick_count::now();
569         oneapi::tbb::parallel_for(
570             blocked_range_with_maps<int>(0, (int)(mapxSize + 1), grain_size, polymap1, polymap2),
571             ApplySplitOverlayCV((*result_map), polymap1, polymap2));
572         t1 = oneapi::tbb::tick_count::now();
573         domainSplitParallelTime = (t1 - t0).seconds() * 1000;
574         std::cout << "Splitting parallel with concurrent_vector and ";
575         if (automatic_threadcount)
576             std::cout << "automatic";
577         else
578             std::cout << nthreads;
579         std::cout << ((nthreads == 1) ? " thread" : " threads");
580         std::cout << " took " << domainSplitParallelTime << " msec : speedup over serial "
581                   << (gSerialTime / domainSplitParallelTime) << "\n";
582         if (gCsvFile.is_open()) {
583             gCsvFile << "," << domainSplitParallelTime;
584         }
585 #if _DEBUG
586         {
587             Polygon_map_t s_result_map;
588             for (concurrent_Polygon_map_t::const_iterator i = (*result_map)->begin();
589                  i != (*result_map)->end();
590                  ++i) {
591                 s_result_map.push_back(*i);
592             }
593             CheckPolygonMap(&s_result_map);
594             ComparePolygonMaps(&s_result_map, gResultMap);
595         }
596 #endif
597         (*result_map)->clear();
598     }
599 
600     if (gCsvFile.is_open()) {
601         gCsvFile << "\n";
602     }
603 }
604 
605 // ------------------------------------------------------
606 
607 class ApplySplitOverlayETS {
608     Polygon_map_t *m_map1, *m_map2;
609     ETS_Polygon_map_t *m_resultMap;
610 
611 public:
612     /*!
613     * @brief functor for columnar parallel version
614     * @param[in] r range of map to be operated on
615     */
616     void operator()(blocked_range_with_maps<int> &r) const {
617         // get yMapSize
618         int r1, g1, b1, r2, g2, b2;
619         int myr = -1;
620         int myg = -1;
621         int myb = -1;
622         int i1, i2, i3, yMapSize;
623         (*m_map1)[0].get(&i1, &i2, &i3, &yMapSize);
624 
625         Flagged_map_t &fmap1 = r.map1();
626         Flagged_map_t &fmap2 = r.map2();
627 
628         // When intersecting polygons from fmap1 and fmap2, if BOTH are flagged
629         // as duplicate, don't add the result to the output map.  We can still
630         // intersect them, because we are keeping track of how much of the polygon
631         // is left over from intersecting, and quitting when the polygon is
632         // used up.
633 
634         for (unsigned int i = 0; i < fmap1.size(); i++) {
635             RPolygon *p1 = fmap1[i].p();
636             bool is_dup = fmap1[i].isDuplicate();
637             int parea = p1->area();
638             p1->getColor(&r1, &g1, &b1);
639             for (unsigned int j = 0; (j < fmap2.size()) && (parea > 0); j++) {
640                 int xl, yl, xh, yh;
641                 RPolygon *p2 = fmap2[j].p();
642                 if (PolygonsOverlap(p1, p2, xl, yl, xh, yh)) {
643                     if (!(is_dup && fmap2[j].isDuplicate())) {
644                         p2->getColor(&r2, &g2, &b2);
645                         myr = r1 + r2;
646                         myg = g1 + g2;
647                         myb = b1 + b2;
648                         (*m_resultMap).local().push_back(RPolygon(xl, yl, xh, yh, myr, myg, myb));
649                     }
650                     parea -= (xh - xl + 1) * (yh - yl + 1);
651                 }
652             }
653         }
654     }
655 
656     ApplySplitOverlayETS(ETS_Polygon_map_t *resultMap, Polygon_map_t *map1, Polygon_map_t *map2)
657             : m_resultMap(resultMap),
658               m_map1(map1),
659               m_map2(map2) {}
660 };
661 
662 /*!
663 * @brief intersects two maps strip-wise, accumulating into an ets variable
664 *
665 * @param[out] resultMap output map (must be allocated)
666 * @param[in] polymap1 map to be intersected
667 * @param[in] polymap2 map to be intersected
668 */
669 void SplitParallelOverlayETS(ETS_Polygon_map_t **result_map,
670                              Polygon_map_t *polymap1,
671                              Polygon_map_t *polymap2) {
672     int nthreads;
673     bool automatic_threadcount = false;
674     double domainSplitParallelTime;
675     oneapi::tbb::tick_count t0, t1;
676     if (gThreadsLow == THREADS_UNSET || gThreadsLow == utility::get_default_num_threads()) {
677         gThreadsLow = gThreadsHigh = utility::get_default_num_threads();
678         automatic_threadcount = true;
679     }
680     *result_map = new ETS_Polygon_map_t;
681 
682     RPolygon *p0 = &((*polymap1)[0]);
683     int mapxSize, mapySize, ignore1, ignore2;
684     p0->get(&ignore1, &ignore2, &mapxSize, &mapySize);
685     // (*result_map)->reserve(mapxSize*mapySize); // can't be any bigger than this
686 
687     int grain_size;
688 #ifdef _DEBUG
689     grain_size = gMapXSize / 4;
690 #else
691     grain_size = gGrainSize;
692 #endif
693     for (nthreads = gThreadsLow; nthreads <= gThreadsHigh; nthreads++) {
694         oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
695                                       nthreads);
696         if (gIsGraphicalVersion) {
697             RPolygon *xp =
698                 new RPolygon(0, 0, gMapXSize - 1, gMapYSize - 1, 0, 0, 0); // Clear the output space
699             delete xp;
700         }
701         // push the map size as the first polygon,
702         // This polygon needs to be first, so we can push it at the start of a combine.
703         // (*result_map)->local.push_back(RPolygon(0,0,mapxSize, mapySize));
704         t0 = oneapi::tbb::tick_count::now();
705         oneapi::tbb::parallel_for(
706             blocked_range_with_maps<int>(0, (int)(mapxSize + 1), grain_size, polymap1, polymap2),
707             ApplySplitOverlayETS((*result_map), polymap1, polymap2));
708         t1 = oneapi::tbb::tick_count::now();
709         domainSplitParallelTime = (t1 - t0).seconds() * 1000;
710         std::cout << "Splitting parallel with ETS and ";
711         if (automatic_threadcount)
712             std::cout << "automatic";
713         else
714             std::cout << nthreads;
715         std::cout << ((nthreads == 1) ? " thread" : " threads");
716         std::cout << " took " << domainSplitParallelTime << " msec : speedup over serial "
717                   << (gSerialTime / domainSplitParallelTime) << "\n";
718         if (gCsvFile.is_open()) {
719             gCsvFile << "," << domainSplitParallelTime;
720         }
721 #if _DEBUG
722         {
723             Polygon_map_t s_result_map;
724             oneapi::tbb::flattened2d<ETS_Polygon_map_t> psv = flatten2d(**result_map);
725             s_result_map.push_back(RPolygon(0, 0, mapxSize, mapySize));
726             for (oneapi::tbb::flattened2d<ETS_Polygon_map_t>::const_iterator ci = psv.begin();
727                  ci != psv.end();
728                  ++ci) {
729                 s_result_map.push_back(*ci);
730             }
731             CheckPolygonMap(&s_result_map);
732             ComparePolygonMaps(&s_result_map, gResultMap);
733         }
734 #endif
735         (*result_map)->clear();
736     }
737 
738     if (gCsvFile.is_open()) {
739         gCsvFile << "\n";
740     }
741 }
742