refs #10530. Add extents calculation.

Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/CutMD.py

@@ -1,6 +1,7 @@
 from mantid.kernel import *
 from mantid.api import *
 from mantid.simpleapi import *
+import numpy as np
 import os.path
 
 
@@ -26,11 +27,110 @@ def PyInit(self):
 
         self.declareProperty(FloatArrayProperty(name='P3Bin', values=[]),
                              doc='Projection 3 binning')
+
+        self.declareProperty(FloatArrayProperty(name='P4Bin', values=[]),
+                             doc='Projection 4 binning')
 
         self.declareProperty(IMDWorkspaceProperty('OutputWorkspace', '',
                              direction=Direction.Output),
                              doc='Output cut workspace')
 
+
+    def __to_mantid_slicing_binning(self, horace_binning, to_cut, dimension_index):
+        dim = to_cut.getDimension(dimension_index)
+        min = dim.getMinimum()
+        max = dim.getMaximum()
+        n_arguments = len(horace_binning)
+        if n_arguments == 0:
+            raise ValueError("binning parameter cannot be empty")
+        elif n_arguments == 1:
+            step_size = horace_binning[0]
+            n_bins = (max - min) / step_size
+        elif n_arguments == 2:
+            min = horace_binning[0]
+            max = horace_binning[1]
+            n_bins = 1
+        elif n_arguments == 3:
+            step_size = horace_binning[1]
+            n_bins = (max - min) / step_size
+            min = horace_binning[0]
+            max = horace_binning[2]
+            pass
+        else:
+            raise ValueError("Too many arguments given to the binning parameter")
+        if min >= max:
+            raise ValueError("Dimension Min >= Max value. Min %.2f Max %.2f", min, max)
+        if n_bins < 1:
+            raise ValueError("Number of bins calculated to be < 1")
+        return (min, max, n_bins)
+
+
+    def __calculate_extents(self, v, u, w, ws):
+        M=np.array([u,v,w])
+        Minv=np.linalg.inv(M)
+
+        # We are assuming that the workspace has dimensions H, K, L in that order. The workspace MUST have 3 dimensions at least for the following to work.
+        Hrange=[ws.getDimension(0).getMinimum(), ws.getDimension(0).getMaximum()]
+        Krange=[ws.getDimension(1).getMinimum(), ws.getDimension(1).getMaximum()]
+        Lrange=[ws.getDimension(2).getMinimum(), ws.getDimension(2).getMaximum()]
+
+        # Create a numpy 2D array. Makes finding minimums and maximums for each transformed coordinates over every corner easier.
+        new_coords = np.empty([8, 3])
+        counter = 0
+        for h in Hrange:
+           for k in Krange:
+               for l in Lrange:
+                   original_corner=np.array([h,k,l])
+                   new_coords[counter]=original_corner.dot(Minv)
+                   counter += 1
+
+        # Get the min max extents
+        extents = list()
+        print new_coords
+        for i in range(0,3):
+            # Vertical slice down each corner for each dimension, then determine the max, min and use as extents
+            extents.append(np.amin(new_coords[:,i]))
+            extents.append(np.amax(new_coords[:,i]))
+
+            # Copy extents for non crystallographic dimensions    
+        non_crystallographic_dimensions = ws.getNumDims() - 3
+        if non_crystallographic_dimensions > 0:
+            for i in range(0, non_crystallographic_dimensions):
+                extents.append(ws.getDimension(i + 3).getMinimum())
+                extents.append(ws.getDimension(i + 3).getMaximum())
+
+        return extents
+
+    def __uvw_from_projection_table(self, projection_table):
+        if not isinstance(projection_table, ITableWorkspace):
+            I = np.identity(3)
+            return (I[0, :], I[1, :], I[2, :])
+        u = np.array(projection_table.column('u'))
+        v = np.array(projection_table.column('v'))
+        w = np.cross(v,u)
+        return (u, v, w)
+
+    def __make_labels(self, projection):
+
+        class Mapping:
+
+            def __init__(self, replace):
+                self.__replace = replace
+
+            def replace(self, entry):
+                if entry > 0:
+                    return self.__replace 
+                return entry
+
+        crystallographic_names = ['zeta', 'eta', 'xi' ]  
+        labels = list()
+        for i in range(len(projection)):
+            cmapping = Mapping(crystallographic_names[i])
+            labels.append( [cmapping.replace(x) for x in projection[i]  ] )
+
+        return labels
+
+
     def PyExec(self):
         to_cut = self.getProperty("InputWorkspace").value
 
@@ -39,25 +139,67 @@ def PyExec(self):
             raise ValueError("Input Workspace must be in reciprocal lattice dimensions (HKL)")
 
 
-        projection = self.getProperty("Projection").value
-        if isinstance(projection, ITableWorkspace):
-            column_names = set(projection.getColumnNames())
+        projection_table = self.getProperty("Projection").value
+        if isinstance(projection_table, ITableWorkspace):
+            column_names = set(projection_table.getColumnNames())
             logger.warning(str(column_names))
-            if not column_names == set(['u', 'v']):
-                if not column_names == set(['u', 'v', 'offsets']):
-                    if not column_names == set(['u', 'v', 'w', 'offsets']):
+            if not column_names == set(['u', 'v', 'type']):
+                if not column_names == set(['u', 'v', 'offsets', 'type']):
+                    if not column_names == set(['u', 'v', 'w', 'offsets', 'type']):
                         raise ValueError("Projection table schema is wrong")
-            if projection.rowCount() != 3:
+            if projection_table.rowCount() != 3:
                 raise ValueError("Projection table expects 3 rows")
+
+        p1_bins = self.getProperty("P1Bin").value
+        p2_bins = self.getProperty("P2Bin").value
+        p3_bins = self.getProperty("P3Bin").value
+        p4_bins = self.getProperty("P4Bin").value # TODO handle 3D only slicing.
+
+        # TODO. THESE ARE WRONG. Need to consider the acutal transformed extents as part of this.
+        xbins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 0);
+        ybins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 1);
+        zbins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 2);
+        #ebins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 3); # TODO. cannot guarantee this one is here
+
+        # TODO. check that workspace is x=H, y=K, z=L before using the projections and slicing.
+        projection = self.__uvw_from_projection_table(projection_table)
+        u,v,w = projection
+
+        extents = self.__calculate_extents(v, u, w, to_cut)
 
+        nopix = True # TODO from input parameter.
+
+        cut_alg_name = "SliceMD"
+        if nopix:
+            cut_alg_name= "BinMD"
+
+        projection_labels = self.__make_labels(projection)
 
-        clone_alg = AlgorithmManager.Instance().create("CloneMDWorkspace")
-        clone_alg.setChild(True)
-        clone_alg.initialize()
-        clone_alg.setProperty("InputWorkspace", to_cut)
-        clone_alg.setPropertyValue("OutputWorkspace", "cloned")
-        clone_alg.execute()
-        cloned = clone_alg.getProperty("OutputWorkspace").value
-        self.setProperty("OutputWorkspace", cloned)
+        cut_alg = AlgorithmManager.Instance().create(cut_alg_name)
+        cut_alg.setChild(True)
+        cut_alg.initialize()
+        cut_alg.setProperty("InputWorkspace", to_cut)
+        cut_alg.setPropertyValue("OutputWorkspace", "sliced")
+        cut_alg.setProperty("NormalizeBasisVectors", False)
+        cut_alg.setProperty("AxisAligned", False)
+        # Now for the basis vectors.
+        for i in range(0, to_cut.getNumDims()):
+            logger.warning("WARNING....")
+            if i <= 2:
+                label = projection_labels[i]
+                unit = "TODO" # Haven't figured out how to do this yet.
+                vec = projection[i]
+                value = "%s, %s, %s" % ( label, unit, ",".join(map(str, vec))) 
+                cut_alg.setPropertyValue("BasisVector{0}".format(i) , value)
+                logger.warning("BasisVector{0}".format(i))
+            if i > 2:
+                raise RuntimeError("Not implmented yet for non-crystallographic basis vector generation.")
+        cut_alg.setProperty("OutputExtents", extents)
+        bins = [ int(xbins[2]), int(ybins[2]), int(zbins[2]) ] # Again this is a hack for 3 dimensional data only.
+        cut_alg.setProperty("OutputBins", bins)
+
+        cut_alg.execute()
+        slice = cut_alg.getProperty("OutputWorkspace").value
+        self.setProperty("OutputWorkspace", slice)
 
 AlgorithmFactory.subscribe(CutMD)

Code/Mantid/Framework/PythonInterface/test/python/mantid/api/CutMDTest.py

@@ -2,6 +2,7 @@
 import testhelpers
 
 from mantid.simpleapi import *
+from mantid.api import IMDHistoWorkspace
 
 
 class CutMDTest(unittest.TestCase):
@@ -23,18 +24,27 @@ def test_exec_throws_if_not_a_hkl_workspace(self):
         test_md = CreateMDWorkspace(Dimensions=3, Extents=[-10,10,-10,10,-10,10], Names="A,B,C", Units="U,U,U")
         # Explicitly set the coordinate system to lab Q.
         SetSpecialCoordinates(InputWorkspace=test_md, SpecialCoordinates='Q (lab frame)')
-        self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws")
+        self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1])
 
     def test_slice_to_original(self):
-        out_md = CutMD(self.__in_md)
-        equals, result = CompareMDWorkspaces(Workspace1=self.__in_md, Workspace2=out_md)
-        self.assertTrue(equals, "Input and output workspaces should be identical" )
+        out_md = CutMD(self.__in_md, P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1])
+        self.assertTrue(isinstance(out_md, IMDHistoWorkspace))
+        # No rotation. Basis vectors should have been left the same, so no extent changes.
+        self.assertEquals(self.__in_md.getDimension(0).getMinimum(), out_md.getDimension(0).getMinimum())
+        self.assertEquals(self.__in_md.getDimension(0).getMaximum(), out_md.getDimension(0).getMaximum())
+        self.assertEquals(self.__in_md.getDimension(1).getMinimum(), out_md.getDimension(1).getMinimum())
+        self.assertEquals(self.__in_md.getDimension(1).getMaximum(), out_md.getDimension(1).getMaximum())
+        self.assertEquals(self.__in_md.getDimension(2).getMinimum(), out_md.getDimension(2).getMinimum())
+        self.assertEquals(self.__in_md.getDimension(2).getMaximum(), out_md.getDimension(2).getMaximum())
+        self.assertEquals("['zeta', 0.0, 0.0]",  out_md.getDimension(0).getName() )
+        self.assertEquals("[0.0, 'eta', 0.0]",  out_md.getDimension(1).getName() )
+        self.assertEquals("[0.0, 0.0, 'xi']",  out_md.getDimension(2).getName() )
 
     def test_wrong_projection_workspace_format_wrong_column_numbers(self):
         projection = CreateEmptyTableWorkspace()
         projection.addColumn("double", "u")
         # missing other columns
-        self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws")
+        self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1])
 
     def test_wrong_table_workspace_format_wrong_row_numbers(self):
         projection = CreateEmptyTableWorkspace()
@@ -43,8 +53,9 @@ def test_wrong_table_workspace_format_wrong_row_numbers(self):
         projection.addColumn("double", "v")
         projection.addColumn("double", "w")
         projection.addColumn("double", "offset")
+        projection.addColumn("string", "type")
         # Incorrect number of rows i.e. zero in this case as none added.
-        self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws")
+        self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1])
 
     def test_run_it(self):
         pass