Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/UnionOfRectangles.cpp @ 1875:b896f20d24ca
added UnionOfRectangles algorithm
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 11 Jan 2022 19:59:40 +0100 |
parents | |
children | b1f510e601d2 |
comparison
equal
deleted
inserted
replaced
1874:08f2476e8f5e | 1875:b896f20d24ca |
---|---|
1 /** | |
2 * Stone of Orthanc | |
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics | |
4 * Department, University Hospital of Liege, Belgium | |
5 * Copyright (C) 2017-2022 Osimis S.A., Belgium | |
6 * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium | |
7 * | |
8 * This program is free software: you can redistribute it and/or | |
9 * modify it under the terms of the GNU Lesser General Public License | |
10 * as published by the Free Software Foundation, either version 3 of | |
11 * the License, or (at your option) any later version. | |
12 * | |
13 * This program is distributed in the hope that it will be useful, but | |
14 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
16 * Lesser General Public License for more details. | |
17 * | |
18 * You should have received a copy of the GNU Lesser General Public | |
19 * License along with this program. If not, see | |
20 * <http://www.gnu.org/licenses/>. | |
21 **/ | |
22 | |
23 | |
24 #include "UnionOfRectangles.h" | |
25 | |
26 #include "Internals/OrientedIntegerLine2D.h" | |
27 #include "Internals/RectanglesIntegerProjection.h" | |
28 #include "SegmentTree.h" | |
29 | |
30 #include <OrthancException.h> | |
31 | |
32 #include <stack> | |
33 | |
34 | |
35 namespace OrthancStone | |
36 { | |
37 class UnionOfRectangles::Payload : public Orthanc::IDynamicObject | |
38 { | |
39 private: | |
40 int counter_; | |
41 Status status_; | |
42 | |
43 public: | |
44 Payload() : | |
45 counter_(0), | |
46 status_(Status_Empty) | |
47 { | |
48 } | |
49 | |
50 int GetCounter() const | |
51 { | |
52 return counter_; | |
53 } | |
54 | |
55 Status GetStatus() const | |
56 { | |
57 return status_; | |
58 } | |
59 | |
60 void SetStatus(Status status) | |
61 { | |
62 status_ = status; | |
63 } | |
64 | |
65 void Increment() | |
66 { | |
67 counter_ ++; | |
68 } | |
69 | |
70 void Decrement() | |
71 { | |
72 if (counter_ == 0) | |
73 { | |
74 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
75 } | |
76 else | |
77 { | |
78 counter_ --; | |
79 } | |
80 } | |
81 }; | |
82 | |
83 | |
84 class UnionOfRectangles::Factory : public SegmentTree::IPayloadFactory | |
85 { | |
86 public: | |
87 virtual Orthanc::IDynamicObject* Create() | |
88 { | |
89 return new Payload; | |
90 } | |
91 }; | |
92 | |
93 | |
94 class UnionOfRectangles::Visitor : public SegmentTree::IVisitor | |
95 { | |
96 private: | |
97 Operation operation_; | |
98 | |
99 public: | |
100 Visitor(Operation operation) : | |
101 operation_(operation) | |
102 { | |
103 } | |
104 | |
105 virtual void Visit(const SegmentTree& node, | |
106 bool fullyInside) | |
107 { | |
108 Payload& payload = node.GetTypedPayload<Payload>(); | |
109 | |
110 if (fullyInside) | |
111 { | |
112 switch (operation_) | |
113 { | |
114 case Operation_Insert: | |
115 payload.Increment(); | |
116 break; | |
117 | |
118 case Operation_Delete: | |
119 payload.Decrement(); | |
120 break; | |
121 | |
122 default: | |
123 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); | |
124 } | |
125 } | |
126 | |
127 if (payload.GetCounter() > 0) | |
128 { | |
129 payload.SetStatus(Status_Full); | |
130 } | |
131 else if (node.IsLeaf()) | |
132 { | |
133 payload.SetStatus(Status_Empty); | |
134 } | |
135 else if (node.GetLeftChild().GetTypedPayload<Payload>().GetStatus() == Status_Empty && | |
136 node.GetRightChild().GetTypedPayload<Payload>().GetStatus() == Status_Empty) | |
137 { | |
138 payload.SetStatus(Status_Empty); | |
139 } | |
140 else | |
141 { | |
142 payload.SetStatus(Status_Partial); | |
143 } | |
144 } | |
145 | |
146 | |
147 // This is the "CONTR()" function from the textbook | |
148 static void IntersectComplement(std::stack<size_t>& stack, | |
149 size_t low, | |
150 size_t high, | |
151 const SegmentTree& node) | |
152 { | |
153 if (low >= high) | |
154 { | |
155 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
156 } | |
157 | |
158 Status status = node.GetTypedPayload<Payload>().GetStatus(); | |
159 | |
160 if (status != Status_Full) | |
161 { | |
162 assert(status == Status_Partial || | |
163 status == Status_Empty); | |
164 | |
165 // Aliases to use the same variable names as in the textbook | |
166 const size_t& b = low; | |
167 const size_t& e = high; | |
168 const size_t& bv = node.GetLowBound(); | |
169 const size_t& ev = node.GetHighBound(); | |
170 | |
171 if (b <= bv && | |
172 ev <= e && | |
173 status == Status_Empty) | |
174 { | |
175 // [B[v], E[v]] is contributed | |
176 if (!stack.empty() && | |
177 stack.top() == bv) | |
178 { | |
179 stack.pop(); // Merge continuous segments | |
180 } | |
181 else | |
182 { | |
183 stack.push(bv); // Beginning of edge | |
184 } | |
185 | |
186 stack.push(ev); // Current termination of edge | |
187 } | |
188 else | |
189 { | |
190 size_t middle = (bv + ev) / 2; | |
191 | |
192 if (b < middle) | |
193 { | |
194 IntersectComplement(stack, b, e, node.GetLeftChild()); | |
195 } | |
196 | |
197 if (middle < e) | |
198 { | |
199 IntersectComplement(stack, b, e, node.GetRightChild()); | |
200 } | |
201 } | |
202 } | |
203 } | |
204 }; | |
205 | |
206 | |
207 static void AddVerticalEdges(std::list<Internals::OrientedIntegerLine2D>& edges, | |
208 std::stack<size_t>& stack, | |
209 size_t x, | |
210 bool isLeft) | |
211 { | |
212 if (stack.size() % 2 != 0) | |
213 { | |
214 throw std::runtime_error("internal error"); | |
215 } | |
216 | |
217 typedef std::pair<size_t, size_t> Interval; | |
218 | |
219 std::list<Interval> intervals; | |
220 | |
221 while (!stack.empty()) | |
222 { | |
223 size_t high = stack.top(); | |
224 stack.pop(); | |
225 size_t low = stack.top(); | |
226 stack.pop(); | |
227 | |
228 if (!intervals.empty() && | |
229 intervals.back().second == low) | |
230 { | |
231 // Extend the previous interval | |
232 intervals.back().second = high; | |
233 } | |
234 else | |
235 { | |
236 intervals.push_back(std::make_pair(low, high)); | |
237 } | |
238 } | |
239 | |
240 for (std::list<Interval>::const_iterator | |
241 it = intervals.begin(); it != intervals.end(); ++it) | |
242 { | |
243 if (it->first == it->second) | |
244 { | |
245 throw std::runtime_error("parameter out of range"); | |
246 } | |
247 | |
248 // By convention, the left sides go downward, and the right go upward | |
249 if (isLeft) | |
250 { | |
251 edges.push_back(Internals::OrientedIntegerLine2D(x, it->first, x, it->second)); | |
252 } | |
253 else | |
254 { | |
255 edges.push_back(Internals::OrientedIntegerLine2D(x, it->second, x, it->first)); | |
256 } | |
257 } | |
258 } | |
259 | |
260 | |
261 class UnionOfRectangles::VerticalSide | |
262 { | |
263 private: | |
264 size_t x_; | |
265 bool isLeft_; | |
266 size_t y1_; | |
267 size_t y2_; | |
268 | |
269 public: | |
270 VerticalSide(size_t x, | |
271 bool isLeft, | |
272 size_t y1, | |
273 size_t y2) : | |
274 x_(x), | |
275 isLeft_(isLeft), | |
276 y1_(y1), | |
277 y2_(y2) | |
278 { | |
279 assert(y1 < y2); | |
280 } | |
281 | |
282 size_t GetX() const | |
283 { | |
284 return x_; | |
285 } | |
286 | |
287 bool IsLeft() const | |
288 { | |
289 return isLeft_; | |
290 } | |
291 | |
292 size_t GetY1() const | |
293 { | |
294 return y1_; | |
295 } | |
296 | |
297 size_t GetY2() const | |
298 { | |
299 return y2_; | |
300 } | |
301 | |
302 bool operator< (const VerticalSide& other) const | |
303 { | |
304 if (x_ < other.x_) | |
305 { | |
306 return true; | |
307 } | |
308 else if (x_ == other.x_) | |
309 { | |
310 return isLeft_; | |
311 } | |
312 else | |
313 { | |
314 return false; | |
315 } | |
316 } | |
317 | |
318 bool Equals(const VerticalSide& other) const | |
319 { | |
320 return (x_ == other.x_ && | |
321 isLeft_ == other.isLeft_); | |
322 } | |
323 }; | |
324 | |
325 | |
326 class UnionOfRectangles::HorizontalJunction | |
327 { | |
328 private: | |
329 size_t x_; | |
330 size_t y_; | |
331 size_t ybis_; | |
332 bool downward_; | |
333 | |
334 public: | |
335 HorizontalJunction(size_t x, | |
336 size_t y, | |
337 size_t ybis, | |
338 bool downward) : | |
339 x_(x), | |
340 y_(y), | |
341 ybis_(ybis), | |
342 downward_(downward) | |
343 { | |
344 } | |
345 | |
346 size_t GetX() const | |
347 { | |
348 return x_; | |
349 } | |
350 | |
351 size_t GetY() const | |
352 { | |
353 return y_; | |
354 } | |
355 | |
356 size_t GetYBis() const | |
357 { | |
358 return ybis_; | |
359 } | |
360 | |
361 bool IsDownward() const | |
362 { | |
363 return downward_; | |
364 } | |
365 | |
366 bool operator< (const HorizontalJunction& other) const | |
367 { | |
368 if (y_ > other.y_) | |
369 { | |
370 return true; | |
371 } | |
372 else if (y_ == other.y_) | |
373 { | |
374 return x_ < other.x_; | |
375 } | |
376 else | |
377 { | |
378 return false; | |
379 } | |
380 } | |
381 }; | |
382 | |
383 | |
384 void UnionOfRectangles::Apply(std::list< std::vector<ScenePoint2D> >& contours, | |
385 const std::list<Extent2D>& rectangles) | |
386 { | |
387 /** | |
388 * STEP 1 | |
389 **/ | |
390 Internals::RectanglesIntegerProjection horizontalProjection(rectangles, true); | |
391 Internals::RectanglesIntegerProjection verticalProjection(rectangles, false); | |
392 | |
393 assert(horizontalProjection.GetProjectedRectanglesCount() == verticalProjection.GetProjectedRectanglesCount()); | |
394 | |
395 | |
396 /** | |
397 * STEP 2 | |
398 **/ | |
399 Factory factory; | |
400 SegmentTree tree(0, verticalProjection.GetEndpointsCount() - 1, factory); | |
401 | |
402 | |
403 /** | |
404 * STEP 3 | |
405 **/ | |
406 std::vector<VerticalSide> verticalSides; | |
407 | |
408 const size_t countRectangles = horizontalProjection.GetProjectedRectanglesCount(); | |
409 | |
410 verticalSides.reserve(2 * countRectangles); | |
411 | |
412 for (size_t i = 0; i < countRectangles; i++) | |
413 { | |
414 size_t horizontalLow = horizontalProjection.GetProjectedRectangleLow(i); | |
415 size_t horizontalHigh = horizontalProjection.GetProjectedRectangleHigh(i); | |
416 size_t verticalLow = verticalProjection.GetProjectedRectangleLow(i); | |
417 size_t verticalHigh = verticalProjection.GetProjectedRectangleHigh(i); | |
418 | |
419 verticalSides.push_back(VerticalSide(horizontalLow, true, verticalLow, verticalHigh)); | |
420 verticalSides.push_back(VerticalSide(horizontalHigh, false, verticalLow, verticalHigh)); | |
421 } | |
422 | |
423 assert(verticalSides.size() == 2 * countRectangles); | |
424 | |
425 std::sort(verticalSides.begin(), verticalSides.end()); | |
426 | |
427 | |
428 /** | |
429 * STEP 4 | |
430 **/ | |
431 std::list<Internals::OrientedIntegerLine2D> verticalEdges; | |
432 std::stack<size_t> stack; | |
433 | |
434 for (size_t i = 0; i < verticalSides.size(); i++) | |
435 { | |
436 if (i > 0 && | |
437 !verticalSides[i].Equals(verticalSides[i - 1])) | |
438 { | |
439 // Output the stack | |
440 AddVerticalEdges(verticalEdges, stack, verticalSides[i - 1].GetX(), | |
441 verticalSides[i - 1].IsLeft()); | |
442 } | |
443 | |
444 size_t y1 = verticalSides[i].GetY1(); | |
445 size_t y2 = verticalSides[i].GetY2(); | |
446 | |
447 if (verticalSides[i].IsLeft()) | |
448 { | |
449 Visitor::IntersectComplement(stack, y1, y2, tree); | |
450 | |
451 Visitor visitor(Operation_Insert); | |
452 tree.VisitSegment(y1, y2, visitor); | |
453 } | |
454 else | |
455 { | |
456 Visitor visitor(Operation_Delete); | |
457 tree.VisitSegment(y1, y2, visitor); | |
458 | |
459 Visitor::IntersectComplement(stack, y1, y2, tree); | |
460 } | |
461 } | |
462 | |
463 if (!verticalSides.empty() && | |
464 !stack.empty()) | |
465 { | |
466 AddVerticalEdges(verticalEdges, stack, verticalSides.back().GetX(), | |
467 verticalSides.back().IsLeft()); | |
468 } | |
469 | |
470 | |
471 /** | |
472 * STEP 5: Horizontal edges | |
473 **/ | |
474 std::vector<HorizontalJunction> horizontalJunctions; | |
475 horizontalJunctions.reserve(2 * verticalEdges.size()); | |
476 | |
477 for (std::list<Internals::OrientedIntegerLine2D>::const_iterator | |
478 it = verticalEdges.begin(); it != verticalEdges.end(); it++) | |
479 { | |
480 assert(it->GetX1() == it->GetX2()); | |
481 horizontalJunctions.push_back(HorizontalJunction(it->GetX1(), it->GetY1(), it->GetY2(), it->IsDownward())); | |
482 horizontalJunctions.push_back(HorizontalJunction(it->GetX1(), it->GetY2(), it->GetY1(), it->IsDownward())); | |
483 } | |
484 | |
485 assert(horizontalJunctions.size() == 2 * verticalEdges.size()); | |
486 std::sort(horizontalJunctions.begin(), horizontalJunctions.end()); | |
487 | |
488 std::list<Internals::OrientedIntegerLine2D> horizontalEdges; | |
489 for (size_t i = 0; i < horizontalJunctions.size(); i += 2) | |
490 { | |
491 size_t x1 = horizontalJunctions[i].GetX(); | |
492 size_t x2 = horizontalJunctions[i + 1].GetX(); | |
493 size_t y = horizontalJunctions[i].GetY(); | |
494 | |
495 if ((horizontalJunctions[i].IsDownward() && y > horizontalJunctions[i].GetYBis()) || | |
496 (!horizontalJunctions[i].IsDownward() && y < horizontalJunctions[i].GetYBis())) | |
497 { | |
498 horizontalEdges.push_back(Internals::OrientedIntegerLine2D(x1, y, x2, y)); | |
499 } | |
500 else | |
501 { | |
502 horizontalEdges.push_back(Internals::OrientedIntegerLine2D(x2, y, x1, y)); | |
503 } | |
504 } | |
505 | |
506 | |
507 /** | |
508 * POST-PROCESSING: Combine the separate sets of horizontal and | |
509 * vertical edges, into a set of 2D chains | |
510 **/ | |
511 std::vector<Internals::OrientedIntegerLine2D> allEdges; | |
512 allEdges.reserve(horizontalEdges.size() + verticalEdges.size()); | |
513 | |
514 for (std::list<Internals::OrientedIntegerLine2D>::const_iterator | |
515 it = horizontalEdges.begin(); it != horizontalEdges.end(); it++) | |
516 { | |
517 allEdges.push_back(*it); | |
518 } | |
519 | |
520 for (std::list<Internals::OrientedIntegerLine2D>::const_iterator | |
521 it = verticalEdges.begin(); it != verticalEdges.end(); it++) | |
522 { | |
523 allEdges.push_back(*it); | |
524 } | |
525 | |
526 assert(allEdges.size() == horizontalEdges.size() + verticalEdges.size()); | |
527 | |
528 std::list<Internals::OrientedIntegerLine2D::Chain> chains; | |
529 Internals::OrientedIntegerLine2D::ExtractChains(chains, allEdges); | |
530 | |
531 contours.clear(); | |
532 | |
533 for (std::list<Internals::OrientedIntegerLine2D::Chain>::const_iterator | |
534 it = chains.begin(); it != chains.end(); ++it) | |
535 { | |
536 assert(!it->empty()); | |
537 | |
538 std::vector<ScenePoint2D> chain; | |
539 chain.reserve(it->size()); | |
540 | |
541 for (Internals::OrientedIntegerLine2D::Chain::const_iterator | |
542 it2 = it->begin(); it2 != it->end(); ++it2) | |
543 { | |
544 chain.push_back(ScenePoint2D(horizontalProjection.GetEndpointCoordinate(it2->first), | |
545 verticalProjection.GetEndpointCoordinate(it2->second))); | |
546 } | |
547 | |
548 assert(!chain.empty()); | |
549 contours.push_back(chain); | |
550 } | |
551 } | |
552 } |