Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructureSetUtils.cpp @ 998:38b6bb0bdd72
added a new set of classes that correctly handle non-convex polygons (not
used yet because of limitations in coordinates computing): DicomStructure2,
DicomStructureSet2, DicomStructurePolygon2, DicomStructureSetSlicer2. Too
many shortcuts have been taken when computing the actual position.
author | Benjamin Golinvaux <bgo@osimis.io> |
---|---|
date | Fri, 20 Sep 2019 11:58:00 +0200 |
parents | |
children | 2d8ab34c8c91 |
comparison
equal
deleted
inserted
replaced
995:9893fa8cd7a6 | 998:38b6bb0bdd72 |
---|---|
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-2019 Osimis S.A., Belgium | |
6 * | |
7 * This program is free software: you can redistribute it and/or | |
8 * modify it under the terms of the GNU Affero General Public License | |
9 * as published by the Free Software Foundation, either version 3 of | |
10 * the License, or (at your option) any later version. | |
11 * | |
12 * This program is distributed in the hope that it will be useful, but | |
13 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 * Affero General Public License for more details. | |
16 * | |
17 * You should have received a copy of the GNU Affero General Public License | |
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. | |
19 **/ | |
20 | |
21 #include "DicomStructureSetUtils.h" | |
22 | |
23 namespace OrthancStone | |
24 { | |
25 | |
26 #if 0 | |
27 void DicomStructure2::PartitionRectangleList(std::vector< std::vector<size_t> > & sets, const std::vector<RtStructRectanglesInSlab> slabCuts) | |
28 { | |
29 // map position ( <slabIndex,rectIndex> )--> disjoint set index | |
30 std::map<std::pair<size_t, size_t>, size_t> posToIndex; | |
31 | |
32 // disjoint set index --> position | |
33 std::map<size_t, std::pair<size_t, size_t> > indexToPos; | |
34 | |
35 size_t nextIndex = 0; | |
36 for (size_t i = 0; i < slabCuts.size(); ++i) | |
37 { | |
38 for (size_t j = 0; j < slabCuts[i].size(); ++j) | |
39 { | |
40 std::pair<size_t, size_t> pos(i, j); | |
41 posToIndex<pos> = nextIndex; | |
42 indexToPos<nextIndex> = pos; | |
43 } | |
44 } | |
45 // nextIndex is now the total rectangle count | |
46 DisjointDataSet ds(nextIndex); | |
47 | |
48 // we loop on all slabs (except the last one) and we connect all rectangles | |
49 if (slabCuts.size() < 2) | |
50 { | |
51 #error write special case | |
52 } | |
53 else | |
54 { | |
55 for (size_t i = 0; i < slabCuts.size() - 1; ++i) | |
56 { | |
57 for (size_t j = 0; j < slabCuts[i].size(); ++j) | |
58 { | |
59 const RtStructRectangleInSlab& r1 = slabCuts[i][j]; | |
60 const size_t r1i = posToIndex(std::pair<size_t, size_t>(i, j)); | |
61 for (size_t k = 0; k < slabCuts[i + 1].size(); ++k) | |
62 { | |
63 const RtStructRectangleInSlab& r2 = slabCuts[i + 1][k]; | |
64 const size_t r2i = posToIndex(std::pair<size_t, size_t>(i, j)); | |
65 // rect.xmin <= rectBottom.xmax && rectBottom.xmin <= rect.xmax | |
66 if ((r1.xmin <= r2.xmax) && (r2.xmin <= r1.xmax)) | |
67 { | |
68 #error now go! | |
69 } | |
70 | |
71 } | |
72 } | |
73 } | |
74 } | |
75 #endif | |
76 | |
77 /* | |
78 | |
79 compute list of segments : | |
80 | |
81 numberOfRectsFromHereOn = 0 | |
82 possibleNext = {in_k,in_kplus1} | |
83 | |
84 for all boundaries: | |
85 - we create a vertical segment and we push it | |
86 - if boundary is a start, numberOfRectsFromHereOn += 1. | |
87 - if we switch from 0 to 1, we start a segment | |
88 - if we switch from 1 to 2, we end the current segment and we record it | |
89 - if boundary is an end, numberOfRectsFromHereOn -= 1. | |
90 - if we switch from 1 to 0, we end the current segment and we record it | |
91 - if we switch from 2 to 1, we start a segment | |
92 */ | |
93 | |
94 // static | |
95 void AddSlabBoundaries( | |
96 std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, | |
97 const std::vector<RtStructRectanglesInSlab> & slabCuts, size_t iSlab) | |
98 { | |
99 if (iSlab < slabCuts.size()) | |
100 { | |
101 const RtStructRectanglesInSlab& slab = slabCuts[iSlab]; | |
102 for (size_t iRect = 0; iRect < slab.size(); ++iRect) | |
103 { | |
104 const RtStructRectangleInSlab& rect = slab[iRect]; | |
105 { | |
106 std::pair<double, RectangleBoundaryKind> boundary(rect.xmin, RectangleBoundaryKind_Start); | |
107 boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); | |
108 } | |
109 { | |
110 std::pair<double, RectangleBoundaryKind> boundary(rect.xmax, RectangleBoundaryKind_End); | |
111 boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); | |
112 } | |
113 } | |
114 } | |
115 } | |
116 | |
117 // static | |
118 void ProcessBoundaryList( | |
119 std::vector< std::pair<Point2D, Point2D> > & segments, | |
120 const std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, | |
121 double y) | |
122 { | |
123 Point2D start; | |
124 Point2D end; | |
125 int curNumberOfSegments = 0; // we count the number of segments. we only draw if it is 1 (not 0 or 2) | |
126 for (size_t i = 0; i < boundaries.size(); ++i) | |
127 { | |
128 switch (boundaries[i].second) | |
129 { | |
130 case RectangleBoundaryKind_Start: | |
131 curNumberOfSegments += 1; | |
132 switch (curNumberOfSegments) | |
133 { | |
134 case 0: | |
135 assert(false); | |
136 break; | |
137 case 1: | |
138 // a new segment has begun! | |
139 start.x = boundaries[i].first; | |
140 start.y = y; | |
141 break; | |
142 case 2: | |
143 // an extra segment has begun : stop the current one (we don't draw overlaps) | |
144 end.x = boundaries[i].first; | |
145 end.y = y; | |
146 segments.push_back(std::pair<Point2D, Point2D>(start, end)); | |
147 break; | |
148 default: | |
149 //assert(false); // seen IRL ! | |
150 break; | |
151 } | |
152 break; | |
153 case RectangleBoundaryKind_End: | |
154 curNumberOfSegments -= 1; | |
155 switch (curNumberOfSegments) | |
156 { | |
157 case 0: | |
158 // a lone (thus active) segment has ended. | |
159 end.x = boundaries[i].first; | |
160 end.y = y; | |
161 segments.push_back(std::pair<Point2D, Point2D>(start, end)); | |
162 break; | |
163 case 1: | |
164 // an extra segment has ended : start a new one one | |
165 start.x = boundaries[i].first; | |
166 start.y = y; | |
167 break; | |
168 default: | |
169 // this should not happen! | |
170 //assert(false); | |
171 break; | |
172 } | |
173 break; | |
174 default: | |
175 assert(false); | |
176 break; | |
177 } | |
178 } | |
179 } | |
180 | |
181 #if 0 | |
182 void ConvertListOfSlabsToSegments( | |
183 std::vector< std::pair<Point2D, Point2D> >& segments, | |
184 const std::vector<RtStructRectanglesInSlab>& slabCuts, | |
185 const size_t totalRectCount) | |
186 { | |
187 #error to delete | |
188 } | |
189 #else | |
190 // See https://www.dropbox.com/s/bllco6q8aazxk44/2019-09-18-rtstruct-cut-algorithm-rect-merge.png | |
191 void ConvertListOfSlabsToSegments( | |
192 std::vector< std::pair<Point2D, Point2D> > & segments, | |
193 const std::vector<RtStructRectanglesInSlab> & slabCuts, | |
194 const size_t totalRectCount) | |
195 { | |
196 if (slabCuts.size() == 0) | |
197 return; | |
198 | |
199 if (totalRectCount > 0) | |
200 segments.reserve(4 * totalRectCount); // worst case, but common. | |
201 | |
202 /* | |
203 VERTICAL | |
204 */ | |
205 for (size_t iSlab = 0; iSlab < slabCuts.size(); ++iSlab) | |
206 { | |
207 for (size_t iRect = 0; iRect < slabCuts[iSlab].size(); ++iRect) | |
208 { | |
209 const RtStructRectangleInSlab& rect = slabCuts[iSlab][iRect]; | |
210 { | |
211 Point2D p1(rect.xmin, rect.ymin); | |
212 Point2D p2(rect.xmin, rect.ymax); | |
213 segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); | |
214 } | |
215 { | |
216 Point2D p1(rect.xmax, rect.ymin); | |
217 Point2D p2(rect.xmax, rect.ymax); | |
218 segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); | |
219 } | |
220 } | |
221 } | |
222 | |
223 /* | |
224 HORIZONTAL | |
225 */ | |
226 | |
227 // if we have N slabs, we have N+1 potential vertical positions for horizontal segments | |
228 // - one for top of slab 0 | |
229 // - N-1 for all positions between two slabs | |
230 // - one for bottom of slab N-1 | |
231 | |
232 // this adds all the horizontal segments for the tops of 3the rectangles | |
233 // in row 0 | |
234 if (slabCuts[0].size() > 0) | |
235 { | |
236 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; | |
237 AddSlabBoundaries(boundaries, slabCuts, 0); | |
238 | |
239 ProcessBoundaryList(segments, boundaries, slabCuts[0][0].ymin); | |
240 } | |
241 | |
242 // this adds all the horizontal segments belonging to two slabs | |
243 for (size_t iSlab = 0; iSlab < slabCuts.size() - 1; ++iSlab) | |
244 { | |
245 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; | |
246 AddSlabBoundaries(boundaries, slabCuts, iSlab); | |
247 AddSlabBoundaries(boundaries, slabCuts, iSlab + 1); | |
248 double curY = 0; | |
249 if (slabCuts[iSlab].size() > 0) | |
250 { | |
251 curY = slabCuts[iSlab][0].ymax; | |
252 ProcessBoundaryList(segments, boundaries, curY); | |
253 } | |
254 else if (slabCuts[iSlab + 1].size() > 0) | |
255 { | |
256 curY = slabCuts[iSlab + 1][0].ymin; | |
257 ProcessBoundaryList(segments, boundaries, curY); | |
258 } | |
259 else | |
260 { | |
261 // nothing to do!! : both slab lists are empty! | |
262 } | |
263 } | |
264 | |
265 // this adds all the horizontal segments for the BOTTOM of the rectangles | |
266 // on last row | |
267 if (slabCuts[slabCuts.size() - 1].size() > 0) | |
268 { | |
269 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; | |
270 AddSlabBoundaries(boundaries, slabCuts, slabCuts.size() - 1); | |
271 | |
272 ProcessBoundaryList(segments, boundaries, slabCuts[slabCuts.size() - 1][0].ymax); | |
273 } | |
274 } | |
275 #endif | |
276 } |