Merge pull request #1033 from vasole/nnma

[ROI Imaging] Support NNMA on stacks of images

PyMca5/PyMcaMath/mva/NNMAModule.py

@@ -63,7 +63,7 @@
                              from numpy or sparse from scipy.sparse
                              package
 
-            k           --   number of componnets to estimate
+            k           --   number of components to estimate
 
             Astart
             Xstart      --   matrices to start iterations. Maybe None
@@ -185,17 +185,6 @@
 """
 import numpy
 import logging
-try:
-    import os
-    os.environ["MDP_DISABLE_SKLEARN"] = "yes"
-    import mdp
-    if mdp.__version__ >= '2.6':
-        MDP = True
-    else:
-        MDP = False
-except Exception:
-    MDP = False
-
 from . import py_nnma
 
 
@@ -223,8 +212,8 @@ def nnma(stack, ncomponents, binning=None,
          mask=None, spectral_mask=None,
          function=None, eps=5e-5, verbose=VERBOSE,
          maxcount=1000, kmeans=False):
-    if kmeans and (not MDP):
-        raise ValueError("K Means not supported")
+    if kmeans:
+        raise ValueError("K Means not supported by this module")
     #I take the defaults for the other parameters
     param = dict(alpha=.1, tau=2, regul=1e-2, sparse_par=1e-1, psi=1e-3)
     if function is None:
@@ -384,13 +373,9 @@ def nnma(stack, ncomponents, binning=None,
     original_intensity = numpy.sum(data)
 
     #final values
-    if kmeans:
-        n_more = 1
-    else:
-        n_more = 0
-    new_images  = numpy.zeros((ncomponents + n_more, r*c), numpy.float32)
-    new_vectors = numpy.zeros((X.shape[0]+n_more, X.shape[1]), numpy.float32)
-    values      = numpy.zeros((ncomponents+n_more,), numpy.float32)
+    new_images  = numpy.zeros((ncomponents, r*c), numpy.float32)
+    new_vectors = numpy.zeros((X.shape[0], X.shape[1]), numpy.float32)
+    values      = numpy.zeros((ncomponents,), numpy.float32)
     for i in range(ncomponents):
         idx = sorted_idx[i]
         if 1:
@@ -407,17 +392,7 @@ def nnma(stack, ncomponents, binning=None,
             new_images[i, maskview] = numpy.sum(numpy.dot(Atmp, Xtmp), axis=1)
         new_vectors[i,:] = X[idx,:]
         values[i] = 100.*total_nnma_intensity[idx][0]/original_intensity
-    new_images.shape = ncomponents + n_more, r, c
-    if kmeans:
-        classifier = mdp.nodes.KMeansClassifier(ncomponents)
-        for i in range(ncomponents):
-            classifier.train(new_vectors[i:i+1])
-        k = 0
-        for i in range(r):
-            for j in range(c):
-                spectrum = data[k:k+1,:]
-                new_images[-1, i,j] = classifier.label(spectrum)[0]
-                k += 1
+    new_images.shape = ncomponents, r, c
     return new_images, values, new_vectors
 
 if __name__ == "__main__":

PyMca5/PyMcaPlugins/NNMAStackPlugin.py

@@ -138,8 +138,8 @@ def calculate(self):
         mcaIndex = stack.info.get('McaIndex')
         shape = stack.data.shape
         stack = None
-        if mcaIndex not in [-1, len(shape) - 1]:
-            raise IndexError("NNMA does not support stacks of images yet")
+        if mcaIndex not in [0, -1, len(shape) - 1]:
+            raise IndexError("NNMA only support stacks of images or spectra")
             return
         if self.configurationWidget is None:
             self.configurationWidget = NNMAParametersDialog(None, regions=True)
@@ -225,9 +225,23 @@ def actualCalculation(self):
                 self._status.setText(text)
 
         oldShape = stack.data.shape
-        result = function(stack, **ddict)
-        if stack.data.shape != oldShape:
-            stack.data.shape = oldShape
+        mcaIndex = stack.info.get('McaIndex')
+        if mcaIndex == 0:
+            # image stack. We need a copy
+            _logger.info("NNMAStackPlugin converting to stack of spectra")
+            data = numpy.zeros(oldShape[1:] + oldShape[0:1], dtype=numpy.float32)
+            data.shape = -1, oldShape[0]
+            for i in range(oldShape[0]):
+                tmpData = stack.data[i]
+                tmpData.shape = -1
+                data[:, i] = tmpData
+            data.shape = oldShape[1:] + oldShape[0:1]
+            result = function(data, **ddict)
+            data = None
+        else:
+            result = function(stack, **ddict)
+            if stack.data.shape != oldShape:
+                stack.data.shape = oldShape
         return result
 
     def threadFinished(self):