Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructurePolygon2.cpp @ 1512:244ad1e4e76a

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