Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/DicomStructurePolygon2.cpp @ 998:38b6bb0bdd72

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