Mercurial > hg > orthanc-stl
changeset 1:0f03a8a0bd6f
encoding of RT-STRUCT as STL
| author | Sebastien Jodogne <s.jodogne@gmail.com> |
|---|---|
| date | Mon, 17 Jul 2023 18:54:31 +0200 |
| parents | 4e889a8e8be2 |
| children | 2bdb9acb7dcf |
| files | CMakeLists.txt Resources/CMake/LibVtkConfiguration.cmake Sources/OrthancExplorer.js Sources/Plugin.cpp Sources/viewer.js |
| diffstat | 5 files changed, 1355 insertions(+), 11 deletions(-) [+] |
line wrap: on
line diff
--- a/CMakeLists.txt Tue Jul 04 19:15:27 2023 +0200 +++ b/CMakeLists.txt Mon Jul 17 18:54:31 2023 +0200 @@ -19,7 +19,20 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. -cmake_minimum_required(VERSION 2.8) +# Native static compilation: +# $ cmake .. -DCMAKE_BUILD_TYPE=Debug -DSTATIC_BUILD=ON && make -j4 + +# LSB compilation: +# $ LSB_CC=gcc-4.8 LSB_CXX=g++-4.8 cmake .. -DCMAKE_BUILD_TYPE=Release -DSTATIC_BUILD=ON -DUSE_LEGACY_JSONCPP=ON -DUSE_LEGACY_BOOST=ON -DCMAKE_TOOLCHAIN_FILE=${HOME}/Subversion/orthanc/OrthancFramework/Resources/Toolchains/LinuxStandardBaseToolchain.cmake && make -j4 + +# MinGW compilation: +# $ cmake .. -DCMAKE_BUILD_TYPE=Release -DSTATIC_BUILD=ON -DUSE_LEGACY_JSONCPP=ON -DUSE_LEGACY_BOOST=ON -DCMAKE_TOOLCHAIN_FILE=${HOME}/Subversion/orthanc/OrthancFramework/Resources/Toolchains/MinGW-W64-Toolchain32.cmake && make -j4 + + +# We request version 3.12 of CMake, as it is the first version to +# support "--parallel" in "cmake --build": +# https://cmake.org/cmake/help/latest/manual/cmake.1.html#build-a-project +cmake_minimum_required(VERSION 3.12) project(OrthancSTL) set(ORTHANC_STL_VERSION "mainline") @@ -48,6 +61,9 @@ # Advanced parameters to fine-tune linking against system libraries SET(USE_SYSTEM_ORTHANC_SDK ON CACHE BOOL "Use the system version of the Orthanc plugin SDK") +# Parameters to fine-tune linking against system libraries +set(USE_SYSTEM_VTK ON CACHE BOOL "Use the system version of VTK") + ##################################################################### ## Download and setup the Orthanc framework @@ -72,11 +88,10 @@ else() include(${ORTHANC_FRAMEWORK_ROOT}/../Resources/CMake/OrthancFrameworkParameters.cmake) - set(ENABLE_LOCALE OFF) # Disable support for locales (notably in Boost) + set(ENABLE_LOCALE ON) set(ENABLE_ZLIB ON) + set(ENABLE_DCMTK ON) - set(ENABLE_MODULE_DICOM OFF CACHE INTERNAL "") - set(ENABLE_MODULE_IMAGES OFF CACHE INTERNAL "") set(ENABLE_MODULE_JOBS OFF CACHE INTERNAL "") include(${ORTHANC_FRAMEWORK_ROOT}/../Resources/CMake/OrthancFrameworkConfiguration.cmake) @@ -104,6 +119,8 @@ ## Platform-specific configuration ##################################################################### +include(${CMAKE_SOURCE_DIR}/Resources/CMake/LibVtkConfiguration.cmake) + add_definitions( -DHAS_ORTHANC_EXCEPTION=1 -DORTHANC_ENABLE_LOGGING_PLUGIN=1 @@ -183,10 +200,14 @@ ${CMAKE_SOURCE_DIR}/Resources/Orthanc/Plugins/OrthancPluginCppWrapper.cpp ${ORTHANC_CORE_SOURCES_DEPENDENCIES} ${ORTHANC_CORE_SOURCES_INTERNAL} + ${ORTHANC_DICOM_SOURCES_DEPENDENCIES} + ${ORTHANC_DICOM_SOURCES_INTERNAL} ${ORTHANC_STL_SOURCES} ) -add_dependencies(OrthancSTL AutogeneratedTarget) +DefineSourceBasenameForTarget(OrthancSTL) + +add_dependencies(OrthancSTL AutogeneratedTarget VTK) message("Setting the version of the library to ${ORTHANC_STL_VERSION}") @@ -194,6 +215,8 @@ VERSION ${ORTHANC_STL_VERSION} SOVERSION ${ORTHANC_STL_VERSION}) +target_link_libraries(OrthancSTL ${DCMTK_LIBRARIES} ${VTK_LIBRARIES}) + install( TARGETS OrthancSTL RUNTIME DESTINATION lib # Destination for Windows
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Resources/CMake/LibVtkConfiguration.cmake Mon Jul 17 18:54:31 2023 +0200 @@ -0,0 +1,190 @@ +# SPDX-FileCopyrightText: 2023 Sebastien Jodogne, UCLouvain, Belgium +# SPDX-License-Identifier: GPL-3.0-or-later + + +# STL plugin for Orthanc +# Copyright (C) 2023 Sebastien Jodogne, UCLouvain, Belgium +# +# This program is free software: you can redistribute it and/or +# modify it under the terms of the GNU 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 +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. + + +if (NOT STATIC_BUILD AND USE_SYSTEM_VTK) + find_package(VTK REQUIRED COMPONENTS + vtkCommonDataModel + vtkFiltersCore + vtkImagingCore + ) + + include_directories(${VTK_INCLUDE_DIRS}) + + # Add a void target + add_custom_target(VTK) + +else() + set(VTK_MAJOR 7) + set(VTK_MINOR 1) + set(VTK_REVISION 1) + set(VTK_MD5 "daee43460f4e95547f0635240ffbc9cb") + + set(VTK_SOURCES_DIR "${CMAKE_BINARY_DIR}/VTK-${VTK_MAJOR}.${VTK_MINOR}.${VTK_REVISION}") + + DownloadPackage( + "${VTK_MD5}" + "https://orthanc.uclouvain.be/third-party-downloads/VTK-${VTK_MAJOR}.${VTK_MINOR}.${VTK_REVISION}.tar.gz" + "${VTK_SOURCES_DIR}") + + if (CMAKE_TOOLCHAIN_FILE) + # Take absolute path to the toolchain + get_filename_component(TMP ${CMAKE_TOOLCHAIN_FILE} REALPATH BASE ${CMAKE_SOURCE_DIR}/..) + list(APPEND VTKCMakeFlags + -DCMAKE_TOOLCHAIN_FILE=${TMP} + -DLSB_CC=${LSB_CC} + -DLSB_CXX=${LSB_CXX} + ) + endif() + + cmake_host_system_information(RESULT NumberOfPhysicalCores QUERY NUMBER_OF_PHYSICAL_CORES) + + include(ExternalProject) + + if (CMAKE_CROSSCOMPILING) + # First, generate the "VTKCompileToolsConfig.cmake" file on the host + # computer. This is necessary before cross-compiling. Explanations: + # https://cmake.org/cmake/help/book/mastering-cmake/chapter/Cross%20Compiling%20With%20CMake.html#cross-compiling-a-complex-project-vtk + externalproject_add(VTKCompileTools + SOURCE_DIR "${VTK_SOURCES_DIR}" + + CMAKE_ARGS + -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE} + -DCMAKE_CXX_COMPILER=g++ + -DCMAKE_C_COMPILER=gcc + + BUILD_COMMAND + ${CMAKE_COMMAND} --build <BINARY_DIR> --config <CONFIG> --parallel ${NumberOfPhysicalCores} -t + VTKCompileToolsConfig.cmake + + INSTALL_COMMAND "" # Skip the install step + ) + + ExternalProject_Get_Property(VTKCompileTools binary_dir) + list(APPEND VTKCMakeFlags + -DVTKCompileTools_DIR=${binary_dir} + ) + + list(APPEND VTKCMakeFlags + -DCMAKE_REQUIRE_LARGE_FILE_SUPPORT=ON + -DCMAKE_REQUIRE_LARGE_FILE_SUPPORT__TRYRUN_OUTPUT="" + -DKWSYS_LFS_WORKS=ON + -DKWSYS_LFS_WORKS__TRYRUN_OUTPUT="" + -DFILE_OFFSET_BITS=64 + + # This simply disables HDF5, which is not used in Orthanc + -DVTK_USE_SYSTEM_HDF5=ON + ) + else() + add_custom_target(VTKCompileTools) # Empty target + endif() + + if (CMAKE_COMPILER_IS_GNUCXX OR + CMAKE_CXX_COMPILER_ID STREQUAL "Clang") + # The following flag is necessary to statically link VTK within + # the shared library containing the plugin + set(Flags "-fPIC") + else() + set(Flags "") + endif() + + externalproject_add(VTK + SOURCE_DIR "${VTK_SOURCES_DIR}" + + CMAKE_ARGS + ${VTKCMakeFlags} + -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE} + -DBUILD_SHARED_LIBS=OFF + -DVTK_Group_Rendering=OFF + + -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} + "-DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS} ${Flags}" + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + "-DCMAKE_C_FLAGS=${CMAKE_C_FLAGS} ${Flags}" + -DCMAKE_OSX_DEPLOYMENT_TARGET=${CMAKE_OSX_DEPLOYMENT_TARGET} + -DCMAKE_OSX_ARCHITECTURES=${CMAKE_OSX_ARCHITECTURES} + + BUILD_COMMAND + ${CMAKE_COMMAND} --build <BINARY_DIR> --config <CONFIG> --parallel ${NumberOfPhysicalCores} -t + vtkCommonCore + vtkCommonDataModel + vtkCommonMath + vtkCommonMisc + vtkCommonSystem + vtkCommonTransforms + vtkFiltersCore + vtkFiltersPoints + vtkImagingCore + + EXCLUDE_FROM_ALL TRUE + + INSTALL_COMMAND "" # Skip the install step + + DEPENDS VTKCompileTools + ) + + if(MSVC) + set(Suffix ".lib") + set(Prefix "") + else() + set(Suffix ".a") + list(GET CMAKE_FIND_LIBRARY_PREFIXES 0 Prefix) + endif() + + foreach(module IN ITEMS + # WARNING: The order of the modules below *is* important! + ImagingCore + FiltersCore + CommonExecutionModel + CommonDataModel + CommonTransforms + CommonMath + CommonMisc + CommonSystem + CommonCore + sys + ) + list(APPEND VTK_LIBRARIES + ${Prefix}vtk${module}-${VTK_MAJOR}.${VTK_MINOR}${Suffix} + ) + endforeach() + + ExternalProject_Get_Property(VTK binary_dir) + ExternalProject_Get_Property(VTK source_dir) + + link_directories(${binary_dir}/lib) + + foreach(dir IN ITEMS + Common/Core + Common/DataModel + Common/ExecutionModel + Common/Misc + Filters/Core + IO/Image + Imaging/Core + Utilities/KWIML + ) + include_directories( + ${source_dir}/${dir} + ${binary_dir}/${dir} + ) + endforeach() + +endif()
--- a/Sources/OrthancExplorer.js Tue Jul 04 19:15:27 2023 +0200 +++ b/Sources/OrthancExplorer.js Mon Jul 17 18:54:31 2023 +0200 @@ -50,7 +50,7 @@ b.insertAfter($('#series-info')); b.click(function() { if ($.mobile.pageData) { - window.open('../stl/viewer.html?instance=' + instanceId); + window.open('../stl/app/viewer.html?instance=' + instanceId); } }); }
--- a/Sources/Plugin.cpp Tue Jul 04 19:15:27 2023 +0200 +++ b/Sources/Plugin.cpp Mon Jul 17 18:54:31 2023 +0200 @@ -24,11 +24,26 @@ #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" +#include <ChunkedBuffer.h> +#include <DicomParsing/ParsedDicomFile.h> +#include <Images/ImageProcessing.h> #include <Logging.h> +#include <OrthancFramework.h> +#include <SerializationToolbox.h> #include <SystemToolbox.h> #include <EmbeddedResources.h> +#include <vtkImageData.h> +#include <vtkImageResize.h> +#include <vtkMarchingCubes.h> +#include <vtkNew.h> +#include <vtkPolyData.h> +#include <vtkTriangle.h> +#include <vtkSmoothPolyDataFilter.h> +#include <vtkPolyDataNormals.h> +#include <vtkImageConstantPad.h> + #include <boost/thread/shared_mutex.hpp> // Forward declaration @@ -125,6 +140,1107 @@ } + + +#include <dcmtk/dcmdata/dcdeftag.h> +#include <dcmtk/dcmdata/dcfilefo.h> +#include <dcmtk/dcmdata/dcitem.h> +#include <dcmtk/dcmdata/dcsequen.h> + +class Extent2D : public boost::noncopyable +{ +private: + bool isEmpty_; + double x1_; + double y1_; + double x2_; + double y2_; + + void CheckNotEmpty() const + { + if (isEmpty_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + +public: + Extent2D() : + isEmpty_(true), + x1_(0), + y1_(0), + x2_(0), + y2_(0) + { + } + + bool IsEmpty() const + { + return isEmpty_; + } + + double GetMinX() const + { + CheckNotEmpty(); + return x1_; + } + + double GetMaxX() const + { + CheckNotEmpty(); + return x2_; + } + + double GetMinY() const + { + CheckNotEmpty(); + return y1_; + } + + double GetMaxY() const + { + CheckNotEmpty(); + return y2_; + } + + double GetWidth() const + { + CheckNotEmpty(); + return x2_ - x1_; + } + + double GetHeight() const + { + CheckNotEmpty(); + return y2_ - y1_; + } + + void Add(double x, + double y) + { + if (isEmpty_) + { + x1_ = x2_ = x; + y1_ = y2_ = y; + isEmpty_ = false; + } + else + { + x1_ = std::min(x1_, x); + x2_ = std::max(x2_, x); + y1_ = std::min(y1_, y); + y2_ = std::max(y2_, y); + } + } +}; + +static void GetStringValue(std::string& value, + DcmItem& item, + const DcmTagKey& key) +{ + const char* s = NULL; + if (!item.findAndGetString(key, s).good() || + s == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + value.assign(s); + } +} + + +static void ListStructuresNames(std::set<std::string>& target, + Orthanc::ParsedDicomFile& source) +{ + target.clear(); + + DcmSequenceOfItems* sequence = NULL; + if (!source.GetDcmtkObject().getDataset()->findAndGetSequence(DCM_StructureSetROISequence, sequence).good() || + sequence == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + for (unsigned long i = 0; i < sequence->card(); i++) + { + DcmItem* item = sequence->getItem(i); + if (item == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + std::string value; + GetStringValue(value, *item, DCM_ROIName); + target.insert(value); + } + } +} + + +static bool IsNear(double a, + double b) +{ + return std::abs(a - b) < 10.0 * std::numeric_limits<double>::epsilon(); +} + + +class Vector3D +{ +private: + double x_; + double y_; + double z_; + +public: + Vector3D() : + x_(0), + y_(0), + z_(0) + { + } + + Vector3D(double x, + double y, + double z) : + x_(x), + y_(y), + z_(z) + { + } + + Vector3D(const Vector3D& from, + const Vector3D& to) : + x_(to.x_ - from.x_), + y_(to.y_ - from.y_), + z_(to.z_ - from.z_) + { + } + + double GetX() const + { + return x_; + } + + double GetY() const + { + return y_; + } + + double GetZ() const + { + return z_; + } + + double ComputeNorm() const + { + return sqrt(x_ * x_ + y_ * y_ + z_ * z_); + } + + void Normalize() + { + double norm = ComputeNorm(); + if (!IsNear(norm, 0)) + { + x_ /= norm; + y_ /= norm; + z_ /= norm; + } + } + + static Vector3D CrossProduct(const Vector3D& u, + const Vector3D& v) + { + return Vector3D(u.GetY() * v.GetZ() - u.GetZ() * v.GetY(), + u.GetZ() * v.GetX() - u.GetX() * v.GetZ(), + u.GetX() * v.GetY() - u.GetY() * v.GetX()); + } + + static double DotProduct(const Vector3D& a, + const Vector3D& b) + { + return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ(); + } +}; + + + +static bool MyParseDouble(double& value, + const std::string& s) +{ +#if 1 + char* end = NULL; + value = strtod(s.c_str(), &end); + return (end == s.c_str() + s.size()); +#else + return Orthanc::SerializationToolbox::ParseDouble(value, s); +#endif +} + + +class StructurePolygon : public boost::noncopyable +{ +private: + std::string roiName_; + std::string referencedSopInstanceUid_; + uint8_t red_; + uint8_t green_; + uint8_t blue_; + std::vector<Vector3D> points_; + +public: + StructurePolygon(Orthanc::ParsedDicomFile& dicom, + unsigned long roiIndex, + unsigned long contourIndex) + { + DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); + + DcmItem* structure = NULL; + DcmItem* roi = NULL; + DcmItem* contour = NULL; + DcmSequenceOfItems* referenced = NULL; + + if (!dataset.findAndGetSequenceItem(DCM_StructureSetROISequence, structure, roiIndex).good() || + structure == NULL || + !dataset.findAndGetSequenceItem(DCM_ROIContourSequence, roi, roiIndex).good() || + roi == NULL || + !roi->findAndGetSequenceItem(DCM_ContourSequence, contour, contourIndex).good() || + contour == NULL || + !contour->findAndGetSequence(DCM_ContourImageSequence, referenced).good() || + referenced == NULL || + referenced->card() != 1) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + GetStringValue(roiName_, *structure, DCM_ROIName); + GetStringValue(referencedSopInstanceUid_, *referenced->getItem(0), DCM_ReferencedSOPInstanceUID); + + std::string s; + GetStringValue(s, *contour, DCM_ContourGeometricType); + if (s != "CLOSED_PLANAR") + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + { + std::string color; + GetStringValue(color, *roi, DCM_ROIDisplayColor); + + std::vector<std::string> tokens; + Orthanc::Toolbox::TokenizeString(tokens, color, '\\'); + + uint32_t r, g, b; + if (tokens.size() != 3 || + !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(r, tokens[0]) || + !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(g, tokens[1]) || + !Orthanc::SerializationToolbox::ParseFirstUnsignedInteger32(b, tokens[2]) || + r > 255 || + g > 255 || + b > 255) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + red_ = r; + green_ = g; + blue_ = b; + } + + { + GetStringValue(s, *contour, DCM_ContourData); + + std::vector<std::string> tokens; + Orthanc::Toolbox::TokenizeString(tokens, s, '\\'); + + GetStringValue(s, *contour, DCM_NumberOfContourPoints); + + uint32_t countPoints; + if (!Orthanc::SerializationToolbox::ParseUnsignedInteger32(countPoints, s) || + tokens.size() != 3 * countPoints) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + points_.reserve(countPoints); + + for (size_t i = 0; i < tokens.size(); i += 3) + { + double x, y, z; + if (!MyParseDouble(x, tokens[i]) || + !MyParseDouble(y, tokens[i + 1]) || + !MyParseDouble(z, tokens[i + 2])) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + points_.push_back(Vector3D(x, y, z)); + } + + assert(points_.size() == countPoints); + } + } + + const std::string& GetRoiName() const + { + return roiName_; + } + + const std::string& GetReferencedSopInstanceUid() const + { + return referencedSopInstanceUid_; + } + + size_t GetPointsCount() const + { + return points_.size(); + } + + const Vector3D& GetPoint(size_t i) const + { + if (i >= points_.size()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + else + { + return points_[i]; + } + } + + bool IsCoplanar(Vector3D& normal) const + { + if (points_.size() < 3) + { + return false; + } + + bool hasNormal = false; + + for (size_t i = 0; i < points_.size(); i++) + { + normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]), + Vector3D(points_[2], points_[0])); + if (!IsNear(normal.ComputeNorm(), 0)) + { + normal.Normalize(); + hasNormal = true; + } + } + + if (!hasNormal) + { + return false; + } + + double a = Vector3D::DotProduct(points_[0], normal); + + for (size_t i = 1; i < points_.size(); i++) + { + double b = Vector3D::DotProduct(points_[i], normal); + if (!IsNear(a, b)) + { + return false; + } + } + + return true; + } + + void Add(Extent2D& extent, + const Vector3D& axisX, + const Vector3D& axisY) const + { + assert(IsNear(1, axisX.ComputeNorm())); + assert(IsNear(1, axisY.ComputeNorm())); + + for (size_t i = 0; i < points_.size(); i++) + { + extent.Add(Vector3D::DotProduct(axisX, points_[i]), + Vector3D::DotProduct(axisY, points_[i])); + } + } +}; + + + +struct IsNearPredicate +{ + bool operator() (const double& a, + const double& b) + { + return IsNear(a, b); + } +}; + +static void RemoveDuplicateValues(std::vector<double>& v) +{ + IsNearPredicate predicate; + std::vector<double>::iterator last = std::unique(v.begin(), v.end(), predicate); + v.erase(last, v.end()); +} + + +class StructureSet : public boost::noncopyable +{ +private: + std::vector<StructurePolygon*> polygons_; + bool hasGeometry_; + Vector3D slicesNormal_; + double slicesSpacing_; + double minProjectionAlongNormal_; + double maxProjectionAlongNormal_; + std::string patientId_; + std::string studyInstanceUid_; + std::string seriesInstanceUid_; + std::string sopInstanceUid_; + + void ComputeGeometry() + { + std::list<double> positionsList; + hasGeometry_ = false; + + for (size_t i = 0; i < polygons_.size(); i++) + { + assert(polygons_[i] != NULL); + + Vector3D n; + if (polygons_[i]->IsCoplanar(n)) + { + const Vector3D& point = polygons_[i]->GetPoint(0); + double z = Vector3D::DotProduct(point, n); + + if (!hasGeometry_) + { + hasGeometry_ = true; + slicesNormal_ = n; + minProjectionAlongNormal_ = z; + maxProjectionAlongNormal_ = z; + } + else if (!IsNear(std::abs(Vector3D::DotProduct(n, slicesNormal_)), 1)) + { + hasGeometry_ = false; + + // RT-STRUCT with non-parallel slices + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + else + { + minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, z); + maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, z); + } + + positionsList.push_back(Vector3D::DotProduct(n, point)); + } + } + + if (hasGeometry_) + { + std::vector<double> positions(positionsList.begin(), positionsList.end()); + assert(!positions.empty()); + + std::sort(positions.begin(), positions.end()); + RemoveDuplicateValues(positions); + assert(!positions.empty()); + + if (positions.size() == 1) + { + hasGeometry_ = false; + return; + } + + std::vector<double> offsets; + offsets.resize(positions.size() - 1); + + for (size_t i = 0; i < offsets.size(); i++) + { + offsets[i] = positions[i + 1] - positions[i]; + assert(offsets[i] > 0); + } + + std::sort(offsets.begin(), offsets.end()); + RemoveDuplicateValues(offsets); + + slicesSpacing_ = offsets[0]; + + for (size_t i = 1; i < offsets.size(); i++) + { + double d = offsets[i] / slicesSpacing_; + if (!IsNear(d, round(d))) + { + // Irregular spacing between the slices + hasGeometry_ = false; + break; + } + } + } + } + + void CheckHasGeometry() const + { + if (!hasGeometry_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + +public: + explicit StructureSet(Orthanc::ParsedDicomFile& dicom) : + hasGeometry_(false), + slicesSpacing_(0), + minProjectionAlongNormal_(0), + maxProjectionAlongNormal_(0) + { + DcmDataset& dataset = *dicom.GetDcmtkObject().getDataset(); + GetStringValue(patientId_, dataset, DCM_PatientID); + GetStringValue(studyInstanceUid_, dataset, DCM_StudyInstanceUID); + GetStringValue(seriesInstanceUid_, dataset, DCM_SeriesInstanceUID); + GetStringValue(sopInstanceUid_, dataset, DCM_SOPInstanceUID); + + DcmSequenceOfItems* rois = NULL; + if (!dataset.findAndGetSequence(DCM_ROIContourSequence, rois).good() || + rois == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + std::vector<DcmSequenceOfItems*> contours(rois->card()); + size_t countPolygons = 0; + + for (unsigned long i = 0; i < rois->card(); i++) + { + DcmSequenceOfItems* contour = NULL; + if (!rois->getItem(i)->findAndGetSequence(DCM_ContourSequence, contour).good() || + contour == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + contours[i] = contour; + countPolygons += contour->card(); + } + } + + polygons_.resize(countPolygons); + + size_t pos = 0; + for (unsigned long i = 0; i < contours.size(); i++) + { + for (unsigned long j = 0; j < contours[i]->card(); j++, pos++) + { + polygons_[pos] = new StructurePolygon(dicom, i, j); + } + } + + assert(pos == countPolygons); + + ComputeGeometry(); + } + + ~StructureSet() + { + for (size_t i = 0; i < polygons_.size(); i++) + { + assert(polygons_[i] != NULL); + delete polygons_[i]; + } + } + + const std::string& GetPatient() const + { + return patientId_; + } + + const std::string& GetStudyInstanceUid() const + { + return studyInstanceUid_; + } + + const std::string& GetSeriesInstanceUid() const + { + return seriesInstanceUid_; + } + + const std::string& GetSopInstanceUid() const + { + return sopInstanceUid_; + } + + std::string HashStudy() const + { + std::string s; + Orthanc::Toolbox::ComputeSHA1(s, patientId_ + "|" + studyInstanceUid_); + return s; + } + + size_t GetPolygonsCount() const + { + return polygons_.size(); + } + + const StructurePolygon& GetPolygon(size_t i) const + { + if (i >= polygons_.size()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + else + { + assert(polygons_[i] != NULL); + return *polygons_[i]; + } + } + + bool HasGeometry() const + { + return hasGeometry_; + } + + const Vector3D& GetSlicesNormal() const + { + CheckHasGeometry(); + return slicesNormal_; + } + + double GetSlicesSpacing() const + { + CheckHasGeometry(); + return slicesSpacing_; + } + + double GetMinProjectionAlongNormal() const + { + CheckHasGeometry(); + return minProjectionAlongNormal_; + } + + double GetMaxProjectionAlongNormal() const + { + CheckHasGeometry(); + return maxProjectionAlongNormal_; + } + + double ProjectAlongNormal(const StructurePolygon& polygon) const + { + CheckHasGeometry(); + return Vector3D::DotProduct(slicesNormal_, polygon.GetPoint(0)); + } + + size_t GetSlicesCount() const + { + CheckHasGeometry(); + double c = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; + assert(c >= 0); + + if (!IsNear(c, round(c))) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + else + { + return static_cast<size_t>(round(c)) + 1; + } + } + + bool LookupSliceIndex(size_t& slice, + const StructurePolygon& polygon) const + { + CheckHasGeometry(); + + double z = ProjectAlongNormal(polygon); + + if (z < minProjectionAlongNormal_ || + z > maxProjectionAlongNormal_) + { + return false; + } + else + { + double c = (z - minProjectionAlongNormal_) / slicesSpacing_; + + if (IsNear(c, round(c))) + { + slice = static_cast<size_t>(round(c)); + return true; + } + else + { + return false; + } + } + } + + bool LookupReferencedSopInstanceUid(std::string& sopInstanceUid) const + { + if (HasGeometry()) + { + for (size_t i = 0; i < polygons_.size(); i++) + { + assert(polygons_[i] != NULL); + + Vector3D n; + if (polygons_[i]->IsCoplanar(n) && + Vector3D::DotProduct(n, slicesNormal_)) + { + sopInstanceUid = polygons_[i]->GetReferencedSopInstanceUid(); + return true; + } + } + } + + return false; + } +}; + + + + +class XorFiller : public Orthanc::ImageProcessing::IPolygonFiller +{ +private: + Orthanc::ImageAccessor& target_; + +public: + XorFiller(Orthanc::ImageAccessor& target) : + target_(target) + { + } + + virtual void Fill(int y, + int x1, + int x2) ORTHANC_OVERRIDE + { + assert(x1 <= x2); + + if (y >= 0 && + y < static_cast<int>(target_.GetHeight())) + { + x1 = std::max(x1, 0); + x2 = std::min(x2, static_cast<int>(target_.GetWidth()) - 1); + + uint8_t* p = reinterpret_cast<uint8_t*>(target_.GetRow(y)) + x1; + + for (int i = x1; i <= x2; i++, p++) + { + *p = (*p ^ 0xff); + } + } + } +}; + + +static void EncodeSTL(std::string& target /* out */, + vtkPolyData& mesh /* in */) +{ + // TODO - Conversion to little endian on big endian + + Orthanc::ChunkedBuffer buffer; + + uint8_t header[80]; + memset(header, 0, sizeof(header)); + buffer.AddChunk(header, sizeof(header)); + + uint32_t n = mesh.GetNumberOfCells(); + buffer.AddChunk(&n, sizeof(n)); + + for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++) + { + vtkCell* cell = mesh.GetCell(i); + vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell); + + double p0[3]; + double p1[3]; + double p2[3]; + triangle->GetPoints()->GetPoint(0, p0); + triangle->GetPoints()->GetPoint(1, p1); + triangle->GetPoints()->GetPoint(2, p2); + + double normal[3]; + vtkTriangle::ComputeNormal(p0, p1, p2, normal); + + float d[4 * 3] = { + static_cast<float>(normal[0]), static_cast<float>(normal[1]), static_cast<float>(normal[2]), + static_cast<float>(p0[0]), static_cast<float>(p0[1]), static_cast<float>(p0[2]), + static_cast<float>(p1[0]), static_cast<float>(p1[1]), static_cast<float>(p1[2]), + static_cast<float>(p2[0]), static_cast<float>(p2[1]), static_cast<float>(p2[2]) }; + buffer.AddChunk(d, sizeof(d)); + + uint16_t a = 0; + buffer.AddChunk(&a, sizeof(a)); + } + + buffer.Flatten(target); +} + + +bool EncodeStructureSetMesh(std::string& stl, + const StructureSet& structureSet, + const std::set<std::string>& roiNames, + unsigned int resolution, + bool smooth) +{ + if (!structureSet.HasGeometry()) + { + return false; + } + + if (resolution < 1) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + if (!IsNear(1, structureSet.GetSlicesNormal().ComputeNorm())) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + // TODO - Axes could be retrieved from the referenced CT volume + Vector3D axisX(1, 0, 0); + Vector3D axisY = Vector3D::CrossProduct(structureSet.GetSlicesNormal(), axisX); + + if (!IsNear(1, axisX.ComputeNorm()) || + !IsNear(1, axisY.ComputeNorm())) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + Extent2D extent; + for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++) + { + structureSet.GetPolygon(i).Add(extent, axisX, axisY); + } + + const int depth = structureSet.GetSlicesCount(); + + vtkNew<vtkImageData> volume; + volume->SetDimensions(resolution, resolution, depth); + volume->AllocateScalars(VTK_UNSIGNED_CHAR, 1); + + assert(sizeof(unsigned char) == 1); + memset(volume->GetScalarPointer(), 0, resolution * resolution * depth); + + for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++) + { + const StructurePolygon& polygon = structureSet.GetPolygon(i); + if (roiNames.find(polygon.GetRoiName()) == roiNames.end()) + { + // This polygon doesn't correspond to a ROI of interest + continue; + } + + size_t j; + if (!structureSet.LookupSliceIndex(j, polygon)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + std::vector<Orthanc::ImageProcessing::ImagePoint> points; + points.reserve(polygon.GetPointsCount()); + for (size_t j = 0; j < polygon.GetPointsCount(); j++) + { + const Vector3D& point = polygon.GetPoint(j); + double x = (Vector3D::DotProduct(point, axisX) - extent.GetMinX()) / extent.GetWidth() * static_cast<double>(resolution); + double y = (Vector3D::DotProduct(point, axisY) - extent.GetMinY()) / extent.GetHeight() * static_cast<double>(resolution); + points.push_back(Orthanc::ImageProcessing::ImagePoint(static_cast<int32_t>(std::floor(x)), + static_cast<int32_t>(std::floor(y)))); + } + + Orthanc::ImageAccessor slice; + slice.AssignWritable(Orthanc::PixelFormat_Grayscale8, resolution, resolution, resolution /* pitch */, + reinterpret_cast<uint8_t*>(volume->GetScalarPointer()) + j * resolution * resolution); + + XorFiller filler(slice); + Orthanc::ImageProcessing::FillPolygon(filler, points); + } + + vtkNew<vtkImageResize> resize; + resize->SetOutputDimensions(resolution, resolution, resolution); + resize->SetInputData(volume.Get()); + resize->Update(); + + resize->GetOutput()->SetSpacing( + extent.GetWidth() / static_cast<double>(resolution), + extent.GetHeight() / static_cast<double>(resolution), + (structureSet.GetMaxProjectionAlongNormal() - structureSet.GetMinProjectionAlongNormal()) / static_cast<double>(resolution)); + + // TODO + // resize->GetOutput()->SetOrigin() + + vtkNew<vtkImageConstantPad> padding; + padding->SetConstant(0); + padding->SetOutputNumberOfScalarComponents(1); + padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution); + padding->SetInputData(resize->GetOutput()); + padding->Update(); + + vtkNew<vtkMarchingCubes> surface; + surface->SetInputData(padding->GetOutput()); + surface->ComputeNormalsOn(); + surface->SetValue(0, 128 /*isoValue*/); + surface->Update(); + + if (smooth) + { + vtkNew<vtkSmoothPolyDataFilter> smoothFilter; + // Apply volume smoothing + // https://examples.vtk.org/site/Cxx/PolyData/SmoothPolyDataFilter/ + smoothFilter->SetInputConnection(surface->GetOutputPort()); + smoothFilter->SetNumberOfIterations(15); + smoothFilter->SetRelaxationFactor(0.1); + smoothFilter->FeatureEdgeSmoothingOff(); + smoothFilter->BoundarySmoothingOn(); + smoothFilter->Update(); + + vtkNew<vtkPolyDataNormals> normalGenerator; + normalGenerator->SetInputConnection(smoothFilter->GetOutputPort()); + normalGenerator->ComputePointNormalsOn(); + normalGenerator->ComputeCellNormalsOn(); + normalGenerator->Update(); + + EncodeSTL(stl, *normalGenerator->GetOutput()); + } + else + { + EncodeSTL(stl, *surface->GetOutput()); + } + + return true; +} + + +void ListStructures(OrthancPluginRestOutput* output, + const char* url, + const OrthancPluginHttpRequest* request) +{ + const std::string instanceId(request->groups[0]); + + if (request->method != OrthancPluginHttpMethod_Get) + { + OrthancPluginSendMethodNotAllowed(OrthancPlugins::GetGlobalContext(), output, "GET"); + return; + } + + std::string dicom; + if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource); + } + else + { + Orthanc::ParsedDicomFile parsed(dicom); + + std::set<std::string> names; + ListStructuresNames(names, parsed); + + Json::Value answer = Json::arrayValue; + + for (std::set<std::string>::const_iterator it = names.begin(); it != names.end(); ++it) + { + answer.append(*it); + } + + std::string s = answer.toStyledString(); + OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::MIME_JSON); + } +} + + +void Encode(OrthancPluginRestOutput* output, + const char* url, + const OrthancPluginHttpRequest* request) +{ + static const char* const KEY_INSTANCE = "Instance"; + static const char* const KEY_RESOLUTION = "Resolution"; + static const char* const KEY_ROI_NAMES = "RoiNames"; + static const char* const KEY_SMOOTH = "Smooth"; + static const char* const KEY_TAGS = "Tags"; + + if (request->method != OrthancPluginHttpMethod_Post) + { + OrthancPluginSendMethodNotAllowed(OrthancPlugins::GetGlobalContext(), output, "POST"); + return; + } + + Json::Value body; + if (!Orthanc::Toolbox::ReadJson(body, request->body, request->bodySize)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest); + } + + const std::string instanceId = Orthanc::SerializationToolbox::ReadString(body, KEY_INSTANCE); + const bool smooth = (body.isMember(KEY_SMOOTH) ? + Orthanc::SerializationToolbox::ReadBoolean(body, KEY_SMOOTH) : + true /* smooth by default */); + const unsigned int resolution = (body.isMember(KEY_RESOLUTION) ? + Orthanc::SerializationToolbox::ReadUnsignedInteger(body, KEY_RESOLUTION) : + 256 /* default value */); + + std::set<std::string> roiNames; + Orthanc::SerializationToolbox::ReadSetOfStrings(roiNames, body, KEY_ROI_NAMES); + + std::string dicom; + if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource); + } + else + { + Orthanc::ParsedDicomFile parsed(dicom); + StructureSet structureSet(parsed); + + std::string stl; + if (!EncodeStructureSetMesh(stl, structureSet, roiNames, resolution, smooth)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, "Cannot encode STL"); + } + else + { + std::string content; + Orthanc::Toolbox::EncodeDataUriScheme(content, "model/stl", stl); + + Json::Value create; + create["Content"] = content; + create["Parent"] = structureSet.HashStudy(); + + if (body.isMember(KEY_TAGS)) + { + create[KEY_TAGS] = body[KEY_TAGS]; + } + else + { + std::string description; + for (std::set<std::string>::const_iterator it = roiNames.begin(); it != roiNames.end(); ++it) + { + if (!description.empty()) + { + description += ", "; + } + + description += *it; + } + + create[KEY_TAGS] = Json::objectValue; + create[KEY_TAGS]["SeriesDescription"] = description; + } + + Json::Value answer; + if (OrthancPlugins::RestApiPost(answer, "/tools/create-dicom", create.toStyledString(), false)) + { + std::string s = answer.toStyledString(); + OrthancPluginAnswerBuffer(OrthancPlugins::GetGlobalContext(), output, s.c_str(), s.size(), Orthanc::MIME_JSON); + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "Cannot create DICOM from STL"); + } + } + } +} + + extern "C" { ORTHANC_PLUGINS_API int32_t OrthancPluginInitialize(OrthancPluginContext* context) @@ -150,6 +1266,8 @@ Orthanc::Logging::Initialize(context); #endif + Orthanc::InitializeFramework("", false); + hasCreateDicomStl_ = OrthancPlugins::CheckMinimalOrthancVersion(1, 12, 1); if (!hasCreateDicomStl_) @@ -160,19 +1278,27 @@ OrthancPluginSetDescription(context, "STL plugin for Orthanc."); - OrthancPlugins::RegisterRestCallback<ServeFile>("/stl/(.*)", true); + OrthancPlugins::RegisterRestCallback<ServeFile>("/stl/app/(.*)", true); // Extend the default Orthanc Explorer with custom JavaScript for STL std::string explorer; Orthanc::EmbeddedResources::GetFileResource(explorer, Orthanc::EmbeddedResources::ORTHANC_EXPLORER); OrthancPluginExtendOrthancExplorer(OrthancPlugins::GetGlobalContext(), explorer.c_str()); + OrthancPlugins::RegisterRestCallback<ListStructures>("/stl/rt-struct/([0-9a-f-]+)", true); + + if (hasCreateDicomStl_) + { + OrthancPlugins::RegisterRestCallback<Encode>("/stl/encode", true); + } + return 0; } ORTHANC_PLUGINS_API void OrthancPluginFinalize() { + Orthanc::FinalizeFramework(); }
--- a/Sources/viewer.js Tue Jul 04 19:15:27 2023 +0200 +++ b/Sources/viewer.js Mon Jul 17 18:54:31 2023 +0200 @@ -65,7 +65,7 @@ const loader = new STLLoader() loader.load( - '../instances/' + instanceId + '/content/0042-0011', + '../../instances/' + instanceId + '/content/0042-0011', function (geometry) { const frustumSize = 200; @@ -73,9 +73,14 @@ geometry.translate(-(geometry.boundingBox.min.x + geometry.boundingBox.max.x) / 2.0, -(geometry.boundingBox.min.y + geometry.boundingBox.max.y) / 2.0, -(geometry.boundingBox.min.z + geometry.boundingBox.max.z) / 2.0); - geometry.scale((frustumSize / 2.0) / (geometry.boundingBox.max.x - geometry.boundingBox.min.x), - (frustumSize / 2.0) / (geometry.boundingBox.max.y - geometry.boundingBox.min.y), - (frustumSize / 2.0) / (geometry.boundingBox.max.z - geometry.boundingBox.min.z)); + + var maxSize = Math.max(geometry.boundingBox.max.x - geometry.boundingBox.min.x, + geometry.boundingBox.max.y - geometry.boundingBox.min.y, + geometry.boundingBox.max.z - geometry.boundingBox.min.z); + + geometry.scale((frustumSize / 2.0) / maxSize, + (frustumSize / 2.0) / maxSize, + (frustumSize / 2.0) / maxSize); const mesh = new THREE.Mesh(geometry, material); scene.add(mesh);