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