Mercurial > hg > orthanc-stone

comparison OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp @ 1891:3716d72161d2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization
author Sebastien Jodogne <s.jodogne@gmail.com>
date Wed, 19 Jan 2022 12:32:15 +0100
parents 6ce81914f7e4
children 925aaf49150c
comparison
equal deleted inserted replaced
1890:6ce81914f7e4 1891:3716d72161d2
21 **/ 21 **/
22 22
23 23
24 #include <gtest/gtest.h> 24 #include <gtest/gtest.h>
25 25
26 #include "../Sources/Toolbox/BucketAccumulator1D.h"
27 #include "../Sources/Toolbox/BucketAccumulator2D.h"
26 #include "../Sources/Toolbox/Internals/OrientedIntegerLine2D.h" 28 #include "../Sources/Toolbox/Internals/OrientedIntegerLine2D.h"
27 #include "../Sources/Toolbox/Internals/RectanglesIntegerProjection.h" 29 #include "../Sources/Toolbox/Internals/RectanglesIntegerProjection.h"
28 #include "../Sources/Toolbox/LinearAlgebra.h" 30 #include "../Sources/Toolbox/LinearAlgebra.h"
29 #include "../Sources/Toolbox/SegmentTree.h" 31 #include "../Sources/Toolbox/SegmentTree.h"
30 #include "../Sources/Toolbox/UnionOfRectangles.h" 32 #include "../Sources/Toolbox/UnionOfRectangles.h"
1005 ASSERT_DOUBLE_EQ(4.5, OrthancStone::LinearAlgebra::ComputeMedian(v)); 1007 ASSERT_DOUBLE_EQ(4.5, OrthancStone::LinearAlgebra::ComputeMedian(v));
1006 } 1008 }
1007 } 1009 }
1008 1010
1009 1011
1010
1011
1012
1013
1014 namespace OrthancStone
1015 {
1016 namespace Internals
1017 {
1018 class BucketMapper : public boost::noncopyable
1019 {
1020 private:
1021 double minValue_;
1022 double maxValue_;
1023 size_t bucketsCount_;
1024
1025 public:
1026 BucketMapper(double minValue,
1027 double maxValue,
1028 size_t bucketsCount) :
1029 minValue_(minValue),
1030 maxValue_(maxValue),
1031 bucketsCount_(bucketsCount)
1032 {
1033 if (minValue >= maxValue ||
1034 bucketsCount <= 0)
1035 {
1036 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
1037 }
1038 }
1039
1040 size_t GetSize() const
1041 {
1042 return bucketsCount_;
1043 }
1044
1045 void CheckIndex(size_t i) const
1046 {
1047 if (i >= bucketsCount_)
1048 {
1049 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
1050 }
1051 }
1052
1053 double GetBucketLow(size_t i) const
1054 {
1055 CheckIndex(i);
1056 double alpha = static_cast<double>(i) / static_cast<double>(bucketsCount_);
1057 return (1.0 - alpha) * minValue_ + alpha * maxValue_;
1058 }
1059
1060 double GetBucketHigh(size_t i) const
1061 {
1062 CheckIndex(i);
1063 double alpha = static_cast<double>(i + 1) / static_cast<double>(bucketsCount_);
1064 return (1.0 - alpha) * minValue_ + alpha * maxValue_;
1065 }
1066
1067 double GetBucketCenter(size_t i) const
1068 {
1069 CheckIndex(i);
1070 double alpha = (static_cast<double>(i) + 0.5) / static_cast<double>(bucketsCount_);
1071 return (1.0 - alpha) * minValue_ + alpha * maxValue_;
1072 }
1073
1074 size_t GetBucketIndex(double value) const
1075 {
1076 if (value < minValue_ ||
1077 value > maxValue_)
1078 {
1079 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
1080 }
1081 else
1082 {
1083 double tmp = (value - minValue_) / (maxValue_ - minValue_) * static_cast<double>(bucketsCount_);
1084 assert(tmp >= 0 && tmp <= static_cast<double>(bucketsCount_));
1085
1086 size_t bucket = static_cast<unsigned int>(std::floor(tmp));
1087 if (bucket == bucketsCount_) // This is the case if "value == maxValue_"
1088 {
1089 return bucketsCount_ - 1;
1090 }
1091 else
1092 {
1093 return bucket;
1094 }
1095 }
1096 }
1097 };
1098 }
1099
1100
1101 class BucketAccumulator1D : public boost::noncopyable
1102 {
1103 private:
1104 struct Bucket
1105 {
1106 size_t count_;
1107 std::list<double> values_;
1108 };
1109
1110 Internals::BucketMapper mapper_;
1111 std::vector<Bucket> buckets_;
1112 bool storeValues_;
1113
1114 public:
1115 BucketAccumulator1D(double minValue,
1116 double maxValue,
1117 size_t countBuckets,
1118 bool storeValues) :
1119 mapper_(minValue, maxValue, countBuckets),
1120 buckets_(countBuckets),
1121 storeValues_(storeValues)
1122 {
1123 }
1124
1125 size_t GetSize() const
1126 {
1127 return mapper_.GetSize();
1128 }
1129
1130 double GetBucketLow(size_t i) const
1131 {
1132 return mapper_.GetBucketLow(i);
1133 }
1134
1135 double GetBucketHigh(size_t i) const
1136 {
1137 return mapper_.GetBucketHigh(i);
1138 }
1139
1140 double GetBucketCenter(size_t i) const
1141 {
1142 return mapper_.GetBucketCenter(i);
1143 }
1144
1145 size_t GetBucketContentSize(size_t i) const
1146 {
1147 mapper_.CheckIndex(i);
1148 return buckets_[i].count_;
1149 }
1150
1151 void AddValue(double value)
1152 {
1153 Bucket& bucket = buckets_[mapper_.GetBucketIndex(value)];
1154
1155 bucket.count_++;
1156
1157 if (storeValues_)
1158 {
1159 bucket.values_.push_back(value);
1160 }
1161 }
1162
1163 size_t FindBestBucket() const
1164 {
1165 size_t best = 0;
1166
1167 for (size_t i = 0; i < buckets_.size(); i++)
1168 {
1169 if (buckets_[i].count_ > buckets_[best].count_)
1170 {
1171 best = i;
1172 }
1173 }
1174
1175 return best;
1176 }
1177
1178 double ComputeBestCenter() const
1179 {
1180 return GetBucketCenter(FindBestBucket());
1181 }
1182
1183 double ComputeBestMedian() const
1184 {
1185 if (!storeValues_)
1186 {
1187 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
1188 }
1189
1190 const std::list<double>& values = buckets_[FindBestBucket()].values_;
1191
1192 std::vector<double> v;
1193 v.reserve(values.size());
1194 for (std::list<double>::const_iterator it = values.begin(); it != values.end(); ++it)
1195 {
1196 v.push_back(*it);
1197 }
1198
1199 return LinearAlgebra::ComputeMedian(v);
1200 }
1201 };
1202
1203
1204 class BucketAccumulator2D : public boost::noncopyable
1205 {
1206 private:
1207 struct Bucket
1208 {
1209 size_t count_;
1210 std::list<double> valuesX_;
1211 std::list<double> valuesY_;
1212 };
1213
1214 Internals::BucketMapper mapperX_;
1215 Internals::BucketMapper mapperY_;
1216 std::vector<Bucket> buckets_;
1217 bool storeValues_;
1218
1219 size_t FindBestInternal() const
1220 {
1221 size_t best = 0;
1222
1223 for (size_t i = 0; i < buckets_.size(); i++)
1224 {
1225 if (buckets_[i].count_ > buckets_[best].count_)
1226 {
1227 best = i;
1228 }
1229 }
1230
1231 return best;
1232 }
1233
1234 public:
1235 BucketAccumulator2D(double minValueX,
1236 double maxValueX,
1237 size_t countBucketsX,
1238 double minValueY,
1239 double maxValueY,
1240 size_t countBucketsY,
1241 bool storeValues) :
1242 mapperX_(minValueX, maxValueX, countBucketsX),
1243 mapperY_(minValueY, maxValueY, countBucketsY),
1244 buckets_(countBucketsX * countBucketsY),
1245 storeValues_(storeValues)
1246 {
1247 }
1248
1249 void AddValue(double valueX,
1250 double valueY)
1251 {
1252 size_t x = mapperX_.GetBucketIndex(valueX);
1253 size_t y = mapperY_.GetBucketIndex(valueY);
1254
1255 Bucket& bucket = buckets_[x + y * mapperX_.GetSize()];
1256
1257 bucket.count_++;
1258
1259 if (storeValues_)
1260 {
1261 bucket.valuesX_.push_back(valueX);
1262 bucket.valuesY_.push_back(valueY);
1263 }
1264 }
1265
1266 void FindBestBucket(size_t& bucketX,
1267 size_t& bucketY) const
1268 {
1269 size_t best = FindBestInternal();
1270 bucketX = best % mapperX_.GetSize();
1271 bucketY = best / mapperX_.GetSize();
1272 }
1273
1274 void ComputeBestCenter(double& x,
1275 double& y) const
1276 {
1277 size_t bucketX, bucketY;
1278 FindBestBucket(bucketX, bucketY);
1279 x = mapperX_.GetBucketCenter(bucketX);
1280 y = mapperY_.GetBucketCenter(bucketY);
1281 }
1282
1283 void ComputeBestMedian(double& x,
1284 double& y) const
1285 {
1286 if (!storeValues_)
1287 {
1288 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls);
1289 }
1290
1291 const std::list<double>& valuesX = buckets_[FindBestInternal()].valuesX_;
1292 const std::list<double>& valuesY = buckets_[FindBestInternal()].valuesY_;
1293
1294 std::vector<double> v;
1295 v.reserve(valuesX.size());
1296 for (std::list<double>::const_iterator it = valuesX.begin(); it != valuesX.end(); ++it)
1297 {
1298 v.push_back(*it);
1299 }
1300
1301 x = LinearAlgebra::ComputeMedian(v);
1302
1303 v.clear();
1304 v.reserve(valuesY.size());
1305 for (std::list<double>::const_iterator it = valuesY.begin(); it != valuesY.end(); ++it)
1306 {
1307 v.push_back(*it);
1308 }
1309
1310 y = LinearAlgebra::ComputeMedian(v);
1311 }
1312 };
1313 }
1314
1315
1316 TEST(BucketAccumulator1D, Basic) 1012 TEST(BucketAccumulator1D, Basic)
1317 { 1013 {
1318 for (int storeValues = 0; storeValues <= 1; storeValues++) 1014 for (int storeValues = 0; storeValues <= 1; storeValues++)
1319 { 1015 {
1320 OrthancStone::BucketAccumulator1D b(-10, 30, 4, storeValues != 0); 1016 OrthancStone::BucketAccumulator1D b(-10, 30, 4, storeValues != 0);
1338 ASSERT_EQ(0u, b.GetBucketContentSize(0)); 1034 ASSERT_EQ(0u, b.GetBucketContentSize(0));
1339 ASSERT_EQ(0u, b.GetBucketContentSize(1)); 1035 ASSERT_EQ(0u, b.GetBucketContentSize(1));
1340 ASSERT_EQ(0u, b.GetBucketContentSize(2)); 1036 ASSERT_EQ(0u, b.GetBucketContentSize(2));
1341 ASSERT_EQ(0u, b.GetBucketContentSize(3)); 1037 ASSERT_EQ(0u, b.GetBucketContentSize(3));
1342 1038
1039 ASSERT_THROW(b.GetBucketIndex(-10.0001), Orthanc::OrthancException);
1040 ASSERT_EQ(0u, b.GetBucketIndex(-10));
1041 ASSERT_EQ(0u, b.GetBucketIndex(-0.0001));
1042 ASSERT_EQ(1u, b.GetBucketIndex(0));
1043 ASSERT_EQ(1u, b.GetBucketIndex(9.9999));
1044 ASSERT_EQ(2u, b.GetBucketIndex(10));
1045 ASSERT_EQ(2u, b.GetBucketIndex(19.9999));
1046 ASSERT_EQ(3u, b.GetBucketIndex(20));
1047 ASSERT_EQ(3u, b.GetBucketIndex(30));
1048 ASSERT_THROW(b.GetBucketIndex(30.0001), Orthanc::OrthancException);
1049
1050 ASSERT_EQ(0u, b.FindBestBucket());
1051 ASSERT_DOUBLE_EQ(-5.0, b.ComputeBestCenter());
1343 ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); // No data point 1052 ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); // No data point
1344 1053
1345 b.AddValue(-10.0); 1054 b.AddValue(-10.0);
1346 b.AddValue(0.0); 1055 b.AddValue(0.0);
1347 b.AddValue(9.9999); 1056 b.AddValue(9.9999);
1348 b.AddValue(10.0); 1057 b.AddValue(10.0);
1349 b.AddValue(20.0); 1058 b.AddValue(20.0);
1350 b.AddValue(29.9999); 1059 b.AddValue(29.9999);
1351 b.AddValue(30.0); 1060 b.AddValue(30.0);
1061 ASSERT_THROW(b.AddValue(-10.00001), Orthanc::OrthancException);
1352 ASSERT_THROW(b.AddValue(30.00001), Orthanc::OrthancException); 1062 ASSERT_THROW(b.AddValue(30.00001), Orthanc::OrthancException);
1353 1063
1354 ASSERT_EQ(3u, b.FindBestBucket()); 1064 ASSERT_EQ(3u, b.FindBestBucket());
1355 ASSERT_EQ(25.0, b.ComputeBestCenter()); 1065 ASSERT_EQ(3u, b.GetBucketContentSize(b.FindBestBucket()));
1066 ASSERT_DOUBLE_EQ(25.0, b.ComputeBestCenter());
1356 1067
1357 ASSERT_EQ(1u, b.GetBucketContentSize(0)); 1068 ASSERT_EQ(1u, b.GetBucketContentSize(0));
1358 ASSERT_EQ(2u, b.GetBucketContentSize(1)); 1069 ASSERT_EQ(2u, b.GetBucketContentSize(1));
1359 ASSERT_EQ(1u, b.GetBucketContentSize(2)); 1070 ASSERT_EQ(1u, b.GetBucketContentSize(2));
1360 ASSERT_EQ(3u, b.GetBucketContentSize(3)); 1071 ASSERT_EQ(3u, b.GetBucketContentSize(3));
1363 { 1074 {
1364 ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); 1075 ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException);
1365 } 1076 }
1366 else 1077 else
1367 { 1078 {
1368 ASSERT_EQ(29.9999, b.ComputeBestMedian()); 1079 ASSERT_DOUBLE_EQ(29.9999, b.ComputeBestMedian());
1369 } 1080 }
1370 } 1081 }
1371 } 1082 }
1372 1083
1373 TEST(BucketAccumulator2D, Basic) 1084 TEST(BucketAccumulator2D, Basic)
1374 { 1085 {
1375 for (int storeValues = 0; storeValues <= 1; storeValues++) 1086 for (int storeValues = 0; storeValues <= 1; storeValues++)
1376 { 1087 {
1377 OrthancStone::BucketAccumulator2D b(-10, 30, 4, 1088 OrthancStone::BucketAccumulator2D b(-10, 30, 4,
1378 0, 3, 3, storeValues != 0); 1089 0, 3, 3, storeValues != 0);
1379 1090
1091 size_t bx, by;
1092 b.FindBestBucket(bx, by);
1093 ASSERT_EQ(0u, bx);
1094 ASSERT_EQ(0u, by);
1095
1096 for (by = 0; by < 3; by++)
1097 {
1098 for (bx = 0; bx < 4; bx++)
1099 {
1100 ASSERT_EQ(0u, b.GetBucketContentSize(bx, by));
1101 }
1102 }
1103
1104 b.GetSize(bx, by);
1105 ASSERT_EQ(4u, bx);
1106 ASSERT_EQ(3u, by);
1107
1108 ASSERT_DOUBLE_EQ(-10.0, b.GetBucketLowX(0));
1109 ASSERT_DOUBLE_EQ(0.0, b.GetBucketLowX(1));
1110 ASSERT_DOUBLE_EQ(10.0, b.GetBucketLowX(2));
1111 ASSERT_DOUBLE_EQ(20.0, b.GetBucketLowX(3));
1112 ASSERT_THROW(b.GetBucketLowX(4), Orthanc::OrthancException);
1113
1114 ASSERT_DOUBLE_EQ(0.0, b.GetBucketLowY(0));
1115 ASSERT_DOUBLE_EQ(1.0, b.GetBucketLowY(1));
1116 ASSERT_DOUBLE_EQ(2.0, b.GetBucketLowY(2));
1117 ASSERT_THROW(b.GetBucketLowY(3), Orthanc::OrthancException);
1118
1119 ASSERT_DOUBLE_EQ(0.0, b.GetBucketHighX(0));
1120 ASSERT_DOUBLE_EQ(10.0, b.GetBucketHighX(1));
1121 ASSERT_DOUBLE_EQ(20.0, b.GetBucketHighX(2));
1122 ASSERT_DOUBLE_EQ(30.0, b.GetBucketHighX(3));
1123 ASSERT_THROW(b.GetBucketHighX(4), Orthanc::OrthancException);
1124
1125 ASSERT_DOUBLE_EQ(1.0, b.GetBucketHighY(0));
1126 ASSERT_DOUBLE_EQ(2.0, b.GetBucketHighY(1));
1127 ASSERT_DOUBLE_EQ(3.0, b.GetBucketHighY(2));
1128 ASSERT_THROW(b.GetBucketHighY(3), Orthanc::OrthancException);
1129
1130 ASSERT_DOUBLE_EQ(-5.0, b.GetBucketCenterX(0));
1131 ASSERT_DOUBLE_EQ(5.0, b.GetBucketCenterX(1));
1132 ASSERT_DOUBLE_EQ(15.0, b.GetBucketCenterX(2));
1133 ASSERT_DOUBLE_EQ(25.0, b.GetBucketCenterX(3));
1134 ASSERT_THROW(b.GetBucketCenterX(4), Orthanc::OrthancException);
1135
1136 ASSERT_DOUBLE_EQ(0.5, b.GetBucketCenterY(0));
1137 ASSERT_DOUBLE_EQ(1.5, b.GetBucketCenterY(1));
1138 ASSERT_DOUBLE_EQ(2.5, b.GetBucketCenterY(2));
1139 ASSERT_THROW(b.GetBucketCenterY(3), Orthanc::OrthancException);
1140
1141 b.GetBucketIndex(bx, by, 5, 2.5);
1142 ASSERT_EQ(1u, bx);
1143 ASSERT_EQ(2u, by);
1144 b.AddValue(4.5, 2.2);
1145 ASSERT_THROW(b.AddValue(-10.001, 2), Orthanc::OrthancException);
1146 ASSERT_THROW(b.AddValue(30.001, 2), Orthanc::OrthancException);
1147 ASSERT_THROW(b.AddValue(0, -0.0001), Orthanc::OrthancException);
1148 ASSERT_THROW(b.AddValue(0, 3.0001), Orthanc::OrthancException);
1149
1150 b.FindBestBucket(bx, by);
1151 ASSERT_EQ(1u, bx);
1152 ASSERT_EQ(2u, by);
1153
1154 for (by = 0; by < 3; by++)
1155 {
1156 for (bx = 0; bx < 4; bx++)
1157 {
1158 ASSERT_EQ((bx == 1u && by == 2u) ? 1u : 0u, b.GetBucketContentSize(bx, by));
1159 }
1160 }
1161
1162 double dx, dy;
1163 b.ComputeBestCenter(dx, dy);
1164 ASSERT_DOUBLE_EQ(5.0, dx);
1165 ASSERT_DOUBLE_EQ(2.5, dy);
1166
1380 if (storeValues == 0) 1167 if (storeValues == 0)
1381 { 1168 {
1382 //ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); 1169 ASSERT_THROW(b.ComputeBestMedian(dx, dy), Orthanc::OrthancException);
1383 } 1170 }
1384 else 1171 else
1385 { 1172 {
1386 //ASSERT_EQ(29.9999, b.ComputeBestMedian()); 1173 b.ComputeBestMedian(dx, dy);
1387 } 1174 ASSERT_DOUBLE_EQ(4.5, dx);
1388 } 1175 ASSERT_DOUBLE_EQ(2.2, dy);
1389 } 1176 }
1177 }
1178 }