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 }