Mercurial > hg > orthanc-stone

comparison UnitTestsSources/GeometryToolboxTests.cpp @ 1765:177f93c8371c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
unit test SlicesSorter.HFP
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 11 May 2021 11:01:22 +0200
parents 9ac2a65d4172
children ae18b5cfbd3e
comparison
equal deleted inserted replaced
1764:adb59bcb0fc5 1765:177f93c8371c
19 **/ 19 **/
20 20
21 21
22 #include <gtest/gtest.h> 22 #include <gtest/gtest.h>
23 23
24 #include "../OrthancStone/Sources/Toolbox/DicomInstanceParameters.h"
24 #include "../OrthancStone/Sources/Toolbox/FiniteProjectiveCamera.h" 25 #include "../OrthancStone/Sources/Toolbox/FiniteProjectiveCamera.h"
25 #include "../OrthancStone/Sources/Toolbox/GenericToolbox.h" 26 #include "../OrthancStone/Sources/Toolbox/GenericToolbox.h"
26 #include "../OrthancStone/Sources/Toolbox/GeometryToolbox.h" 27 #include "../OrthancStone/Sources/Toolbox/GeometryToolbox.h"
28 #include "../OrthancStone/Sources/Toolbox/SlicesSorter.h"
27 29
28 #include <Logging.h> 30 #include <Logging.h>
29 #include <OrthancException.h> 31 #include <OrthancException.h>
30 32
31 #include <boost/math/special_functions/round.hpp> 33 #include <boost/math/special_functions/round.hpp>
929 double r = d.GetNormal() [0] / 2.0; 931 double r = d.GetNormal() [0] / 2.0;
930 ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]); 932 ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]);
931 ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]); 933 ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]);
932 } 934 }
933 } 935 }
936
937
938 TEST(SlicesSorter, HFP)
939 {
940 // 2021-04-27-repro-bug-HFP-HFS-cartman
941
942 {
943 // This is the last instance in the CT series ("InstanceNumber" is 368):
944 // CT1.2.752.243.1.1.20210202150623868.3730.61448.dcm
945 const OrthancStone::CoordinateSystem3D system("300\\302.5\\323.11", "-1\\0\\0\\0\\-1\\0");
946
947 /**
948 * The first instance in the series ("InstanceNumber" is 1) is
949 * CT1.2.752.243.1.1.20210202150623381.2000.76318.dcm, and its
950 * "ImagePositionPatient" is "300\\302.5\\690.11". It cannot be
951 * taken as the origin of the volume, otherwise the Z axis is
952 * shifted by the depth of the volume.
953 **/
954
955 const double spacingXY = 1.171875;
956 const double width = 512;
957 const double height = 512;
958 const double depth = 368; // Number of instances in the series
959
960 OrthancStone::VolumeImageGeometry geometry;
961 geometry.SetAxialGeometry(system);
962 geometry.SetSizeInVoxels(width, height, depth);
963 geometry.SetVoxelDimensions(spacingXY, spacingXY, 1 /* pixel spacing Z */);
964
965 OrthancStone::Vector p, q;
966
967 p = OrthancStone::LinearAlgebra::CreateVector(0.5 / double(width),
968 0.5 / double(height),
969 0.5 / double(depth), 1);
970 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
971 ASSERT_FLOAT_EQ(300, q[0]);
972 ASSERT_FLOAT_EQ(302.5, q[1]);
973 ASSERT_FLOAT_EQ(323.11, q[2]);
974 ASSERT_FLOAT_EQ(1, q[3]);
975
976 p = OrthancStone::LinearAlgebra::CreateVector((width - 0.5) / double(width),
977 (height - 0.5) / double(height),
978 (depth - 0.5) / double(depth), 1);
979 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
980
981 ASSERT_FLOAT_EQ(300.0 - (width - 1.0) * spacingXY, q[0]); // "X" is swapped
982 ASSERT_FLOAT_EQ(302.5 - (height - 1.0) * spacingXY, q[1]); // "Y" is swapped
983 ASSERT_FLOAT_EQ(690.11, q[2]);
984 ASSERT_FLOAT_EQ(1, q[3]);
985 }
986
987 {
988 // DOSE instance: RD1.2.752.243.1.1.20210202150624529.3790.85357_DoseTPS.dcm
989 OrthancStone::CoordinateSystem3D system("-217.0492\\-161.4141\\376.61", "1\\0\\0\\0\\1\\0");
990 const double spacingXY = 3;
991 const double width = 146; // Columns
992 const double height = 84; // Row
993 const double depth = 86; // Number of frames, same as the length of "GridFrameOffsetVector"
994
995 OrthancStone::VolumeImageGeometry geometry;
996 geometry.SetAxialGeometry(system);
997 geometry.SetSizeInVoxels(width, height, depth);
998 geometry.SetVoxelDimensions(spacingXY, spacingXY, 3 /* pixel spacing Z, cf. "GridFrameOffsetVector" */);
999
1000 OrthancStone::Vector p, q;
1001
1002 p = OrthancStone::LinearAlgebra::CreateVector(0.5 / double(width), 0.5 / double(height), 0.5 / double(depth), 1);
1003 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1004 ASSERT_FLOAT_EQ(-217.0492, q[0]);
1005 ASSERT_FLOAT_EQ(-161.4141, q[1]);
1006 ASSERT_FLOAT_EQ(376.61, q[2]);
1007 ASSERT_FLOAT_EQ(1, q[3]);
1008
1009 p = OrthancStone::LinearAlgebra::CreateVector((width - 0.5) / double(width),
1010 (height - 0.5) / double(height),
1011 (depth - 0.5) / double(depth), 1);
1012 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1013
1014 ASSERT_FLOAT_EQ(-217.0492 + (width - 1.0) * spacingXY, q[0]);
1015 ASSERT_FLOAT_EQ(-161.4141 + (height - 1.0) * spacingXY, q[1]);
1016 ASSERT_FLOAT_EQ(376.61 + 255.0 /* last item in "GridFrameOffsetVector" */, q[2]);
1017 ASSERT_FLOAT_EQ(1, q[3]);
1018 }
1019
1020 for (unsigned int upward = 0; upward < 2; upward++)
1021 {
1022 OrthancStone::SlicesSorter slices;
1023
1024 for (unsigned int i = 0; i < 368; i++)
1025 {
1026 unsigned int z = (upward ? 323 + i : 690 - i);
1027 OrthancStone::CoordinateSystem3D p("300\\302.5\\" + boost::lexical_cast<std::string>(z) + ".11",
1028 "-1\\0\\0\\0\\-1\\0");
1029 slices.AddSlice(p, new Orthanc::SingleValueObject<unsigned int>(z));
1030 }
1031
1032 slices.Sort();
1033
1034 double spacingZ;
1035 ASSERT_TRUE(slices.ComputeSpacingBetweenSlices(spacingZ));
1036 ASSERT_FLOAT_EQ(1, spacingZ);
1037
1038 ASSERT_FLOAT_EQ(300.0, slices.GetSliceGeometry(0).GetOrigin() [0]);
1039 ASSERT_FLOAT_EQ(302.5, slices.GetSliceGeometry(0).GetOrigin() [1]);
1040 ASSERT_FLOAT_EQ(323.11, slices.GetSliceGeometry(0).GetOrigin() [2]);
1041
1042 ASSERT_FLOAT_EQ(300.0, slices.GetSliceGeometry(367).GetOrigin() [0]);
1043 ASSERT_FLOAT_EQ(302.5, slices.GetSliceGeometry(367).GetOrigin() [1]);
1044 ASSERT_FLOAT_EQ(690.11, slices.GetSliceGeometry(367).GetOrigin() [2]);
1045
1046 ASSERT_EQ(323u, dynamic_cast<const Orthanc::SingleValueObject<unsigned int>&>(slices.GetSlicePayload(0)).GetValue());
1047 ASSERT_EQ(690u, dynamic_cast<const Orthanc::SingleValueObject<unsigned int>&>(slices.GetSlicePayload(367)).GetValue());
1048
1049 const unsigned int width = 512;
1050 const unsigned int height = 512;
1051 const unsigned int depth = 368;
1052 const double spacingXY = 1.171875;
1053
1054 OrthancStone::VolumeImageGeometry geometry;
1055 geometry.SetSizeInVoxels(width, height, depth);
1056 geometry.SetVoxelDimensions(spacingXY, spacingXY, 1);
1057 geometry.SetAxialGeometry(slices.GetSliceGeometry(0));
1058
1059 OrthancStone::Vector p, q;
1060
1061 p = OrthancStone::LinearAlgebra::CreateVector(0, 0, 0, 1);
1062 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1063 ASSERT_FLOAT_EQ(300 + spacingXY / 2.0, q[0]);
1064 ASSERT_FLOAT_EQ(302.5 + spacingXY / 2.0, q[1]);
1065 ASSERT_FLOAT_EQ(323.11 - 0.5, q[2]);
1066 ASSERT_FLOAT_EQ(1, q[3]);
1067
1068 p = OrthancStone::LinearAlgebra::CreateVector(1, 1, 1, 1);
1069 q = OrthancStone::LinearAlgebra::Product(geometry.GetTransform(), p);
1070 ASSERT_FLOAT_EQ(300 - double(width) * spacingXY + spacingXY / 2.0, q[0]);
1071 ASSERT_FLOAT_EQ(302.5 - double(height) * spacingXY + spacingXY / 2.0, q[1]);
1072 ASSERT_FLOAT_EQ(323.11 + 368.0 - 0.5, q[2]);
1073 ASSERT_FLOAT_EQ(1, q[3]);
1074 }
1075 }