Mercurial > hg > orthanc-stone

comparison OrthancStone/Resources/Graveyard/RTStructTentativeReimplementation-BGO/DicomStructurePolygon2.cpp @ 1908:affde38b84de

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