Mercurial > hg > orthanc-stl
changeset 32:976da5476810
reorganization
| author | Sebastien Jodogne <s.jodogne@gmail.com> |
|---|---|
| date | Thu, 04 Apr 2024 18:35:54 +0200 (13 months ago) |
| parents | ab231760799d |
| children | 2460b376d3f7 |
| files | .hgignore CMakeLists.txt Sources/Extent2D.cpp Sources/Extent2D.h Sources/Plugin.cpp Sources/Toolbox.cpp Sources/Toolbox.h Sources/VTKToolbox.cpp Sources/VTKToolbox.h Sources/Vector3D.cpp Sources/Vector3D.h |
| diffstat | 11 files changed, 870 insertions(+), 537 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgignore Thu Apr 04 18:05:13 2024 +0200 +++ b/.hgignore Thu Apr 04 18:35:54 2024 +0200 @@ -7,3 +7,4 @@ Three/three.js-*.tar.gz Three/dist i/ +*.orig
--- a/CMakeLists.txt Thu Apr 04 18:05:13 2024 +0200 +++ b/CMakeLists.txt Thu Apr 04 18:35:54 2024 +0200 @@ -203,7 +203,12 @@ ##################################################################### add_library(OrthancSTL SHARED + Sources/Extent2D.cpp Sources/Plugin.cpp + Sources/Toolbox.cpp + Sources/VTKToolbox.cpp + Sources/Vector3D.cpp + ${AUTOGENERATED_SOURCES} ${CMAKE_SOURCE_DIR}/Resources/Orthanc/Plugins/OrthancPluginCppWrapper.cpp ${ORTHANC_CORE_SOURCES_DEPENDENCIES}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Extent2D.cpp Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,107 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#include "Extent2D.h" + +#include <OrthancException.h> + + +void Extent2D::CheckNotEmpty() const +{ + if (isEmpty_) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } +} + + +Extent2D::Extent2D() : + isEmpty_(true), + x1_(0), + y1_(0), + x2_(0), + y2_(0) +{ +} + + +double Extent2D::GetMinX() const +{ + CheckNotEmpty(); + return x1_; +} + + +double Extent2D::GetMaxX() const +{ + CheckNotEmpty(); + return x2_; +} + + +double Extent2D::GetMinY() const +{ + CheckNotEmpty(); + return y1_; +} + + +double Extent2D::GetMaxY() const +{ + CheckNotEmpty(); + return y2_; +} + + +double Extent2D::GetWidth() const +{ + CheckNotEmpty(); + return x2_ - x1_; +} + + +double Extent2D::GetHeight() const +{ + CheckNotEmpty(); + return y2_ - y1_; +} + + +void Extent2D::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); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Extent2D.h Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,63 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#pragma once + +#include <boost/noncopyable.hpp> + + +class Extent2D : public boost::noncopyable +{ +private: + bool isEmpty_; + double x1_; + double y1_; + double x2_; + double y2_; + + void CheckNotEmpty() const; + +public: + Extent2D(); + + bool IsEmpty() const + { + return isEmpty_; + } + + double GetMinX() const; + + double GetMaxX() const; + + double GetMinY() const; + + double GetMaxY() const; + + double GetWidth() const; + + double GetHeight() const; + + void Add(double x, + double y); +};
--- a/Sources/Plugin.cpp Thu Apr 04 18:05:13 2024 +0200 +++ b/Sources/Plugin.cpp Thu Apr 04 18:35:54 2024 +0200 @@ -22,13 +22,16 @@ **/ +#include "VTKToolbox.h" +#include "Vector3D.h" +#include "Toolbox.h" +#include "Extent2D.h" + #include "../Resources/Orthanc/Plugins/OrthancPluginCppWrapper.h" #include <EmbeddedResources.h> #include <DicomFormat/DicomInstanceHasher.h> -#include <Compression/GzipCompressor.h> -#include <ChunkedBuffer.h> #include <DicomParsing/FromDcmtkBridge.h> #include <DicomParsing/ParsedDicomFile.h> #include <Images/ImageProcessing.h> @@ -37,20 +40,8 @@ #include <SerializationToolbox.h> #include <SystemToolbox.h> -#include <vtkImageConstantPad.h> -#include <vtkImageData.h> -#include <vtkImageResize.h> -#include <vtkMarchingCubes.h> -#include <vtkNew.h> -#include <vtkPolyData.h> -#include <vtkPolyDataNormals.h> -#include <vtkSmoothPolyDataFilter.h> -#include <vtkTriangle.h> - #include <boost/thread/shared_mutex.hpp> -#include <nifti1_io.h> - #define ORTHANC_PLUGIN_NAME "stl" @@ -162,92 +153,6 @@ #include <dcmtk/dcmdata/dcsequen.h> #include <dcmtk/dcmdata/dcuid.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 std::string GetStringValue(DcmItem& item, const DcmTagKey& key) @@ -292,128 +197,6 @@ } -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 ComputeSquaredNorm() const - { - return x_ * x_ + y_ * y_ + z_ * z_; - } - - double ComputeNorm() const - { - return sqrt(ComputeSquaredNorm()); - } - - 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 AreParallel(const Vector3D& a, - const Vector3D& b) - { - if (IsNear(a.ComputeSquaredNorm(), 1) && - IsNear(b.ComputeSquaredNorm(), 1)) - { - return IsNear(std::abs(Vector3D::DotProduct(a, b)), 1); - } - else - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange, - "Only applicable to normalized vectors"); - } - } -}; - - - -static bool MyParseDouble(double& value, - const std::string& s) -{ -#if 1 - // This version is much faster, as "ParseDouble()" internally uses - // "boost::lexical_cast<double>()" - 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: @@ -496,9 +279,9 @@ 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])) + if (!Toolbox::MyParseDouble(x, tokens[i]) || + !Toolbox::MyParseDouble(y, tokens[i + 1]) || + !Toolbox::MyParseDouble(z, tokens[i + 2])) { throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); } @@ -550,7 +333,7 @@ { normal = Vector3D::CrossProduct(Vector3D(points_[1], points_[0]), Vector3D(points_[2], points_[0])); - if (!IsNear(normal.ComputeNorm(), 0)) + if (!Toolbox::IsNear(normal.ComputeNorm(), 0)) { normal.Normalize(); hasNormal = true; @@ -567,7 +350,7 @@ for (size_t i = 1; i < points_.size(); i++) { double b = Vector3D::DotProduct(points_[i], normal); - if (!IsNear(a, b)) + if (!Toolbox::IsNear(a, b)) { return false; } @@ -580,8 +363,8 @@ const Vector3D& axisX, const Vector3D& axisY) const { - assert(IsNear(1, axisX.ComputeNorm())); - assert(IsNear(1, axisY.ComputeNorm())); + assert(Toolbox::IsNear(1, axisX.ComputeNorm())); + assert(Toolbox::IsNear(1, axisY.ComputeNorm())); for (size_t i = 0; i < points_.size(); i++) { @@ -593,23 +376,6 @@ -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: @@ -789,7 +555,7 @@ double d = (z - minProjectionAlongNormal_) / slicesSpacing_; - if (IsNear(d, round(d))) + if (Toolbox::IsNear(d, round(d))) { if (d < 0.0 || d > static_cast<double>(slicesCount_) - 1.0) @@ -857,7 +623,7 @@ // Only keep unique projections std::sort(projections.begin(), projections.end()); - RemoveDuplicateValues(projections); + Toolbox::RemoveDuplicateValues(projections); assert(!projections.empty()); if (projections.size() == 1) @@ -884,7 +650,7 @@ } std::sort(spacings.begin(), spacings.end()); - RemoveDuplicateValues(spacings); + Toolbox::RemoveDuplicateValues(spacings); if (spacings.empty()) { @@ -921,7 +687,7 @@ while (it != candidates.end()) { double d = (projections[*it] - projections[reference]) / slicesSpacing_; - if (IsNear(d, round(d))) + if (Toolbox::IsNear(d, round(d))) { countSupport ++; } @@ -960,7 +726,7 @@ for (size_t i = 0; i < projections.size(); i++) { double d = (projections[i] - bestProjection) / slicesSpacing_; - if (IsNear(d, round(d))) + if (Toolbox::IsNear(d, round(d))) { minProjectionAlongNormal_ = std::min(minProjectionAlongNormal_, projections[i]); maxProjectionAlongNormal_ = std::max(maxProjectionAlongNormal_, projections[i]); @@ -968,7 +734,7 @@ } double d = (maxProjectionAlongNormal_ - minProjectionAlongNormal_) / slicesSpacing_; - if (IsNear(d, round(d))) + if (Toolbox::IsNear(d, round(d))) { slicesCount_ = static_cast<size_t>(round(d)) + 1; } @@ -1077,171 +843,6 @@ }; -static void WriteFloat(Orthanc::ChunkedBuffer& buffer, - float value, - Orthanc::Endianness endianness) -{ - switch (endianness) - { - case Orthanc::Endianness_Little: - buffer.AddChunk(&value, sizeof(float)); - break; - - case Orthanc::Endianness_Big: - { - uint8_t tmp[4]; - tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; - tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; - tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; - tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; - buffer.AddChunk(&tmp, sizeof(float)); - break; - } - - default: - throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); - } -} - - -static void WriteInteger(Orthanc::ChunkedBuffer& buffer, - uint32_t value, - Orthanc::Endianness endianness) -{ - switch (endianness) - { - case Orthanc::Endianness_Little: - buffer.AddChunk(&value, sizeof(uint32_t)); - break; - - case Orthanc::Endianness_Big: - { - uint8_t tmp[4]; - tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; - tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; - tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; - tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; - buffer.AddChunk(&tmp, sizeof(uint32_t)); - break; - } - - default: - throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); - } -} - - -static void EncodeSTL(std::string& target /* out */, - vtkPolyData& mesh /* in */) -{ - const Orthanc::Endianness endianness = Orthanc::Toolbox::DetectEndianness(); - - Orthanc::ChunkedBuffer buffer; - - uint8_t header[80]; - memset(header, 0, sizeof(header)); - buffer.AddChunk(header, sizeof(header)); // This doesn't depend on endianness - - WriteInteger(buffer, mesh.GetNumberOfCells(), endianness); - - 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); - - WriteFloat(buffer, normal[0], endianness); - WriteFloat(buffer, normal[1], endianness); - WriteFloat(buffer, normal[2], endianness); - WriteFloat(buffer, p0[0], endianness); - WriteFloat(buffer, p0[1], endianness); - WriteFloat(buffer, p0[2], endianness); - WriteFloat(buffer, p1[0], endianness); - WriteFloat(buffer, p1[1], endianness); - WriteFloat(buffer, p1[2], endianness); - WriteFloat(buffer, p2[0], endianness); - WriteFloat(buffer, p2[1], endianness); - WriteFloat(buffer, p2[2], endianness); - - uint16_t a = 0; - buffer.AddChunk(&a, sizeof(a)); // This doesn't depend on endianness - } - - buffer.Flatten(target); -} - - -bool EncodeVolume(std::string& stl, - vtkImageData* volume, - unsigned int resolution, - bool smooth) -{ - if (volume == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); - } - - vtkNew<vtkImageResize> resize; - resize->SetOutputDimensions(resolution, resolution, resolution); - resize->SetInputData(volume); - resize->Update(); - - vtkNew<vtkImageConstantPad> padding; - padding->SetConstant(0); - padding->SetOutputNumberOfScalarComponents(1); - padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution); - padding->SetInputData(resize->GetOutput()); - padding->Update(); - - double range[2]; - padding->GetOutput()->GetScalarRange(range); - - const double isoValue = (range[0] + range[1]) / 2.0; - - vtkNew<vtkMarchingCubes> surface; - surface->SetInputData(padding->GetOutput()); - surface->ComputeNormalsOn(); - surface->SetValue(0, 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; -} - - static Orthanc::ParsedDicomFile* LoadInstance(const std::string& instanceId) { std::string dicom; @@ -1301,12 +902,12 @@ double x1, x2, x3, y1, y2, y3; if (items.size() == 6 && - MyParseDouble(x1, items[0]) && - MyParseDouble(x2, items[1]) && - MyParseDouble(x3, items[2]) && - MyParseDouble(y1, items[3]) && - MyParseDouble(y2, items[4]) && - MyParseDouble(y3, items[5])) + Toolbox::MyParseDouble(x1, items[0]) && + Toolbox::MyParseDouble(x2, items[1]) && + Toolbox::MyParseDouble(x3, items[2]) && + Toolbox::MyParseDouble(y1, items[3]) && + Toolbox::MyParseDouble(y2, items[4]) && + Toolbox::MyParseDouble(y3, items[5])) { axisX = Vector3D(x1, x2, x3); axisY = Vector3D(y1, y2, y3); @@ -1337,7 +938,7 @@ throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); } - if (!IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) + if (!Toolbox::IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) { throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); } @@ -1347,8 +948,8 @@ Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY); - if (!IsNear(1, axisX.ComputeNorm()) || - !IsNear(1, axisY.ComputeNorm()) || + if (!Toolbox::IsNear(1, axisX.ComputeNorm()) || + !Toolbox::IsNear(1, axisY.ComputeNorm()) || !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal())) { throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); @@ -1417,7 +1018,7 @@ volume->SetSpacing(pixelSpacingX, pixelSpacingY, pixelSpacingZ); - return EncodeVolume(stl, volume.Get(), resolution, smooth); + return VTKToolbox::EncodeVolume(stl, volume.Get(), resolution, smooth); } @@ -1674,114 +1275,6 @@ -class NiftiHeader : public boost::noncopyable -{ -private: - nifti_image* image_; - -public: - NiftiHeader(const std::string& nifti) - { - nifti_1_header header; - if (nifti.size() < sizeof(header)) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - memcpy(&header, nifti.c_str(), sizeof(header)); - if (!nifti_hdr_looks_good(&header)) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - - image_ = nifti_convert_nhdr2nim(header, "dummy_filename"); - if (image_ == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); - } - } - - ~NiftiHeader() - { - nifti_image_free(image_); - } - - const nifti_image& GetInfo() const - { - assert(image_ != NULL); - return *image_; - } -}; - - -static void LoadNifti(vtkImageData* volume, - std::string& nifti) -{ - if (volume == NULL) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); - } - - const uint8_t* p = reinterpret_cast<const uint8_t*>(nifti.c_str()); - - if (nifti.size() >= 2 && - p[0] == 0x1f && - p[1] == 0x8b) - { - Orthanc::GzipCompressor compressor; - std::string uncompressed; - Orthanc::IBufferCompressor::Uncompress(uncompressed, compressor, nifti); - nifti.swap(uncompressed); - } - - NiftiHeader header(nifti); - - if (header.GetInfo().ndim != 3) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, - "Only 3D NIfTI volumes are allowed"); - } - - size_t itemSize; - int vtkType; - - switch (header.GetInfo().datatype) - { - case DT_UNSIGNED_CHAR: - itemSize = 1; - vtkType = VTK_UNSIGNED_CHAR; - break; - - case DT_FLOAT: - itemSize = sizeof(float); - vtkType = VTK_FLOAT; - break; - - case DT_DOUBLE: - itemSize = sizeof(double); - vtkType = VTK_DOUBLE; - break; - - default: - throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); - } - - assert(static_cast<int>(header.GetInfo().nvox) == header.GetInfo().nx * header.GetInfo().ny * header.GetInfo().nz); - - const size_t pixelDataOffset = sizeof(nifti_1_header) + 4 /* extension */; - - if (nifti.size() != pixelDataOffset + header.GetInfo().nvox * itemSize) - { - throw Orthanc::OrthancException(Orthanc::ErrorCode_CorruptedFile); - } - - volume->SetDimensions(header.GetInfo().nx, header.GetInfo().ny, header.GetInfo().nz); - volume->AllocateScalars(vtkType, 1); - volume->SetSpacing(header.GetInfo().dx, header.GetInfo().dy, header.GetInfo().dz); - memcpy(volume->GetScalarPointer(), &nifti[pixelDataOffset], header.GetInfo().nvox * itemSize); -} - - void EncodeNifti(OrthancPluginRestOutput* output, const char* url, const OrthancPluginHttpRequest* request) @@ -1818,10 +1311,10 @@ 256 /* default value */); vtkNew<vtkImageData> volume; - LoadNifti(volume.Get(), nifti); + VTKToolbox::LoadNifti(volume.Get(), nifti); std::string stl; - if (!EncodeVolume(stl, volume.Get(), resolution, smooth)) + if (!VTKToolbox::EncodeVolume(stl, volume.Get(), resolution, smooth)) { throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, "Cannot encode STL from NIfTI"); }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Toolbox.cpp Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,72 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#include "Toolbox.h" + +#include <algorithm> +#include <cmath> + + +namespace Toolbox +{ + bool IsNear(double a, + double b) + { + return std::abs(a - b) < 10.0 * std::numeric_limits<double>::epsilon(); + } + + + 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 + } + + + namespace + { + struct IsNearPredicate + { + bool operator() (const double& a, + const double& b) + { + return Toolbox::IsNear(a, b); + } + }; + } + + + 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()); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Toolbox.h Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,42 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#pragma once + +#include <string> +#include <vector> + + +namespace Toolbox +{ + bool IsNear(double a, + double b); + + // This version is much faster, as "ParseDouble()" internally uses + // "boost::lexical_cast<double>()" + bool MyParseDouble(double& value, + const std::string& s); + + void RemoveDuplicateValues(std::vector<double>& v); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/VTKToolbox.cpp Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,319 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#include "VTKToolbox.h" + +#include <ChunkedBuffer.h> +#include <Compression/GzipCompressor.h> +#include <Endianness.h> +#include <OrthancException.h> +#include <Toolbox.h> + +#include <vtkTriangle.h> +#include <vtkImageResize.h> +#include <vtkImageConstantPad.h> +#include <vtkMarchingCubes.h> +#include <vtkSmoothPolyDataFilter.h> +#include <vtkPolyDataNormals.h> + +#include <nifti1_io.h> + + +namespace VTKToolbox +{ + static void WriteFloat(Orthanc::ChunkedBuffer& buffer, + float value, + Orthanc::Endianness endianness) + { + switch (endianness) + { + case Orthanc::Endianness_Little: + buffer.AddChunk(&value, sizeof(float)); + break; + + case Orthanc::Endianness_Big: + { + uint8_t tmp[4]; + tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; + tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; + tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; + tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; + buffer.AddChunk(&tmp, sizeof(float)); + break; + } + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + + + static void WriteInteger(Orthanc::ChunkedBuffer& buffer, + uint32_t value, + Orthanc::Endianness endianness) + { + switch (endianness) + { + case Orthanc::Endianness_Little: + buffer.AddChunk(&value, sizeof(uint32_t)); + break; + + case Orthanc::Endianness_Big: + { + uint8_t tmp[4]; + tmp[0] = reinterpret_cast<uint8_t*>(&value) [3]; + tmp[1] = reinterpret_cast<uint8_t*>(&value) [2]; + tmp[2] = reinterpret_cast<uint8_t*>(&value) [1]; + tmp[3] = reinterpret_cast<uint8_t*>(&value) [0]; + buffer.AddChunk(&tmp, sizeof(uint32_t)); + break; + } + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + + + void EncodeSTL(std::string& target /* out */, + vtkPolyData& mesh /* in */) + { + const Orthanc::Endianness endianness = Orthanc::Toolbox::DetectEndianness(); + + Orthanc::ChunkedBuffer buffer; + + uint8_t header[80]; + memset(header, 0, sizeof(header)); + buffer.AddChunk(header, sizeof(header)); // This doesn't depend on endianness + + WriteInteger(buffer, mesh.GetNumberOfCells(), endianness); + + 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); + + WriteFloat(buffer, normal[0], endianness); + WriteFloat(buffer, normal[1], endianness); + WriteFloat(buffer, normal[2], endianness); + WriteFloat(buffer, p0[0], endianness); + WriteFloat(buffer, p0[1], endianness); + WriteFloat(buffer, p0[2], endianness); + WriteFloat(buffer, p1[0], endianness); + WriteFloat(buffer, p1[1], endianness); + WriteFloat(buffer, p1[2], endianness); + WriteFloat(buffer, p2[0], endianness); + WriteFloat(buffer, p2[1], endianness); + WriteFloat(buffer, p2[2], endianness); + + uint16_t a = 0; + buffer.AddChunk(&a, sizeof(a)); // This doesn't depend on endianness + } + + buffer.Flatten(target); + } + + + bool EncodeVolume(std::string& stl, + vtkImageData* volume, + unsigned int resolution, + bool smooth) + { + if (volume == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); + } + + vtkNew<vtkImageResize> resize; + resize->SetOutputDimensions(resolution, resolution, resolution); + resize->SetInputData(volume); + resize->Update(); + + vtkNew<vtkImageConstantPad> padding; + padding->SetConstant(0); + padding->SetOutputNumberOfScalarComponents(1); + padding->SetOutputWholeExtent(-1, resolution, -1, resolution, -1, resolution); + padding->SetInputData(resize->GetOutput()); + padding->Update(); + + double range[2]; + padding->GetOutput()->GetScalarRange(range); + + const double isoValue = (range[0] + range[1]) / 2.0; + + vtkNew<vtkMarchingCubes> surface; + surface->SetInputData(padding->GetOutput()); + surface->ComputeNormalsOn(); + surface->SetValue(0, 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; + } + + + namespace + { + class NiftiHeader : public boost::noncopyable + { + private: + nifti_image* image_; + + public: + NiftiHeader(const std::string& nifti) + { + nifti_1_header header; + if (nifti.size() < sizeof(header)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + memcpy(&header, nifti.c_str(), sizeof(header)); + if (!nifti_hdr_looks_good(&header)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + image_ = nifti_convert_nhdr2nim(header, "dummy_filename"); + if (image_ == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + } + + ~NiftiHeader() + { + nifti_image_free(image_); + } + + const nifti_image& GetInfo() const + { + assert(image_ != NULL); + return *image_; + } + }; + } + + + void LoadNifti(vtkImageData* volume /* out */, + std::string& nifti /* in */) + { + if (volume == NULL) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_NullPointer); + } + + const uint8_t* p = reinterpret_cast<const uint8_t*>(nifti.c_str()); + + if (nifti.size() >= 2 && + p[0] == 0x1f && + p[1] == 0x8b) + { + Orthanc::GzipCompressor compressor; + std::string uncompressed; + Orthanc::IBufferCompressor::Uncompress(uncompressed, compressor, nifti); + nifti.swap(uncompressed); + } + + NiftiHeader header(nifti); + + if (header.GetInfo().ndim != 3) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat, + "Only 3D NIfTI volumes are allowed"); + } + + size_t itemSize; + int vtkType; + + switch (header.GetInfo().datatype) + { + case DT_UNSIGNED_CHAR: + itemSize = 1; + vtkType = VTK_UNSIGNED_CHAR; + break; + + case DT_FLOAT: + itemSize = sizeof(float); + vtkType = VTK_FLOAT; + break; + + case DT_DOUBLE: + itemSize = sizeof(double); + vtkType = VTK_DOUBLE; + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + + assert(static_cast<int>(header.GetInfo().nvox) == header.GetInfo().nx * header.GetInfo().ny * header.GetInfo().nz); + + const size_t pixelDataOffset = sizeof(nifti_1_header) + 4 /* extension */; + + if (nifti.size() != pixelDataOffset + header.GetInfo().nvox * itemSize) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_CorruptedFile); + } + + volume->SetDimensions(header.GetInfo().nx, header.GetInfo().ny, header.GetInfo().nz); + volume->AllocateScalars(vtkType, 1); + volume->SetSpacing(header.GetInfo().dx, header.GetInfo().dy, header.GetInfo().dz); + memcpy(volume->GetScalarPointer(), &nifti[pixelDataOffset], header.GetInfo().nvox * itemSize); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/VTKToolbox.h Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,43 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#pragma once + +#include <vtkPolyData.h> +#include <vtkImageData.h> + + +namespace VTKToolbox +{ + void EncodeSTL(std::string& target /* out */, + vtkPolyData& mesh /* in */); + + bool EncodeVolume(std::string& stl, + vtkImageData* volume, + unsigned int resolution, + bool smooth); + + void LoadNifti(vtkImageData* volume /* out */, + std::string& nifti /* in */); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Vector3D.cpp Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,114 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#include "Vector3D.h" + +#include "Toolbox.h" + +#include <OrthancException.h> + +#include <cmath> + + +Vector3D::Vector3D() : + x_(0), + y_(0), + z_(0) +{ +} + + +Vector3D::Vector3D(double x, + double y, + double z) : + x_(x), + y_(y), + z_(z) +{ +} + + +Vector3D::Vector3D(const Vector3D& from, + const Vector3D& to) : + x_(to.x_ - from.x_), + y_(to.y_ - from.y_), + z_(to.z_ - from.z_) +{ +} + + +double Vector3D::ComputeSquaredNorm() const +{ + return x_ * x_ + y_ * y_ + z_ * z_; +} + + +double Vector3D::ComputeNorm() const +{ + return sqrt(ComputeSquaredNorm()); +} + + +void Vector3D::Normalize() +{ + double norm = ComputeNorm(); + if (!Toolbox::IsNear(norm, 0)) + { + x_ /= norm; + y_ /= norm; + z_ /= norm; + } +} + + +Vector3D 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()); +} + + +double Vector3D::DotProduct(const Vector3D& a, + const Vector3D& b) +{ + return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ(); +} + + +bool Vector3D::AreParallel(const Vector3D& a, + const Vector3D& b) +{ + if (Toolbox::IsNear(a.ComputeSquaredNorm(), 1) && + Toolbox::IsNear(b.ComputeSquaredNorm(), 1)) + { + return Toolbox::IsNear(std::abs(Vector3D::DotProduct(a, b)), 1); + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange, + "Only applicable to normalized vectors"); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Sources/Vector3D.h Thu Apr 04 18:35:54 2024 +0200 @@ -0,0 +1,74 @@ +/** + * SPDX-FileCopyrightText: 2023-2024 Sebastien Jodogne, UCLouvain, Belgium + * SPDX-License-Identifier: GPL-3.0-or-later + */ + +/** + * STL plugin for Orthanc + * Copyright (C) 2023-2024 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/>. + **/ + + +#pragma once + + +class Vector3D +{ +private: + double x_; + double y_; + double z_; + +public: + Vector3D(); + + Vector3D(double x, + double y, + double z); + + Vector3D(const Vector3D& from, + const Vector3D& to); + + double GetX() const + { + return x_; + } + + double GetY() const + { + return y_; + } + + double GetZ() const + { + return z_; + } + + double ComputeSquaredNorm() const; + + double ComputeNorm() const; + + void Normalize(); + + static Vector3D CrossProduct(const Vector3D& u, + const Vector3D& v); + + static double DotProduct(const Vector3D& a, + const Vector3D& b); + + static bool AreParallel(const Vector3D& a, + const Vector3D& b); +};