Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/DicomStructureSet.cpp @ 1909:782ba9eb6f22

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improved variable names
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 01 Feb 2022 08:56:36 +0100
parents affde38b84de
children f81cdf283859
comparison
equal deleted inserted replaced
1908:affde38b84de 1909:782ba9eb6f22
241 return true; 241 return true;
242 } 242 }
243 } 243 }
244 } 244 }
245 245
246 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& slice, 246 bool DicomStructureSet::Polygon::IsOnSlice(const CoordinateSystem3D& cuttingPlane,
247 const Vector& estimatedNormal, 247 const Vector& estimatedNormal,
248 double estimatedSliceThickness) const 248 double estimatedSliceThickness) const
249 { 249 {
250 bool isOpposite = false; 250 bool isOpposite = false;
251 251
254 return false; 254 return false;
255 } 255 }
256 else if (hasSlice_) 256 else if (hasSlice_)
257 { 257 {
258 // Use the actual geometry of this specific slice 258 // Use the actual geometry of this specific slice
259 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), geometry_.GetNormal())) 259 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, cuttingPlane.GetNormal(), geometry_.GetNormal()))
260 { 260 {
261 return false; 261 return false;
262 } 262 }
263 else 263 else
264 { 264 {
265 double d = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), geometry_.GetNormal()); 265 double d = GeometryToolbox::ProjectAlongNormal(cuttingPlane.GetOrigin(), geometry_.GetNormal());
266 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, sliceThickness_ / 2.0)); 266 return (LinearAlgebra::IsNear(d, projectionAlongNormal_, sliceThickness_ / 2.0));
267 } 267 }
268 } 268 }
269 else 269 else
270 { 270 {
271 // Use the estimated geometry for the global RT-STRUCT volume 271 // Use the estimated geometry for the global RT-STRUCT volume
272 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, slice.GetNormal(), estimatedNormal)) 272 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, cuttingPlane.GetNormal(), estimatedNormal))
273 { 273 {
274 return false; 274 return false;
275 } 275 }
276 else 276 else
277 { 277 {
278 double d1 = GeometryToolbox::ProjectAlongNormal(slice.GetOrigin(), estimatedNormal); 278 double d1 = GeometryToolbox::ProjectAlongNormal(cuttingPlane.GetOrigin(), estimatedNormal);
279 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal); 279 double d2 = GeometryToolbox::ProjectAlongNormal(points_.front(), estimatedNormal);
280 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0)); 280 return (LinearAlgebra::IsNear(d1, d2, estimatedSliceThickness / 2.0));
281 } 281 }
282 } 282 }
283 } 283 }
284 284
285 bool DicomStructureSet::Polygon::Project(double& x1, 285 bool DicomStructureSet::Polygon::Project(double& x1,
286 double& y1, 286 double& y1,
287 double& x2, 287 double& x2,
288 double& y2, 288 double& y2,
289 const CoordinateSystem3D& slice, 289 const CoordinateSystem3D& cuttingPlane,
290 const Vector& estimatedNormal, 290 const Vector& estimatedNormal,
291 double estimatedSliceThickness) const 291 double estimatedSliceThickness) const
292 { 292 {
293 #if 0
294 if (points_.size() <= 1)
295 {
296 return false;
297 }
298
299 Vector normal = estimatedNormal;
300 double thickness = estimatedSliceThickness;
301 if (hasSlice_)
302 {
303 normal = geometry_.GetNormal();
304 thickness = sliceThickness_;
305 }
306
307 bool isOpposite;
308 if (!GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisX()) &&
309 !GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cuttingPlane.GetAxisY()))
310 {
311 printf("UUUU\n");
312 return false;
313 }
314
315 return false;
316
317 #else
293 if (!hasSlice_ || 318 if (!hasSlice_ ||
294 points_.size() <= 1) 319 points_.size() <= 1)
295 { 320 {
296 return false; 321 return false;
297 } 322 }
298 323
299 double x, y; 324 double x, y;
300 geometry_.ProjectPoint2(x, y, slice.GetOrigin()); 325 geometry_.ProjectPoint2(x, y, cuttingPlane.GetOrigin());
301 326
302 bool isOpposite; 327 bool isOpposite;
303 if (GeometryToolbox::IsParallelOrOpposite 328 if (GeometryToolbox::IsParallelOrOpposite
304 (isOpposite, slice.GetNormal(), geometry_.GetAxisY())) 329 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisY()))
305 { 330 {
306 // plane is constant Y 331 // plane is constant Y
307 332
308 if (y < extent_.GetY1() || 333 if (y < extent_.GetY1() ||
309 y > extent_.GetY2()) 334 y > extent_.GetY2())
359 sliceThickness_ / 2.0 * geometry_.GetNormal()); 384 sliceThickness_ / 2.0 * geometry_.GetNormal());
360 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) - 385 Vector p2 = (geometry_.MapSliceToWorldCoordinates(xmax, y) -
361 sliceThickness_ / 2.0 * geometry_.GetNormal()); 386 sliceThickness_ / 2.0 * geometry_.GetNormal());
362 387
363 // then to the cutting plane geometry... 388 // then to the cutting plane geometry...
364 slice.ProjectPoint2(x1, y1, p1); 389 cuttingPlane.ProjectPoint2(x1, y1, p1);
365 slice.ProjectPoint2(x2, y2, p2); 390 cuttingPlane.ProjectPoint2(x2, y2, p2);
366 return true; 391 return true;
367 } 392 }
368 } 393 }
369 else if (GeometryToolbox::IsParallelOrOpposite 394 else if (GeometryToolbox::IsParallelOrOpposite
370 (isOpposite, slice.GetNormal(), geometry_.GetAxisX())) 395 (isOpposite, cuttingPlane.GetNormal(), geometry_.GetAxisX()))
371 { 396 {
372 // plane is constant X => Sagittal view (remember that in the 397 // plane is constant X => Sagittal view (remember that in the
373 // sagittal projection, the normal must be swapped) 398 // sagittal projection, the normal must be swapped)
374 399
375 400
420 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) + 445 Vector p1 = (geometry_.MapSliceToWorldCoordinates(x, ymin) +
421 sliceThickness_ / 2.0 * geometry_.GetNormal()); 446 sliceThickness_ / 2.0 * geometry_.GetNormal());
422 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) - 447 Vector p2 = (geometry_.MapSliceToWorldCoordinates(x, ymax) -
423 sliceThickness_ / 2.0 * geometry_.GetNormal()); 448 sliceThickness_ / 2.0 * geometry_.GetNormal());
424 449
425 slice.ProjectPoint2(x1, y1, p1); 450 cuttingPlane.ProjectPoint2(x1, y1, p1);
426 slice.ProjectPoint2(x2, y2, p2); 451 cuttingPlane.ProjectPoint2(x2, y2, p2);
427 452
428 return true; 453 return true;
429 } 454 }
430 } 455 }
431 else 456 else
432 { 457 {
433 // Should not happen 458 // Should not happen
434 return false; 459 return false;
435 } 460 }
461 #endif
436 } 462 }
437 463
438 464
439 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const 465 const DicomStructureSet::Structure& DicomStructureSet::GetStructure(size_t index) const
440 { 466 {
812 } 838 }
813 } 839 }
814 840
815 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains, 841 bool DicomStructureSet::ProjectStructure(std::vector< std::vector<ScenePoint2D> >& chains,
816 const Structure& structure, 842 const Structure& structure,
817 const CoordinateSystem3D& sourceSlice) const 843 const CoordinateSystem3D& cuttingPlane) const
818 { 844 {
819 const CoordinateSystem3D slice = CoordinateSystem3D::NormalizeCuttingPlane(sourceSlice); 845 const CoordinateSystem3D cutting = CoordinateSystem3D::NormalizeCuttingPlane(cuttingPlane);
820 846
821 chains.clear(); 847 chains.clear();
822 848
823 Vector normal = GetNormal(); 849 Vector normal = GetNormal();
824 850
825 bool isOpposite; 851 bool isOpposite;
826 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetNormal())) 852 if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetNormal()))
827 { 853 {
828 // This is an axial projection 854 // This is an axial projection
829 855
830 chains.reserve(structure.polygons_.size()); 856 chains.reserve(structure.polygons_.size());
831 857
832 for (Polygons::const_iterator polygon = structure.polygons_.begin(); 858 for (Polygons::const_iterator polygon = structure.polygons_.begin();
833 polygon != structure.polygons_.end(); ++polygon) 859 polygon != structure.polygons_.end(); ++polygon)
834 { 860 {
835 const Points& points = polygon->GetPoints(); 861 const Points& points = polygon->GetPoints();
836 862
837 if (polygon->IsOnSlice(slice, GetEstimatedNormal(), GetEstimatedSliceThickness()) && 863 if (polygon->IsOnSlice(cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()) &&
838 !points.empty()) 864 !points.empty())
839 { 865 {
840 chains.push_back(std::vector<ScenePoint2D>()); 866 chains.push_back(std::vector<ScenePoint2D>());
841 chains.back().reserve(points.size() + 1); 867 chains.back().reserve(points.size() + 1);
842 868
843 for (Points::const_iterator p = points.begin(); 869 for (Points::const_iterator p = points.begin();
844 p != points.end(); ++p) 870 p != points.end(); ++p)
845 { 871 {
846 double x, y; 872 double x, y;
847 slice.ProjectPoint2(x, y, *p); 873 cutting.ProjectPoint2(x, y, *p);
848 chains.back().push_back(ScenePoint2D(x, y)); 874 chains.back().push_back(ScenePoint2D(x, y));
849 } 875 }
850 876
851 double x0, y0; 877 double x0, y0;
852 slice.ProjectPoint2(x0, y0, points.front()); 878 cutting.ProjectPoint2(x0, y0, points.front());
853 chains.back().push_back(ScenePoint2D(x0, y0)); 879 chains.back().push_back(ScenePoint2D(x0, y0));
854 } 880 }
855 } 881 }
856 882
857 return true; 883 return true;
858 } 884 }
859 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || 885 else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetAxisX()) ||
860 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) 886 GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, cutting.GetAxisY()))
861 { 887 {
862 // Sagittal or coronal projection 888 // Sagittal or coronal projection
863 889
864 #if USE_BOOST_UNION_FOR_POLYGONS == 1 890 #if USE_BOOST_UNION_FOR_POLYGONS == 1
865 std::vector<BoostPolygon> projected; 891 std::vector<BoostPolygon> projected;
867 for (Polygons::const_iterator polygon = structure.polygons_.begin(); 893 for (Polygons::const_iterator polygon = structure.polygons_.begin();
868 polygon != structure.polygons_.end(); ++polygon) 894 polygon != structure.polygons_.end(); ++polygon)
869 { 895 {
870 double x1, y1, x2, y2; 896 double x1, y1, x2, y2;
871 897
872 if (polygon->Project(x1, y1, x2, y2, slice, GetEstimatedNormal(), GetEstimatedSliceThickness())) 898 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()))
873 { 899 {
874 projected.push_back(CreateRectangle(x1, y1, x2, y2)); 900 projected.push_back(CreateRectangle(x1, y1, x2, y2));
875 } 901 }
876 } 902 }
877 903
897 for (Polygons::const_iterator polygon = structure.polygons_.begin(); 923 for (Polygons::const_iterator polygon = structure.polygons_.begin();
898 polygon != structure.polygons_.end(); ++polygon) 924 polygon != structure.polygons_.end(); ++polygon)
899 { 925 {
900 double x1, y1, x2, y2; 926 double x1, y1, x2, y2;
901 927
902 if (polygon->Project(x1, y1, x2, y2, slice, GetEstimatedNormal(), GetEstimatedSliceThickness())) 928 if (polygon->Project(x1, y1, x2, y2, cutting, GetEstimatedNormal(), GetEstimatedSliceThickness()))
903 { 929 {
904 rectangles.push_back(Extent2D(x1, y1, x2, y2)); 930 rectangles.push_back(Extent2D(x1, y1, x2, y2));
905 } 931 }
906 } 932 }
907 933
927 } 953 }
928 } 954 }
929 955
930 956
931 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer, 957 void DicomStructureSet::ProjectOntoLayer(PolylineSceneLayer& layer,
932 const CoordinateSystem3D& plane, 958 const CoordinateSystem3D& cuttingPlane,
933 size_t structureIndex, 959 size_t structureIndex,
934 const Color& color) const 960 const Color& color) const
935 { 961 {
936 std::vector< std::vector<ScenePoint2D> > chains; 962 std::vector< std::vector<ScenePoint2D> > chains;
937 963
938 if (ProjectStructure(chains, structureIndex, plane)) 964 if (ProjectStructure(chains, structureIndex, cuttingPlane))
939 { 965 {
940 for (size_t j = 0; j < chains.size(); j++) 966 for (size_t j = 0; j < chains.size(); j++)
941 { 967 {
942 layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue()); 968 layer.AddChain(chains[j], false, color.GetRed(), color.GetGreen(), color.GetBlue());
943 } 969 }