Mercurial > hg > orthanc-stone
view OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp @ 2166:f14260c15151
fix
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Mon, 07 Oct 2024 12:35:07 +0200 |
parents | 16c01cc201e7 |
children |
line wrap: on
line source
/** * Stone of Orthanc * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics * Department, University Hospital of Liege, Belgium * Copyright (C) 2017-2023 Osimis S.A., Belgium * Copyright (C) 2021-2024 Sebastien Jodogne, ICTEAM UCLouvain, Belgium * * This program is free software: you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this program. If not, see * <http://www.gnu.org/licenses/>. **/ #include "CoordinateSystem3D.h" #include "LinearAlgebra.h" #include "GeometryToolbox.h" #include <Logging.h> #include <Toolbox.h> #include <OrthancException.h> namespace OrthancStone { void CoordinateSystem3D::CheckAndComputeNormal() { /** * DICOM expects normal vectors to define the axes: "The row and * column direction cosine vectors shall be normal, i.e., the dot * product of each direction cosine vector with itself shall be * unity." * http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html **/ if (!LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) || !LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0)) { LOG(WARNING) << "Invalid 3D geometry: Axes are not normal vectors"; SetupCanonical(); } /** * The vectors within "Image Orientation Patient" must be * orthogonal, according to the DICOM specification: "The row and * column direction cosine vectors shall be orthogonal, i.e., * their dot product shall be zero." * http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html * * The "0.00001" threshold is needed for KNIX (on this sample * image, the inner product equals "0.000003", which is rejected * by "LinearAlgebra::IsCloseToZero()"). **/ else if (!LinearAlgebra::IsNear(0, boost::numeric::ublas::inner_prod(axisX_, axisY_), 0.00001)) { LOG(WARNING) << "Invalid 3D geometry: Image orientation patient is not orthogonal"; SetupCanonical(); } else { LinearAlgebra::CrossProduct(normal_, axisX_, axisY_); d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]); // Just a sanity check, it should be useless by construction (*) assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0)); } } void CoordinateSystem3D::SetupCanonical() { valid_ = false; LinearAlgebra::AssignVector(origin_, 0, 0, 0); LinearAlgebra::AssignVector(axisX_, 1, 0, 0); LinearAlgebra::AssignVector(axisY_, 0, 1, 0); LinearAlgebra::AssignVector(normal_, 0, 0, 1); d_ = 0; } CoordinateSystem3D::CoordinateSystem3D(const Vector& origin, const Vector& axisX, const Vector& axisY) : valid_(true), origin_(origin), axisX_(axisX), axisY_(axisY) { CheckAndComputeNormal(); } void CoordinateSystem3D::Setup(const std::string& imagePositionPatient, const std::string& imageOrientationPatient) { valid_ = true; std::string tmpPosition = Orthanc::Toolbox::StripSpaces(imagePositionPatient); std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient); Vector orientation; if (!LinearAlgebra::ParseVector(origin_, tmpPosition) || !LinearAlgebra::ParseVector(orientation, tmpOrientation) || origin_.size() != 3 || orientation.size() != 6) { LOG(WARNING) << "Bad 3D geometry: image position/orientation patient: \"" << tmpPosition << "\" / \"" << tmpOrientation << "\""; SetupCanonical(); } else { axisX_.resize(3); axisX_[0] = orientation[0]; axisX_[1] = orientation[1]; axisX_[2] = orientation[2]; axisY_.resize(3); axisY_[0] = orientation[3]; axisY_[1] = orientation[4]; axisY_[2] = orientation[5]; CheckAndComputeNormal(); } } CoordinateSystem3D::CoordinateSystem3D(const IDicomDataset& dicom) { std::string a, b; if (dicom.GetStringValue(a, Orthanc::DicomPath(Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT)) && dicom.GetStringValue(b, Orthanc::DicomPath(Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT))) { Setup(a, b); } else { SetupCanonical(); } } CoordinateSystem3D::CoordinateSystem3D(const Orthanc::DicomMap& dicom) { std::string a, b; if (dicom.LookupStringValue(a, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, false) && dicom.LookupStringValue(b, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, false)) { Setup(a, b); } else { SetupCanonical(); } } void CoordinateSystem3D::SetOrigin(const Vector& origin) { if (origin.size() != 3) { throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); } else { origin_ = origin; } } Vector CoordinateSystem3D::MapSliceToWorldCoordinates(double x, double y) const { return origin_ + x * axisX_ + y * axisY_; } double CoordinateSystem3D::ProjectAlongNormal(const Vector& point) const { return boost::numeric::ublas::inner_prod(point, normal_); } void CoordinateSystem3D::ProjectPoint2(double& offsetX, double& offsetY, const Vector& point) const { // Project the point onto the slice double projectionX,projectionY,projectionZ; GeometryToolbox::ProjectPointOntoPlane2(projectionX, projectionY, projectionZ, point, normal_, origin_); // As the axes are orthonormal vectors thanks to // CheckAndComputeNormal(), the following dot products give the // offset of the origin of the slice wrt. the origin of the // reference plane https://en.wikipedia.org/wiki/Vector_projection offsetX = axisX_[0] * (projectionX - origin_[0]) + axisX_[1] * (projectionY - origin_[1]) + axisX_[2] * (projectionZ - origin_[2]); offsetY = axisY_[0] * (projectionX - origin_[0]) + axisY_[1] * (projectionY - origin_[1]) + axisY_[2] * (projectionZ - origin_[2]); } void CoordinateSystem3D::ProjectPoint(double& offsetX, double& offsetY, const Vector& point) const { // Project the point onto the slice Vector projection; GeometryToolbox::ProjectPointOntoPlane(projection, point, normal_, origin_); // As the axes are orthonormal vectors thanks to // CheckAndComputeNormal(), the following dot products give the // offset of the origin of the slice wrt. the origin of the // reference plane https://en.wikipedia.org/wiki/Vector_projection offsetX = boost::numeric::ublas::inner_prod(axisX_, projection - origin_); offsetY = boost::numeric::ublas::inner_prod(axisY_, projection - origin_); } bool CoordinateSystem3D::IntersectSegment(Vector& p, const Vector& edgeFrom, const Vector& edgeTo) const { return GeometryToolbox::IntersectPlaneAndSegment(p, normal_, d_, edgeFrom, edgeTo); } bool CoordinateSystem3D::IntersectLine(Vector& p, const Vector& origin, const Vector& direction) const { return GeometryToolbox::IntersectPlaneAndLine(p, normal_, d_, origin, direction); } double CoordinateSystem3D::ComputeDistance(const Vector& p) const { /** * "normal_" is an unit vector (*) => sqrt(a_1^2+a_2^2+a_3^2) = 1, * and the denominator equals 1 by construction. * https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane#Closest_point_and_distance_for_a_hyperplane_and_arbitrary_point **/ return std::abs(boost::numeric::ublas::inner_prod(p, normal_) + d_); } bool CoordinateSystem3D::ComputeDistance(double& distance, const CoordinateSystem3D& a, const CoordinateSystem3D& b) { bool opposite = false; // Ignored if (GeometryToolbox::IsParallelOrOpposite(opposite, a.GetNormal(), b.GetNormal())) { distance = std::abs(a.ProjectAlongNormal(a.GetOrigin()) - a.ProjectAlongNormal(b.GetOrigin())); return true; } else { return false; } } std::ostream& operator<< (std::ostream& s, const CoordinateSystem3D& that) { s << "origin: " << that.origin_ << " normal: " << that.normal_ << " axisX: " << that.axisX_ << " axisY: " << that.axisY_ << " D: " << that.d_; return s; } CoordinateSystem3D CoordinateSystem3D::NormalizeCuttingPlane(const CoordinateSystem3D& plane) { double ox, oy; plane.ProjectPoint(ox, oy, LinearAlgebra::CreateVector(0, 0, 0)); CoordinateSystem3D normalized(plane); normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy)); return normalized; } CoordinateSystem3D CoordinateSystem3D::CreateFromPlaneGeneralForm(double a, double b, double c, double d) { /** * "a*x + b*y + c*z + d = 0" => The un-normalized normal is vector * (a,b,c). **/ Vector normal; LinearAlgebra::AssignVector(normal, a, b, c); /** * Choose the origin of plane, as the point that is closest to the * origin of the axes (0,0,0). * https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane#Restatement_using_linear_algebra **/ double squaredNorm = a * a + b * b + c * c; if (LinearAlgebra::IsCloseToZero(squaredNorm)) { throw Orthanc::OrthancException(Orthanc::ErrorCode_BadGeometry, "Singular matrix"); } Vector origin = -d * normal / squaredNorm; /** * Select the X axis by computing a vector that is perpendicular * to the normal. * * "Exactly 1 and only 1 of the bools get set; b0/b1/b2 gets set * if dimension "i" has magnitude strictly less than all * subsequent dimensions and not greater than all previous * dimensions. We then have a unit vector with a single non-zero * dimension that corresponds to a dimension of minimum magnitude * in "normal". The cross product of this with "normal" is * orthogonal to "normal" by definition of cross product. Consider * now that the cross product is numerically unstable only when * the two vectors are very closely aligned. Consider that our * unit vector is large in only a single dimension and that that * dimension corresponds to the dimension where "normal" was * small. It's thus guaranteed to be loosely orthogonal to * "normal" before taking the cross product, with least * orthogonality in the case where all dimensions of "normal" are * equal. In this least-orthogonal case, we're still quite * orthogonal given that our unit vector has all but one dimension * 0 whereas "normal" has all equal. We thus avoid the unstable * case of taking the cross product of two nearly-aligned * vectors." https://stackoverflow.com/a/43454629/881731 **/ bool b0 = (normal[0] < normal[1]) && (normal[0] < normal[2]); bool b1 = (normal[1] <= normal[0]) && (normal[1] < normal[2]); bool b2 = (normal[2] <= normal[0]) && (normal[2] <= normal[1]); Vector swap = LinearAlgebra::CreateVector(b0 ? 1 : 0, b1 ? 1 : 0, b2 ? 1 : 0); Vector axisX; LinearAlgebra::CrossProduct(axisX, normal, swap); LinearAlgebra::NormalizeVector(axisX); /** * The Y axis follows as the cross-product of the normal and the X * axis. **/ Vector axisY; LinearAlgebra::CrossProduct(axisY, axisX, normal); LinearAlgebra::NormalizeVector(axisY); return CoordinateSystem3D(origin, axisX, axisY); } CoordinateSystem3D CoordinateSystem3D::CreateFromThreePoints(const Vector& a, const Vector& b, const Vector& c) { Vector axisX = b - a; LinearAlgebra::NormalizeVector(axisX); Vector normal; LinearAlgebra::CrossProduct(normal, axisX, c - a); Vector axisY; LinearAlgebra::CrossProduct(axisY, axisX, normal); LinearAlgebra::NormalizeVector(axisY); return CoordinateSystem3D(a, axisX, axisY); } static std::string GetOrientationString(const Vector& v) { /** * This function directly comes from David Clunie: * https://sites.google.com/site/dicomnotes/ * * It is also a C++ reimplementation of function * "getOrientationString()" from Cornerstone: * https://bitbucket.org/osimis/osimis-webviewer-plugin/src/master/frontend/local_dependencies/cornerstoneTools/cornerstoneTools.js **/ const char orientationX = v[0] < 0 ? 'R' : 'L'; const char orientationY = v[1] < 0 ? 'A' : 'P'; const char orientationZ = v[2] < 0 ? 'F' : 'H'; double absX = std::abs(v[0]); double absY = std::abs(v[1]); double absZ = std::abs(v[2]); std::string result; static const double THRESHOLD = 0.0001; for (unsigned int i = 0; i < 3; i++) { if (absX > THRESHOLD && absX > absY && absX > absZ) { result.push_back(orientationX); absX = 0; } else if (absY > THRESHOLD && absY > absX && absY > absZ) { result.push_back(orientationY); absY = 0; } else if (absZ > THRESHOLD && absZ > absX && absZ > absY) { result.push_back(orientationZ); absZ = 0; } else { break; } } return result; } static std::string InvertOrientationString(const std::string& source) { std::string target; target.resize(source.length()); for (unsigned int i = 0; i < source.length(); i++) { switch (source[i]) { case 'H': target[i] = 'F'; break; case 'F': target[i] = 'H'; break; case 'R': target[i] = 'L'; break; case 'L': target[i] = 'R'; break; case 'A': target[i] = 'P'; break; case 'P': target[i] = 'A'; break; default: throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); } } return target; } void CoordinateSystem3D::GetOrientationMarkers(std::string& top /* out */, std::string& bottom /* out */, std::string& left /* out */, std::string& right /* out */) const { if (IsValid()) { /** * This is a C++ reimplementation of function * "getOrientationMarkers()" from Cornerstone: * https://bitbucket.org/osimis/osimis-webviewer-plugin/src/master/frontend/local_dependencies/cornerstoneTools/cornerstoneTools.js * * ImageOrientationPatient is row cosines then column cosines * (i.e. horizontal then vertical). * https://dicom.nema.org/medical/dicom/current/output/html/part03.html#sect_C.7.6.2 **/ const Vector& rowCosines = axisX_; // horizontal const Vector& columnCosines = axisY_; // vertical bottom = GetOrientationString(columnCosines); right = GetOrientationString(rowCosines); top = InvertOrientationString(bottom); left = InvertOrientationString(right); } else { throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); } } }