Mercurial > hg > orthanc-stl

comparison Sources/StructureSetGeometry.cpp @ 35:ee3bc8f7df5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 04 Apr 2024 20:37:08 +0200
parents
children 8a1daa321afe
comparison
equal deleted inserted replaced
34:bee2017f3088 35:ee3bc8f7df5b
1 /**
2 * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
3 * SPDX-License-Identifier: GPL-3.0-or-later
4 */
5
6 /**
7 * STL plugin for Orthanc
8 * Copyright (C) 2023-2024 Sebastien Jodogne, UCLouvain, Belgium
9 *
10 * This program is free software: you can redistribute it and/or
11 * modify it under the terms of the GNU General Public License as
12 * published by the Free Software Foundation, either version 3 of the
13 * License, or (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful, but
16 * WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with this program. If not, see <http://www.gnu.org/licenses/>.
22 **/
23
24
25 #include "StructureSetGeometry.h"
26
27 #include "STLToolbox.h"
28
29 #include <OrthancException.h>
30
31 #include <list>
32
33
34 bool StructureSetGeometry::LookupProjectionIndex(size_t& index,
35 double z) const
36 {
37 if (slicesCount_ == 0)
38 {
39 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
40 }
41
42 if (z < minProjectionAlongNormal_ ||
43 z > maxProjectionAlongNormal_)
44 {
45 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
46 }
47
48 assert(slicesSpacing_ > 0 &&
49 minProjectionAlongNormal_ < maxProjectionAlongNormal_);
50
51 double d = (z - minProjectionAlongNormal_) / slicesSpacing_;
52
53 if (STLToolbox::IsNear(d, round(d)))
54 {
55 if (d < 0.0 ||
56 d > static_cast<double>(slicesCount_) - 1.0)
57 {
58 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
59 }
60 else
61 {
62 index = static_cast<size_t>(round(d));
63 return true;
64 }
65 }
66 else
67 {
68 return false;
69 }
70 }
71
72
73 StructureSetGeometry::StructureSetGeometry(const StructureSet& structures,
74 bool strict)
75 {
76 bool isValid = false;
77
78 std::vector<double> projections;
79 projections.reserve(structures.GetPolygonsCount());
80
81 for (size_t i = 0; i < structures.GetPolygonsCount(); i++)
82 {
83 Vector3D normal;
84 if (structures.GetPolygon(i).IsCoplanar(normal))
85 {
86 // Initialize the normal of the whole volume, if need be
87 if (!isValid)
88 {
89 isValid = true;
90 slicesNormal_ = normal;
91 }
92
93 if (Vector3D::AreParallel(normal, slicesNormal_))
94 {
95 // This is a valid slice (it is parallel to the normal)
96 const Vector3D& point = structures.GetPolygon(i).GetPoint(0);
97 projections.push_back(Vector3D::DotProduct(point, slicesNormal_));
98 }
99 else
100 {
101 // RT-STRUCT with non-parallel slices
102 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
103 }
104 }
105 else
106 {
107 // Ignore slices that are not coplanar
108 }
109 }
110
111 if (projections.empty())
112 {
113 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented,
114 "Structure set without a valid geometry");
115 }
116
117 // Only keep unique projections
118
119 std::sort(projections.begin(), projections.end());
120 STLToolbox::RemoveDuplicateValues(projections);
121 assert(!projections.empty());
122
123 if (projections.size() == 1)
124 {
125 // Volume with one single slice
126 minProjectionAlongNormal_ = projections[0];
127 maxProjectionAlongNormal_ = projections[0];
128 slicesSpacing_ = 1; // Arbitrary value
129 slicesCount_ = 1;
130 return;
131 }
132
133
134 // Compute the most probable spacing between the slices
135
136 {
137 std::vector<double> spacings;
138 spacings.resize(projections.size() - 1);
139
140 for (size_t i = 0; i < spacings.size(); i++)
141 {
142 spacings[i] = projections[i + 1] - projections[i];
143 assert(spacings[i] > 0);
144 }
145
146 std::sort(spacings.begin(), spacings.end());
147 STLToolbox::RemoveDuplicateValues(spacings);
148
149 if (spacings.empty())
150 {
151 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
152 }
153
154 slicesSpacing_ = spacings[spacings.size() / 10]; // Take the 90% percentile of smallest spacings
155 assert(slicesSpacing_ > 0);
156 }
157
158
159 // Find the projection along the normal with the largest support
160
161 bool first = true;
162 size_t bestSupport;
163 double bestProjection;
164
165 std::list<size_t> candidates;
166 for (size_t i = 0; i < projections.size(); i++)
167 {
168 candidates.push_back(i);
169 }
170
171 while (!candidates.empty())
172 {
173 std::list<size_t> next;
174
175 size_t countSupport = 0;
176
177 std::list<size_t>::const_iterator it = candidates.begin();
178 size_t reference = *it;
179 it++;
180
181 while (it != candidates.end())
182 {
183 double d = (projections[*it] - projections[reference]) / slicesSpacing_;
184 if (STLToolbox::IsNear(d, round(d)))
185 {
186 countSupport ++;
187 }
188 else
189 {
190 next.push_back(*it);
191 }
192
193 it++;
194 }
195
196 if (first ||
197 countSupport > bestSupport)
198 {
199 first = false;
200 bestSupport = countSupport;
201 bestProjection = projections[reference];
202 }
203
204 if (strict &&
205 !next.empty())
206 {
207 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat,
208 "Structure set with multiple support, which is not allowed in Strict mode");
209 }
210
211 candidates.swap(next);
212 }
213
214
215 // Compute the range of the projections
216
217 minProjectionAlongNormal_ = bestProjection;
218 maxProjectionAlongNormal_ = bestProjection;
219
220 for (size_t i = 0; i < projections.size(); i++)
221 {
222 double d = (projections[i] - bestProjection) / slicesSpacing_;
223 if (STLToolbox::IsNear(d, round(d)))
224 {
225 minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]);
226 maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]);
227 }
228 }
229
230 double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_;
231 if (STLToolbox::IsNear(d, round(d)))
232 {
233 slicesCount_ = static_cast<size_t>(round(d)) + 1;
234 }
235 else
236 {
237 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
238 }
239
240
241 // Sanity check
242
243 size_t a, b;
244 if (!LookupProjectionIndex(a, minProjectionAlongNormal_) ||
245 !LookupProjectionIndex(b, maxProjectionAlongNormal_) ||
246 a != 0 ||
247 b + 1 != slicesCount_)
248 {
249 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
250 }
251 }
252
253
254 bool StructureSetGeometry::ProjectAlongNormal(double& z,
255 const StructurePolygon& polygon) const
256 {
257 Vector3D normal;
258 if (polygon.IsCoplanar(normal) &&
259 Vector3D::AreParallel(normal, slicesNormal_))
260 {
261 z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_);
262 return true;
263 }
264 else
265 {
266 return false;
267 }
268 }
269
270
271 bool StructureSetGeometry::LookupSliceIndex(size_t& slice,
272 const StructurePolygon& polygon) const
273 {
274 double z;
275 return (ProjectAlongNormal(z, polygon) &&
276 LookupProjectionIndex(slice, z));
277 }