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