Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructure2.cpp @ 998:38b6bb0bdd72
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 } |