# HG changeset patch # User Sebastien Jodogne # Date 1606219150 -3600 # Node ID 5b8b88e5bfd6451d668519914562bb5e304a1d44 # Parent 1393e3393a0b2fa911b68a34959da62fc0e76394 successfully running unit tests in WebAssembly diff -r 1393e3393a0b -r 5b8b88e5bfd6 Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp --- a/Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -1102,21 +1102,26 @@ { OrthancStone::DicomInstanceParameters params(dicom); - // TODO - WINDOWING - if (params.GetPresetWindowingsCount() > 0) + for (size_t i = 0; i < params.GetWindowingPresetsCount(); i++) { - GetViewport().presetWindowingCenter_ = params.GetPresetWindowingCenter(0); - GetViewport().presetWindowingWidth_ = params.GetPresetWindowingWidth(0); - LOG(INFO) << "Preset windowing: " << params.GetPresetWindowingCenter(0) - << "," << params.GetPresetWindowingWidth(0); - - GetViewport().windowingCenter_ = params.GetPresetWindowingCenter(0); - GetViewport().windowingWidth_ = params.GetPresetWindowingWidth(0); + LOG(INFO) << "Preset windowing " << i << "/" << params.GetWindowingPresetsCount() + << ": " << params.GetWindowingPresetCenter(i) + << "," << params.GetWindowingPresetWidth(i); + } + + // TODO - WINDOWING + if (params.GetWindowingPresetsCount() > 0) + { + GetViewport().presetWindowingCenter_ = params.GetWindowingPresetCenter(0); + GetViewport().presetWindowingWidth_ = params.GetWindowingPresetWidth(0); + + GetViewport().windowingCenter_ = params.GetWindowingPresetCenter(0); + GetViewport().windowingWidth_ = params.GetWindowingPresetWidth(0); } else { LOG(INFO) << "No preset windowing"; - GetViewport().ResetPresetWindowing(); + GetViewport().ResetWindowingPreset(); } } @@ -1411,7 +1416,7 @@ } - void ResetPresetWindowing() + void ResetWindowingPreset() { presetWindowingCenter_ = 128; presetWindowingWidth_ = 256; @@ -1790,7 +1795,7 @@ emscripten_set_keydown_callback(EMSCRIPTEN_EVENT_TARGET_WINDOW, this, false, OnKey); emscripten_set_keyup_callback(EMSCRIPTEN_EVENT_TARGET_WINDOW, this, false, OnKey); - ResetPresetWindowing(); + ResetWindowingPreset(); } static EM_BOOL OnKey(int eventType, @@ -1898,7 +1903,7 @@ LOG(INFO) << "Number of frames in series: " << frames_->GetFramesCount(); - ResetPresetWindowing(); + ResetWindowingPreset(); ClearViewport(); prefetchQueue_.clear(); @@ -2127,7 +2132,7 @@ } - void SetPresetWindowing() + void SetWindowingPreset() { SetWindowing(presetWindowingCenter_, presetWindowingWidth_); } @@ -2328,6 +2333,11 @@ { return cineRate_; } + + void FormatWindowingPresets(Json::Value& target) const + { + target = Json::arrayValue; + } }; @@ -2818,11 +2828,11 @@ EMSCRIPTEN_KEEPALIVE - void SetPresetWindowing(const char* canvas) + void SetWindowingPreset(const char* canvas) { try { - GetViewport(canvas)->SetPresetWindowing(); + GetViewport(canvas)->SetWindowingPreset(); } EXTERN_CATCH_EXCEPTIONS; } @@ -3017,5 +3027,18 @@ } EXTERN_CATCH_EXCEPTIONS; return 0; - } + } + + + EMSCRIPTEN_KEEPALIVE + void LoadWindowingPresets(const char* canvas) + { + try + { + Json::Value v; + GetViewport(canvas)->FormatWindowingPresets(v); + stringBuffer_ = v.toStyledString(); + } + EXTERN_CATCH_EXCEPTIONS; + } } diff -r 1393e3393a0b -r 5b8b88e5bfd6 OrthancStone/Sources/Loaders/SeriesFramesLoader.cpp --- a/OrthancStone/Sources/Loaders/SeriesFramesLoader.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/OrthancStone/Sources/Loaders/SeriesFramesLoader.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -307,77 +307,6 @@ } - static void GetBounds(float& low, - float& high, - double center, // in - double width) // in - { - low = static_cast(center - width / 2.0); - high = static_cast(center + width / 2.0); - } - - - void SeriesFramesLoader::GetPreviewWindowing(float& center, - float& width, - size_t index) const - { - const Orthanc::DicomMap& instance = frames_.GetInstance(index); - const DicomInstanceParameters& parameters = frames_.GetInstanceParameters(index); - - size_t s = parameters.GetPresetWindowingsCount(); - - if (s > 0) - { - // Use the largest windowing given all the preset windowings - // that are available in the DICOM tags - float low, high; - GetBounds(low, high, parameters.GetPresetWindowingCenter(0), - parameters.GetPresetWindowingWidth(0)); - - for (size_t i = 1; i < s; i++) - { - float a, b; - GetBounds(a, b, parameters.GetPresetWindowingCenter(i), - parameters.GetPresetWindowingWidth(i)); - low = std::min(low, a); - high = std::max(high, b); - } - - assert(low <= high); - - if (LinearAlgebra::IsNear(low, high)) - { - // Cannot infer a suitable windowing from the available tags - center = 128.0f; - width = 256.0f; - } - else - { - center = (low + high) / 2.0f; - width = (high - low); - printf(">> %f %f\n", center, width); - } - } - else - { - float a, b; - if (instance.ParseFloat(a, Orthanc::DICOM_TAG_SMALLEST_IMAGE_PIXEL_VALUE) && - instance.ParseFloat(b, Orthanc::DICOM_TAG_LARGEST_IMAGE_PIXEL_VALUE) && - a < b) - { - center = (a + b) / 2.0f; - width = (b - a); - } - else - { - // Cannot infer a suitable windowing from the available tags - center = 128.0f; - width = 256.0f; - } - } - } - - Orthanc::IDynamicObject& SeriesFramesLoader::FrameLoadedMessage::GetUserPayload() const { if (userPayload_) @@ -485,8 +414,10 @@ if (source.HasDicomWebRendered() && quality == 0) { + const DicomInstanceParameters& parameters = frames_.GetInstanceParameters(index); + float c, w; - GetPreviewWindowing(c, w, index); + parameters.GetWindowingPresetsUnion(c, w); std::map arguments, headers; arguments["window"] = (boost::lexical_cast(c) + "," + diff -r 1393e3393a0b -r 5b8b88e5bfd6 OrthancStone/Sources/Loaders/SeriesFramesLoader.h --- a/OrthancStone/Sources/Loaders/SeriesFramesLoader.h Tue Nov 24 07:40:19 2020 +0100 +++ b/OrthancStone/Sources/Loaders/SeriesFramesLoader.h Tue Nov 24 12:59:10 2020 +0100 @@ -76,10 +76,6 @@ void Handle(const HttpCommand::SuccessMessage& message); - void GetPreviewWindowing(float& center, - float& width, - size_t index) const; - public: class FrameLoadedMessage : public OriginMessage { diff -r 1393e3393a0b -r 5b8b88e5bfd6 OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp --- a/OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -182,10 +182,10 @@ bool ok = false; - if (LinearAlgebra::ParseVector(presetWindowingCenters_, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) && - LinearAlgebra::ParseVector(presetWindowingWidths_, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH)) + if (LinearAlgebra::ParseVector(windowingPresetCenters_, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) && + LinearAlgebra::ParseVector(windowingPresetWidths_, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH)) { - if (presetWindowingCenters_.size() == presetWindowingWidths_.size()) + if (windowingPresetCenters_.size() == windowingPresetWidths_.size()) { ok = true; } @@ -199,8 +199,8 @@ if (!ok) { // Don't use "Vector::clear()", as it has not the same meaning as "std::vector::clear()" - presetWindowingCenters_.resize(0); - presetWindowingWidths_.resize(0); + windowingPresetCenters_.resize(0); + windowingPresetWidths_.resize(0); } // This computes the "IndexInSeries" metadata from Orthanc (check @@ -397,18 +397,31 @@ } - size_t DicomInstanceParameters::GetPresetWindowingsCount() const + size_t DicomInstanceParameters::GetWindowingPresetsCount() const { - assert(data_.presetWindowingCenters_.size() == data_.presetWindowingWidths_.size()); - return data_.presetWindowingCenters_.size(); + assert(data_.windowingPresetCenters_.size() == data_.windowingPresetWidths_.size()); + return data_.windowingPresetCenters_.size(); } - float DicomInstanceParameters::GetPresetWindowingCenter(size_t i) const + float DicomInstanceParameters::GetWindowingPresetCenter(size_t i) const { - if (i < GetPresetWindowingsCount()) + if (i < GetWindowingPresetsCount()) + { + return static_cast(data_.windowingPresetCenters_[i]); + } + else { - return static_cast(data_.presetWindowingCenters_[i]); + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + } + + + float DicomInstanceParameters::GetWindowingPresetWidth(size_t i) const + { + if (i < GetWindowingPresetsCount()) + { + return static_cast(data_.windowingPresetWidths_[i]); } else { @@ -417,19 +430,71 @@ } - float DicomInstanceParameters::GetPresetWindowingWidth(size_t i) const + static void GetWindowingBounds(float& low, + float& high, + double center, // in + double width) // in + { + low = static_cast(center - width / 2.0); + high = static_cast(center + width / 2.0); + } + + + void DicomInstanceParameters::GetWindowingPresetsUnion(float& center, + float& width) const { - if (i < GetPresetWindowingsCount()) + assert(tags_.get() != NULL); + size_t s = GetWindowingPresetsCount(); + + if (s > 0) { - return static_cast(data_.presetWindowingWidths_[i]); + // Use the largest windowing given all the preset windowings + // that are available in the DICOM tags + float low, high; + GetWindowingBounds(low, high, GetWindowingPresetCenter(0), GetWindowingPresetWidth(0)); + + for (size_t i = 1; i < s; i++) + { + float a, b; + GetWindowingBounds(a, b, GetWindowingPresetCenter(i), GetWindowingPresetWidth(i)); + low = std::min(low, a); + high = std::max(high, b); + } + + assert(low <= high); + + if (LinearAlgebra::IsNear(low, high)) + { + // Cannot infer a suitable windowing from the available tags + center = 128.0f; + width = 256.0f; + } + else + { + center = (low + high) / 2.0f; + width = (high - low); + } } else { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + float a, b; + if (tags_->ParseFloat(a, Orthanc::DICOM_TAG_SMALLEST_IMAGE_PIXEL_VALUE) && + tags_->ParseFloat(b, Orthanc::DICOM_TAG_LARGEST_IMAGE_PIXEL_VALUE) && + a < b) + { + center = (a + b) / 2.0f; + width = (b - a); + } + else + { + // Cannot infer a suitable windowing from the available tags + center = 128.0f; + width = 256.0f; + } } } - + Orthanc::ImageAccessor* DicomInstanceParameters::ConvertToFloat(const Orthanc::ImageAccessor& pixelData) const { std::unique_ptr converted(new Orthanc::Image(Orthanc::PixelFormat_Float32, @@ -494,9 +559,9 @@ texture.reset(new FloatTextureSceneLayer(*converted)); } - if (GetPresetWindowingsCount() > 0) + if (GetWindowingPresetsCount() > 0) { - texture->SetCustomWindowing(GetPresetWindowingCenter(0), GetPresetWindowingWidth(0)); + texture->SetCustomWindowing(GetWindowingPresetCenter(0), GetWindowingPresetWidth(0)); } switch (GetImageInformation().GetPhotometricInterpretation()) diff -r 1393e3393a0b -r 5b8b88e5bfd6 OrthancStone/Sources/Toolbox/DicomInstanceParameters.h --- a/OrthancStone/Sources/Toolbox/DicomInstanceParameters.h Tue Nov 24 07:40:19 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/DicomInstanceParameters.h Tue Nov 24 12:59:10 2020 +0100 @@ -55,8 +55,8 @@ bool hasRescale_; double rescaleIntercept_; double rescaleSlope_; - Vector presetWindowingCenters_; - Vector presetWindowingWidths_; + Vector windowingPresetCenters_; + Vector windowingPresetWidths_; bool hasIndexInSeries_; unsigned int indexInSeries_; std::string doseUnits_; @@ -181,11 +181,14 @@ double GetRescaleSlope() const; - size_t GetPresetWindowingsCount() const; + size_t GetWindowingPresetsCount() const; + + float GetWindowingPresetCenter(size_t i) const; - float GetPresetWindowingCenter(size_t i) const; + float GetWindowingPresetWidth(size_t i) const; - float GetPresetWindowingWidth(size_t i) const; + void GetWindowingPresetsUnion(float& center, + float& width) const; Orthanc::PixelFormat GetExpectedPixelFormat() const; diff -r 1393e3393a0b -r 5b8b88e5bfd6 OrthancStone/Sources/Toolbox/GenericToolbox.h --- a/OrthancStone/Sources/Toolbox/GenericToolbox.h Tue Nov 24 07:40:19 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/GenericToolbox.h Tue Nov 24 12:59:10 2020 +0100 @@ -255,7 +255,7 @@ } while (*p >= '0' && *p <= '9') { - r = (r * 10) + (*p - '0'); // 1 12 123 123 12345 + r = (r * 10) + static_cast(*p - '0'); // 1 12 123 123 12345 ++p; } r *= neg; diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/CMakeLists.txt --- a/UnitTestsSources/CMakeLists.txt Tue Nov 24 07:40:19 2020 +0100 +++ b/UnitTestsSources/CMakeLists.txt Tue Nov 24 12:59:10 2020 +0100 @@ -40,24 +40,9 @@ set(ENABLE_PUGIXML ON) include(${ORTHANC_STONE_ROOT}/../Resources/CMake/OrthancStoneConfiguration.cmake) +include(${CMAKE_CURRENT_SOURCE_DIR}/UnitTestsSources.cmake) -add_executable(UnitTests - ${CMAKE_SOURCE_DIR}/GenericToolboxTests.cpp - ${CMAKE_SOURCE_DIR}/ImageToolboxTests.cpp - ${CMAKE_SOURCE_DIR}/PixelTestPatternsTests.cpp - ${CMAKE_SOURCE_DIR}/SortedFramesTests.cpp - ${CMAKE_SOURCE_DIR}/TestMessageBroker.cpp - ${CMAKE_SOURCE_DIR}/TestStrategy.cpp - ${CMAKE_SOURCE_DIR}/TestStructureSet.cpp - ${CMAKE_SOURCE_DIR}/UnitTestsMain.cpp - - ${AUTOGENERATED_SOURCES} - ${BOOST_EXTENDED_SOURCES} - ${GOOGLE_TEST_SOURCES} - ${ORTHANC_STONE_SOURCES} - ) - -target_link_libraries(UnitTests ${DCMTK_LIBRARIES}) +add_executable(UnitTests ${UNIT_TESTS_SOURCES}) ##################################################################### diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/DicomTests.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/DicomTests.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,109 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2020 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + **/ + + +#include + +#include "../OrthancStone/Sources/Toolbox/DicomInstanceParameters.h" + +#include + + +static void SetupUids(Orthanc::DicomMap& m) +{ + m.SetValue(Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, "my_study", false); + m.SetValue(Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, "my_series", false); + m.SetValue(Orthanc::DICOM_TAG_SOP_INSTANCE_UID, "my_sop", false); +} + + +TEST(DicomInstanceParameters, Basic) +{ + Orthanc::DicomMap m; + SetupUids(m); + + std::unique_ptr p; + p.reset(OrthancStone::DicomInstanceParameters(m).Clone()); + + ASSERT_TRUE(p->GetOrthancInstanceIdentifier().empty()); + ASSERT_EQ(3u, p->GetTags().GetSize()); + ASSERT_EQ("my_study", p->GetStudyInstanceUid()); + ASSERT_EQ("my_series", p->GetSeriesInstanceUid()); + ASSERT_EQ("my_sop", p->GetSopInstanceUid()); + ASSERT_EQ(OrthancStone::SopClassUid_Other, p->GetSopClassUid()); + ASSERT_EQ(1u, p->GetNumberOfFrames()); + ASSERT_EQ(0u, p->GetWidth()); + ASSERT_EQ(0u, p->GetHeight()); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsCloseToZero(p->GetSliceThickness())); + ASSERT_FLOAT_EQ(1, p->GetPixelSpacingX()); + ASSERT_FLOAT_EQ(1, p->GetPixelSpacingY()); + ASSERT_FALSE(p->GetGeometry().IsValid()); + ASSERT_THROW(p->GetImageInformation(), Orthanc::OrthancException); + ASSERT_FALSE(p->GetFrameGeometry(0).IsValid()); + ASSERT_THROW(p->IsColor(), Orthanc::OrthancException); // Accesses DicomImageInformation + ASSERT_FALSE(p->HasRescale()); + ASSERT_THROW(p->GetRescaleIntercept(), Orthanc::OrthancException); + ASSERT_THROW(p->GetRescaleSlope(), Orthanc::OrthancException); + ASSERT_EQ(0u, p->GetWindowingPresetsCount()); + ASSERT_THROW(p->GetWindowingPresetCenter(0), Orthanc::OrthancException); + ASSERT_THROW(p->GetWindowingPresetWidth(0), Orthanc::OrthancException); + + float c, w; + p->GetWindowingPresetsUnion(c, w); + ASSERT_FLOAT_EQ(128.0f, c); + ASSERT_FLOAT_EQ(256.0f, w); + + ASSERT_THROW(p->GetExpectedPixelFormat(), Orthanc::OrthancException); + ASSERT_FALSE(p->HasIndexInSeries()); + ASSERT_THROW(p->GetIndexInSeries(), Orthanc::OrthancException); + ASSERT_TRUE(p->GetDoseUnits().empty()); + ASSERT_DOUBLE_EQ(1.0, p->GetDoseGridScaling()); + ASSERT_DOUBLE_EQ(1.0, p->ApplyRescale(1.0)); + + double s; + ASSERT_FALSE(p->ComputeRegularSpacing(s)); + ASSERT_TRUE(p->GetFrameOfReferenceUid().empty()); +} + + +TEST(DicomInstanceParameters, Windowing) +{ + Orthanc::DicomMap m; + SetupUids(m); + m.SetValue(Orthanc::DICOM_TAG_WINDOW_CENTER, "10\\100\\1000", false); + m.SetValue(Orthanc::DICOM_TAG_WINDOW_WIDTH, "50\\60\\70", false); + + OrthancStone::DicomInstanceParameters p(m); + ASSERT_EQ(3u, p.GetWindowingPresetsCount()); + ASSERT_FLOAT_EQ(10, p.GetWindowingPresetCenter(0)); + ASSERT_FLOAT_EQ(100, p.GetWindowingPresetCenter(1)); + ASSERT_FLOAT_EQ(1000, p.GetWindowingPresetCenter(2)); + ASSERT_FLOAT_EQ(50, p.GetWindowingPresetWidth(0)); + ASSERT_FLOAT_EQ(60, p.GetWindowingPresetWidth(1)); + ASSERT_FLOAT_EQ(70, p.GetWindowingPresetWidth(2)); + + const float a = 10.0f - 50.0f / 2.0f; + const float b = 1000.0f + 70.0f / 2.0f; + + float c, w; + p.GetWindowingPresetsUnion(c, w); + ASSERT_FLOAT_EQ((a + b) / 2.0f, c); + ASSERT_FLOAT_EQ(b - a, w); +} diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/GenericToolboxTests.cpp --- a/UnitTestsSources/GenericToolboxTests.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/UnitTestsSources/GenericToolboxTests.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -26,6 +26,7 @@ #include #include +#include // For PRId64 #include TEST(GenericToolbox, TestLegitDoubleString) @@ -115,15 +116,27 @@ EXPECT_FALSE(LegitIntegerString("\n0.")); } -TEST(GenericToolbox, TestStringToDouble) + +/** + * The very long "TestStringToDouble" was split in 4 parts. Otherwise, + * while running in WebAssembly (at least in "Debug" CMAKE_BUILD_TYPE + * with Emscripten 2.0.0), one get error "failed to asynchronously + * prepare wasm: CompileError: WebAssembly.instantiate(): Compiling + * function [...] failed: local count too large". This is because the + * function is too long. + **/ +TEST(GenericToolbox, TestStringToDouble1) { using OrthancStone::GenericToolbox::StringToDouble; const double TOLERANCE = 0.00000000000001; double r = 0.0; - bool ok = StringToDouble(r, "0.0001"); - EXPECT_TRUE(ok); - EXPECT_NEAR(0.0001, r, TOLERANCE); + + { + bool ok = StringToDouble(r, "0.0001"); + EXPECT_TRUE(ok); + EXPECT_NEAR(0.0001, r, TOLERANCE); + } { bool ok = StringToDouble(r, "0.0001"); @@ -964,6 +977,14 @@ EXPECT_TRUE(ok); EXPECT_NEAR(-0.83191012798055900000, r, TOLERANCE); } +} + +TEST(GenericToolbox, TestStringToDouble2) +{ + using OrthancStone::GenericToolbox::StringToDouble; + + const double TOLERANCE = 0.00000000000001; + double r = 0.0; { bool ok = StringToDouble(r, "2.58090819604341000000"); @@ -1816,6 +1837,14 @@ EXPECT_TRUE(ok); EXPECT_NEAR(0.89377268220461400000, r, TOLERANCE); } +} + +TEST(GenericToolbox, TestStringToDouble3) +{ + using OrthancStone::GenericToolbox::StringToDouble; + + const double TOLERANCE = 0.00000000000001; + double r = 0.0; { bool ok = StringToDouble(r, "1.45674815825147000000"); @@ -2782,6 +2811,14 @@ EXPECT_TRUE(ok); EXPECT_NEAR(-1.49620375847259000000, r, TOLERANCE); } +} + +TEST(GenericToolbox, TestStringToDouble4) +{ + using OrthancStone::GenericToolbox::StringToDouble; + + const double TOLERANCE = 0.00000000000001; + double r = 0.0; { bool ok = StringToDouble(r, "-2.04336416841016000000"); @@ -4177,6 +4214,7 @@ for (double b = DBL_EPSILON; b < DBL_MAX && i < COUNT; ++i, b *= FACTOR) { int64_t bi = static_cast(b); + char txt[1024]; #if defined(_MSC_VER) # if (_MSC_VER > 1800) @@ -4185,8 +4223,9 @@ sprintf_s(txt, "%I64d", bi); # endif #else - snprintf(txt, sizeof(txt) - 1, "%ld", bi); + snprintf(txt, sizeof(txt) - 1, "%" PRId64, bi); // https://stackoverflow.com/a/9225648/881731 #endif + int64_t r = 0; bool ok = StringToInteger(r, txt); EXPECT_TRUE(ok); @@ -4194,11 +4233,11 @@ #if 0 if (ok) { - printf("OK for b = %.17f bi = %lld txt = \"%s\" and r = %lld\n", b, bi, txt, r); + printf("OK for b = %.17f bi = %" PRId64 " txt = \"%s\" and r = %" PRId64 "\n", b, bi, txt, r); } else { - printf("NOK for b = %.17f bi = %lld txt = \"%s\" and r = %lld\n", b, bi, txt,r); + printf("NOK for b = %.17f bi = %" PRId64 " txt = \"%s\" and r = %" PRId64 "\n", b, bi, txt, r); ok = StringToInteger(r, txt); } #endif diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/GeometryToolboxTests.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/GeometryToolboxTests.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,933 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2020 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + **/ + + +#include + +#include "../OrthancStone/Sources/Toolbox/FiniteProjectiveCamera.h" +#include "../OrthancStone/Sources/Toolbox/GenericToolbox.h" +#include "../OrthancStone/Sources/Toolbox/GeometryToolbox.h" + +#include +#include + +#include + + +TEST(GeometryToolbox, Interpolation) +{ + using namespace OrthancStone::GeometryToolbox; + + // https://en.wikipedia.org/wiki/Bilinear_interpolation#Application_in_image_processing + ASSERT_DOUBLE_EQ(146.1f, ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95)); + + ASSERT_DOUBLE_EQ(91, ComputeBilinearInterpolationUnitSquare(0, 0, 91, 210, 162, 95)); + ASSERT_DOUBLE_EQ(210, ComputeBilinearInterpolationUnitSquare(1, 0, 91, 210, 162, 95)); + ASSERT_DOUBLE_EQ(162, ComputeBilinearInterpolationUnitSquare(0, 1, 91, 210, 162, 95)); + ASSERT_DOUBLE_EQ(95, ComputeBilinearInterpolationUnitSquare(1, 1, 91, 210, 162, 95)); + + ASSERT_DOUBLE_EQ(123.35f, ComputeTrilinearInterpolationUnitSquare + (0.5f, 0.2f, 0.7f, + 91, 210, 162, 95, + 51, 190, 80, 92)); + + ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95), + ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 0, + 91, 210, 162, 95, + 51, 190, 80, 92)); + + ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 51, 190, 80, 92), + ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 1, + 91, 210, 162, 95, + 51, 190, 80, 92)); +} + + +static bool CompareMatrix(const OrthancStone::Matrix& a, + const OrthancStone::Matrix& b, + double threshold = 0.00000001) +{ + if (a.size1() != b.size1() || + a.size2() != b.size2()) + { + return false; + } + + for (size_t i = 0; i < a.size1(); i++) + { + for (size_t j = 0; j < a.size2(); j++) + { + if (fabs(a(i, j) - b(i, j)) > threshold) + { + LOG(ERROR) << "Too large difference in component (" + << i << "," << j << "): " << a(i,j) << " != " << b(i,j); + return false; + } + } + } + + return true; +} + + +static bool CompareVector(const OrthancStone::Vector& a, + const OrthancStone::Vector& b, + double threshold = 0.00000001) +{ + if (a.size() != b.size()) + { + return false; + } + + for (size_t i = 0; i < a.size(); i++) + { + if (fabs(a(i) - b(i)) > threshold) + { + LOG(ERROR) << "Too large difference in component " + << i << ": " << a(i) << " != " << b(i); + return false; + } + } + + return true; +} + + + +TEST(FiniteProjectiveCamera, Decomposition1) +{ + // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd + // edition" (page 163) + const double p[12] = { + 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6, + -1.03528e+2, 2.33212e+1, 4.59607e+2, -6.32525e+5, + 7.07107e-1, -3.53553e-1, 6.12372e-1, -9.18559e+2 + }; + + OrthancStone::FiniteProjectiveCamera camera(p); + ASSERT_EQ(3u, camera.GetMatrix().size1()); + ASSERT_EQ(4u, camera.GetMatrix().size2()); + ASSERT_EQ(3u, camera.GetIntrinsicParameters().size1()); + ASSERT_EQ(3u, camera.GetIntrinsicParameters().size2()); + ASSERT_EQ(3u, camera.GetRotation().size1()); + ASSERT_EQ(3u, camera.GetRotation().size2()); + ASSERT_EQ(3u, camera.GetCenter().size()); + + ASSERT_NEAR(1000.0, camera.GetCenter()[0], 0.01); + ASSERT_NEAR(2000.0, camera.GetCenter()[1], 0.01); + ASSERT_NEAR(1500.0, camera.GetCenter()[2], 0.01); + + ASSERT_NEAR(468.2, camera.GetIntrinsicParameters() (0, 0), 0.1); + ASSERT_NEAR(91.2, camera.GetIntrinsicParameters() (0, 1), 0.1); + ASSERT_NEAR(300.0, camera.GetIntrinsicParameters() (0, 2), 0.1); + ASSERT_NEAR(427.2, camera.GetIntrinsicParameters() (1, 1), 0.1); + ASSERT_NEAR(200.0, camera.GetIntrinsicParameters() (1, 2), 0.1); + ASSERT_NEAR(1.0, camera.GetIntrinsicParameters() (2, 2), 0.1); + + ASSERT_NEAR(0, camera.GetIntrinsicParameters() (1, 0), 0.0000001); + ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 0), 0.0000001); + ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 1), 0.0000001); + + ASSERT_NEAR(0.41380, camera.GetRotation() (0, 0), 0.00001); + ASSERT_NEAR(0.90915, camera.GetRotation() (0, 1), 0.00001); + ASSERT_NEAR(0.04708, camera.GetRotation() (0, 2), 0.00001); + ASSERT_NEAR(-0.57338, camera.GetRotation() (1, 0), 0.00001); + ASSERT_NEAR(0.22011, camera.GetRotation() (1, 1), 0.00001); + ASSERT_NEAR(0.78917, camera.GetRotation() (1, 2), 0.00001); + ASSERT_NEAR(0.70711, camera.GetRotation() (2, 0), 0.00001); + ASSERT_NEAR(-0.35355, camera.GetRotation() (2, 1), 0.00001); + ASSERT_NEAR(0.61237, camera.GetRotation() (2, 2), 0.00001); + + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); + + OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), + camera.GetRotation(), + camera.GetCenter()); + + ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); + ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); + ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); + ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); +} + + +TEST(FiniteProjectiveCamera, Decomposition2) +{ + const double p[] = { 1188.111986, 580.205341, -808.445330, 128000.000000, -366.466264, 1446.510501, 418.499736, 128000.000000, -0.487118, 0.291726, -0.823172, 500.000000 }; + const double k[] = { -1528.494743, 0.000000, 256.000000, 0.000000, 1528.494743, 256.000000, 0.000000, 0.000000, 1.000000 }; + const double r[] = { -0.858893, -0.330733, 0.391047, -0.158171, 0.897503, 0.411668, -0.487118, 0.291726, -0.823172 }; + const double c[] = { 243.558936, -145.863085, 411.585964 }; + + OrthancStone::FiniteProjectiveCamera camera(p); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); + + OrthancStone::FiniteProjectiveCamera camera2(k, r, c); + ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix(), 1)); + ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters(), 0.001)); + ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation(), 0.000001)); + ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter(), 0.0001)); +} + + +TEST(FiniteProjectiveCamera, Decomposition3) +{ + const double p[] = { 10, 0, 0, 0, + 0, 20, 0, 0, + 0, 0, 30, 0 }; + + OrthancStone::FiniteProjectiveCamera camera(p); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); + ASSERT_DOUBLE_EQ(10, camera.GetIntrinsicParameters() (0, 0)); + ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); + ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); +} + + +TEST(FiniteProjectiveCamera, Decomposition4) +{ + const double p[] = { 1, 0, 0, 10, + 0, 1, 0, 20, + 0, 0, 1, 30 }; + + OrthancStone::FiniteProjectiveCamera camera(p); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); + ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (0, 0)); + ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (1, 1)); + ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (2, 2)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); + ASSERT_DOUBLE_EQ(-10, camera.GetCenter() (0)); + ASSERT_DOUBLE_EQ(-20, camera.GetCenter() (1)); + ASSERT_DOUBLE_EQ(-30, camera.GetCenter() (2)); +} + + +TEST(FiniteProjectiveCamera, Decomposition5) +{ + const double p[] = { 0, 0, 10, 0, + 0, 20, 0, 0, + 30, 0, 0, 0 }; + + OrthancStone::FiniteProjectiveCamera camera(p); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); + ASSERT_DOUBLE_EQ(-10, camera.GetIntrinsicParameters() (0, 0)); + ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); + ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); + ASSERT_DOUBLE_EQ(-1, camera.GetRotation() (0, 2)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); + ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 0)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); + ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); + + OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), + camera.GetRotation(), + camera.GetCenter()); + ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); + ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); + ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); + ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); +} + + +static double GetCosAngle(const OrthancStone::Vector& a, + const OrthancStone::Vector& b) +{ + // Returns the cosine of the angle between two vectors + // https://en.wikipedia.org/wiki/Dot_product#Geometric_definition + return boost::numeric::ublas::inner_prod(a, b) / + (boost::numeric::ublas::norm_2(a) * boost::numeric::ublas::norm_2(b)); +} + + +TEST(FiniteProjectiveCamera, Ray) +{ + const double pp[] = { -1499.650894, 2954.618773, -259.737419, 637891.819097, + -2951.517707, -1501.019129, -285.785281, 637891.819097, + 0.008528, 0.003067, -0.999959, 2491.764918 }; + + const OrthancStone::FiniteProjectiveCamera camera(pp); + + ASSERT_NEAR(-21.2492, camera.GetCenter() (0), 0.0001); + ASSERT_NEAR(-7.64234, camera.GetCenter() (1), 0.00001); + ASSERT_NEAR(2491.66, camera.GetCenter() (2), 0.01); + + // Image plane that led to these parameters, with principal point at + // (256,256). The image has dimensions 512x512. + OrthancStone::Vector o = + OrthancStone::LinearAlgebra::CreateVector(7.009620, 2.521030, -821.942000); + OrthancStone::Vector ax = + OrthancStone::LinearAlgebra::CreateVector(-0.453219, 0.891399, -0.001131); + OrthancStone::Vector ay = + OrthancStone::LinearAlgebra::CreateVector(-0.891359, -0.453210, -0.008992); + + OrthancStone::CoordinateSystem3D imagePlane(o, ax, ay); + + // Back-projection of the principal point + { + OrthancStone::Vector ray = camera.GetRayDirection(256, 256); + + // The principal axis vector is orthogonal to the image plane + // (i.e. parallel to the plane normal), in the opposite direction + // ("-1" corresponds to "cos(pi)"). + ASSERT_NEAR(-1, GetCosAngle(ray, imagePlane.GetNormal()), 0.0000001); + + // Forward projection of principal axis, resulting in the principal point + double x, y; + camera.ApplyFinite(x, y, camera.GetCenter() - ray); + + ASSERT_NEAR(256, x, 0.00001); + ASSERT_NEAR(256, y, 0.00001); + } + + // Back-projection of the 4 corners of the image + std::vector cx, cy; + cx.push_back(0); + cy.push_back(0); + cx.push_back(512); + cy.push_back(0); + cx.push_back(512); + cy.push_back(512); + cx.push_back(0); + cy.push_back(512); + + bool first = true; + double angle; + + for (size_t i = 0; i < cx.size(); i++) + { + OrthancStone::Vector ray = camera.GetRayDirection(cx[i], cy[i]); + + // Check that the angle wrt. principal axis is the same for all + // the 4 corners + double a = GetCosAngle(ray, imagePlane.GetNormal()); + if (first) + { + first = false; + angle = a; + } + else + { + ASSERT_NEAR(angle, a, 0.000001); + } + + // Forward projection of the ray, going back to the original point + double x, y; + camera.ApplyFinite(x, y, camera.GetCenter() - ray); + + ASSERT_NEAR(cx[i], x, 0.00001); + ASSERT_NEAR(cy[i], y, 0.00001); + + // Alternative construction, by computing the intersection of the + // ray with the image plane + OrthancStone::Vector p; + ASSERT_TRUE(imagePlane.IntersectLine(p, camera.GetCenter(), -ray)); + imagePlane.ProjectPoint(x, y, p); + ASSERT_NEAR(cx[i], x + 256, 0.01); + ASSERT_NEAR(cy[i], y + 256, 0.01); + } +} + + +TEST(Matrix, Inverse1) +{ + OrthancStone::Matrix a, b; + + a.resize(0, 0); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(0u, b.size1()); + ASSERT_EQ(0u, b.size2()); + + a.resize(2, 3); + ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); + + a.resize(1, 1); + a(0, 0) = 45.0; + + ASSERT_DOUBLE_EQ(45, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(1u, b.size1()); + ASSERT_EQ(1u, b.size2()); + ASSERT_DOUBLE_EQ(1.0 / 45.0, b(0, 0)); + + a(0, 0) = 0; + ASSERT_DOUBLE_EQ(0, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); +} + + +TEST(Matrix, Inverse2) +{ + OrthancStone::Matrix a, b; + a.resize(2, 2); + a(0, 0) = 4; + a(0, 1) = 3; + a(1, 0) = 3; + a(1, 1) = 2; + + ASSERT_DOUBLE_EQ(-1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(2u, b.size1()); + ASSERT_EQ(2u, b.size2()); + + ASSERT_DOUBLE_EQ(-2, b(0, 0)); + ASSERT_DOUBLE_EQ(3, b(0, 1)); + ASSERT_DOUBLE_EQ(3, b(1, 0)); + ASSERT_DOUBLE_EQ(-4, b(1, 1)); + + a(0, 0) = 1; + a(0, 1) = 2; + a(1, 0) = 3; + a(1, 1) = 4; + + ASSERT_DOUBLE_EQ(-2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + + ASSERT_DOUBLE_EQ(-2, b(0, 0)); + ASSERT_DOUBLE_EQ(1, b(0, 1)); + ASSERT_DOUBLE_EQ(1.5, b(1, 0)); + ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); +} + + +TEST(Matrix, Inverse3) +{ + OrthancStone::Matrix a, b; + a.resize(3, 3); + a(0, 0) = 7; + a(0, 1) = 2; + a(0, 2) = 1; + a(1, 0) = 0; + a(1, 1) = 3; + a(1, 2) = -1; + a(2, 0) = -3; + a(2, 1) = 4; + a(2, 2) = -2; + + ASSERT_DOUBLE_EQ(1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(3u, b.size1()); + ASSERT_EQ(3u, b.size2()); + + ASSERT_DOUBLE_EQ(-2, b(0, 0)); + ASSERT_DOUBLE_EQ(8, b(0, 1)); + ASSERT_DOUBLE_EQ(-5, b(0, 2)); + ASSERT_DOUBLE_EQ(3, b(1, 0)); + ASSERT_DOUBLE_EQ(-11, b(1, 1)); + ASSERT_DOUBLE_EQ(7, b(1, 2)); + ASSERT_DOUBLE_EQ(9, b(2, 0)); + ASSERT_DOUBLE_EQ(-34, b(2, 1)); + ASSERT_DOUBLE_EQ(21, b(2, 2)); + + + a(0, 0) = 1; + a(0, 1) = 2; + a(0, 2) = 2; + a(1, 0) = 1; + a(1, 1) = 0; + a(1, 2) = 1; + a(2, 0) = 1; + a(2, 1) = 2; + a(2, 2) = 1; + + ASSERT_DOUBLE_EQ(2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(3u, b.size1()); + ASSERT_EQ(3u, b.size2()); + + ASSERT_DOUBLE_EQ(-1, b(0, 0)); + ASSERT_DOUBLE_EQ(1, b(0, 1)); + ASSERT_DOUBLE_EQ(1, b(0, 2)); + ASSERT_DOUBLE_EQ(0, b(1, 0)); + ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); + ASSERT_DOUBLE_EQ(0.5, b(1, 2)); + ASSERT_DOUBLE_EQ(1, b(2, 0)); + ASSERT_DOUBLE_EQ(0, b(2, 1)); + ASSERT_DOUBLE_EQ(-1, b(2, 2)); +} + + +TEST(Matrix, Inverse4) +{ + OrthancStone::Matrix a, b; + a.resize(4, 4); + a(0, 0) = 2; + a(0, 1) = 1; + a(0, 2) = 2; + a(0, 3) = -3; + a(1, 0) = -2; + a(1, 1) = 2; + a(1, 2) = -1; + a(1, 3) = -1; + a(2, 0) = 2; + a(2, 1) = 2; + a(2, 2) = -3; + a(2, 3) = -1; + a(3, 0) = 3; + a(3, 1) = -2; + a(3, 2) = -3; + a(3, 3) = -1; + + OrthancStone::LinearAlgebra::InvertMatrix(b, a); + ASSERT_EQ(4u, b.size1()); + ASSERT_EQ(4u, b.size2()); + + b *= 134.0; // This is the determinant + + ASSERT_DOUBLE_EQ(8, b(0, 0)); + ASSERT_DOUBLE_EQ(-44, b(0, 1)); + ASSERT_DOUBLE_EQ(30, b(0, 2)); + ASSERT_DOUBLE_EQ(-10, b(0, 3)); + ASSERT_DOUBLE_EQ(2, b(1, 0)); + ASSERT_DOUBLE_EQ(-11, b(1, 1)); + ASSERT_DOUBLE_EQ(41, b(1, 2)); + ASSERT_DOUBLE_EQ(-36, b(1, 3)); + ASSERT_DOUBLE_EQ(16, b(2, 0)); + ASSERT_DOUBLE_EQ(-21, b(2, 1)); + ASSERT_DOUBLE_EQ(-7, b(2, 2)); + ASSERT_DOUBLE_EQ(-20, b(2, 3)); + ASSERT_DOUBLE_EQ(-28, b(3, 0)); + ASSERT_DOUBLE_EQ(-47, b(3, 1)); + ASSERT_DOUBLE_EQ(29, b(3, 2)); + ASSERT_DOUBLE_EQ(-32, b(3, 3)); +} + + +TEST(FiniteProjectiveCamera, Calibration) +{ + unsigned int volumeWidth = 512; + unsigned int volumeHeight = 512; + unsigned int volumeDepth = 110; + + OrthancStone::Vector camera = OrthancStone::LinearAlgebra::CreateVector + (-1000, -5000, -static_cast(volumeDepth) * 32); + + OrthancStone::Vector principalPoint = OrthancStone::LinearAlgebra::CreateVector + (volumeWidth/2, volumeHeight/2, volumeDepth * 2); + + OrthancStone::FiniteProjectiveCamera c(camera, principalPoint, 0, 512, 512, 1, 1); + + double swapv[9] = { 1, 0, 0, + 0, -1, 512, + 0, 0, 1 }; + OrthancStone::Matrix swap; + OrthancStone::LinearAlgebra::FillMatrix(swap, 3, 3, swapv); + + OrthancStone::Matrix p = OrthancStone::LinearAlgebra::Product(swap, c.GetMatrix()); + p /= p(2,3); + + ASSERT_NEAR( 1.04437, p(0,0), 0.00001); + ASSERT_NEAR(-0.0703111, p(0,1), 0.00000001); + ASSERT_NEAR(-0.179283, p(0,2), 0.000001); + ASSERT_NEAR( 61.7431, p(0,3), 0.0001); + ASSERT_NEAR( 0.11127, p(1,0), 0.000001); + ASSERT_NEAR(-0.595541, p(1,1), 0.000001); + ASSERT_NEAR( 0.872211, p(1,2), 0.000001); + ASSERT_NEAR( 203.748, p(1,3), 0.001); + ASSERT_NEAR( 3.08593e-05, p(2,0), 0.0000000001); + ASSERT_NEAR( 0.000129138, p(2,1), 0.000000001); + ASSERT_NEAR( 9.18901e-05, p(2,2), 0.0000000001); + ASSERT_NEAR( 1, p(2,3), 0.0000001); +} + + +static bool IsEqualRotationVector(OrthancStone::Vector a, + OrthancStone::Vector b) +{ + if (a.size() != b.size() || + a.size() != 3) + { + return false; + } + else + { + OrthancStone::LinearAlgebra::NormalizeVector(a); + OrthancStone::LinearAlgebra::NormalizeVector(b); + return OrthancStone::LinearAlgebra::IsCloseToZero(boost::numeric::ublas::norm_2(a - b)); + } +} + + +TEST(GeometryToolbox, AlignVectorsWithRotation) +{ + OrthancStone::Vector a, b; + OrthancStone::Matrix r; + + OrthancStone::LinearAlgebra::AssignVector(a, -200, 200, -846.63); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); + ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); + + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, b, a); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); + ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, b), a)); + + OrthancStone::LinearAlgebra::AssignVector(a, 1, 0, 0); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); + ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); + + OrthancStone::LinearAlgebra::AssignVector(a, 0, 1, 0); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); + ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); + + OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 1); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); + ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); + ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); + + OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 0); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + ASSERT_THROW(OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b), Orthanc::OrthancException); + + // TODO: Deal with opposite vectors + + /* + OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, -1); + OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); + OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); + OrthancStone::LinearAlgebra::Print(r); + OrthancStone::LinearAlgebra::Print(boost::numeric::ublas::prod(r, a)); + */ +} + + +static bool IsEqualVectorL1(OrthancStone::Vector a, + OrthancStone::Vector b) +{ + if (a.size() != b.size()) + { + return false; + } + else + { + for (size_t i = 0; i < a.size(); i++) + { + if (!OrthancStone::LinearAlgebra::IsNear(a[i], b[i], 0.0001)) + { + return false; + } + } + + return true; + } +} + + +TEST(VolumeImageGeometry, Basic) +{ + using namespace OrthancStone; + + VolumeImageGeometry g; + g.SetSizeInVoxels(10, 20, 30); + g.SetVoxelDimensions(1, 2, 3); + + Vector p = g.GetCoordinates(0, 0, 0); + ASSERT_EQ(3u, p.size()); + ASSERT_DOUBLE_EQ(-1.0 / 2.0, p[0]); + ASSERT_DOUBLE_EQ(-2.0 / 2.0, p[1]); + ASSERT_DOUBLE_EQ(-3.0 / 2.0, p[2]); + + p = g.GetCoordinates(1, 1, 1); + ASSERT_DOUBLE_EQ(-1.0 / 2.0 + 10.0 * 1.0, p[0]); + ASSERT_DOUBLE_EQ(-2.0 / 2.0 + 20.0 * 2.0, p[1]); + ASSERT_DOUBLE_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]); + + VolumeProjection proj; + ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry().GetNormal())); + ASSERT_EQ(VolumeProjection_Axial, proj); + ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry().GetNormal())); + ASSERT_EQ(VolumeProjection_Coronal, proj); + ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry().GetNormal())); + ASSERT_EQ(VolumeProjection_Sagittal, proj); + + ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Axial)); + ASSERT_EQ(20u, g.GetProjectionHeight(VolumeProjection_Axial)); + ASSERT_EQ(30u, g.GetProjectionDepth(VolumeProjection_Axial)); + ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Coronal)); + ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Coronal)); + ASSERT_EQ(20u, g.GetProjectionDepth(VolumeProjection_Coronal)); + ASSERT_EQ(20u, g.GetProjectionWidth(VolumeProjection_Sagittal)); + ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Sagittal)); + ASSERT_EQ(10u, g.GetProjectionDepth(VolumeProjection_Sagittal)); + + p = g.GetVoxelDimensions(VolumeProjection_Axial); + ASSERT_EQ(3u, p.size()); + ASSERT_DOUBLE_EQ(1, p[0]); + ASSERT_DOUBLE_EQ(2, p[1]); + ASSERT_DOUBLE_EQ(3, p[2]); + p = g.GetVoxelDimensions(VolumeProjection_Coronal); + ASSERT_EQ(3u, p.size()); + ASSERT_DOUBLE_EQ(1, p[0]); + ASSERT_DOUBLE_EQ(3, p[1]); + ASSERT_DOUBLE_EQ(2, p[2]); + p = g.GetVoxelDimensions(VolumeProjection_Sagittal); + ASSERT_EQ(3u, p.size()); + ASSERT_DOUBLE_EQ(2, p[0]); + ASSERT_DOUBLE_EQ(3, p[1]); + ASSERT_DOUBLE_EQ(1, p[2]); + + // Loop over all the voxels of the volume + for (unsigned int z = 0; z < g.GetDepth(); z++) + { + const float zz = (0.5f + static_cast(z)) / static_cast(g.GetDepth()); // Z-center of the voxel + + for (unsigned int y = 0; y < g.GetHeight(); y++) + { + const float yy = (0.5f + static_cast(y)) / static_cast(g.GetHeight()); // Y-center of the voxel + + for (unsigned int x = 0; x < g.GetWidth(); x++) + { + const float xx = (0.5f + static_cast(x)) / static_cast(g.GetWidth()); // X-center of the voxel + + const float sx = 1.0f; + const float sy = 2.0f; + const float sz = 3.0f; + + Vector p = g.GetCoordinates(xx, yy, zz); + + Vector q = (g.GetAxialGeometry().MapSliceToWorldCoordinates( + static_cast(x) * sx, + static_cast(y) * sy) + + z * sz * g.GetAxialGeometry().GetNormal()); + ASSERT_TRUE(IsEqualVectorL1(p, q)); + + q = (g.GetCoronalGeometry().MapSliceToWorldCoordinates( + static_cast(x) * sx, + static_cast(g.GetDepth() - 1 - z) * sz) + + y * sy * g.GetCoronalGeometry().GetNormal()); + ASSERT_TRUE(IsEqualVectorL1(p, q)); + + /** + * WARNING: In sagittal geometry, the normal points to + * REDUCING X-axis in the 3D world. This is necessary to keep + * the right-hand coordinate system. Hence the "-". + **/ + q = (g.GetSagittalGeometry().MapSliceToWorldCoordinates( + static_cast(y) * sy, + static_cast(g.GetDepth() - 1 - z) * sz) + + x * sx * (-g.GetSagittalGeometry().GetNormal())); + ASSERT_TRUE(IsEqualVectorL1(p, q)); + } + } + } + + ASSERT_EQ(0, (int) VolumeProjection_Axial); + ASSERT_EQ(1, (int) VolumeProjection_Coronal); + ASSERT_EQ(2, (int) VolumeProjection_Sagittal); + + for (int p = 0; p < 3; p++) + { + VolumeProjection projection = (VolumeProjection) p; + const CoordinateSystem3D& s = g.GetProjectionGeometry(projection); + + ASSERT_THROW(g.GetProjectionSlice(projection, g.GetProjectionDepth(projection)), Orthanc::OrthancException); + + for (unsigned int i = 0; i < g.GetProjectionDepth(projection); i++) + { + CoordinateSystem3D plane = g.GetProjectionSlice(projection, i); + + if (projection == VolumeProjection_Sagittal) + { + ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast(i) * + (-s.GetNormal()) * g.GetVoxelDimensions(projection)[2])); + } + else + { + ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast(i) * + s.GetNormal() * g.GetVoxelDimensions(projection)[2])); + } + + ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisX(), s.GetAxisX())); + ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisY(), s.GetAxisY())); + + unsigned int slice; + VolumeProjection q; + ASSERT_TRUE(g.DetectSlice(q, slice, plane)); + ASSERT_EQ(projection, q); + ASSERT_EQ(i, slice); + } + } +} + + +TEST(LinearAlgebra, ParseVectorLocale) +{ + OrthancStone::Vector v; + + ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.2")); + ASSERT_EQ(1u, v.size()); + ASSERT_DOUBLE_EQ(1.2, v[0]); + + ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1.2e+2")); + ASSERT_EQ(1u, v.size()); + ASSERT_DOUBLE_EQ(-120.0, v[0]); + + ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1e-2\\2")); + ASSERT_EQ(2u, v.size()); + ASSERT_DOUBLE_EQ(-0.01, v[0]); + ASSERT_DOUBLE_EQ(2.0, v[1]); + + ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.3671875\\1.3671875")); + ASSERT_EQ(2u, v.size()); + ASSERT_DOUBLE_EQ(1.3671875, v[0]); + ASSERT_DOUBLE_EQ(1.3671875, v[1]); +} + +TEST(GenericToolbox, NormalizeUuid) +{ + std::string ref("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); + + { + std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space left + { + std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space right + { + std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space l & r + { + std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space left + case + { + std::string u(" 44CA5051-14ef-4d2f-8bd7-dB20bfb61fbb"); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space right + case + { + std::string u("44ca5051-14EF-4D2f-8bd7-db20bfb61fbB "); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // space l & r + case + { + std::string u(" 44cA5051-14Ef-4d2f-8bD7-db20bfb61fbb "); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_EQ(ref, u); + } + + // no + { + std::string u(" 44ca5051-14ef-4d2f-8bd7- db20bfb61fbb"); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_NE(ref, u); + } + + // no + { + std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fb"); + OrthancStone::GenericToolbox::NormalizeUuid(u); + ASSERT_NE(ref, u); + } +} + + +TEST(CoordinateSystem3D, Basic) +{ + using namespace OrthancStone; + + { + CoordinateSystem3D c; + ASSERT_FALSE(c.IsValid()); + ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); + ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); + ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); + + ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 0))); + ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(5, 0, 0))); + ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 5, 0))); + ASSERT_DOUBLE_EQ(5, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 5))); + } + + { + CoordinateSystem3D c("nope1", "nope2"); + ASSERT_FALSE(c.IsValid()); + ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); + ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); + ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); + } + + { + // https://www.vedantu.com/maths/perpendicular-distance-of-a-point-from-a-plane + CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, 4, -4, -6); + ASSERT_DOUBLE_EQ(3, c.ComputeDistance(LinearAlgebra::CreateVector(0, 3, 6))); + } + + { + // https://mathinsight.org/distance_point_plane_examples + CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, -2, 5, 8); + ASSERT_DOUBLE_EQ(39.0 / sqrt(33.0), c.ComputeDistance(LinearAlgebra::CreateVector(4, -4, 3))); + } + + { + // https://www.ck12.org/calculus/distance-between-a-point-and-a-plane/lesson/Distance-Between-a-Point-and-a-Plane-MAT-ALY/ + const Vector a = LinearAlgebra::CreateVector(3, 6, 9); + const Vector b = LinearAlgebra::CreateVector(9, 6, 3); + const Vector c = LinearAlgebra::CreateVector(6, -9, 9); + CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); + ASSERT_DOUBLE_EQ(0, d.ComputeDistance(a)); + ASSERT_DOUBLE_EQ(0, d.ComputeDistance(b)); + ASSERT_DOUBLE_EQ(0, d.ComputeDistance(c)); + } + + { + // https://tutorial.math.lamar.edu/classes/calcii/eqnsofplanes.aspx + const Vector a = LinearAlgebra::CreateVector(1, -2, 0); + const Vector b = LinearAlgebra::CreateVector(3, 1, 4); + const Vector c = LinearAlgebra::CreateVector(0, -1, 2); + CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); + double r = d.GetNormal() [0] / 2.0; + ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]); + ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]); + } +} diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/TestStructureSet.cpp --- a/UnitTestsSources/TestStructureSet.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/UnitTestsSources/TestStructureSet.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -29,7 +29,6 @@ #endif #include "../OrthancStone/Sources/Loaders/DicomStructureSetLoader.h" -#include "../OrthancStone/Sources/Loaders/GenericLoadersContext.h" #include "../OrthancStone/Sources/Loaders/OrthancSeriesVolumeProgressiveLoader.h" #include "../OrthancStone/Sources/Toolbox/DicomStructureSet2.h" #include "../OrthancStone/Sources/Toolbox/DicomStructureSetUtils.h" @@ -38,8 +37,6 @@ #include #include -#include - #include #include diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/UnitTestsMain.cpp --- a/UnitTestsSources/UnitTestsMain.cpp Tue Nov 24 07:40:19 2020 +0100 +++ b/UnitTestsSources/UnitTestsMain.cpp Tue Nov 24 12:59:10 2020 +0100 @@ -21,923 +21,15 @@ #include +#include "../OrthancStone/Sources/StoneEnumerations.h" #include "../OrthancStone/Sources/StoneInitialization.h" -#include "../OrthancStone/Sources/Toolbox/FiniteProjectiveCamera.h" -#include "../OrthancStone/Sources/Toolbox/GenericToolbox.h" -#include "../OrthancStone/Sources/Toolbox/GeometryToolbox.h" -#include "../OrthancStone/Sources/Volumes/ImageBuffer3D.h" - -#include -#include -#include -#include - -#include -#include -#include -#include - - -TEST(GeometryToolbox, Interpolation) -{ - using namespace OrthancStone::GeometryToolbox; - - // https://en.wikipedia.org/wiki/Bilinear_interpolation#Application_in_image_processing - ASSERT_DOUBLE_EQ(146.1f, ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95)); - - ASSERT_DOUBLE_EQ(91, ComputeBilinearInterpolationUnitSquare(0, 0, 91, 210, 162, 95)); - ASSERT_DOUBLE_EQ(210, ComputeBilinearInterpolationUnitSquare(1, 0, 91, 210, 162, 95)); - ASSERT_DOUBLE_EQ(162, ComputeBilinearInterpolationUnitSquare(0, 1, 91, 210, 162, 95)); - ASSERT_DOUBLE_EQ(95, ComputeBilinearInterpolationUnitSquare(1, 1, 91, 210, 162, 95)); - - ASSERT_DOUBLE_EQ(123.35f, ComputeTrilinearInterpolationUnitSquare - (0.5f, 0.2f, 0.7f, - 91, 210, 162, 95, - 51, 190, 80, 92)); - - ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 91, 210, 162, 95), - ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 0, - 91, 210, 162, 95, - 51, 190, 80, 92)); - - ASSERT_DOUBLE_EQ(ComputeBilinearInterpolationUnitSquare(0.5f, 0.2f, 51, 190, 80, 92), - ComputeTrilinearInterpolationUnitSquare(0.5f, 0.2f, 1, - 91, 210, 162, 95, - 51, 190, 80, 92)); -} - - -static bool CompareMatrix(const OrthancStone::Matrix& a, - const OrthancStone::Matrix& b, - double threshold = 0.00000001) -{ - if (a.size1() != b.size1() || - a.size2() != b.size2()) - { - return false; - } - - for (size_t i = 0; i < a.size1(); i++) - { - for (size_t j = 0; j < a.size2(); j++) - { - if (fabs(a(i, j) - b(i, j)) > threshold) - { - LOG(ERROR) << "Too large difference in component (" - << i << "," << j << "): " << a(i,j) << " != " << b(i,j); - return false; - } - } - } - - return true; -} - - -static bool CompareVector(const OrthancStone::Vector& a, - const OrthancStone::Vector& b, - double threshold = 0.00000001) -{ - if (a.size() != b.size()) - { - return false; - } - - for (size_t i = 0; i < a.size(); i++) - { - if (fabs(a(i) - b(i)) > threshold) - { - LOG(ERROR) << "Too large difference in component " - << i << ": " << a(i) << " != " << b(i); - return false; - } - } - - return true; -} - - - -TEST(FiniteProjectiveCamera, Decomposition1) -{ - // Example 6.2 of "Multiple View Geometry in Computer Vision - 2nd - // edition" (page 163) - const double p[12] = { - 3.53553e+2, 3.39645e+2, 2.77744e+2, -1.44946e+6, - -1.03528e+2, 2.33212e+1, 4.59607e+2, -6.32525e+5, - 7.07107e-1, -3.53553e-1, 6.12372e-1, -9.18559e+2 - }; - - OrthancStone::FiniteProjectiveCamera camera(p); - ASSERT_EQ(3u, camera.GetMatrix().size1()); - ASSERT_EQ(4u, camera.GetMatrix().size2()); - ASSERT_EQ(3u, camera.GetIntrinsicParameters().size1()); - ASSERT_EQ(3u, camera.GetIntrinsicParameters().size2()); - ASSERT_EQ(3u, camera.GetRotation().size1()); - ASSERT_EQ(3u, camera.GetRotation().size2()); - ASSERT_EQ(3u, camera.GetCenter().size()); - - ASSERT_NEAR(1000.0, camera.GetCenter()[0], 0.01); - ASSERT_NEAR(2000.0, camera.GetCenter()[1], 0.01); - ASSERT_NEAR(1500.0, camera.GetCenter()[2], 0.01); - - ASSERT_NEAR(468.2, camera.GetIntrinsicParameters() (0, 0), 0.1); - ASSERT_NEAR(91.2, camera.GetIntrinsicParameters() (0, 1), 0.1); - ASSERT_NEAR(300.0, camera.GetIntrinsicParameters() (0, 2), 0.1); - ASSERT_NEAR(427.2, camera.GetIntrinsicParameters() (1, 1), 0.1); - ASSERT_NEAR(200.0, camera.GetIntrinsicParameters() (1, 2), 0.1); - ASSERT_NEAR(1.0, camera.GetIntrinsicParameters() (2, 2), 0.1); - - ASSERT_NEAR(0, camera.GetIntrinsicParameters() (1, 0), 0.0000001); - ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 0), 0.0000001); - ASSERT_NEAR(0, camera.GetIntrinsicParameters() (2, 1), 0.0000001); - - ASSERT_NEAR(0.41380, camera.GetRotation() (0, 0), 0.00001); - ASSERT_NEAR(0.90915, camera.GetRotation() (0, 1), 0.00001); - ASSERT_NEAR(0.04708, camera.GetRotation() (0, 2), 0.00001); - ASSERT_NEAR(-0.57338, camera.GetRotation() (1, 0), 0.00001); - ASSERT_NEAR(0.22011, camera.GetRotation() (1, 1), 0.00001); - ASSERT_NEAR(0.78917, camera.GetRotation() (1, 2), 0.00001); - ASSERT_NEAR(0.70711, camera.GetRotation() (2, 0), 0.00001); - ASSERT_NEAR(-0.35355, camera.GetRotation() (2, 1), 0.00001); - ASSERT_NEAR(0.61237, camera.GetRotation() (2, 2), 0.00001); - - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); - - OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), - camera.GetRotation(), - camera.GetCenter()); - - ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); - ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); - ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); - ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); -} - - -TEST(FiniteProjectiveCamera, Decomposition2) -{ - const double p[] = { 1188.111986, 580.205341, -808.445330, 128000.000000, -366.466264, 1446.510501, 418.499736, 128000.000000, -0.487118, 0.291726, -0.823172, 500.000000 }; - const double k[] = { -1528.494743, 0.000000, 256.000000, 0.000000, 1528.494743, 256.000000, 0.000000, 0.000000, 1.000000 }; - const double r[] = { -0.858893, -0.330733, 0.391047, -0.158171, 0.897503, 0.411668, -0.487118, 0.291726, -0.823172 }; - const double c[] = { 243.558936, -145.863085, 411.585964 }; - - OrthancStone::FiniteProjectiveCamera camera(p); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); - - OrthancStone::FiniteProjectiveCamera camera2(k, r, c); - ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix(), 1)); - ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters(), 0.001)); - ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation(), 0.000001)); - ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter(), 0.0001)); -} - - -TEST(FiniteProjectiveCamera, Decomposition3) -{ - const double p[] = { 10, 0, 0, 0, - 0, 20, 0, 0, - 0, 0, 30, 0 }; - - OrthancStone::FiniteProjectiveCamera camera(p); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); - ASSERT_DOUBLE_EQ(10, camera.GetIntrinsicParameters() (0, 0)); - ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); - ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); -} - - -TEST(FiniteProjectiveCamera, Decomposition4) -{ - const double p[] = { 1, 0, 0, 10, - 0, 1, 0, 20, - 0, 0, 1, 30 }; - - OrthancStone::FiniteProjectiveCamera camera(p); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); - ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (0, 0)); - ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (1, 1)); - ASSERT_DOUBLE_EQ(1, camera.GetIntrinsicParameters() (2, 2)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (0, 0)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 2)); - ASSERT_DOUBLE_EQ(-10, camera.GetCenter() (0)); - ASSERT_DOUBLE_EQ(-20, camera.GetCenter() (1)); - ASSERT_DOUBLE_EQ(-30, camera.GetCenter() (2)); -} - - -TEST(FiniteProjectiveCamera, Decomposition5) -{ - const double p[] = { 0, 0, 10, 0, - 0, 20, 0, 0, - 30, 0, 0, 0 }; - OrthancStone::FiniteProjectiveCamera camera(p); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(camera.GetRotation())); - ASSERT_DOUBLE_EQ(-10, camera.GetIntrinsicParameters() (0, 0)); - ASSERT_DOUBLE_EQ(20, camera.GetIntrinsicParameters() (1, 1)); - ASSERT_DOUBLE_EQ(30, camera.GetIntrinsicParameters() (2, 2)); - ASSERT_DOUBLE_EQ(-1, camera.GetRotation() (0, 2)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (1, 1)); - ASSERT_DOUBLE_EQ(1, camera.GetRotation() (2, 0)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (0)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (1)); - ASSERT_DOUBLE_EQ(0, camera.GetCenter() (2)); - - OrthancStone::FiniteProjectiveCamera camera2(camera.GetIntrinsicParameters(), - camera.GetRotation(), - camera.GetCenter()); - ASSERT_TRUE(CompareMatrix(camera.GetMatrix(), camera2.GetMatrix())); - ASSERT_TRUE(CompareMatrix(camera.GetIntrinsicParameters(), camera2.GetIntrinsicParameters())); - ASSERT_TRUE(CompareMatrix(camera.GetRotation(), camera2.GetRotation())); - ASSERT_TRUE(CompareVector(camera.GetCenter(), camera2.GetCenter())); -} - - -static double GetCosAngle(const OrthancStone::Vector& a, - const OrthancStone::Vector& b) -{ - // Returns the cosine of the angle between two vectors - // https://en.wikipedia.org/wiki/Dot_product#Geometric_definition - return boost::numeric::ublas::inner_prod(a, b) / - (boost::numeric::ublas::norm_2(a) * boost::numeric::ublas::norm_2(b)); -} - - -TEST(FiniteProjectiveCamera, Ray) -{ - const double pp[] = { -1499.650894, 2954.618773, -259.737419, 637891.819097, - -2951.517707, -1501.019129, -285.785281, 637891.819097, - 0.008528, 0.003067, -0.999959, 2491.764918 }; - - const OrthancStone::FiniteProjectiveCamera camera(pp); - - ASSERT_NEAR(-21.2492, camera.GetCenter() (0), 0.0001); - ASSERT_NEAR(-7.64234, camera.GetCenter() (1), 0.00001); - ASSERT_NEAR(2491.66, camera.GetCenter() (2), 0.01); - - // Image plane that led to these parameters, with principal point at - // (256,256). The image has dimensions 512x512. - OrthancStone::Vector o = - OrthancStone::LinearAlgebra::CreateVector(7.009620, 2.521030, -821.942000); - OrthancStone::Vector ax = - OrthancStone::LinearAlgebra::CreateVector(-0.453219, 0.891399, -0.001131); - OrthancStone::Vector ay = - OrthancStone::LinearAlgebra::CreateVector(-0.891359, -0.453210, -0.008992); - - OrthancStone::CoordinateSystem3D imagePlane(o, ax, ay); - - // Back-projection of the principal point - { - OrthancStone::Vector ray = camera.GetRayDirection(256, 256); - - // The principal axis vector is orthogonal to the image plane - // (i.e. parallel to the plane normal), in the opposite direction - // ("-1" corresponds to "cos(pi)"). - ASSERT_NEAR(-1, GetCosAngle(ray, imagePlane.GetNormal()), 0.0000001); - - // Forward projection of principal axis, resulting in the principal point - double x, y; - camera.ApplyFinite(x, y, camera.GetCenter() - ray); - - ASSERT_NEAR(256, x, 0.00001); - ASSERT_NEAR(256, y, 0.00001); - } - - // Back-projection of the 4 corners of the image - std::vector cx, cy; - cx.push_back(0); - cy.push_back(0); - cx.push_back(512); - cy.push_back(0); - cx.push_back(512); - cy.push_back(512); - cx.push_back(0); - cy.push_back(512); - - bool first = true; - double angle; - - for (size_t i = 0; i < cx.size(); i++) - { - OrthancStone::Vector ray = camera.GetRayDirection(cx[i], cy[i]); - - // Check that the angle wrt. principal axis is the same for all - // the 4 corners - double a = GetCosAngle(ray, imagePlane.GetNormal()); - if (first) - { - first = false; - angle = a; - } - else - { - ASSERT_NEAR(angle, a, 0.000001); - } - - // Forward projection of the ray, going back to the original point - double x, y; - camera.ApplyFinite(x, y, camera.GetCenter() - ray); - - ASSERT_NEAR(cx[i], x, 0.00001); - ASSERT_NEAR(cy[i], y, 0.00001); - - // Alternative construction, by computing the intersection of the - // ray with the image plane - OrthancStone::Vector p; - ASSERT_TRUE(imagePlane.IntersectLine(p, camera.GetCenter(), -ray)); - imagePlane.ProjectPoint(x, y, p); - ASSERT_NEAR(cx[i], x + 256, 0.01); - ASSERT_NEAR(cy[i], y + 256, 0.01); - } -} - - -TEST(Matrix, Inverse1) -{ - OrthancStone::Matrix a, b; - - a.resize(0, 0); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(0u, b.size1()); - ASSERT_EQ(0u, b.size2()); - - a.resize(2, 3); - ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); - - a.resize(1, 1); - a(0, 0) = 45.0; - - ASSERT_DOUBLE_EQ(45, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(1u, b.size1()); - ASSERT_EQ(1u, b.size2()); - ASSERT_DOUBLE_EQ(1.0 / 45.0, b(0, 0)); - - a(0, 0) = 0; - ASSERT_DOUBLE_EQ(0, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - ASSERT_THROW(OrthancStone::LinearAlgebra::InvertMatrix(b, a), Orthanc::OrthancException); -} - - -TEST(Matrix, Inverse2) -{ - OrthancStone::Matrix a, b; - a.resize(2, 2); - a(0, 0) = 4; - a(0, 1) = 3; - a(1, 0) = 3; - a(1, 1) = 2; - - ASSERT_DOUBLE_EQ(-1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(2u, b.size1()); - ASSERT_EQ(2u, b.size2()); - - ASSERT_DOUBLE_EQ(-2, b(0, 0)); - ASSERT_DOUBLE_EQ(3, b(0, 1)); - ASSERT_DOUBLE_EQ(3, b(1, 0)); - ASSERT_DOUBLE_EQ(-4, b(1, 1)); - - a(0, 0) = 1; - a(0, 1) = 2; - a(1, 0) = 3; - a(1, 1) = 4; - - ASSERT_DOUBLE_EQ(-2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - - ASSERT_DOUBLE_EQ(-2, b(0, 0)); - ASSERT_DOUBLE_EQ(1, b(0, 1)); - ASSERT_DOUBLE_EQ(1.5, b(1, 0)); - ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); -} - - -TEST(Matrix, Inverse3) -{ - OrthancStone::Matrix a, b; - a.resize(3, 3); - a(0, 0) = 7; - a(0, 1) = 2; - a(0, 2) = 1; - a(1, 0) = 0; - a(1, 1) = 3; - a(1, 2) = -1; - a(2, 0) = -3; - a(2, 1) = 4; - a(2, 2) = -2; - - ASSERT_DOUBLE_EQ(1, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(3u, b.size1()); - ASSERT_EQ(3u, b.size2()); - - ASSERT_DOUBLE_EQ(-2, b(0, 0)); - ASSERT_DOUBLE_EQ(8, b(0, 1)); - ASSERT_DOUBLE_EQ(-5, b(0, 2)); - ASSERT_DOUBLE_EQ(3, b(1, 0)); - ASSERT_DOUBLE_EQ(-11, b(1, 1)); - ASSERT_DOUBLE_EQ(7, b(1, 2)); - ASSERT_DOUBLE_EQ(9, b(2, 0)); - ASSERT_DOUBLE_EQ(-34, b(2, 1)); - ASSERT_DOUBLE_EQ(21, b(2, 2)); - - - a(0, 0) = 1; - a(0, 1) = 2; - a(0, 2) = 2; - a(1, 0) = 1; - a(1, 1) = 0; - a(1, 2) = 1; - a(2, 0) = 1; - a(2, 1) = 2; - a(2, 2) = 1; - - ASSERT_DOUBLE_EQ(2, OrthancStone::LinearAlgebra::ComputeDeterminant(a)); - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(3u, b.size1()); - ASSERT_EQ(3u, b.size2()); - - ASSERT_DOUBLE_EQ(-1, b(0, 0)); - ASSERT_DOUBLE_EQ(1, b(0, 1)); - ASSERT_DOUBLE_EQ(1, b(0, 2)); - ASSERT_DOUBLE_EQ(0, b(1, 0)); - ASSERT_DOUBLE_EQ(-0.5, b(1, 1)); - ASSERT_DOUBLE_EQ(0.5, b(1, 2)); - ASSERT_DOUBLE_EQ(1, b(2, 0)); - ASSERT_DOUBLE_EQ(0, b(2, 1)); - ASSERT_DOUBLE_EQ(-1, b(2, 2)); -} - +#include -TEST(Matrix, Inverse4) -{ - OrthancStone::Matrix a, b; - a.resize(4, 4); - a(0, 0) = 2; - a(0, 1) = 1; - a(0, 2) = 2; - a(0, 3) = -3; - a(1, 0) = -2; - a(1, 1) = 2; - a(1, 2) = -1; - a(1, 3) = -1; - a(2, 0) = 2; - a(2, 1) = 2; - a(2, 2) = -3; - a(2, 3) = -1; - a(3, 0) = 3; - a(3, 1) = -2; - a(3, 2) = -3; - a(3, 3) = -1; - - OrthancStone::LinearAlgebra::InvertMatrix(b, a); - ASSERT_EQ(4u, b.size1()); - ASSERT_EQ(4u, b.size2()); - - b *= 134.0; // This is the determinant - - ASSERT_DOUBLE_EQ(8, b(0, 0)); - ASSERT_DOUBLE_EQ(-44, b(0, 1)); - ASSERT_DOUBLE_EQ(30, b(0, 2)); - ASSERT_DOUBLE_EQ(-10, b(0, 3)); - ASSERT_DOUBLE_EQ(2, b(1, 0)); - ASSERT_DOUBLE_EQ(-11, b(1, 1)); - ASSERT_DOUBLE_EQ(41, b(1, 2)); - ASSERT_DOUBLE_EQ(-36, b(1, 3)); - ASSERT_DOUBLE_EQ(16, b(2, 0)); - ASSERT_DOUBLE_EQ(-21, b(2, 1)); - ASSERT_DOUBLE_EQ(-7, b(2, 2)); - ASSERT_DOUBLE_EQ(-20, b(2, 3)); - ASSERT_DOUBLE_EQ(-28, b(3, 0)); - ASSERT_DOUBLE_EQ(-47, b(3, 1)); - ASSERT_DOUBLE_EQ(29, b(3, 2)); - ASSERT_DOUBLE_EQ(-32, b(3, 3)); -} - - -TEST(FiniteProjectiveCamera, Calibration) -{ - unsigned int volumeWidth = 512; - unsigned int volumeHeight = 512; - unsigned int volumeDepth = 110; - - OrthancStone::Vector camera = OrthancStone::LinearAlgebra::CreateVector - (-1000, -5000, -static_cast(volumeDepth) * 32); - - OrthancStone::Vector principalPoint = OrthancStone::LinearAlgebra::CreateVector - (volumeWidth/2, volumeHeight/2, volumeDepth * 2); - - OrthancStone::FiniteProjectiveCamera c(camera, principalPoint, 0, 512, 512, 1, 1); - - double swapv[9] = { 1, 0, 0, - 0, -1, 512, - 0, 0, 1 }; - OrthancStone::Matrix swap; - OrthancStone::LinearAlgebra::FillMatrix(swap, 3, 3, swapv); - - OrthancStone::Matrix p = OrthancStone::LinearAlgebra::Product(swap, c.GetMatrix()); - p /= p(2,3); - - ASSERT_NEAR( 1.04437, p(0,0), 0.00001); - ASSERT_NEAR(-0.0703111, p(0,1), 0.00000001); - ASSERT_NEAR(-0.179283, p(0,2), 0.000001); - ASSERT_NEAR( 61.7431, p(0,3), 0.0001); - ASSERT_NEAR( 0.11127, p(1,0), 0.000001); - ASSERT_NEAR(-0.595541, p(1,1), 0.000001); - ASSERT_NEAR( 0.872211, p(1,2), 0.000001); - ASSERT_NEAR( 203.748, p(1,3), 0.001); - ASSERT_NEAR( 3.08593e-05, p(2,0), 0.0000000001); - ASSERT_NEAR( 0.000129138, p(2,1), 0.000000001); - ASSERT_NEAR( 9.18901e-05, p(2,2), 0.0000000001); - ASSERT_NEAR( 1, p(2,3), 0.0000001); -} - - -static bool IsEqualRotationVector(OrthancStone::Vector a, - OrthancStone::Vector b) -{ - if (a.size() != b.size() || - a.size() != 3) - { - return false; - } - else - { - OrthancStone::LinearAlgebra::NormalizeVector(a); - OrthancStone::LinearAlgebra::NormalizeVector(b); - return OrthancStone::LinearAlgebra::IsCloseToZero(boost::numeric::ublas::norm_2(a - b)); - } -} - - -TEST(GeometryToolbox, AlignVectorsWithRotation) -{ - OrthancStone::Vector a, b; - OrthancStone::Matrix r; - - OrthancStone::LinearAlgebra::AssignVector(a, -200, 200, -846.63); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); - ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); - - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, b, a); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); - ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, b), a)); - - OrthancStone::LinearAlgebra::AssignVector(a, 1, 0, 0); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); - ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); - - OrthancStone::LinearAlgebra::AssignVector(a, 0, 1, 0); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); - ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); - - OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 1); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); - ASSERT_TRUE(OrthancStone::LinearAlgebra::IsRotationMatrix(r)); - ASSERT_TRUE(IsEqualRotationVector(OrthancStone::LinearAlgebra::Product(r, a), b)); - - OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, 0); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - ASSERT_THROW(OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b), Orthanc::OrthancException); - - // TODO: Deal with opposite vectors - - /* - OrthancStone::LinearAlgebra::AssignVector(a, 0, 0, -1); - OrthancStone::LinearAlgebra::AssignVector(b, 0, 0, 1); - OrthancStone::GeometryToolbox::AlignVectorsWithRotation(r, a, b); - OrthancStone::LinearAlgebra::Print(r); - OrthancStone::LinearAlgebra::Print(boost::numeric::ublas::prod(r, a)); - */ -} - - -static bool IsEqualVectorL1(OrthancStone::Vector a, - OrthancStone::Vector b) -{ - if (a.size() != b.size()) - { - return false; - } - else - { - for (size_t i = 0; i < a.size(); i++) - { - if (!OrthancStone::LinearAlgebra::IsNear(a[i], b[i], 0.0001)) - { - return false; - } - } - - return true; - } -} - - -TEST(VolumeImageGeometry, Basic) -{ - using namespace OrthancStone; - - VolumeImageGeometry g; - g.SetSizeInVoxels(10, 20, 30); - g.SetVoxelDimensions(1, 2, 3); - - Vector p = g.GetCoordinates(0, 0, 0); - ASSERT_EQ(3u, p.size()); - ASSERT_DOUBLE_EQ(-1.0 / 2.0, p[0]); - ASSERT_DOUBLE_EQ(-2.0 / 2.0, p[1]); - ASSERT_DOUBLE_EQ(-3.0 / 2.0, p[2]); - - p = g.GetCoordinates(1, 1, 1); - ASSERT_DOUBLE_EQ(-1.0 / 2.0 + 10.0 * 1.0, p[0]); - ASSERT_DOUBLE_EQ(-2.0 / 2.0 + 20.0 * 2.0, p[1]); - ASSERT_DOUBLE_EQ(-3.0 / 2.0 + 30.0 * 3.0, p[2]); - - VolumeProjection proj; - ASSERT_TRUE(g.DetectProjection(proj, g.GetAxialGeometry().GetNormal())); - ASSERT_EQ(VolumeProjection_Axial, proj); - ASSERT_TRUE(g.DetectProjection(proj, g.GetCoronalGeometry().GetNormal())); - ASSERT_EQ(VolumeProjection_Coronal, proj); - ASSERT_TRUE(g.DetectProjection(proj, g.GetSagittalGeometry().GetNormal())); - ASSERT_EQ(VolumeProjection_Sagittal, proj); - - ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Axial)); - ASSERT_EQ(20u, g.GetProjectionHeight(VolumeProjection_Axial)); - ASSERT_EQ(30u, g.GetProjectionDepth(VolumeProjection_Axial)); - ASSERT_EQ(10u, g.GetProjectionWidth(VolumeProjection_Coronal)); - ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Coronal)); - ASSERT_EQ(20u, g.GetProjectionDepth(VolumeProjection_Coronal)); - ASSERT_EQ(20u, g.GetProjectionWidth(VolumeProjection_Sagittal)); - ASSERT_EQ(30u, g.GetProjectionHeight(VolumeProjection_Sagittal)); - ASSERT_EQ(10u, g.GetProjectionDepth(VolumeProjection_Sagittal)); - - p = g.GetVoxelDimensions(VolumeProjection_Axial); - ASSERT_EQ(3u, p.size()); - ASSERT_DOUBLE_EQ(1, p[0]); - ASSERT_DOUBLE_EQ(2, p[1]); - ASSERT_DOUBLE_EQ(3, p[2]); - p = g.GetVoxelDimensions(VolumeProjection_Coronal); - ASSERT_EQ(3u, p.size()); - ASSERT_DOUBLE_EQ(1, p[0]); - ASSERT_DOUBLE_EQ(3, p[1]); - ASSERT_DOUBLE_EQ(2, p[2]); - p = g.GetVoxelDimensions(VolumeProjection_Sagittal); - ASSERT_EQ(3u, p.size()); - ASSERT_DOUBLE_EQ(2, p[0]); - ASSERT_DOUBLE_EQ(3, p[1]); - ASSERT_DOUBLE_EQ(1, p[2]); +#if defined(__EMSCRIPTEN__) +# include +#endif - // Loop over all the voxels of the volume - for (unsigned int z = 0; z < g.GetDepth(); z++) - { - const float zz = (0.5f + static_cast(z)) / static_cast(g.GetDepth()); // Z-center of the voxel - - for (unsigned int y = 0; y < g.GetHeight(); y++) - { - const float yy = (0.5f + static_cast(y)) / static_cast(g.GetHeight()); // Y-center of the voxel - - for (unsigned int x = 0; x < g.GetWidth(); x++) - { - const float xx = (0.5f + static_cast(x)) / static_cast(g.GetWidth()); // X-center of the voxel - - const float sx = 1.0f; - const float sy = 2.0f; - const float sz = 3.0f; - - Vector p = g.GetCoordinates(xx, yy, zz); - - Vector q = (g.GetAxialGeometry().MapSliceToWorldCoordinates( - static_cast(x) * sx, - static_cast(y) * sy) + - z * sz * g.GetAxialGeometry().GetNormal()); - ASSERT_TRUE(IsEqualVectorL1(p, q)); - - q = (g.GetCoronalGeometry().MapSliceToWorldCoordinates( - static_cast(x) * sx, - static_cast(g.GetDepth() - 1 - z) * sz) + - y * sy * g.GetCoronalGeometry().GetNormal()); - ASSERT_TRUE(IsEqualVectorL1(p, q)); - - /** - * WARNING: In sagittal geometry, the normal points to - * REDUCING X-axis in the 3D world. This is necessary to keep - * the right-hand coordinate system. Hence the "-". - **/ - q = (g.GetSagittalGeometry().MapSliceToWorldCoordinates( - static_cast(y) * sy, - static_cast(g.GetDepth() - 1 - z) * sz) + - x * sx * (-g.GetSagittalGeometry().GetNormal())); - ASSERT_TRUE(IsEqualVectorL1(p, q)); - } - } - } - - ASSERT_EQ(0, (int) VolumeProjection_Axial); - ASSERT_EQ(1, (int) VolumeProjection_Coronal); - ASSERT_EQ(2, (int) VolumeProjection_Sagittal); - - for (int p = 0; p < 3; p++) - { - VolumeProjection projection = (VolumeProjection) p; - const CoordinateSystem3D& s = g.GetProjectionGeometry(projection); - - ASSERT_THROW(g.GetProjectionSlice(projection, g.GetProjectionDepth(projection)), Orthanc::OrthancException); - - for (unsigned int i = 0; i < g.GetProjectionDepth(projection); i++) - { - CoordinateSystem3D plane = g.GetProjectionSlice(projection, i); - - if (projection == VolumeProjection_Sagittal) - { - ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast(i) * - (-s.GetNormal()) * g.GetVoxelDimensions(projection)[2])); - } - else - { - ASSERT_TRUE(IsEqualVectorL1(plane.GetOrigin(), s.GetOrigin() + static_cast(i) * - s.GetNormal() * g.GetVoxelDimensions(projection)[2])); - } - - ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisX(), s.GetAxisX())); - ASSERT_TRUE(IsEqualVectorL1(plane.GetAxisY(), s.GetAxisY())); - - unsigned int slice; - VolumeProjection q; - ASSERT_TRUE(g.DetectSlice(q, slice, plane)); - ASSERT_EQ(projection, q); - ASSERT_EQ(i, slice); - } - } -} - - -TEST(LinearAlgebra, ParseVectorLocale) -{ - OrthancStone::Vector v; - - ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.2")); - ASSERT_EQ(1u, v.size()); - ASSERT_DOUBLE_EQ(1.2, v[0]); - - ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1.2e+2")); - ASSERT_EQ(1u, v.size()); - ASSERT_DOUBLE_EQ(-120.0, v[0]); - - ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "-1e-2\\2")); - ASSERT_EQ(2u, v.size()); - ASSERT_DOUBLE_EQ(-0.01, v[0]); - ASSERT_DOUBLE_EQ(2.0, v[1]); - - ASSERT_TRUE(OrthancStone::LinearAlgebra::ParseVector(v, "1.3671875\\1.3671875")); - ASSERT_EQ(2u, v.size()); - ASSERT_DOUBLE_EQ(1.3671875, v[0]); - ASSERT_DOUBLE_EQ(1.3671875, v[1]); -} - -TEST(GenericToolbox, NormalizeUuid) -{ - std::string ref("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); - - { - std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space left - { - std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb"); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space right - { - std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space l & r - { - std::string u(" 44ca5051-14ef-4d2f-8bd7-db20bfb61fbb "); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space left + case - { - std::string u(" 44CA5051-14ef-4d2f-8bd7-dB20bfb61fbb"); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space right + case - { - std::string u("44ca5051-14EF-4D2f-8bd7-db20bfb61fbB "); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // space l & r + case - { - std::string u(" 44cA5051-14Ef-4d2f-8bD7-db20bfb61fbb "); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_EQ(ref, u); - } - - // no - { - std::string u(" 44ca5051-14ef-4d2f-8bd7- db20bfb61fbb"); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_NE(ref, u); - } - - // no - { - std::string u("44ca5051-14ef-4d2f-8bd7-db20bfb61fb"); - OrthancStone::GenericToolbox::NormalizeUuid(u); - ASSERT_NE(ref, u); - } -} - - -TEST(CoordinateSystem3D, Basic) -{ - using namespace OrthancStone; - - { - CoordinateSystem3D c; - ASSERT_FALSE(c.IsValid()); - ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); - ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); - ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); - - ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 0))); - ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(5, 0, 0))); - ASSERT_DOUBLE_EQ(0, c.ComputeDistance(LinearAlgebra::CreateVector(0, 5, 0))); - ASSERT_DOUBLE_EQ(5, c.ComputeDistance(LinearAlgebra::CreateVector(0, 0, 5))); - } - - { - CoordinateSystem3D c("nope1", "nope2"); - ASSERT_FALSE(c.IsValid()); - ASSERT_DOUBLE_EQ(c.GetNormal()[0], 0); - ASSERT_DOUBLE_EQ(c.GetNormal()[1], 0); - ASSERT_DOUBLE_EQ(c.GetNormal()[2], 1); - } - - { - // https://www.vedantu.com/maths/perpendicular-distance-of-a-point-from-a-plane - CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, 4, -4, -6); - ASSERT_DOUBLE_EQ(3, c.ComputeDistance(LinearAlgebra::CreateVector(0, 3, 6))); - } - - { - // https://mathinsight.org/distance_point_plane_examples - CoordinateSystem3D c = CoordinateSystem3D::CreateFromPlaneGeneralForm(2, -2, 5, 8); - ASSERT_DOUBLE_EQ(39.0 / sqrt(33.0), c.ComputeDistance(LinearAlgebra::CreateVector(4, -4, 3))); - } - - { - // https://www.ck12.org/calculus/distance-between-a-point-and-a-plane/lesson/Distance-Between-a-Point-and-a-Plane-MAT-ALY/ - const Vector a = LinearAlgebra::CreateVector(3, 6, 9); - const Vector b = LinearAlgebra::CreateVector(9, 6, 3); - const Vector c = LinearAlgebra::CreateVector(6, -9, 9); - CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); - ASSERT_DOUBLE_EQ(0, d.ComputeDistance(a)); - ASSERT_DOUBLE_EQ(0, d.ComputeDistance(b)); - ASSERT_DOUBLE_EQ(0, d.ComputeDistance(c)); - } - - { - // https://tutorial.math.lamar.edu/classes/calcii/eqnsofplanes.aspx - const Vector a = LinearAlgebra::CreateVector(1, -2, 0); - const Vector b = LinearAlgebra::CreateVector(3, 1, 4); - const Vector c = LinearAlgebra::CreateVector(0, -1, 2); - CoordinateSystem3D d = CoordinateSystem3D::CreateFromThreePoints(a, b, c); - double r = d.GetNormal() [0] / 2.0; - ASSERT_DOUBLE_EQ(-8 * r, d.GetNormal() [1]); - ASSERT_DOUBLE_EQ(5 * r, d.GetNormal() [2]); - } -} TEST(Enumerations, Basic) @@ -965,13 +57,37 @@ int main(int argc, char **argv) { - Orthanc::Logging::Initialize(); - Orthanc::Logging::EnableInfoLevel(true); +#if defined(__EMSCRIPTEN__) + std::string output; +#endif + + int result; + + { + Orthanc::Logging::Initialize(); + Orthanc::Logging::EnableInfoLevel(true); - ::testing::InitGoogleTest(&argc, argv); - int result = RUN_ALL_TESTS(); + ::testing::InitGoogleTest(&argc, argv); + +#if defined(__EMSCRIPTEN__) + ::testing::internal::CaptureStdout(); +#endif + + result = RUN_ALL_TESTS(); + + Orthanc::Logging::Finalize(); - Orthanc::Logging::Finalize(); +#if defined(__EMSCRIPTEN__) + output = testing::internal::GetCapturedStdout(); +#endif + } + +#if defined(__EMSCRIPTEN__) + EM_ASM({ + document.getElementById("output").innerHTML = UTF8ToString($0); + }, + output.c_str()); +#endif return result; } diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/UnitTestsSources.cmake --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/UnitTestsSources.cmake Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,36 @@ +# Stone of Orthanc +# Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics +# Department, University Hospital of Liege, Belgium +# Copyright (C) 2017-2020 Osimis S.A., Belgium +# +# This program is free software: you can redistribute it and/or +# modify it under the terms of the GNU Affero General Public License +# as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + + +set(UNIT_TESTS_SOURCES + ${CMAKE_CURRENT_LIST_DIR}/DicomTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/GenericToolboxTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/GeometryToolboxTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/ImageToolboxTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/PixelTestPatternsTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/SortedFramesTests.cpp + ${CMAKE_CURRENT_LIST_DIR}/TestMessageBroker.cpp + ${CMAKE_CURRENT_LIST_DIR}/TestStrategy.cpp + ${CMAKE_CURRENT_LIST_DIR}/TestStructureSet.cpp + ${CMAKE_CURRENT_LIST_DIR}/UnitTestsMain.cpp + + ${AUTOGENERATED_SOURCES} + ${BOOST_EXTENDED_SOURCES} + ${GOOGLE_TEST_SOURCES} + ${ORTHANC_STONE_SOURCES} + ) diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/WebAssembly/CMakeLists.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/WebAssembly/CMakeLists.txt Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,99 @@ +# Stone of Orthanc +# Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics +# Department, University Hospital of Liege, Belgium +# Copyright (C) 2017-2020 Osimis S.A., Belgium +# +# This program is free software: you can redistribute it and/or +# modify it under the terms of the GNU Affero General Public License +# as published by the Free Software Foundation, either version 3 of +# the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . + + +cmake_minimum_required(VERSION 2.8.3) +cmake_policy(SET CMP0058 NEW) + +project(OrthancStone) + +set(ORTHANC_STONE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/../../wasm-binaries/" CACHE PATH "Where to put the WebAssembly binaries") + + +# Configuration of the Emscripten compiler for WebAssembly target +# --------------------------------------------------------------- +set(USE_WASM ON CACHE BOOL "") + +set(WASM_FLAGS "-s WASM=1 -s FETCH=1 -s ASSERTIONS=1 -s DISABLE_EXCEPTION_CATCHING=0") +if (CMAKE_BUILD_TYPE STREQUAL "Debug") + set(WASM_FLAGS "${WASM_FLAGS} -s SAFE_HEAP=1") +endif() + +set(WASM_LINKER_FLAGS "${WASM_LINKER_FLAGS} -s EXTRA_EXPORTED_RUNTIME_METHODS='[\"ccall\", \"cwrap\"]'") +set(WASM_LINKER_FLAGS "${WASM_LINKER_FLAGS} -s ERROR_ON_UNDEFINED_SYMBOLS=1") +set(WASM_LINKER_FLAGS "${WASM_LINKER_FLAGS} -s ALLOW_MEMORY_GROWTH=1 -s TOTAL_MEMORY=268435456") # 256MB + resize +set(WASM_LINKER_FLAGS "${WASM_LINKER_FLAGS} -s DISABLE_DEPRECATED_FIND_EVENT_TARGET_BEHAVIOR=1") +add_definitions( + -DDISABLE_DEPRECATED_FIND_EVENT_TARGET_BEHAVIOR=1 +) + + +# Stone of Orthanc configuration +# --------------------------------------------------------------- +set(ALLOW_DOWNLOADS ON) + +include(${CMAKE_SOURCE_DIR}/../../Applications/Platforms/WebAssembly/OrthancStoneWebAssemblyParameters.cmake) + +SET(ENABLE_DCMTK OFF) # Not necessary +SET(ENABLE_GOOGLE_TEST OFF) +SET(ENABLE_LOCALE ON) # Necessary for text rendering +SET(ORTHANC_SANDBOXED ON) + +# Needed to redirect std::cout to a
+set(ORTHANC_ENABLE_LOGGING_STDIO OFF CACHE INTERNAL "") + +# this will set up the build system for Stone of Orthanc and will +# populate the ORTHANC_STONE_SOURCES CMake variable +include(${CMAKE_SOURCE_DIR}/../../Applications/Platforms/WebAssembly/OrthancStoneWebAssemblyConfiguration.cmake) + + +################################################################################ + +# Define the WASM module +# --------------------------------------------------------------- + +set(USE_SYSTEM_GOOGLE_TEST OFF CACHE BOOL "Use the system version of Google Test") +set(USE_GOOGLE_TEST_DEBIAN_PACKAGE OFF CACHE BOOL "Use the sources of Google Test shipped with libgtest-dev (Debian only)") +mark_as_advanced(USE_GOOGLE_TEST_DEBIAN_PACKAGE) +include(${ORTHANC_STONE_ROOT}/../Resources/Orthanc/CMake/DownloadPackage.cmake) +include(${ORTHANC_STONE_ROOT}/../Resources/Orthanc/CMake/GoogleTestConfiguration.cmake) + + +include(${CMAKE_SOURCE_DIR}/../UnitTestsSources.cmake) +add_executable(UnitTests ${UNIT_TESTS_SOURCES}) + + +# Declare installation files for the module +# --------------------------------------------------------------- +install( + TARGETS UnitTests + RUNTIME DESTINATION ${ORTHANC_STONE_INSTALL_PREFIX}/UnitTests/ + ) + +# Declare installation files for the companion files (web scaffolding) +# please note that ${CMAKE_CURRENT_BINARY_DIR}/RtViewerWasm.js +# (the generated JS loader for the WASM module) is handled by the `install1` +# section above: it is considered to be the binary output of +# the linker. +# --------------------------------------------------------------- +install( + FILES + ${CMAKE_SOURCE_DIR}/index.html + ${CMAKE_CURRENT_BINARY_DIR}/UnitTests.wasm + DESTINATION ${ORTHANC_STONE_INSTALL_PREFIX}/UnitTests + ) diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/WebAssembly/NOTES.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/WebAssembly/NOTES.txt Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,18 @@ +Native compilation (without Docker) +=================================== + +Install Emscripten: +https://emscripten.org/docs/getting_started/downloads.html + +Then, if the installation path is "~/Downloads/emsdk/": + +# source ~/Downloads/emsdk/emsdk_env.sh +# mkdir Build && cd Build +# cmake .. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_TOOLCHAIN_FILE=${EMSDK}/upstream/emscripten/cmake/Modules/Platform/Emscripten.cmake -DALLOW_DOWNLOADS=ON -G Ninja +# ninja install + +=> The binaries will be put in "../../../wasm-binaries/UnitTests/" + +# cd `pwd`/../../../wasm-binaries/UnitTests +# python3 -m http.server 8000 +# firefox http://localhost:8000/index.html diff -r 1393e3393a0b -r 5b8b88e5bfd6 UnitTestsSources/WebAssembly/index.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnitTestsSources/WebAssembly/index.html Tue Nov 24 12:59:10 2020 +0100 @@ -0,0 +1,14 @@ + + + + + Stone of Orthanc + + +

Stone of Orthanc - Unit tests

+ 
+      Executing the tests...
+
+ + +