# HG changeset patch # User Sebastien Jodogne # Date 1605205628 -3600 # Node ID 4a43106bc122a53c764f7a5796b50b79cdf23bc6 # Parent adff3cd78775d4639aff054904983a68916ed4ba cross-hair starts to work diff -r adff3cd78775 -r 4a43106bc122 Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp --- a/Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp Thu Nov 12 16:57:15 2020 +0100 +++ b/Applications/StoneWebViewer/WebAssembly/StoneWebViewer.cpp Thu Nov 12 19:27:08 2020 +0100 @@ -901,6 +901,8 @@ size_t currentFrame, size_t countFrames, DisplayedFrameQuality quality) = 0; + + virtual void SignalCrosshair(const OrthancStone::Vector& click) = 0; }; private: @@ -1963,9 +1965,10 @@ class Interactor : public OrthancStone::DefaultViewportInteractor { private: - WebViewerAction leftAction_; - WebViewerAction middleAction_; - WebViewerAction rightAction_; + ViewerViewport& viewer_; + WebViewerAction leftAction_; + WebViewerAction middleAction_; + WebViewerAction rightAction_; bool IsAction(const OrthancStone::PointerEvent& event, WebViewerAction action) @@ -1987,9 +1990,11 @@ } public: - Interactor(WebViewerAction leftAction, + Interactor(ViewerViewport& viewer, + WebViewerAction leftAction, WebViewerAction middleAction, WebViewerAction rightAction) : + viewer_(viewer), leftAction_(leftAction), middleAction_(middleAction), rightAction_(rightAction) @@ -2005,9 +2010,28 @@ unsigned int viewportWidth, unsigned int viewportHeight) ORTHANC_OVERRIDE { - if (IsAction(event, WebViewerAction_Crosshair)) + boost::shared_ptr lock1(viewport.lock()); + + if (lock1 && + IsAction(event, WebViewerAction_Crosshair)) { - printf("CROSS-HAIR!\n"); + OrthancStone::CoordinateSystem3D plane; + if (viewer_.GetCurrentPlane(plane)) + { + std::unique_ptr lock2(lock1->Lock()); + + const OrthancStone::ScenePoint2D p = event.GetMainPosition(); + double x = p.GetX(); + double y = p.GetY(); + lock2->GetController().GetCanvasToSceneTransform().Apply(x, y); + + OrthancStone::Vector click = plane.MapSliceToWorldCoordinates(x, y); + if (viewer_.observer_.get() != NULL) + { + viewer_.observer_->SignalCrosshair(click); + } + } + return NULL; // No need for a tracker, this is just a click } else @@ -2024,7 +2048,7 @@ WebViewerAction rightAction) { assert(viewport_ != NULL); - viewport_->AcquireInteractor(new Interactor(leftAction, middleAction, rightAction)); + viewport_->AcquireInteractor(new Interactor(*this, leftAction, middleAction, rightAction)); } void FitForPrint() @@ -2067,6 +2091,25 @@ hasFocusOnInstance_ = false; } } + + void FocusOnPoint(const OrthancStone::Vector& p) + { + static const double MAX_DISTANCE = 0.5; // 0.5 cm => TODO parameter? + + OrthancStone::LinearAlgebra::Print(p); + size_t frameIndex; + if (cursor_.get() != NULL && + frames_.get() != NULL && + frames_->FindClosestFrame(frameIndex, p, MAX_DISTANCE)) + { + cursor_->SetCurrentIndex(frameIndex); + DisplayCurrentFrame(); + } + else + { + printf("nope\n"); + } + } }; @@ -2172,7 +2215,16 @@ UpdateReferenceLines(); - }; + } + + virtual void SignalCrosshair(const OrthancStone::Vector& click) ORTHANC_OVERRIDE + { + for (Viewports::const_iterator it = allViewports_.begin(); it != allViewports_.end(); ++it) + { + assert(it->second != NULL); + it->second->FocusOnPoint(click); + } + } }; diff -r adff3cd78775 -r 4a43106bc122 OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp --- a/OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp Thu Nov 12 16:57:15 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/CoordinateSystem3D.cpp Thu Nov 12 19:27:08 2020 +0100 @@ -29,6 +29,8 @@ #include #include +#include + namespace OrthancStone { void CoordinateSystem3D::CheckAndComputeNormal() @@ -61,7 +63,7 @@ d_ = -(normal_[0] * origin_[0] + normal_[1] * origin_[1] + normal_[2] * origin_[2]); - // Just a sanity check, it should be useless by construction + // Just a sanity check, it should be useless by construction (*) assert(LinearAlgebra::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0)); } } @@ -229,6 +231,18 @@ } + 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) @@ -265,4 +279,81 @@ 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); + } } diff -r adff3cd78775 -r 4a43106bc122 OrthancStone/Sources/Toolbox/CoordinateSystem3D.h --- a/OrthancStone/Sources/Toolbox/CoordinateSystem3D.h Thu Nov 12 16:57:15 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/CoordinateSystem3D.h Thu Nov 12 19:27:08 2020 +0100 @@ -135,6 +135,9 @@ const Vector& origin, const Vector& direction) const; + // Point-to-plane distance + double ComputeDistance(const Vector& p) const; + // Returns "false" is the two planes are not parallel static bool ComputeDistance(double& distance, const CoordinateSystem3D& a, @@ -143,5 +146,14 @@ // Normalize a cutting plane so that the origin (0,0,0) of the 3D // world is mapped to the origin of its (x,y) coordinate system static CoordinateSystem3D NormalizeCuttingPlane(const CoordinateSystem3D& plane); + + // Construct one possible coordinate system from the general form + // of the equation of a plane "a*x+b*y+c*z+d=0". Note that the + // axes are not determined in this case, and so they are chosen + // arbitrarily. + static CoordinateSystem3D CreateFromPlaneGeneralForm(double a, + double b, + double c, + double d); }; } diff -r adff3cd78775 -r 4a43106bc122 OrthancStone/Sources/Toolbox/SortedFrames.cpp --- a/OrthancStone/Sources/Toolbox/SortedFrames.cpp Thu Nov 12 16:57:15 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/SortedFrames.cpp Thu Nov 12 19:27:08 2020 +0100 @@ -42,6 +42,13 @@ } + double SortedFrames::Frame::ComputeDistance(const Vector& p) const + { + const CoordinateSystem3D& plane = instance_->GetFrameGeometry(frameNumber_); + return plane.ComputeDistance(p); + } + + const DicomInstanceParameters& SortedFrames::GetInstance(size_t instanceIndex) const { if (instanceIndex >= instances_.size()) @@ -399,4 +406,42 @@ sorted_ = true; } } + + + bool SortedFrames::FindClosestFrame(size_t& frameIndex, + const Vector& point, + double maximumDistance) const + { + if (sorted_) + { + if (frames_.empty()) + { + return false; + } + else + { + frameIndex = 0; + double closestDistance = frames_[0].ComputeDistance(point); + + for (size_t i = 1; i < frames_.size(); i++) + { + double d = frames_[i].ComputeDistance(point); + printf("%f ", d); + if (d < closestDistance) + { + frameIndex = i; + closestDistance = d; + } + } + + printf("\n>> %f\n", closestDistance); + return (closestDistance <= maximumDistance); + } + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls, + "Sort() has not been called"); + } + } } diff -r adff3cd78775 -r 4a43106bc122 OrthancStone/Sources/Toolbox/SortedFrames.h --- a/OrthancStone/Sources/Toolbox/SortedFrames.h Thu Nov 12 16:57:15 2020 +0100 +++ b/OrthancStone/Sources/Toolbox/SortedFrames.h Thu Nov 12 19:27:08 2020 +0100 @@ -58,6 +58,8 @@ { return frameNumber_; } + + double ComputeDistance(const Vector& p) const; }; @@ -147,5 +149,9 @@ unsigned int frameNumber) const; void Sort(); + + bool FindClosestFrame(size_t& frameIndex, + const Vector& point, + double maximumDistance) const; }; } diff -r adff3cd78775 -r 4a43106bc122 UnitTestsSources/UnitTestsMain.cpp --- a/UnitTestsSources/UnitTestsMain.cpp Thu Nov 12 16:57:15 2020 +0100 +++ b/UnitTestsSources/UnitTestsMain.cpp Thu Nov 12 19:27:08 2020 +0100 @@ -887,6 +887,11 @@ ASSERT_FLOAT_EQ(c.GetNormal()[0], 0); ASSERT_FLOAT_EQ(c.GetNormal()[1], 0); ASSERT_FLOAT_EQ(c.GetNormal()[2], 1); + + ASSERT_FLOAT_EQ(0, c.ComputeDistance(OrthancStone::LinearAlgebra::CreateVector(0, 0, 0))); + ASSERT_FLOAT_EQ(0, c.ComputeDistance(OrthancStone::LinearAlgebra::CreateVector(5, 0, 0))); + ASSERT_FLOAT_EQ(0, c.ComputeDistance(OrthancStone::LinearAlgebra::CreateVector(0, 5, 0))); + ASSERT_FLOAT_EQ(5, c.ComputeDistance(OrthancStone::LinearAlgebra::CreateVector(0, 0, 5))); } { @@ -896,6 +901,13 @@ ASSERT_FLOAT_EQ(c.GetNormal()[1], 0); ASSERT_FLOAT_EQ(c.GetNormal()[2], 1); } + + { + // https://www.vedantu.com/maths/perpendicular-distance-of-a-point-from-a-plane + OrthancStone::CoordinateSystem3D c = + OrthancStone::CoordinateSystem3D::CreateFromPlaneGeneralForm(2, 4, -4, -6); + ASSERT_FLOAT_EQ(3, c.ComputeDistance(OrthancStone::LinearAlgebra::CreateVector(0, 3, 6))); + } }