Mercurial > hg > orthanc-stone
comparison Framework/Toolbox/DicomStructurePolygon2.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 | 58eed6bbcabb |
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 "DicomStructurePolygon2.h" | |
22 | |
23 #include "../Toolbox/LinearAlgebra.h" | |
24 | |
25 namespace OrthancStone | |
26 { | |
27 void DicomStructurePolygon2::ComputeDependentProperties() | |
28 { | |
29 ORTHANC_ASSERT(state_ == Building); | |
30 | |
31 for (size_t j = 0; j < points_.size(); ++j) | |
32 { | |
33 // TODO: move to AddPoint! | |
34 const Point3D& p = points_[j]; | |
35 if (p[0] < minX_) | |
36 minX_ = p[0]; | |
37 if (p[0] > maxX_) | |
38 maxX_ = p[0]; | |
39 | |
40 if (p[1] < minY_) | |
41 minY_ = p[1]; | |
42 if (p[1] > maxY_) | |
43 maxY_ = p[1]; | |
44 | |
45 if (p[2] < minZ_) | |
46 minZ_ = p[2]; | |
47 if (p[2] > maxZ_) | |
48 maxZ_ = p[2]; | |
49 } | |
50 | |
51 if (LinearAlgebra::IsNear(minX_, maxX_)) | |
52 { | |
53 LinearAlgebra::AssignVector(normal_, 1, 0, 0); | |
54 //ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX, maxX)); | |
55 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_)); | |
56 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_)); | |
57 } | |
58 else if (LinearAlgebra::IsNear(minY_, maxY_)) | |
59 { | |
60 LinearAlgebra::AssignVector(normal_, 0, 1, 0); | |
61 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_)); | |
62 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minZ_, maxZ_)); | |
63 } | |
64 else if (LinearAlgebra::IsNear(minZ_, maxZ_)) | |
65 { | |
66 LinearAlgebra::AssignVector(normal_, 0, 0, 1); | |
67 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minX_, maxX_)); | |
68 ORTHANC_ASSERT(!LinearAlgebra::IsNear(minY_, maxY_)); | |
69 } | |
70 else | |
71 { | |
72 LOG(ERROR) << "The contour is not coplanar and not parallel to any axis."; | |
73 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
74 } | |
75 state_ = Valid; | |
76 } | |
77 | |
78 | |
79 void DicomStructurePolygon2::ProjectOnConstantPlane( | |
80 std::vector<Point2D>& intersections, const CoordinateSystem3D& plane) const | |
81 { | |
82 // the plane can either have constant X, or constant Y. | |
83 // - for constant Z planes, use the ProjectOnParallelPlane method | |
84 // - other type of planes are not supported | |
85 | |
86 // V is the coordinate that is constant in the plane | |
87 double planeV = 0.0; | |
88 | |
89 // if true, then "u" in the code is "x" and "v" is "y". | |
90 // (v is constant in the plane) | |
91 bool uvxy = false; | |
92 // if true, then "u" in the code is "y" and "v" is "x" | |
93 // (v is constant in the plane) | |
94 bool uvyx = false; | |
95 | |
96 size_t uindex = static_cast<size_t>(-1); | |
97 size_t vindex = static_cast<size_t>(-1); | |
98 | |
99 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[2], 0.0)); | |
100 | |
101 if (LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)) | |
102 { | |
103 // normal is 1,0,0 (or -1,0,0). | |
104 // plane is constant X | |
105 uindex = 1; | |
106 vindex = 0; | |
107 | |
108 uvxy = false; | |
109 uvyx = true; | |
110 planeV = plane.GetOrigin()[0]; | |
111 if (planeV < minX_) | |
112 return; | |
113 if (planeV > maxX_) | |
114 return; | |
115 } | |
116 else if (LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)) | |
117 { | |
118 // normal is 0,1,0 (or 0,-1,0). | |
119 // plane is constant Y | |
120 uindex = 0; | |
121 vindex = 1; | |
122 | |
123 uvxy = true; | |
124 uvyx = false; | |
125 planeV = plane.GetOrigin()[1]; | |
126 if (planeV < minY_) | |
127 return; | |
128 if (planeV > maxY_) | |
129 return; | |
130 } | |
131 else | |
132 { | |
133 // if the following assertion(s) fail(s), it means the plane is NOT a constant-X or constant-Y plane | |
134 LOG(ERROR) << "Plane normal must be (a,0,0) or (0,a,0), with a == -1 or a == 1"; | |
135 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); | |
136 } | |
137 | |
138 size_t pointCount = GetPointCount(); | |
139 if (pointCount >= 3) | |
140 { | |
141 // this vector will contain the coordinates of the intersection points | |
142 // between the plane and the polygon. | |
143 // these are expressed in the U coordinate, that is either X or Y, | |
144 // depending upon the plane orientation | |
145 std::vector<double> uIntersections; | |
146 | |
147 // we loop on the segments of the polygon (TODO: optimize) | |
148 // and we compute the intersection between each segment and the cut | |
149 // cutting plane (slice) has a constant X | |
150 | |
151 for (size_t iPoint = 0; iPoint < pointCount; ++iPoint) | |
152 { | |
153 double u1 = points_[iPoint][uindex]; | |
154 double v1 = points_[iPoint][vindex]; | |
155 | |
156 double u2 = 0; | |
157 double v2 = 0; | |
158 | |
159 if (iPoint < pointCount - 1) | |
160 { | |
161 u2 = points_[iPoint + 1][uindex]; | |
162 v2 = points_[iPoint + 1][vindex]; | |
163 } | |
164 else | |
165 { | |
166 u2 = points_[0][uindex]; | |
167 v2 = points_[0][vindex]; | |
168 } | |
169 | |
170 // Check if the segment intersects the plane | |
171 if ((std::min(v1, v2) <= planeV) && (std::max(v1, v2) >= planeV)) | |
172 { | |
173 // special case: the segment is parallel to the plane but close to it | |
174 if (LinearAlgebra::IsNear(v1, v2)) | |
175 { | |
176 // in that case, we choose to label both points as an intersection | |
177 double x, y; | |
178 plane.ProjectPoint(x, y, points_[iPoint]); | |
179 intersections.push_back(Point2D(x, y)); | |
180 | |
181 plane.ProjectPoint(x, y, points_[iPoint + 1]); | |
182 intersections.push_back(Point2D(x, y)); | |
183 } | |
184 else | |
185 { | |
186 // we are looking for u so that (u,planeV) belongs to the segment | |
187 // let's define alpha = (u-u2)/(u1-u2) --> u = alpha*(u1-u2) + u2 | |
188 // alpha = (v2-planeV)/(v2-v1) | |
189 // because the following two triangles are similar | |
190 // [ (planeY,x) , (y2,x2), (planeY,x2) ] or | |
191 // [ (planeX,y) , (x2,y2), (planeX,y2) ] | |
192 // and | |
193 // [ (y1 ,x1) , (y2,x2), (y1 ,x2) ] or | |
194 // [ (x1 ,y1) , (x2,y2), (x1 ,y2) ] | |
195 | |
196 /* | |
197 void CoordinateSystem3D::ProjectPoint(double& offsetX, | |
198 double& offsetY, | |
199 const Vector& point) const | |
200 */ | |
201 double alpha = (v2 - planeV) / (v2 - v1); | |
202 | |
203 // get rid of numerical oddities | |
204 if (alpha < 0.0) | |
205 alpha = 0.0; | |
206 if (alpha > 1.0) | |
207 alpha = 1.0; | |
208 double u = alpha * (u1 - u2) + u2; | |
209 | |
210 // here is the intersection in world coordinates | |
211 Vector intersection; | |
212 if(uvxy) | |
213 LinearAlgebra::AssignVector(intersection, u, planeV, minZ_); | |
214 else | |
215 LinearAlgebra::AssignVector(intersection, planeV, u, minZ_); | |
216 | |
217 // and we convert it to plane coordinates | |
218 { | |
219 double xi, yi; | |
220 plane.ProjectPoint(xi, yi, intersection); | |
221 | |
222 // we consider that the x axis is always parallel to the polygons | |
223 // TODO: is this hypothesis safe?????? | |
224 uIntersections.insert(std::lower_bound(uIntersections.begin(), uIntersections.end(), xi), xi); | |
225 } | |
226 } | |
227 } | |
228 } // end of for (size_t iPoint = 0; iPoint < pointCount; ++iPoint) | |
229 | |
230 // now we convert the intersections to plane points | |
231 // we consider that the x axis is always parallel to the polygons | |
232 // TODO: same hypothesis as above: plane is perpendicular to polygons, | |
233 // plane is parallel to the XZ (constant Y) or YZ (constant X) 3D planes | |
234 for (size_t i = 0; i < uIntersections.size(); ++i) | |
235 { | |
236 double x = uIntersections[i]; | |
237 intersections.push_back(Point2D(x, minZ_)); | |
238 } | |
239 } // end of if (pointCount >= 3) | |
240 else | |
241 { | |
242 LOG(ERROR) << "This polygon has " << pointCount << " vertices, which is less than 3 --> skipping"; | |
243 } | |
244 } | |
245 | |
246 void OrthancStone::DicomStructurePolygon2::ProjectOnParallelPlane( | |
247 std::vector< std::pair<Point2D, Point2D> >& segments, | |
248 const CoordinateSystem3D& plane) const | |
249 { | |
250 if (points_.size() < 3) | |
251 return; | |
252 | |
253 // the plane is horizontal | |
254 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[0], 0.0)); | |
255 ORTHANC_ASSERT(LinearAlgebra::IsNear(plane.GetNormal()[1], 0.0)); | |
256 | |
257 size_t pointCount = GetPointCount(); | |
258 | |
259 segments.clear(); | |
260 segments.reserve(points_.size()); | |
261 // since the returned values need to be expressed in the supplied coordinate | |
262 // system, we need to subtract origin_ from the returned points | |
263 | |
264 double planeOriginX = plane.GetOrigin()[0]; | |
265 double planeOriginY = plane.GetOrigin()[1]; | |
266 | |
267 // precondition: points_.size() >= 3 | |
268 for (size_t j = 0; j < points_.size()-1; ++j) | |
269 { | |
270 // segment between point j and j+1 | |
271 | |
272 const Point3D& point0 = GetPoint(j); | |
273 // subtract plane origin x and y | |
274 Point2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY); | |
275 | |
276 const Point3D& point1 = GetPoint(j+1); | |
277 // subtract plane origin x and y | |
278 Point2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY); | |
279 | |
280 segments.push_back(std::pair<Point2D, Point2D>(p0,p1)); | |
281 } | |
282 | |
283 | |
284 // final segment | |
285 | |
286 const Point3D& point0 = GetPoint(points_.size() - 1); | |
287 // subtract plane origin x and y | |
288 Point2D p0(point0[0] - planeOriginX, point0[1] - planeOriginY); | |
289 | |
290 const Point3D& point1 = GetPoint(0); | |
291 // subtract plane origin x and y | |
292 Point2D p1(point1[0] - planeOriginX, point1[1] - planeOriginY); | |
293 | |
294 segments.push_back(std::pair<Point2D, Point2D>(p0, p1)); | |
295 } | |
296 | |
297 double OrthancStone::DicomStructurePolygon2::GetZ() const | |
298 { | |
299 ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[0], 0.0)); | |
300 ORTHANC_ASSERT(LinearAlgebra::IsNear(normal_[1], 0.0)); | |
301 ORTHANC_ASSERT(LinearAlgebra::IsNear(minZ_, maxZ_)); | |
302 return minZ_; | |
303 } | |
304 | |
305 | |
306 } |