--- 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<OrthancStone::IViewport> lock1(viewport.lock());
+      
+      if (lock1 &&
+          IsAction(event, WebViewerAction_Crosshair))
       {
-        printf("CROSS-HAIR!\n");
+        OrthancStone::CoordinateSystem3D plane;
+        if (viewer_.GetCurrentPlane(plane))
+        {
+          std::unique_ptr<OrthancStone::IViewport::ILock> 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);
+    }
+  }
 };
 
 

--- 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 <Toolbox.h>
 #include <OrthancException.h>
 
+#include <boost/math/constants/constants.hpp>
+
 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);
+  }
 }

--- 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);
   };
 }

--- 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");
+    }
+  }
 }

--- 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;
   };
 }

--- 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)));
+  }
 }