Mercurial > hg > orthanc-stone
comparison OrthancStone/Sources/Toolbox/SlicesSorter.cpp @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Toolbox/SlicesSorter.cpp@30deba7bc8e2 |
children | 85e117739eca |
comparison
equal
deleted
inserted
replaced
1511:9dfeee74c1e6 | 1512:244ad1e4e76a |
---|---|
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-2020 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 | |
22 #include "SlicesSorter.h" | |
23 | |
24 #include "GeometryToolbox.h" | |
25 | |
26 #include <OrthancException.h> | |
27 | |
28 namespace OrthancStone | |
29 { | |
30 class SlicesSorter::SliceWithDepth : public boost::noncopyable | |
31 { | |
32 private: | |
33 CoordinateSystem3D geometry_; | |
34 double depth_; | |
35 | |
36 std::unique_ptr<Orthanc::IDynamicObject> payload_; | |
37 | |
38 public: | |
39 SliceWithDepth(const CoordinateSystem3D& geometry, | |
40 Orthanc::IDynamicObject* payload) : | |
41 geometry_(geometry), | |
42 depth_(0), | |
43 payload_(payload) | |
44 { | |
45 } | |
46 | |
47 void SetNormal(const Vector& normal) | |
48 { | |
49 depth_ = boost::numeric::ublas::inner_prod(geometry_.GetOrigin(), normal); | |
50 } | |
51 | |
52 double GetDepth() const | |
53 { | |
54 return depth_; | |
55 } | |
56 | |
57 const CoordinateSystem3D& GetGeometry() const | |
58 { | |
59 return geometry_; | |
60 } | |
61 | |
62 bool HasPayload() const | |
63 { | |
64 return (payload_.get() != NULL); | |
65 } | |
66 | |
67 const Orthanc::IDynamicObject& GetPayload() const | |
68 { | |
69 if (HasPayload()) | |
70 { | |
71 return *payload_; | |
72 } | |
73 else | |
74 { | |
75 LOG(ERROR) << "SlicesSorter::SliceWithDepth::GetPayload(): (!HasPayload())"; | |
76 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
77 } | |
78 } | |
79 }; | |
80 | |
81 | |
82 struct SlicesSorter::Comparator | |
83 { | |
84 bool operator() (const SliceWithDepth* const& a, | |
85 const SliceWithDepth* const& b) const | |
86 { | |
87 return a->GetDepth() < b->GetDepth(); | |
88 } | |
89 }; | |
90 | |
91 | |
92 SlicesSorter::~SlicesSorter() | |
93 { | |
94 for (size_t i = 0; i < slices_.size(); i++) | |
95 { | |
96 assert(slices_[i] != NULL); | |
97 delete slices_[i]; | |
98 } | |
99 } | |
100 | |
101 | |
102 void SlicesSorter::AddSlice(const CoordinateSystem3D& slice, | |
103 Orthanc::IDynamicObject* payload) | |
104 { | |
105 slices_.push_back(new SliceWithDepth(slice, payload)); | |
106 } | |
107 | |
108 | |
109 const SlicesSorter::SliceWithDepth& SlicesSorter::GetSlice(size_t i) const | |
110 { | |
111 if (i >= slices_.size()) | |
112 { | |
113 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
114 } | |
115 else | |
116 { | |
117 assert(slices_[i] != NULL); | |
118 return *slices_[i]; | |
119 } | |
120 } | |
121 | |
122 | |
123 const CoordinateSystem3D& SlicesSorter::GetSliceGeometry(size_t i) const | |
124 { | |
125 return GetSlice(i).GetGeometry(); | |
126 } | |
127 | |
128 | |
129 bool SlicesSorter::HasSlicePayload(size_t i) const | |
130 { | |
131 return GetSlice(i).HasPayload(); | |
132 } | |
133 | |
134 | |
135 const Orthanc::IDynamicObject& SlicesSorter::GetSlicePayload(size_t i) const | |
136 { | |
137 return GetSlice(i).GetPayload(); | |
138 } | |
139 | |
140 | |
141 void SlicesSorter::SetNormal(const Vector& normal) | |
142 { | |
143 for (size_t i = 0; i < slices_.size(); i++) | |
144 { | |
145 slices_[i]->SetNormal(normal); | |
146 } | |
147 | |
148 hasNormal_ = true; | |
149 } | |
150 | |
151 | |
152 void SlicesSorter::SortInternal() | |
153 { | |
154 if (!hasNormal_) | |
155 { | |
156 LOG(ERROR) << "SlicesSorter::SortInternal(): (!hasNormal_)"; | |
157 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
158 } | |
159 | |
160 Comparator comparator; | |
161 std::sort(slices_.begin(), slices_.end(), comparator); | |
162 } | |
163 | |
164 | |
165 void SlicesSorter::FilterNormal(const Vector& normal) | |
166 { | |
167 size_t pos = 0; | |
168 | |
169 for (size_t i = 0; i < slices_.size(); i++) | |
170 { | |
171 if (GeometryToolbox::IsParallel(normal, slices_[i]->GetGeometry().GetNormal())) | |
172 { | |
173 // This slice is compatible with the selected normal | |
174 slices_[pos] = slices_[i]; | |
175 pos += 1; | |
176 } | |
177 else | |
178 { | |
179 delete slices_[i]; | |
180 slices_[i] = NULL; | |
181 } | |
182 } | |
183 | |
184 slices_.resize(pos); | |
185 } | |
186 | |
187 | |
188 bool SlicesSorter::SelectNormal(Vector& normal) const | |
189 { | |
190 std::vector<Vector> normalCandidates; | |
191 std::vector<unsigned int> normalCount; | |
192 | |
193 bool found = false; | |
194 | |
195 for (size_t i = 0; !found && i < GetSlicesCount(); i++) | |
196 { | |
197 const Vector& normal = GetSlice(i).GetGeometry().GetNormal(); | |
198 | |
199 bool add = true; | |
200 for (size_t j = 0; add && j < normalCandidates.size(); j++) // (*) | |
201 { | |
202 if (GeometryToolbox::IsParallel(normal, normalCandidates[j])) | |
203 { | |
204 normalCount[j] += 1; | |
205 add = false; | |
206 } | |
207 } | |
208 | |
209 if (add) | |
210 { | |
211 if (normalCount.size() > 2) | |
212 { | |
213 // To get linear-time complexity in (*). This heuristics | |
214 // allows the series to have one single frame that is | |
215 // not parallel to the others (such a frame could be a | |
216 // generated preview) | |
217 found = false; | |
218 } | |
219 else | |
220 { | |
221 normalCandidates.push_back(normal); | |
222 normalCount.push_back(1); | |
223 } | |
224 } | |
225 } | |
226 | |
227 for (size_t i = 0; !found && i < normalCandidates.size(); i++) | |
228 { | |
229 unsigned int count = normalCount[i]; | |
230 if (count == GetSlicesCount() || | |
231 count + 1 == GetSlicesCount()) | |
232 { | |
233 normal = normalCandidates[i]; | |
234 found = true; | |
235 } | |
236 } | |
237 | |
238 return found; | |
239 } | |
240 | |
241 | |
242 bool SlicesSorter::Sort() | |
243 { | |
244 if (GetSlicesCount() > 0) | |
245 { | |
246 Vector normal; | |
247 if (SelectNormal(normal)) | |
248 { | |
249 FilterNormal(normal); | |
250 SetNormal(normal); | |
251 SortInternal(); | |
252 return true; | |
253 } | |
254 } | |
255 | |
256 return false; | |
257 } | |
258 | |
259 | |
260 bool SlicesSorter::LookupClosestSlice(size_t& index, | |
261 double& distance, | |
262 const CoordinateSystem3D& slice) const | |
263 { | |
264 // TODO Turn this linear-time lookup into a log-time lookup, | |
265 // keeping track of whether the slices are sorted along the normal | |
266 | |
267 bool found = false; | |
268 | |
269 distance = std::numeric_limits<double>::infinity(); | |
270 | |
271 for (size_t i = 0; i < slices_.size(); i++) | |
272 { | |
273 assert(slices_[i] != NULL); | |
274 | |
275 double tmp; | |
276 if (CoordinateSystem3D::ComputeDistance(tmp, slices_[i]->GetGeometry(), slice)) | |
277 { | |
278 if (!found || | |
279 tmp < distance) | |
280 { | |
281 index = i; | |
282 distance = tmp; | |
283 found = true; | |
284 } | |
285 } | |
286 } | |
287 | |
288 return found; | |
289 } | |
290 | |
291 | |
292 bool SlicesSorter::ComputeSpacingBetweenSlices(double& spacing /* out */) const | |
293 { | |
294 if (GetSlicesCount() <= 1) | |
295 { | |
296 // This is a volume that is empty or that contains one single | |
297 // slice: Choose a dummy z-dimension for voxels | |
298 spacing = 1.0; | |
299 return true; | |
300 } | |
301 | |
302 const OrthancStone::CoordinateSystem3D& reference = GetSliceGeometry(0); | |
303 | |
304 double referencePosition = reference.ProjectAlongNormal(reference.GetOrigin()); | |
305 | |
306 double p = reference.ProjectAlongNormal(GetSliceGeometry(1).GetOrigin()); | |
307 spacing = p - referencePosition; | |
308 | |
309 if (spacing <= 0) | |
310 { | |
311 LOG(ERROR) << "SlicesSorter::ComputeSpacingBetweenSlices(): (spacing <= 0)"; | |
312 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls, | |
313 "Please call the Sort() method before"); | |
314 } | |
315 | |
316 for (size_t i = 1; i < GetSlicesCount(); i++) | |
317 { | |
318 OrthancStone::Vector p = reference.GetOrigin() + spacing * static_cast<double>(i) * reference.GetNormal(); | |
319 double d = boost::numeric::ublas::norm_2(p - GetSliceGeometry(i).GetOrigin()); | |
320 | |
321 if (!OrthancStone::LinearAlgebra::IsNear(d, 0, 0.001 /* tolerance expressed in mm */)) | |
322 { | |
323 return false; | |
324 } | |
325 } | |
326 | |
327 return true; | |
328 } | |
329 | |
330 | |
331 bool SlicesSorter::AreAllSlicesDistinct() const | |
332 { | |
333 if (GetSlicesCount() <= 1) | |
334 { | |
335 return true; | |
336 } | |
337 else | |
338 { | |
339 const OrthancStone::CoordinateSystem3D& reference = GetSliceGeometry(0); | |
340 double previousPosition = reference.ProjectAlongNormal(GetSliceGeometry(0).GetOrigin()); | |
341 | |
342 for (size_t i = 1; i < GetSlicesCount(); i++) | |
343 { | |
344 double position = reference.ProjectAlongNormal(GetSliceGeometry(i).GetOrigin()); | |
345 | |
346 if (OrthancStone::LinearAlgebra::IsNear(position, previousPosition, 0.001 /* tolerance expressed in mm */)) | |
347 { | |
348 return false; | |
349 } | |
350 | |
351 previousPosition = position; | |
352 } | |
353 | |
354 return true; | |
355 } | |
356 } | |
357 } |