Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructure2.cpp @ 998:38b6bb0bdd72

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added a new set of classes that correctly handle non-convex polygons (not used yet because of limitations in coordinates computing): DicomStructure2, DicomStructureSet2, DicomStructurePolygon2, DicomStructureSetSlicer2. Too many shortcuts have been taken when computing the actual position.
author Benjamin Golinvaux <bgo@osimis.io>
date Fri, 20 Sep 2019 11:58:00 +0200
parents
children 29f5f2031310
comparison
equal deleted inserted replaced
995:9893fa8cd7a6 998:38b6bb0bdd72
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-2019 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 #include "DicomStructure2.h"
22
23 #include "../Toolbox/GeometryToolbox.h"
24 #include "../Toolbox/DisjointDataSet.h"
25
26 namespace OrthancStone
27 {
28 // see header
29 //void DicomStructure2::ComputeNormal()
30 //{
31 // try
32 // {
33 // if (polygons_.size() > 0)
34 // {
35
36 // // TODO: check all polygons are OK
37 // const DicomStructurePolygon2 polygon = polygons_[0];
38 // $$$$$$$$$$$$$$$$$
39 // state_ = NormalComputed;
40 // }
41 // else
42 // {
43 // // bogus! no polygons. Let's assign a "nothing here" value
44 // LinearAlgebra::AssignVector(normal_, 0, 0, 0);
45 // state_ = Invalid;
46 // }
47 // }
48 // catch (const Orthanc::OrthancException& e)
49 // {
50 // state_ = Invalid;
51 // if (e.HasDetails())
52 // {
53 // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What() << " Details: " << e.GetDetails();
54 // }
55 // else
56 // {
57 // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What();
58 // }
59 // throw;
60 // }
61 // catch (const std::exception& e)
62 // {
63 // state_ = Invalid;
64 // LOG(ERROR) << "std::exception in ComputeNormal: " << e.what();
65 // throw;
66 // }
67 // catch (...)
68 // {
69 // state_ = Invalid;
70 // LOG(ERROR) << "Unknown exception in ComputeNormal";
71 // throw;
72 // }
73 //}
74
75 void DicomStructure2::ComputeSliceThickness()
76 {
77 if (state_ != NormalComputed)
78 {
79 LOG(ERROR) << "DicomStructure2::ComputeSliceThickness - state must be NormalComputed";
80 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
81 }
82 if (polygons_.size() < 2)
83 {
84 // cannot compute thickness if there are not at least 2 slabs (contours)
85 sliceThickness_ = 1.0;
86 state_ = Invalid;
87 }
88 else
89 {
90 // normal can be (1,0,0), (0,1,0) or (0,0,1), nothing else.
91 // these can be compared with == (exact double representation)
92 if (normal_[0] == 1)
93 {
94 // in a single polygon, all the points have the same X
95 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[0] - polygons_[1].GetPoint(0)[0]);
96 }
97 else if (normal_[1] == 1)
98 {
99 // in a single polygon, all the points have the same X
100 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[1] - polygons_[1].GetPoint(0)[1]);
101 }
102 else if (normal_[2] == 1)
103 {
104 // in a single polygon, all the points have the same X
105 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[2] - polygons_[1].GetPoint(0)[2]);
106 }
107 else
108 {
109 ORTHANC_ASSERT(false);
110 state_ = Invalid;
111 }
112 }
113 state_ = Valid;
114 }
115
116 void DicomStructure2::AddPolygon(const DicomStructurePolygon2& polygon)
117 {
118 if (state_ != Building)
119 {
120 LOG(ERROR) << "DicomStructure2::AddPolygon - can only add polygon while building";
121 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
122 }
123 polygons_.push_back(polygon);
124 }
125
126 void DicomStructure2::ComputeDependentProperties()
127 {
128 if (state_ != Building)
129 {
130 LOG(ERROR) << "DicomStructure2::ComputeDependentProperties - can only be called once";
131 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
132 }
133 for (size_t i = 0; i < polygons_.size(); ++i)
134 {
135 // "compute" the polygon normal
136 polygons_[i].ComputeDependentProperties();
137 }
138 if (polygons_.size() > 0)
139 {
140 normal_ = polygons_[0].GetNormal();
141 state_ = NormalComputed;
142 }
143 else
144 {
145 LinearAlgebra::AssignVector(normal_, 0, 0, 0);
146 state_ = Invalid; // THIS MAY HAPPEN !!! (for instance for instance 72c773ac-5059f2c4-2e6a9120-4fd4bca1-45701661 :) )
147 }
148 if (polygons_.size() >= 2)
149 ComputeSliceThickness(); // this will change state_ from NormalComputed to Valid
150 }
151
152 OrthancStone::Vector DicomStructure2::GetNormal() const
153 {
154 if (state_ != Valid && state_ != Invalid)
155 {
156 LOG(ERROR) << "DicomStructure2::GetNormal() -- please call ComputeDependentProperties first.";
157 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
158 }
159 if (state_ == Invalid)
160 {
161 LOG(ERROR) << "DicomStructure2::GetNormal() -- The Dicom structure is invalid. The normal is set to 0,0,0";
162 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
163 }
164 return normal_;
165 }
166
167 const DicomStructurePolygon2* DicomStructure2::GetPolygonClosestToSlice(
168 const CoordinateSystem3D& plane) const
169 {
170 ORTHANC_ASSERT(state_ == Valid);
171
172 // we assume 0,0,1 for now
173 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0));
174 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0));
175
176 for (size_t i = 0; i < polygons_.size(); ++i)
177 {
178 const DicomStructurePolygon2& polygon = polygons_[i];
179
180 // "height" of cutting plane
181 double cutZ = plane.GetOrigin()[2];
182
183 if (LinearAlgebra::IsNear(
184 cutZ, polygon.GetZ(),
185 sliceThickness_ / 2.0 /* in mm */))
186 return &polygon;
187 }
188 return NULL;
189 }
190
191
192 bool DicomStructure2::Project(std::vector< std::pair<Point2D, Point2D> > & segments, const CoordinateSystem3D & plane) const
193 {
194 segments.clear();
195
196 Vector normal = GetNormal();
197
198 size_t totalRectCount = 0;
199
200 // dummy var
201 bool isOpposite = false;
202
203 // This is an axial projection
204 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, plane.GetNormal()))
205 {
206 const DicomStructurePolygon2* polygon = GetPolygonClosestToSlice(plane);
207 if (polygon)
208 {
209 polygon->ProjectOnParallelPlane(segments, plane);
210 }
211 }
212 else
213 {
214 // let's compute the dot product of the plane normal and the polygons
215 // normal.
216 double dot = LinearAlgebra::DotProduct(plane.GetNormal(), normal);
217
218 if (LinearAlgebra::IsNear(dot, 0))
219 {
220 // Coronal or sagittal projection
221
222 // vector of vector of rectangles that will be merged in a single big contour:
223
224 // each polygon slab cut by a perpendicular plane yields 0..* rectangles
225 std::vector< RtStructRectanglesInSlab > rectanglesForEachSlab;
226
227 for (size_t i = 0; i < polygons_.size(); ++i)
228 {
229 // book an entry for this slab
230 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab());
231
232 // let's compute the intersection between the polygon and the plane
233 // intersections are in plane coords
234 std::vector<Point2D> intersections;
235
236 polygons_[i].ProjectOnConstantPlane(intersections, plane);
237
238 // for each pair of intersections, we add a rectangle.
239 if ((intersections.size() % 2) != 0)
240 {
241 LOG(WARNING) << "Odd number of intersections between structure "
242 << name_ << ", polygon # " << i
243 << " and plane where X axis is parallel to polygon normal vector";
244 }
245
246 size_t numRects = intersections.size() / 2;
247
248 // we keep count of the total number of rects for vector pre-allocations
249 totalRectCount += numRects;
250
251 for (size_t iRect = 0; iRect < numRects; ++iRect)
252 {
253 RtStructRectangleInSlab rectangle;
254 ORTHANC_ASSERT(LinearAlgebra::IsNear(intersections[2 * iRect].y, intersections[2 * iRect + 1].y));
255 ORTHANC_ASSERT((2 * iRect + 1) < intersections.size());
256 double x1 = intersections[2 * iRect].x;
257 double x2 = intersections[2 * iRect + 1].x;
258 double y1 = intersections[2 * iRect].y - sliceThickness_ * 0.5;
259 double y2 = intersections[2 * iRect].y + sliceThickness_ * 0.5;
260
261 rectangle.xmin = std::min(x1, x2);
262 rectangle.xmax = std::max(x1, x2);
263 rectangle.ymin = std::min(y1, y2);
264 rectangle.ymax = std::max(y1, y2);
265
266 // TODO: keep them sorted!!!!
267
268 rectanglesForEachSlab.back().push_back(rectangle);
269 }
270 }
271 // now we need to merge all the slabs into a set of polygons (1 or more)
272 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, totalRectCount);
273 }
274 else
275 {
276 // plane is not perpendicular to the polygons
277 // 180.0 / [Math]::Pi = 57.2957795130823
278 double acDot = 57.2957795130823 * acos(dot);
279 LOG(ERROR) << "DicomStructure2::Project -- cutting plane must be "
280 << "perpendicular to the contours, but dot product is: "
281 << dot << " and (180/pi)*acos(dot) = " << acDot;
282 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
283 }
284 }
285 return segments.size() != 0;
286 }
287 }