Mercurial > hg > orthanc-stone
diff OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp @ 1890:6ce81914f7e4
added classes BucketAccumulator1D/2D
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 18 Jan 2022 22:08:55 +0100 |
parents | a2955abe4c2e |
children | 3716d72161d2 |
line wrap: on
line diff
--- a/OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp Tue Jan 18 17:52:43 2022 +0100 +++ b/OrthancStone/UnitTestsSources/ComputationalGeometryTests.cpp Tue Jan 18 22:08:55 2022 +0100 @@ -25,6 +25,7 @@ #include "../Sources/Toolbox/Internals/OrientedIntegerLine2D.h" #include "../Sources/Toolbox/Internals/RectanglesIntegerProjection.h" +#include "../Sources/Toolbox/LinearAlgebra.h" #include "../Sources/Toolbox/SegmentTree.h" #include "../Sources/Toolbox/UnionOfRectangles.h" @@ -973,3 +974,416 @@ ASSERT_TRUE(contours.front()[12].IsEqual(OrthancStone::ScenePoint2D(2, 5))); } } + + +TEST(LinearAlgebra, ComputeMedian) +{ + { + std::vector<double> v; + ASSERT_THROW(OrthancStone::LinearAlgebra::ComputeMedian(v), Orthanc::OrthancException); + + v.push_back(1); + v.push_back(3); + v.push_back(3); + v.push_back(6); + v.push_back(7); + v.push_back(8); + v.push_back(9); + ASSERT_DOUBLE_EQ(6.0, OrthancStone::LinearAlgebra::ComputeMedian(v)); + } + + { + std::vector<double> v; + v.push_back(1); + v.push_back(2); + v.push_back(3); + v.push_back(4); + v.push_back(5); + v.push_back(6); + v.push_back(8); + v.push_back(9); + ASSERT_DOUBLE_EQ(4.5, OrthancStone::LinearAlgebra::ComputeMedian(v)); + } +} + + + + + + +namespace OrthancStone +{ + namespace Internals + { + class BucketMapper : public boost::noncopyable + { + private: + double minValue_; + double maxValue_; + size_t bucketsCount_; + + public: + BucketMapper(double minValue, + double maxValue, + size_t bucketsCount) : + minValue_(minValue), + maxValue_(maxValue), + bucketsCount_(bucketsCount) + { + if (minValue >= maxValue || + bucketsCount <= 0) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + } + + size_t GetSize() const + { + return bucketsCount_; + } + + void CheckIndex(size_t i) const + { + if (i >= bucketsCount_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + } + + double GetBucketLow(size_t i) const + { + CheckIndex(i); + double alpha = static_cast<double>(i) / static_cast<double>(bucketsCount_); + return (1.0 - alpha) * minValue_ + alpha * maxValue_; + } + + double GetBucketHigh(size_t i) const + { + CheckIndex(i); + double alpha = static_cast<double>(i + 1) / static_cast<double>(bucketsCount_); + return (1.0 - alpha) * minValue_ + alpha * maxValue_; + } + + double GetBucketCenter(size_t i) const + { + CheckIndex(i); + double alpha = (static_cast<double>(i) + 0.5) / static_cast<double>(bucketsCount_); + return (1.0 - alpha) * minValue_ + alpha * maxValue_; + } + + size_t GetBucketIndex(double value) const + { + if (value < minValue_ || + value > maxValue_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + else + { + double tmp = (value - minValue_) / (maxValue_ - minValue_) * static_cast<double>(bucketsCount_); + assert(tmp >= 0 && tmp <= static_cast<double>(bucketsCount_)); + + size_t bucket = static_cast<unsigned int>(std::floor(tmp)); + if (bucket == bucketsCount_) // This is the case if "value == maxValue_" + { + return bucketsCount_ - 1; + } + else + { + return bucket; + } + } + } + }; + } + + + class BucketAccumulator1D : public boost::noncopyable + { + private: + struct Bucket + { + size_t count_; + std::list<double> values_; + }; + + Internals::BucketMapper mapper_; + std::vector<Bucket> buckets_; + bool storeValues_; + + public: + BucketAccumulator1D(double minValue, + double maxValue, + size_t countBuckets, + bool storeValues) : + mapper_(minValue, maxValue, countBuckets), + buckets_(countBuckets), + storeValues_(storeValues) + { + } + + size_t GetSize() const + { + return mapper_.GetSize(); + } + + double GetBucketLow(size_t i) const + { + return mapper_.GetBucketLow(i); + } + + double GetBucketHigh(size_t i) const + { + return mapper_.GetBucketHigh(i); + } + + double GetBucketCenter(size_t i) const + { + return mapper_.GetBucketCenter(i); + } + + size_t GetBucketContentSize(size_t i) const + { + mapper_.CheckIndex(i); + return buckets_[i].count_; + } + + void AddValue(double value) + { + Bucket& bucket = buckets_[mapper_.GetBucketIndex(value)]; + + bucket.count_++; + + if (storeValues_) + { + bucket.values_.push_back(value); + } + } + + size_t FindBestBucket() const + { + size_t best = 0; + + for (size_t i = 0; i < buckets_.size(); i++) + { + if (buckets_[i].count_ > buckets_[best].count_) + { + best = i; + } + } + + return best; + } + + double ComputeBestCenter() const + { + return GetBucketCenter(FindBestBucket()); + } + + double ComputeBestMedian() const + { + if (!storeValues_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + + const std::list<double>& values = buckets_[FindBestBucket()].values_; + + std::vector<double> v; + v.reserve(values.size()); + for (std::list<double>::const_iterator it = values.begin(); it != values.end(); ++it) + { + v.push_back(*it); + } + + return LinearAlgebra::ComputeMedian(v); + } + }; + + + class BucketAccumulator2D : public boost::noncopyable + { + private: + struct Bucket + { + size_t count_; + std::list<double> valuesX_; + std::list<double> valuesY_; + }; + + Internals::BucketMapper mapperX_; + Internals::BucketMapper mapperY_; + std::vector<Bucket> buckets_; + bool storeValues_; + + size_t FindBestInternal() const + { + size_t best = 0; + + for (size_t i = 0; i < buckets_.size(); i++) + { + if (buckets_[i].count_ > buckets_[best].count_) + { + best = i; + } + } + + return best; + } + + public: + BucketAccumulator2D(double minValueX, + double maxValueX, + size_t countBucketsX, + double minValueY, + double maxValueY, + size_t countBucketsY, + bool storeValues) : + mapperX_(minValueX, maxValueX, countBucketsX), + mapperY_(minValueY, maxValueY, countBucketsY), + buckets_(countBucketsX * countBucketsY), + storeValues_(storeValues) + { + } + + void AddValue(double valueX, + double valueY) + { + size_t x = mapperX_.GetBucketIndex(valueX); + size_t y = mapperY_.GetBucketIndex(valueY); + + Bucket& bucket = buckets_[x + y * mapperX_.GetSize()]; + + bucket.count_++; + + if (storeValues_) + { + bucket.valuesX_.push_back(valueX); + bucket.valuesY_.push_back(valueY); + } + } + + void FindBestBucket(size_t& bucketX, + size_t& bucketY) const + { + size_t best = FindBestInternal(); + bucketX = best % mapperX_.GetSize(); + bucketY = best / mapperX_.GetSize(); + } + + void ComputeBestCenter(double& x, + double& y) const + { + size_t bucketX, bucketY; + FindBestBucket(bucketX, bucketY); + x = mapperX_.GetBucketCenter(bucketX); + y = mapperY_.GetBucketCenter(bucketY); + } + + void ComputeBestMedian(double& x, + double& y) const + { + if (!storeValues_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + + const std::list<double>& valuesX = buckets_[FindBestInternal()].valuesX_; + const std::list<double>& valuesY = buckets_[FindBestInternal()].valuesY_; + + std::vector<double> v; + v.reserve(valuesX.size()); + for (std::list<double>::const_iterator it = valuesX.begin(); it != valuesX.end(); ++it) + { + v.push_back(*it); + } + + x = LinearAlgebra::ComputeMedian(v); + + v.clear(); + v.reserve(valuesY.size()); + for (std::list<double>::const_iterator it = valuesY.begin(); it != valuesY.end(); ++it) + { + v.push_back(*it); + } + + y = LinearAlgebra::ComputeMedian(v); + } + }; +} + + +TEST(BucketAccumulator1D, Basic) +{ + for (int storeValues = 0; storeValues <= 1; storeValues++) + { + OrthancStone::BucketAccumulator1D b(-10, 30, 4, storeValues != 0); + ASSERT_EQ(4u, b.GetSize()); + + ASSERT_DOUBLE_EQ(-10.0, b.GetBucketLow(0)); + ASSERT_DOUBLE_EQ(0.0, b.GetBucketLow(1)); + ASSERT_DOUBLE_EQ(10.0, b.GetBucketLow(2)); + ASSERT_DOUBLE_EQ(20.0, b.GetBucketLow(3)); + + ASSERT_DOUBLE_EQ(0.0, b.GetBucketHigh(0)); + ASSERT_DOUBLE_EQ(10.0, b.GetBucketHigh(1)); + ASSERT_DOUBLE_EQ(20.0, b.GetBucketHigh(2)); + ASSERT_DOUBLE_EQ(30.0, b.GetBucketHigh(3)); + + ASSERT_DOUBLE_EQ(-5.0, b.GetBucketCenter(0)); + ASSERT_DOUBLE_EQ(5.0, b.GetBucketCenter(1)); + ASSERT_DOUBLE_EQ(15.0, b.GetBucketCenter(2)); + ASSERT_DOUBLE_EQ(25.0, b.GetBucketCenter(3)); + + ASSERT_EQ(0u, b.GetBucketContentSize(0)); + ASSERT_EQ(0u, b.GetBucketContentSize(1)); + ASSERT_EQ(0u, b.GetBucketContentSize(2)); + ASSERT_EQ(0u, b.GetBucketContentSize(3)); + + ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); // No data point + + b.AddValue(-10.0); + b.AddValue(0.0); + b.AddValue(9.9999); + b.AddValue(10.0); + b.AddValue(20.0); + b.AddValue(29.9999); + b.AddValue(30.0); + ASSERT_THROW(b.AddValue(30.00001), Orthanc::OrthancException); + + ASSERT_EQ(3u, b.FindBestBucket()); + ASSERT_EQ(25.0, b.ComputeBestCenter()); + + ASSERT_EQ(1u, b.GetBucketContentSize(0)); + ASSERT_EQ(2u, b.GetBucketContentSize(1)); + ASSERT_EQ(1u, b.GetBucketContentSize(2)); + ASSERT_EQ(3u, b.GetBucketContentSize(3)); + + if (storeValues == 0) + { + ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); + } + else + { + ASSERT_EQ(29.9999, b.ComputeBestMedian()); + } + } +} + +TEST(BucketAccumulator2D, Basic) +{ + for (int storeValues = 0; storeValues <= 1; storeValues++) + { + OrthancStone::BucketAccumulator2D b(-10, 30, 4, + 0, 3, 3, storeValues != 0); + + if (storeValues == 0) + { + //ASSERT_THROW(b.ComputeBestMedian(), Orthanc::OrthancException); + } + else + { + //ASSERT_EQ(29.9999, b.ComputeBestMedian()); + } + } +}