Mercurial > hg > orthanc-stone
comparison OrthancStone/Resources/Graveyard/RTStructTentativeReimplementation-BGO/DicomStructure2.cpp @ 1908:affde38b84de
moved tentative bgo reimplementation of rt-struct into graveyard
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 01 Feb 2022 08:38:32 +0100 |
parents | OrthancStone/Sources/Toolbox/DicomStructure2.cpp@14c8f339d480 |
children | 07964689cb0b |
comparison
equal
deleted
inserted
replaced
1907:0208f99b8bde | 1908:affde38b84de |
---|---|
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-2022 Osimis S.A., Belgium | |
6 * Copyright (C) 2021-2022 Sebastien Jodogne, ICTEAM UCLouvain, Belgium | |
7 * | |
8 * This program is free software: you can redistribute it and/or | |
9 * modify it under the terms of the GNU Lesser General Public License | |
10 * as published by the Free Software Foundation, either version 3 of | |
11 * the License, or (at your option) any later version. | |
12 * | |
13 * This program is distributed in the hope that it will be useful, but | |
14 * WITHOUT ANY WARRANTY; without even the implied warranty of | |
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
16 * Lesser General Public License for more details. | |
17 * | |
18 * You should have received a copy of the GNU Lesser General Public | |
19 * License along with this program. If not, see | |
20 * <http://www.gnu.org/licenses/>. | |
21 **/ | |
22 | |
23 #ifdef BGO_ENABLE_DICOMSTRUCTURESETLOADER2 | |
24 | |
25 #include "DicomStructure2.h" | |
26 | |
27 #include "GeometryToolbox.h" | |
28 #include "DisjointDataSet.h" | |
29 | |
30 #include <Logging.h> | |
31 | |
32 namespace OrthancStone | |
33 { | |
34 // see header | |
35 //void DicomStructure2::ComputeNormal() | |
36 //{ | |
37 // try | |
38 // { | |
39 // if (polygons_.size() > 0) | |
40 // { | |
41 | |
42 // // TODO: check all polygons are OK | |
43 // const DicomStructurePolygon2 polygon = polygons_[0]; | |
44 // $$$$$$$$$$$$$$$$$ | |
45 // state_ = NormalComputed; | |
46 // } | |
47 // else | |
48 // { | |
49 // // bogus! no polygons. Let's assign a "nothing here" value | |
50 // LinearAlgebra::AssignVector(normal_, 0, 0, 0); | |
51 // state_ = Invalid; | |
52 // } | |
53 // } | |
54 // catch (const Orthanc::OrthancException& e) | |
55 // { | |
56 // state_ = Invalid; | |
57 // if (e.HasDetails()) | |
58 // { | |
59 // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What() << " Details: " << e.GetDetails(); | |
60 // } | |
61 // else | |
62 // { | |
63 // LOG(ERROR) << "OrthancException in ComputeNormal: " << e.What(); | |
64 // } | |
65 // throw; | |
66 // } | |
67 // catch (const std::exception& e) | |
68 // { | |
69 // state_ = Invalid; | |
70 // LOG(ERROR) << "std::exception in ComputeNormal: " << e.what(); | |
71 // throw; | |
72 // } | |
73 // catch (...) | |
74 // { | |
75 // state_ = Invalid; | |
76 // LOG(ERROR) << "Unknown exception in ComputeNormal"; | |
77 // throw; | |
78 // } | |
79 //} | |
80 | |
81 void DicomStructure2::ComputeSliceThickness() | |
82 { | |
83 if (state_ != NormalComputed) | |
84 { | |
85 LOG(ERROR) << "DicomStructure2::ComputeSliceThickness - state must be NormalComputed"; | |
86 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
87 } | |
88 if (polygons_.size() < 2) | |
89 { | |
90 // cannot compute thickness if there are not at least 2 slabs (structures) | |
91 sliceThickness_ = 1.0; | |
92 state_ = Invalid; | |
93 } | |
94 else | |
95 { | |
96 // normal can be (1,0,0), (0,1,0) or (0,0,1), nothing else. | |
97 // these can be compared with == (exact double representation) | |
98 if (normal_[0] == 1) | |
99 { | |
100 // in a single polygon, all the points have the same X | |
101 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[0] - polygons_[1].GetPoint(0)[0]); | |
102 } | |
103 else if (normal_[1] == 1) | |
104 { | |
105 // in a single polygon, all the points have the same X | |
106 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[1] - polygons_[1].GetPoint(0)[1]); | |
107 } | |
108 else if (normal_[2] == 1) | |
109 { | |
110 // in a single polygon, all the points have the same X | |
111 sliceThickness_ = fabs(polygons_[0].GetPoint(0)[2] - polygons_[1].GetPoint(0)[2]); | |
112 } | |
113 else | |
114 { | |
115 ORTHANC_ASSERT(false); | |
116 state_ = Invalid; | |
117 } | |
118 } | |
119 state_ = Valid; | |
120 } | |
121 | |
122 void DicomStructure2::AddPolygon(const DicomStructurePolygon2& polygon) | |
123 { | |
124 if (state_ != Building) | |
125 { | |
126 LOG(ERROR) << "DicomStructure2::AddPolygon - can only add polygon while building"; | |
127 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
128 } | |
129 polygons_.push_back(polygon); | |
130 } | |
131 | |
132 void DicomStructure2::ComputeDependentProperties() | |
133 { | |
134 if (state_ != Building) | |
135 { | |
136 LOG(ERROR) << "DicomStructure2::ComputeDependentProperties - can only be called once"; | |
137 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
138 } | |
139 for (size_t i = 0; i < polygons_.size(); ++i) | |
140 { | |
141 // "compute" the polygon normal | |
142 polygons_[i].ComputeDependentProperties(); | |
143 } | |
144 if (polygons_.size() > 0) | |
145 { | |
146 normal_ = polygons_[0].GetNormal(); | |
147 state_ = NormalComputed; | |
148 } | |
149 else | |
150 { | |
151 LinearAlgebra::AssignVector(normal_, 0, 0, 0); | |
152 state_ = Invalid; // THIS MAY HAPPEN !!! (for instance for instance 72c773ac-5059f2c4-2e6a9120-4fd4bca1-45701661 :) ) | |
153 } | |
154 if (polygons_.size() >= 2) | |
155 ComputeSliceThickness(); // this will change state_ from NormalComputed to Valid | |
156 } | |
157 | |
158 Vector DicomStructure2::GetNormal() const | |
159 { | |
160 if (state_ != Valid && state_ != Invalid) | |
161 { | |
162 LOG(ERROR) << "DicomStructure2::GetNormal() -- please call ComputeDependentProperties first."; | |
163 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
164 } | |
165 if (state_ == Invalid) | |
166 { | |
167 LOG(ERROR) << "DicomStructure2::GetNormal() -- The Dicom structure is invalid. The normal is set to 0,0,0"; | |
168 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); | |
169 } | |
170 return normal_; | |
171 } | |
172 | |
173 const DicomStructurePolygon2* DicomStructure2::GetPolygonClosestToSlice( | |
174 const CoordinateSystem3D& plane) const | |
175 { | |
176 ORTHANC_ASSERT(state_ == Valid); | |
177 | |
178 // we assume 0,0,1 for now | |
179 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)); | |
180 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)); | |
181 | |
182 for (size_t i = 0; i < polygons_.size(); ++i) | |
183 { | |
184 const DicomStructurePolygon2& polygon = polygons_[i]; | |
185 | |
186 // "height" of cutting plane | |
187 double cutZ = plane.GetOrigin()[2]; | |
188 | |
189 if (LinearAlgebra::IsNear( | |
190 cutZ, polygon.GetZ(), | |
191 sliceThickness_ / 2.0 /* in mm */)) | |
192 return &polygon; | |
193 } | |
194 return NULL; | |
195 } | |
196 | |
197 | |
198 bool DicomStructure2::Project(std::vector< std::pair<ScenePoint2D, ScenePoint2D> > & segments, const CoordinateSystem3D & plane) const | |
199 { | |
200 segments.clear(); | |
201 | |
202 Vector normal = GetNormal(); | |
203 | |
204 size_t totalRectCount = 0; | |
205 | |
206 // dummy var | |
207 bool isOpposite = false; | |
208 | |
209 // This is an axial projection | |
210 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, plane.GetNormal())) | |
211 { | |
212 const DicomStructurePolygon2* polygon = GetPolygonClosestToSlice(plane); | |
213 if (polygon) | |
214 { | |
215 polygon->ProjectOnParallelPlane(segments, plane); | |
216 } | |
217 } | |
218 else | |
219 { | |
220 // let's compute the dot product of the plane normal and the polygons | |
221 // normal. | |
222 double dot = LinearAlgebra::DotProduct(plane.GetNormal(), normal); | |
223 | |
224 if (LinearAlgebra::IsNear(dot, 0)) | |
225 { | |
226 // Coronal or sagittal projection | |
227 | |
228 // vector of vector of rectangles that will be merged in a single big contour: | |
229 | |
230 // each polygon slab cut by a perpendicular plane yields 0..* rectangles | |
231 std::vector< RtStructRectanglesInSlab > rectanglesForEachSlab; | |
232 | |
233 for (size_t i = 0; i < polygons_.size(); ++i) | |
234 { | |
235 // book an entry for this slab | |
236 rectanglesForEachSlab.push_back(RtStructRectanglesInSlab()); | |
237 | |
238 // let's compute the intersection between the polygon and the plane | |
239 // intersections are in plane coords | |
240 std::vector<ScenePoint2D> intersections; | |
241 | |
242 polygons_[i].ProjectOnConstantPlane(intersections, plane); | |
243 | |
244 // for each pair of intersections, we add a rectangle. | |
245 if ((intersections.size() % 2) != 0) | |
246 { | |
247 LOG(WARNING) << "Odd number of intersections between structure " | |
248 << name_ << ", polygon # " << i | |
249 << " and plane where X axis is parallel to polygon normal vector"; | |
250 } | |
251 | |
252 size_t numRects = intersections.size() / 2; | |
253 | |
254 // we keep count of the total number of rects for vector pre-allocations | |
255 totalRectCount += numRects; | |
256 | |
257 for (size_t iRect = 0; iRect < numRects; ++iRect) | |
258 { | |
259 RtStructRectangleInSlab rectangle; | |
260 ORTHANC_ASSERT(LinearAlgebra::IsNear(intersections[2 * iRect].GetY(), intersections[2 * iRect + 1].GetY())); | |
261 ORTHANC_ASSERT((2 * iRect + 1) < intersections.size()); | |
262 double x1 = intersections[2 * iRect].GetX(); | |
263 double x2 = intersections[2 * iRect + 1].GetX(); | |
264 double y1 = intersections[2 * iRect].GetY() - sliceThickness_ * 0.5; | |
265 double y2 = intersections[2 * iRect].GetY() + sliceThickness_ * 0.5; | |
266 | |
267 rectangle.xmin = std::min(x1, x2); | |
268 rectangle.xmax = std::max(x1, x2); | |
269 rectangle.ymin = std::min(y1, y2); | |
270 rectangle.ymax = std::max(y1, y2); | |
271 | |
272 // TODO: keep them sorted!!!! | |
273 | |
274 rectanglesForEachSlab.back().push_back(rectangle); | |
275 } | |
276 } | |
277 // now we need to merge all the slabs into a set of polygons (1 or more) | |
278 ConvertListOfSlabsToSegments(segments, rectanglesForEachSlab, totalRectCount); | |
279 } | |
280 else | |
281 { | |
282 // plane is not perpendicular to the polygons | |
283 // 180.0 / [Math]::Pi = 57.2957795130823 | |
284 double acDot = 57.2957795130823 * acos(dot); | |
285 LOG(ERROR) << "DicomStructure2::Project -- cutting plane must be " | |
286 << "perpendicular to the structures, but dot product is: " | |
287 << dot << " and (180/pi)*acos(dot) = " << acDot; | |
288 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); | |
289 } | |
290 } | |
291 return segments.size() != 0; | |
292 } | |
293 } | |
294 | |
295 #endif | |
296 // BGO_ENABLE_DICOMSTRUCTURESETLOADER2 | |
297 |