Mercurial > hg > orthanc-stone

comparison OrthancStone/Resources/Graveyard/RTStructTentativeReimplementation-BGO/DicomStructureSetUtils.cpp @ 1908:affde38b84de

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved tentative bgo reimplementation of rt-struct into graveyard
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 01 Feb 2022 08:38:32 +0100
parents OrthancStone/Sources/Toolbox/DicomStructureSetUtils.cpp@14c8f339d480
children 07964689cb0b
comparison
equal deleted inserted replaced
1907:0208f99b8bde 1908:affde38b84de
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 #include "DicomStructureSetUtils.h"
24
25 namespace OrthancStone
26 {
27
28 #if 0
29 void DicomStructure2::PartitionRectangleList(std::vector< std::vector<size_t> > & sets, const std::vector<RtStructRectanglesInSlab> slabCuts)
30 {
31 // map position ( <slabIndex,rectIndex> )--> disjoint set index
32 std::map<std::pair<size_t, size_t>, size_t> posToIndex;
33
34 // disjoint set index --> position
35 std::map<size_t, std::pair<size_t, size_t> > indexToPos;
36
37 size_t nextIndex = 0;
38 for (size_t i = 0; i < slabCuts.size(); ++i)
39 {
40 for (size_t j = 0; j < slabCuts[i].size(); ++j)
41 {
42 std::pair<size_t, size_t> pos(i, j);
43 posToIndex<pos> = nextIndex;
44 indexToPos<nextIndex> = pos;
45 }
46 }
47 // nextIndex is now the total rectangle count
48 DisjointDataSet ds(nextIndex);
49
50 // we loop on all slabs (except the last one) and we connect all rectangles
51 if (slabCuts.size() < 2)
52 {
53 #error write special case
54 }
55 else
56 {
57 for (size_t i = 0; i < slabCuts.size() - 1; ++i)
58 {
59 for (size_t j = 0; j < slabCuts[i].size(); ++j)
60 {
61 const RtStructRectangleInSlab& r1 = slabCuts[i][j];
62 const size_t r1i = posToIndex(std::pair<size_t, size_t>(i, j));
63 for (size_t k = 0; k < slabCuts[i + 1].size(); ++k)
64 {
65 const RtStructRectangleInSlab& r2 = slabCuts[i + 1][k];
66 const size_t r2i = posToIndex(std::pair<size_t, size_t>(i, j));
67 // rect.xmin <= rectBottom.xmax && rectBottom.xmin <= rect.xmax
68 if ((r1.xmin <= r2.xmax) && (r2.xmin <= r1.xmax))
69 {
70 #error now go!
71 }
72
73 }
74 }
75 }
76 }
77 #endif
78
79 /*
80
81 compute list of segments :
82
83 numberOfRectsFromHereOn = 0
84 possibleNext = {in_k,in_kplus1}
85
86 for all boundaries:
87 - we create a vertical segment and we push it
88 - if boundary is a start, numberOfRectsFromHereOn += 1.
89 - if we switch from 0 to 1, we start a segment
90 - if we switch from 1 to 2, we end the current segment and we record it
91 - if boundary is an end, numberOfRectsFromHereOn -= 1.
92 - if we switch from 1 to 0, we end the current segment and we record it
93 - if we switch from 2 to 1, we start a segment
94 */
95
96 // static
97 void AddSlabBoundaries(
98 std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries,
99 const std::vector<RtStructRectanglesInSlab> & slabCuts, size_t iSlab)
100 {
101 if (iSlab < slabCuts.size())
102 {
103 const RtStructRectanglesInSlab& slab = slabCuts[iSlab];
104 for (size_t iRect = 0; iRect < slab.size(); ++iRect)
105 {
106 const RtStructRectangleInSlab& rect = slab[iRect];
107 {
108 std::pair<double, RectangleBoundaryKind> boundary(rect.xmin, RectangleBoundaryKind_Start);
109 boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary);
110 }
111 {
112 std::pair<double, RectangleBoundaryKind> boundary(rect.xmax, RectangleBoundaryKind_End);
113 boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary);
114 }
115 }
116 }
117 }
118
119 // static
120 void ProcessBoundaryList(
121 std::vector< std::pair<ScenePoint2D, ScenePoint2D> > & segments,
122 const std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries,
123 double y)
124 {
125 ScenePoint2D start;
126 ScenePoint2D end;
127 int curNumberOfSegments = 0; // we count the number of segments. we only draw if it is 1 (not 0 or 2)
128 for (size_t i = 0; i < boundaries.size(); ++i)
129 {
130 switch (boundaries[i].second)
131 {
132 case RectangleBoundaryKind_Start:
133 curNumberOfSegments += 1;
134 switch (curNumberOfSegments)
135 {
136 case 0:
137 assert(false);
138 break;
139 case 1:
140 // a new segment has begun!
141 start = ScenePoint2D(boundaries[i].first, y);
142 break;
143 case 2:
144 // an extra segment has begun : stop the current one (we don't draw overlaps)
145 end = ScenePoint2D(boundaries[i].first, y);
146 segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(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 = ScenePoint2D(boundaries[i].first, y);
160 segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(start, end));
161 break;
162 case 1:
163 // an extra segment has ended : start a new one one
164 start = ScenePoint2D(boundaries[i].first, y);
165 break;
166 default:
167 // this should not happen!
168 //assert(false);
169 break;
170 }
171 break;
172 default:
173 assert(false);
174 break;
175 }
176 }
177 }
178
179 #if 0
180 void ConvertListOfSlabsToSegments(
181 std::vector< std::pair<ScenePoint2D, ScenePoint2D> >& segments,
182 const std::vector<RtStructRectanglesInSlab>& slabCuts,
183 const size_t totalRectCount)
184 {
185 #error to delete
186 }
187 #else
188 // See https://www.dropbox.com/s/bllco6q8aazxk44/2019-09-18-rtstruct-cut-algorithm-rect-merge.png
189 void ConvertListOfSlabsToSegments(
190 std::vector< std::pair<ScenePoint2D, ScenePoint2D> > & segments,
191 const std::vector<RtStructRectanglesInSlab> & slabCuts,
192 const size_t totalRectCount)
193 {
194 if (slabCuts.size() == 0)
195 return;
196
197 if (totalRectCount > 0)
198 segments.reserve(4 * totalRectCount); // worst case, but common.
199
200 /*
201 VERTICAL
202 */
203 for (size_t iSlab = 0; iSlab < slabCuts.size(); ++iSlab)
204 {
205 for (size_t iRect = 0; iRect < slabCuts[iSlab].size(); ++iRect)
206 {
207 const RtStructRectangleInSlab& rect = slabCuts[iSlab][iRect];
208 {
209 ScenePoint2D p1(rect.xmin, rect.ymin);
210 ScenePoint2D p2(rect.xmin, rect.ymax);
211 segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(p1, p2));
212 }
213 {
214 ScenePoint2D p1(rect.xmax, rect.ymin);
215 ScenePoint2D p2(rect.xmax, rect.ymax);
216 segments.push_back(std::pair<ScenePoint2D, ScenePoint2D>(p1, p2));
217 }
218 }
219 }
220
221 /*
222 HORIZONTAL
223 */
224
225 // if we have N slabs, we have N+1 potential vertical positions for horizontal segments
226 // - one for top of slab 0
227 // - N-1 for all positions between two slabs
228 // - one for bottom of slab N-1
229
230 // this adds all the horizontal segments for the tops of 3the rectangles
231 // in row 0
232 if (slabCuts[0].size() > 0)
233 {
234 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries;
235 AddSlabBoundaries(boundaries, slabCuts, 0);
236
237 ProcessBoundaryList(segments, boundaries, slabCuts[0][0].ymin);
238 }
239
240 // this adds all the horizontal segments belonging to two slabs
241 for (size_t iSlab = 0; iSlab < slabCuts.size() - 1; ++iSlab)
242 {
243 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries;
244 AddSlabBoundaries(boundaries, slabCuts, iSlab);
245 AddSlabBoundaries(boundaries, slabCuts, iSlab + 1);
246 double curY = 0;
247 if (slabCuts[iSlab].size() > 0)
248 {
249 curY = slabCuts[iSlab][0].ymax;
250 ProcessBoundaryList(segments, boundaries, curY);
251 }
252 else if (slabCuts[iSlab + 1].size() > 0)
253 {
254 curY = slabCuts[iSlab + 1][0].ymin;
255 ProcessBoundaryList(segments, boundaries, curY);
256 }
257 else
258 {
259 // nothing to do!! : both slab lists are empty!
260 }
261 }
262
263 // this adds all the horizontal segments for the BOTTOM of the rectangles
264 // on last row
265 if (slabCuts[slabCuts.size() - 1].size() > 0)
266 {
267 std::vector<std::pair<double, RectangleBoundaryKind> > boundaries;
268 AddSlabBoundaries(boundaries, slabCuts, slabCuts.size() - 1);
269
270 ProcessBoundaryList(segments, boundaries, slabCuts[slabCuts.size() - 1][0].ymax);
271 }
272 }
273 #endif
274 }