Mercurial > hg > orthanc-stone

comparison OrthancStone/Resources/Graveyard/RTStructTentativeReimplementation-BGO/DicomStructure2.cpp @ 1908:affde38b84de

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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