--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Framework/Toolbox/SliceGeometry.cpp	Fri Oct 14 15:34:11 2016 +0200
@@ -0,0 +1,165 @@
+/**
+ * Stone of Orthanc
+ * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
+ * Department, University Hospital of Liege, 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.
+ *
+ * In addition, as a special exception, the copyright holders of this
+ * program give permission to link the code of its release with the
+ * OpenSSL project's "OpenSSL" library (or with modified versions of it
+ * that use the same license as the "OpenSSL" library), and distribute
+ * the linked executables. You must obey the GNU General Public License
+ * in all respects for all of the code used other than "OpenSSL". If you
+ * modify file(s) with this exception, you may extend this exception to
+ * your version of the file(s), but you are not obligated to do so. If
+ * you do not wish to do so, delete this exception statement from your
+ * version. If you delete this exception statement from all source files
+ * in the program, then also delete it here.
+ * 
+ * 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 "SliceGeometry.h"
+
+#include "GeometryToolbox.h"
+
+#include "../Orthanc/Core/Logging.h"
+#include "../Orthanc/Core/Toolbox.h"
+#include "../Orthanc/Core/OrthancException.h"
+
+namespace OrthancStone
+{
+  void SliceGeometry::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 (!GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisX_), 1.0) ||
+        !GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(axisY_), 1.0))
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+    }
+
+    // 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
+    if (!GeometryToolbox::IsCloseToZero(boost::numeric::ublas::inner_prod(axisX_, axisY_)))
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+    }
+
+    GeometryToolbox::CrossProduct(normal_, axisX_, axisY_);
+
+    // Just a sanity check, it should be useless by construction
+    assert(GeometryToolbox::IsNear(boost::numeric::ublas::norm_2(normal_), 1.0));
+  }
+
+
+  void SliceGeometry::SetupCanonical()
+  {
+    GeometryToolbox::AssignVector(origin_, 0, 0, 0);
+    GeometryToolbox::AssignVector(axisX_, 1, 0, 0);
+    GeometryToolbox::AssignVector(axisY_, 0, 1, 0);
+    CheckAndComputeNormal();
+  }
+
+
+  SliceGeometry::SliceGeometry(const Vector& origin,
+                               const Vector& axisX,
+                               const Vector& axisY) :
+    origin_(origin),
+    axisX_(axisX),
+    axisY_(axisY)
+  {
+    CheckAndComputeNormal();
+  }
+
+
+  void SliceGeometry::Setup(const std::string& imagePositionPatient,
+                            const std::string& imageOrientationPatient)
+  {
+    std::string tmpPosition = Orthanc::Toolbox::StripSpaces(imagePositionPatient);
+    std::string tmpOrientation = Orthanc::Toolbox::StripSpaces(imageOrientationPatient);
+
+    Vector orientation;
+    if (!GeometryToolbox::ParseVector(origin_, tmpPosition) ||
+        !GeometryToolbox::ParseVector(orientation, tmpOrientation) ||
+        origin_.size() != 3 ||
+        orientation.size() != 6)
+    {
+      throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat);
+    }
+
+    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();
+  }   
+
+
+  SliceGeometry::SliceGeometry(const DicomDataset& dicom)
+  {
+    if (dicom.HasTag(DICOM_TAG_IMAGE_POSITION_PATIENT) &&
+        dicom.HasTag(DICOM_TAG_IMAGE_ORIENTATION_PATIENT))
+    {
+      Setup(dicom.GetStringValue(DICOM_TAG_IMAGE_POSITION_PATIENT),
+            dicom.GetStringValue(DICOM_TAG_IMAGE_ORIENTATION_PATIENT));
+    }
+    else
+    {
+      SetupCanonical();
+    }
+  }   
+
+
+  Vector SliceGeometry::MapSliceToWorldCoordinates(double x,
+                                                   double y) const
+  {
+    return origin_ + x * axisX_ + y * axisY_;
+  }
+
+
+  double SliceGeometry::ProjectAlongNormal(const Vector& point) const
+  {
+    return boost::numeric::ublas::inner_prod(point, normal_);
+  }
+
+
+  void SliceGeometry::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_);
+  }
+}