Mercurial > hg > orthanc-stl
comparison Sources/StructureSetGeometry.cpp @ 35:ee3bc8f7df5b
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 } |