Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructure2.cpp @ 1512:244ad1e4e76a

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