Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructureSetUtils.cpp @ 998:38b6bb0bdd72

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }