Refs #3930: Method on MDHistoWorkspace to get the signal at a given coordinate

Code/Mantid/Framework/API/inc/MantidAPI/IMDWorkspace.h

@@ -67,6 +67,21 @@ namespace Mantid
       /// Creates a new iterator pointing to the first cell in the workspace
       virtual IMDIterator* createIterator() const;
 
+      /// Returns the (normalized) signal at a given coordinates
+      virtual signal_t getSignalAtCoord(const coord_t * coords) const
+      {
+        UNUSED_ARG(coords);
+        throw std::runtime_error("getSignalAtCoord() not implemented.");
+      }
+
+      //-------------------------------------------------------------------------------------------
+      /// Returns the (normalized) signal at a given coordinates
+      /// @param coords :: coordinate as a VMD vector
+      signal_t getSignalAtCoord(const Mantid::Kernel::VMD & coords) const
+      {
+        return this->getSignalAtCoord(coords.getBareArray());
+      }
+
       virtual ~IMDWorkspace();
 
     };

Code/Mantid/Framework/Kernel/inc/MantidKernel/VMD.h

@@ -256,6 +256,13 @@ namespace Kernel
     double& operator[](const size_t index)
     { return data[index]; }
 
+    //-------------------------------------------------------------------------------------------
+    /** @return the bare data array directly. */
+    const double * getBareArray() const
+    {
+      return data;
+    }
+
     //-------------------------------------------------------------------------------------------
     /** Return a simple string representation of the vector
      * @param separator :: string to place between values, one space is the default

Code/Mantid/Framework/Kernel/test/VMDTest.h

@@ -65,6 +65,16 @@ public:
     TS_ASSERT_DELTA( a[3], 4.0, 1e-5);
   }
 
+  void test_getBareArray()
+  {
+    VMD b(1,2,3,4);
+    const double * a = b.getBareArray();
+    TS_ASSERT_DELTA( a[0], 1.0, 1e-5);
+    TS_ASSERT_DELTA( a[1], 2.0, 1e-5);
+    TS_ASSERT_DELTA( a[2], 3.0, 1e-5);
+    TS_ASSERT_DELTA( a[3], 4.0, 1e-5);
+  }
+
   void test_operators_throw_ifNonMatchingDimensions()
   {
     VMD a(1,2,3);

Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h

@@ -92,6 +92,8 @@ namespace MDEvents
 
     Mantid::Kernel::VMD getCenter(size_t linearIndex) const;
 
+    /// Returns the (normalized) signal at a given coordinates
+    signal_t getSignalAtCoord(const coord_t * coords) const;
 
     /// Sets the signal at the specified index.
     void setSignalAt(size_t index, signal_t value)
@@ -245,6 +247,9 @@ namespace MDEvents
     /// Vector of the length of the box in each dimension
     coord_t * m_boxLength;
 
+    /// Vector of the origin in each dimension
+    coord_t * m_origin;
+
     /// For converting to/from linear index to tdimensions
     size_t * m_indexMaker;
     /// Max index into each dimension

Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspaceIterator.h

@@ -12,7 +12,9 @@ namespace MDEvents
 {
 
   /** An implementation of IMDIterator that iterates through
-    a MDHistoWorkspace.
+    a MDHistoWorkspace. It treats the bin in the workspace as
+    a box containing a single "event", at the center of each bin
+    and with the proper signal/error.
     
     @author Janik Zikovsky
     @date 2011-10-06

Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp

@@ -7,8 +7,7 @@
 #include "MantidMDEvents/MDHistoWorkspaceIterator.h"
 
 using namespace Mantid::Kernel;
-using Mantid::Geometry::IMDDimension_sptr;
-using Mantid::Geometry::IMDDimension;
+using namespace Mantid::Geometry;
 
 namespace Mantid
 {
@@ -55,6 +54,7 @@ namespace MDEvents
     delete [] m_boxLength;
     delete [] m_indexMaker;
     delete [] m_indexMax;
+    delete [] m_origin;
   }
 
   //----------------------------------------------------------------------------------------------
@@ -191,10 +191,14 @@ namespace MDEvents
       } // (for each dimension)
     }
 
-    // Now set the m_boxLength
+    // Now set the m_boxLength and origin
     m_boxLength = new coord_t[nd];
+    m_origin = new coord_t[nd];
     for (size_t d=0; d<nd; d++)
+    {
       m_boxLength[d] = m_dimensions[d]->getX(1) - m_dimensions[d]->getX(0);
+      m_origin[d] = m_dimensions[d]->getX(0);
+    }
 
     // The index calculator
     m_indexMax = new size_t[numDimensions];
@@ -235,6 +239,7 @@ namespace MDEvents
     return out;
   }
 
+
   //----------------------------------------------------------------------------------------------
   /** Return the position of the center of a bin at a given position
    *
@@ -256,6 +261,29 @@ namespace MDEvents
   }
 
 
+  //----------------------------------------------------------------------------------------------
+  /** Get the signal at a particular coordinate in the workspace.
+   *
+   * @param coords :: numDimensions-sized array of the coordinates to look at
+   * @return the (normalized) signal at a given coordinates.
+   *         NaN if outside the range of this workspace
+   */
+  signal_t MDHistoWorkspace::getSignalAtCoord(const coord_t * coords) const
+  {
+    // Build up the linear index, dimension by dimension
+    size_t linearIndex = 0;
+    for (size_t d=0; d<numDimensions; d++)
+    {
+      coord_t x = coords[d] - m_origin[d];
+      size_t ix = size_t(x / m_boxLength[d]);
+      if (ix >= m_indexMax[d] || (x<0))
+        return std::numeric_limits<signal_t>::quiet_NaN();
+      linearIndex += ix * m_indexMaker[d];
+    }
+    return m_signals[linearIndex];
+  }
+
+
   //----------------------------------------------------------------------------------------------
   /// Creates a new iterator pointing to the first cell in the workspace
   Mantid::API::IMDIterator* MDHistoWorkspace::createIterator() const

Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h

 #ifndef MANTID_MDEVENTS_MDHISTOWORKSPACETEST_H_
 #define MANTID_MDEVENTS_MDHISTOWORKSPACETEST_H_
 
+#include "MantidAPI/IMDIterator.h"
 #include "MantidGeometry/MDGeometry/MDHistoDimension.h"
 #include "MantidKernel/System.h"
 #include "MantidKernel/Timer.h"
+#include "MantidKernel/VMD.h"
 #include "MantidMDEvents/MDHistoWorkspace.h"
+#include "MantidMDEvents/MDHistoWorkspaceIterator.h"
+#include "MantidTestHelpers/MDEventsTestHelper.h"
 #include <boost/math/special_functions/fpclassify.hpp>
 #include <cxxtest/TestSuite.h>
 #include <iomanip>
 #include <iostream>
-#include "MantidKernel/VMD.h"
-#include "MantidAPI/IMDIterator.h"
-#include "MantidMDEvents/MDHistoWorkspaceIterator.h"
+#include "MantidAPI/IMDWorkspace.h"
 
 using namespace Mantid::MDEvents;
 using namespace Mantid::Geometry;
+using namespace Mantid::API;
+using namespace Mantid::Kernel;
 using namespace Mantid;
-using Mantid::Kernel::VMD;
-using Mantid::API::IMDIterator;
 
 class MDHistoWorkspaceTest : public CxxTest::TestSuite
 {
@@ -339,6 +341,25 @@ public:
   }
 
 
+  //---------------------------------------------------------------------------------------------------
+  void test_getSignalAtCoord()
+  {
+    // 2D workspace with signal[i] = i (linear index)
+    MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, 2, 10);
+    for (size_t i=0; i<100; i++)
+      ws->setSignalAt(i, double(i));
+    IMDWorkspace_sptr iws(ws);
+    TS_ASSERT_DELTA( iws->getSignalAtCoord(VMD(0.5, 0.5)), 0.0, 1e-6);
+    TS_ASSERT_DELTA( iws->getSignalAtCoord(VMD(1.5, 0.5)), 1.0, 1e-6);
+    TS_ASSERT_DELTA( iws->getSignalAtCoord(VMD(1.5, 1.5)), 11.0, 1e-6);
+    TS_ASSERT_DELTA( iws->getSignalAtCoord(VMD(9.5, 9.5)), 99.0, 1e-6);
+    // Out of range = NaN
+    TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(-0.01, 2.5)) ) );
+    TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(3.5, -0.02)) ) );
+    TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(10.01, 2.5)) ) );
+    TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(3.5, 10.02)) ) );
+  }
+
 };