#/*########################################################################## # # The PyMca X-Ray Fluorescence Toolkit # # Copyright (c) 2004-2022 European Synchrotron Radiation Facility # # This file is part of the PyMca X-ray Fluorescence Toolkit developed at # the ESRF. # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. # #############################################################################*/ __author__ = "V.A. Sole" __contact__ = "sole@esrf.fr" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" """ Base class to handle stacks. """ from PyMca5.PyMcaCore import DataObject import numpy import time import os import sys import glob import logging logger = logging.getLogger(__name__) PLUGINS_DIR = None try: import PyMca5 if os.path.exists(os.path.join(os.path.dirname(__file__), "PyMcaPlugins")): from PyMca5 import PyMcaPlugins PLUGINS_DIR = os.path.dirname(PyMcaPlugins.__file__) else: directory = os.path.dirname(__file__) while True: if os.path.exists(os.path.join(directory, "PyMcaPlugins")): PLUGINS_DIR = os.path.join(directory, "PyMcaPlugins") break directory = os.path.dirname(directory) if len(directory) < 5: break userPluginsDirectory = PyMca5.getDefaultUserPluginsDirectory() PYMCA_PLUGINS_DIR = PLUGINS_DIR if userPluginsDirectory is not None: if PLUGINS_DIR is None: PLUGINS_DIR = userPluginsDirectory else: PLUGINS_DIR = [PLUGINS_DIR, userPluginsDirectory] except: PYMCA_PLUGINS_DIR = None pass class StackBase(object): def __init__(self): self._stack = DataObject.DataObject() self._stack.x = None self._stackImageData = None self._selectionMask = None self._finiteData = True self._ROIDict = {'name': "ICR", 'type': "CHANNEL", 'calibration': [0, 1.0, 0.0], 'from': 0, 'to': -1} self._ROIImageDict = {'ROI': None, 'Maximum': None, 'Minimum': None, 'Left': None, 'Middle': None, 'Right': None, 'Background': None} self.__ROIImageCalculationIsUsingSuppliedEnergyAxis = False self._ROIImageList = [] self._ROIImageNames = [] self.__pluginDirList = [] self.pluginList = [] self.pluginInstanceDict = {} self.getPlugins() # beyond 5 million elements, iterate to calculate the sums # preventing huge intermediate use of memory when calculating # the sums. self._dynamicLimit = 5.0E6 self._tryNumpy = True def setPluginDirectoryList(self, dirlist): for directory in dirlist: if not os.path.exists(directory): raise IOError("Directory:\n%s\ndoes not exist." % directory) self.__pluginDirList = dirlist def getPluginDirectoryList(self): return self.__pluginDirList def getPlugins(self): """ Import or reloads all the available plugins. It returns the number of plugins loaded. """ if PLUGINS_DIR is not None: if self.__pluginDirList == []: if type(PLUGINS_DIR) == type([]): self.__pluginDirList = PLUGINS_DIR else: self.__pluginDirList = [PLUGINS_DIR] self.pluginList = [] for directory in self.__pluginDirList: if directory is None: continue if not os.path.exists(directory): raise IOError("Directory:\n%s\ndoes not exist." % directory) fileList = glob.glob(os.path.join(directory, "*.py")) targetMethod = 'getStackPluginInstance' # prevent unnecessary imports moduleList = [] for fname in fileList: # in Python 3, rb implies bytes and not strings f = open(fname, 'r') lines = f.readlines() f.close() f = None for line in lines: if line.startswith("def"): if line.split(" ")[1].startswith(targetMethod): moduleList.append(fname) break for module in moduleList: try: pluginName = os.path.basename(module)[:-3] if directory == PYMCA_PLUGINS_DIR: plugin = "PyMcaPlugins." + pluginName else: plugin = pluginName if directory not in sys.path: sys.path.insert(0, directory) if pluginName in self.pluginList: idx = self.pluginList.index(pluginName) del self.pluginList[idx] if plugin in self.pluginInstanceDict.keys(): del self.pluginInstanceDict[plugin] if plugin in sys.modules: if hasattr(sys.modules[plugin], targetMethod): if sys.version < '3.0': reload(sys.modules[plugin]) else: import imp imp.reload(sys.modules[plugin]) else: try: __import__(plugin) except: if directory == PYMCA_PLUGINS_DIR: plugin = "PyMca5.PyMcaPlugins." + pluginName __import__(plugin) else: raise if hasattr(sys.modules[plugin], targetMethod): self.pluginInstanceDict[plugin] = \ sys.modules[plugin].getStackPluginInstance(self) self.pluginList.append(plugin) except: logger.debug("Problem importing module %s", plugin) if logger.getEffectiveLevel() == logging.DEBUG: raise return len(self.pluginList) def setStack(self, stack, mcaindex=2, fileindex=None): #unfortunaly python 3 reports #isinstance(stack, DataObject.DataObject) as false #for DataObject derived classes like OmnicMap!!!! if id(stack) == id(self._stack): # just updated pass elif hasattr(stack, "shape") and\ hasattr(stack, "dtype"): # array like self._stack.x = None self._stack.data = stack self._stack.info['SourceName'] = "Data of unknown origin" elif isinstance(stack, DataObject.DataObject) or\ ("DataObject.DataObject" in ("%s" % type(stack))) or\ ("QStack" in ("%s" % type(stack))) or\ ("Map" in ("%s" % type(stack)))or\ ("Stack" in ("%s" % type(stack))) or\ (hasattr(stack, "data") and hasattr(stack, "info")): self._stack = stack self._stack.info['SourceName'] = stack.info.get('SourceName', "Data of unknown origin") else: self._stack.x = None self._stack.data = stack self._stack.info['SourceName'] = "Data of unknown origin" info = self._stack.info mcaIndex = info.get('McaIndex', mcaindex) if (mcaIndex < 0) and (len(self._stack.data.shape) == 3): mcaIndex = len(self._stack.data.shape) + mcaIndex fileIndex = info.get('FileIndex', fileindex) if fileIndex is None: if mcaIndex == 2: fileIndex = 0 elif mcaIndex == 0: fileIndex = 1 else: fileIndex = 0 for i in range(3): if i not in [mcaIndex, fileIndex]: otherIndex = i break self.mcaIndex = mcaIndex self.fileIndex = fileIndex self.otherIndex = otherIndex self._stack.info['McaCalib'] = info.get('McaCalib', [0.0, 1.0, 0.0]) self._stack.info['Channel0'] = info.get('Channel0', 0.0) self._stack.info['McaIndex'] = mcaIndex self._stack.info['FileIndex'] = fileIndex self._stack.info['OtherIndex'] = otherIndex self.stackUpdated(info.get("positioners", None)) def stackUpdated(self, positioners=None): """ Recalculates the different images associated to the stack """ self._tryNumpy = True if hasattr(self._stack.data, "size"): if self._stack.data.size > self._dynamicLimit: self._tryNumpy = False else: # is not a numpy ndarray in any case self._tryNumpy = False previousStackImageSize = None if self._stackImageData is not None: previousStackImageSize = self._stackImageData.size mcaMax = None if self._tryNumpy and isinstance(self._stack.data, numpy.ndarray): self._stackImageData = numpy.sum(self._stack.data, axis=self.mcaIndex, dtype=numpy.float64) #original ICR mca logger.debug("(self.otherIndex, self.fileIndex) = (%d, %d)", self.otherIndex, self.fileIndex) i = max(self.otherIndex, self.fileIndex) j = min(self.otherIndex, self.fileIndex) mcaData0 = numpy.sum(numpy.sum(self._stack.data, axis=i, dtype=numpy.float64), j) # max MCA if self.mcaIndex == -1 or (self.mcaIndex == (len(self._stack.data.shape) - 1)): mcaMax = numpy.nanmax(numpy.nanmax(self._stack.data, axis=0), axis=0) elif self.mcaIndex == 0: mcaMax = numpy.nanmax(numpy.nanmax(self._stack.data, axis=-1), axis=-1) else: _logger.info("Unsupported index for max spectrum calculation") else: t0 = time.time() shape = self._stack.data.shape if self.mcaIndex in [2, -1]: self._stackImageData = numpy.zeros((shape[0], shape[1]), dtype=numpy.float64) mcaData0 = numpy.zeros((shape[2],), numpy.float64) mcaMax = mcaData0 + numpy.NINF step = 1 for i in range(shape[0]): tmpData = self._stack.data[i:i+step,:,:] numpy.add(self._stackImageData[i:i+step,:], numpy.sum(tmpData, 2), self._stackImageData[i:i+step,:]) tmpData.shape = step*shape[1], shape[2] tmpMax = numpy.nanmax(tmpData, axis=0) numpy.add(mcaData0, numpy.sum(tmpData, 0), mcaData0) mcaMax =numpy.max([mcaMax, tmpMax], axis=0) elif self.mcaIndex == 0: self._stackImageData = numpy.zeros((shape[1], shape[2]), dtype=numpy.float64) mcaData0 = numpy.zeros((shape[0],), numpy.float64) mcaMax = mcaData0 + numpy.NINF step = 1 for i in range(shape[0]): tmpData = self._stack.data[i:i+step,:,:] tmpData.shape = tmpData.shape[1:] numpy.add(self._stackImageData, tmpData, self._stackImageData) mcaData0[i] = tmpData.sum() mcaMax[i] = numpy.nanmax(tmpData) else: raise ValueError("Unhandled case 1D index = %d" % self.mcaIndex) logger.debug("Print dynamic loading elapsed = %f", time.time() - t0) logger.debug("__stackImageData.shape = %s", self._stackImageData.shape) if previousStackImageSize: if previousStackImageSize != self._stackImageData.size: self._clearPositioners() calib = self._stack.info.get('McaCalib', [0.0, 1.0, 0.0]) dataObject = DataObject.DataObject() dataObject.info = {"McaCalib": calib, "selectiontype": "1D", "SourceName": "Stack", "Key": "SUM"} if "McaLiveTime" in self._stack.info: if hasattr(self._stack.info["McaLiveTime"], "sum"): dataObject.info["McaLiveTime"] = \ self._stack.info["McaLiveTime"].sum() else: print("Not an array. Skipping time information") #del dataObject.info["McaLiveTime"] if not hasattr(self._stack, 'x'): self._stack.x = None if self._stack.x in [None, []]: self._stack.x = [numpy.arange(len(mcaData0)).astype(numpy.float64)+\ self._stack.info.get('Channel0', 0.0)] # for the time being it can only contain one axis dataObject.x = [self._stack.x[0]] dataObject.y = [mcaData0] #store the original spectrum self._mcaData0 = dataObject #store the original max spectrum self._mcaMax = mcaMax #add the original image self.showOriginalImage() #add the mca goodData = numpy.isfinite(self._mcaData0.y[0].sum()) if goodData: self._finiteData = True self.showOriginalMca() else: self._finiteData = False self.handleNonFiniteData() #calculate the ROIs self._ROIDict = {'name': "ICR", 'type': "CHANNEL", 'calibration': calib, 'from': dataObject.x[0][0], 'to': dataObject.x[0][-1]} self.updateROIImages() if positioners is not None: try: self.setPositioners(positioners) except: logging.error("Error setting positioners. Ignoring them") self._clearPositioners() for key in self.pluginInstanceDict.keys(): self.pluginInstanceDict[key].stackUpdated() def getStackOriginalCurve(self): # TODO: Make sure copies are returned x = self._mcaData0.x[0] y = self._mcaData0.y[0] legend = "Stack SUM" info = self._mcaData0.info return [x, y, legend, info] def isStackFinite(self): """ Returns True if stack does not contain inf or nans Returns False if stack is not finite """ return self._finiteData def handleNonFiniteData(self): text = "Your data contain infinite values or nans.\n" text += "Pixels containing those values will be ignored." logger.info(text) def updateROIImages(self, ddict=None): if ddict is None: updateROIDict = False ddict = self._ROIDict else: updateROIDict = True xw = ddict['calibration'][0] + \ ddict['calibration'][1] * self._mcaData0.x[0] + \ ddict['calibration'][2] * (self._mcaData0.x[0] ** 2) if ddict["name"] == "ICR": i1 = 0 i2 = self._stack.data.shape[self.mcaIndex] imiddle = int(0.5 * (i1 + i2)) pos = 0.5 * (ddict['from'] + ddict['to']) if ddict["type"].upper() != "CHANNEL": imiddle = max(numpy.nonzero(xw <= pos)[0]) elif ddict["type"].upper() != "CHANNEL": #energy like ROI if xw[0] < xw[-1]: i1 = numpy.nonzero(ddict['from'] <= xw)[0] if len(i1): i1 = min(i1) else: logger.debug("updateROIImages: nothing to be made") return i2 = numpy.nonzero(xw <= ddict['to'])[0] if len(i2): i2 = max(i2) + 1 else: logger.debug("updateROIImages: nothing to be made") return pos = 0.5 * (ddict['from'] + ddict['to']) imiddle = max(numpy.nonzero(xw <= pos)[0]) else: i2 = numpy.nonzero(ddict['from'] <= xw)[0] if len(i2): i2 = max(i2) else: logger.debug("updateROIImages: nothing to be made") return i1 = numpy.nonzero(xw <= ddict['to'])[0] if len(i1): i1 = min(i1) + 1 else: logger.debug("updateROIImages: nothing to be made") return pos = 0.5 * (ddict['from'] + ddict['to']) imiddle = min(numpy.nonzero(xw <= pos)[0]) else: i1 = numpy.nonzero(ddict['from'] <= self._mcaData0.x[0])[0] if len(i1): if self._mcaData0.x[0][0] > self._mcaData0.x[0][-1]: i1 = max(i1) else: i1 = min(i1) else: i1 = 0 i1 = max(i1, 0) i2 = numpy.nonzero(self._mcaData0.x[0] <= ddict['to'])[0] if len(i2): if self._mcaData0.x[0][0] > self._mcaData0.x[0][-1]: i2 = min(i2) else: i2 = max(i2) else: i2 = 0 i2 = min(i2 + 1, self._stack.data.shape[self.mcaIndex]) pos = 0.5 * (ddict['from'] + ddict['to']) if self._mcaData0.x[0][0] > self._mcaData0.x[0][-1]: imiddle = min(numpy.nonzero(self._mcaData0.x[0] <= pos)[0]) else: imiddle = max(numpy.nonzero(self._mcaData0.x[0] <= pos)[0]) xw = self._mcaData0.x[0] self._ROIImageDict = self.calculateROIImages(i1, i2, imiddle, energy=xw) if updateROIDict: self._ROIDict.update(ddict) roiKeys = ['ROI', 'Maximum', 'Minimum', 'Left', 'Middle', 'Right', 'Background'] nImages = len(roiKeys) imageList = [None] * nImages for i in range(nImages): key = roiKeys[i] imageList[i] = self._ROIImageDict[key] title = "%s" % ddict["name"] if ddict["name"] == "ICR": cursor = "Energy" if abs(ddict['calibration'][0]) < 1.0e-5: if abs(ddict['calibration'][1] - 1) < 1.0e-5: if abs(ddict['calibration'][2]) < 1.0e-5: cursor = "Channel" elif ddict["type"].upper() == "CHANNEL": cursor = "Channel" else: cursor = ddict["type"] imageNames = [title, '%s Maximum' % title, '%s Minimum' % title, '%s %.6g' % (cursor, xw[i1]), '%s %.6g' % (cursor, xw[imiddle]), '%s %.6g' % (cursor, xw[(i2 - 1)]), '%s Background' % title] if self.__ROIImageCalculationIsUsingSuppliedEnergyAxis: imageNames[1] = "%s %s at Max." % (title, cursor) imageNames[2] = "%s %s at Min." % (title, cursor) self.showROIImageList(imageList, image_names=imageNames) def showOriginalImage(self): logger.debug("showOriginalImage to be implemented") def showOriginalMca(self): logger.debug("showOriginalMca to be implemented") def showROIImageList(self, imageList, image_names=None): self._ROIImageList = imageList self._ROIImageNames = image_names self._stackROIImageListUpdated() def _stackROIImageListUpdated(self): for key in self.pluginInstanceDict.keys(): self.pluginInstanceDict[key].stackROIImageListUpdated() def getStackROIImagesAndNames(self): return self._ROIImageList, self._ROIImageNames def getStackOriginalImage(self): return self._stackImageData def calculateMcaDataObject(self, normalize=False, mask=None): #original ICR mca if self._stackImageData is None: return if mask is None: selectionMask = self._selectionMask else: selectionMask = mask mcaData = None goodData = numpy.isfinite(self._mcaData0.y[0].sum()) logger.debug("Stack data is not finite") if (selectionMask is None) and goodData: if normalize: logger.debug("Case 1") npixels = self._stackImageData.shape[0] *\ self._stackImageData.shape[1] * 1.0 dataObject = DataObject.DataObject() dataObject.info.update(self._mcaData0.info) dataObject.x = [self._mcaData0.x[0]] dataObject.y = [self._mcaData0.y[0] / float(npixels)] if "McaLiveTime" in dataObject.info: dataObject.info["McaLiveTime"] /= float(npixels) else: logger.debug("Case 2") dataObject = self._mcaData0 return dataObject #deal with NaN and inf values if selectionMask is None: if (self._ROIImageDict["ROI"] is not None) and\ (self.mcaIndex != 0): actualSelectionMask = numpy.isfinite(self._ROIImageDict["ROI"]) else: actualSelectionMask = numpy.isfinite(self._stackImageData) else: if (self._ROIImageDict["ROI"] is not None) and\ (self.mcaIndex != 0): actualSelectionMask = selectionMask * numpy.isfinite(self._ROIImageDict["ROI"]) else: actualSelectionMask = selectionMask * numpy.isfinite(self._stackImageData) npixels = actualSelectionMask.sum() if (npixels == 0) and goodData: if normalize: logger.debug("Case 3") npixels = self._stackImageData.shape[0] * self._stackImageData.shape[1] * 1.0 dataObject = DataObject.DataObject() dataObject.info.update(self._mcaData0.info) dataObject.x = [self._mcaData0.x[0]] dataObject.y = [self._mcaData0.y[0] / float(npixels)] if "McaLiveTime" in dataObject.info: dataObject.info["McaLiveTime"] /= float(npixels) else: logger.debug("Case 4") dataObject = self._mcaData0 return dataObject mcaData = numpy.zeros(self._mcaData0.y[0].shape, numpy.float64) n_nonselected = self._stackImageData.shape[0] *\ self._stackImageData.shape[1] - npixels if goodData: if n_nonselected < npixels: arrayMask = (actualSelectionMask == 0) else: arrayMask = (actualSelectionMask > 0) else: arrayMask = (actualSelectionMask > 0) logger.debug("Reached MCA calculation") cleanMask = numpy.nonzero(arrayMask) logger.debug("self.fileIndex, self.mcaIndex = %d , %d", self.fileIndex, self.mcaIndex) t0 = time.time() if len(cleanMask[0]) and len(cleanMask[1]): logger.debug("USING MASK") cleanMask = numpy.array(cleanMask).transpose() if self.fileIndex == 2: if self.mcaIndex == 0: if isinstance(self._stack.data, numpy.ndarray): logger.debug("In memory case 0") for r, c in cleanMask: mcaData += self._stack.data[:, r, c] else: logger.debug("Dynamic loading case 0") #no other choice than to read all images #for the time being, one by one rMin = cleanMask[0][0] rMax = cleanMask[-1][0] cMin = cleanMask[:, 1].min() cMax = cleanMask[:, 1].max() #rMin, cMin = cleanMask.min(axis=0) #rMax, cMax = cleanMask.max(axis=0) tmpMask = arrayMask[rMin:(rMax+1),cMin:(cMax+1)] tmpData = numpy.zeros((1, rMax-rMin+1,cMax-cMin+1)) for i in range(self._stack.data.shape[0]): tmpData[0:1,:,:] = self._stack.data[i:i+1,rMin:(rMax+1),cMin:(cMax+1)] #multiplication is faster than selection mcaData[i] = (tmpData[0]*tmpMask).sum(dtype=numpy.float64) elif self.mcaIndex == 1: if isinstance(self._stack.data, numpy.ndarray): for r, c in cleanMask: mcaData += self._stack.data[r,:,c] else: raise IndexError("Dynamic loading case 1") else: raise IndexError("Wrong combination of indices. Case 0") elif self.fileIndex == 1: if self.mcaIndex == 0: if isinstance(self._stack.data, numpy.ndarray): logger.debug("In memory case 2") for r, c in cleanMask: mcaData += self._stack.data[:, r, c] else: logger.debug("Dynamic loading case 2") #no other choice than to read all images #for the time being, one by one if 1: rMin = cleanMask[0][0] rMax = cleanMask[-1][0] cMin = cleanMask[:, 1].min() cMax = cleanMask[:, 1].max() #rMin, cMin = cleanMask.min(axis=0) #rMax, cMax = cleanMask.max(axis=0) tmpMask = arrayMask[rMin:(rMax + 1), cMin:(cMax + 1)] tmpData = numpy.zeros((1, rMax - rMin + 1, cMax - cMin + 1)) for i in range(self._stack.data.shape[0]): tmpData[0:1, :, :] = self._stack.data[i:i + 1, rMin:(rMax + 1), cMin:(cMax + 1)] #multiplication is faster than selection mcaData[i] = (tmpData[0] * tmpMask).sum(dtype=numpy.float64) if 0: tmpData = numpy.zeros((1, self._stack.data.shape[1], self._stack.data.shape[2])) for i in range(self._stack.data.shape[0]): tmpData[0:1, :, :] = self._stack.data[i:i + 1,:,:] #multiplication is faster than selection #tmpData[arrayMask].sum() in my machine mcaData[i] = (tmpData[0] * arrayMask).sum(dtype=numpy.float64) elif self.mcaIndex == 2: if isinstance(self._stack.data, numpy.ndarray): logger.debug("In memory case 3") for r, c in cleanMask: mcaData += self._stack.data[r, c, :] else: logger.debug("Dynamic loading case 3") #try to minimize access to the file row_list = [] row_dict = {} for r, c in cleanMask: if r not in row_list: row_list.append(r) row_dict[r] = [] row_dict[r].append(c) for r in row_list: tmpMcaData = self._stack.data[r:r + 1, row_dict[r], :] tmpMcaData.shape = -1, mcaData.shape[0] mcaData += numpy.sum(tmpMcaData, axis=0, dtype=numpy.float64) else: raise IndexError("Wrong combination of indices. Case 1") elif self.fileIndex == 0: if self.mcaIndex == 1: if isinstance(self._stack.data, numpy.ndarray): logger.debug("In memory case 4") for r, c in cleanMask: mcaData += self._stack.data[r, :, c] else: raise IndexError("Dynamic loading case 4") elif self.mcaIndex in [2, -1]: if isinstance(self._stack.data, numpy.ndarray): logger.debug("In memory case 5") for r, c in cleanMask: mcaData += self._stack.data[r, c, :] else: logger.debug("Dynamic loading case 5") #try to minimize access to the file row_list = [] row_dict = {} for r, c in cleanMask: if r not in row_list: row_list.append(r) row_dict[r] = [] row_dict[r].append(c) for r in row_list: tmpMcaData = self._stack.data[r:r + 1, row_dict[r], :] tmpMcaData.shape = -1, mcaData.shape[0] mcaData += tmpMcaData.sum(axis=0, dtype=numpy.float64) else: raise IndexError("Wrong combination of indices. Case 2") else: raise IndexError("File index undefined") else: logger.debug("NOT USING MASK !") logger.debug("Mca sum elapsed = %f", time.time() - t0) if goodData: if n_nonselected < npixels: mcaData = self._mcaData0.y[0] - mcaData if normalize: mcaData = mcaData / float(npixels) calib = self._stack.info['McaCalib'] dataObject = DataObject.DataObject() dataObject.info = {"McaCalib": calib, "selectiontype": "1D", "SourceName": "Stack", "Key": "Selection"} if "McaLiveTime" in self._stack.info: selectedPixels = actualSelectionMask > 0 liveTime = self._stack.info["McaLiveTime"][:] liveTime.shape = actualSelectionMask.shape liveTime = liveTime[selectedPixels].sum() if normalize: liveTime = liveTime / float(npixels) dataObject.info["McaLiveTime"] = liveTime dataObject.x = [self._mcaData0.x[0]] dataObject.y = [mcaData] return dataObject def calculateROIImages(self, index1, index2, imiddle=None, energy=None): logger.debug("Calculating ROI images") i1 = min(index1, index2) i2 = max(index1, index2) if imiddle is None: imiddle = int(0.5 * (i1 + i2)) if energy is None: energy = self._mcaData0.x[0] if i1 == i2: dummy = numpy.zeros(self._stackImageData.shape, numpy.float64) imageDict = {'ROI': dummy, 'Maximum': dummy, 'Minimum': dummy, 'Left': dummy, 'Middle': dummy, 'Right': dummy, 'Background': dummy} return imageDict isUsingSuppliedEnergyAxis = False if self.fileIndex == 0: if self.mcaIndex == 1: leftImage = self._stack.data[:, i1, :] middleImage = self._stack.data[:, imiddle, :] rightImage = self._stack.data[:, i2 - 1, :] dataImage = self._stack.data[:, i1:i2, :] background = 0.5 * (i2 - i1) * (leftImage + rightImage) roiImage = numpy.sum(dataImage, axis=1, dtype=numpy.float64) maxImage = energy[(numpy.argmax(dataImage, axis=1) + i1)] minImage = energy[(numpy.argmin(dataImage, axis=1) + i1)] isUsingSuppliedEnergyAxis = True else: t0 = time.time() if self._tryNumpy and\ isinstance(self._stack.data, numpy.ndarray): leftImage = self._stack.data[:, :, i1] middleImage = self._stack.data[:, :, imiddle] rightImage = self._stack.data[:, :, i2 - 1] dataImage = self._stack.data[:, :, i1:i2] background = 0.5 * (i2 - i1) * (leftImage + rightImage) roiImage = numpy.sum(dataImage, axis=2, dtype=numpy.float64) maxImage = energy[numpy.argmax(dataImage, axis=2) + i1] minImage = energy[numpy.argmin(dataImage, axis=2) + i1] isUsingSuppliedEnergyAxis = True logger.debug("Case 1 ROI image calculation elapsed = %f ", time.time() - t0) else: shape = self._stack.data.shape roiImage = numpy.zeros(self._stackImageData.shape, numpy.float64) background = roiImage * 1 leftImage = roiImage * 1 middleImage = roiImage * 1 rightImage = roiImage * 1 maxImage = numpy.zeros(self._stackImageData.shape, numpy.uint) minImage = numpy.zeros(self._stackImageData.shape, numpy.uint) step = 1 for i in range(shape[0]): tmpData = self._stack.data[i:i+step,:, i1:i2] * 1 numpy.add(roiImage[i:i+step,:], numpy.sum(tmpData, axis=2,dtype=numpy.float64), roiImage[i:i+step,:]) minImage[i:i + step,:] = i1 + numpy.argmin(tmpData, axis=2) maxImage[i:i + step, :] = i1 + numpy.argmax(tmpData, axis=2) leftImage[i:i + step, :] += tmpData[:, :, 0] middleImage[i:i + step, :] += tmpData[:, :, imiddle - i1] rightImage[i:i + step, :] += tmpData[:, :, -1] background = 0.5 * (i2 - i1) * (leftImage + rightImage) isUsingSuppliedEnergyAxis = True minImage = energy[minImage] maxImage = energy[maxImage] logger.debug("2 Dynamic ROI image calculation elapsed = %f ", time.time() - t0) elif self.fileIndex == 1: if self.mcaIndex == 0: t0 = time.time() if isinstance(self._stack.data, numpy.ndarray) and\ self._tryNumpy: leftImage = self._stack.data[i1, :, :] middleImage= self._stack.data[imiddle, :, :] rightImage = self._stack.data[i2 - 1, :, :] dataImage = self._stack.data[i1:i2, :, :] # this calculation is very slow but it is extremely useful # for XANES studies if 1: maxImage = energy[numpy.argmax(dataImage, axis=0) + i1] minImage = energy[numpy.argmin(dataImage, axis=0) + i1] else: # this is slower, but uses less memory maxImage = numpy.zeros(leftImage.shape, numpy.int32) minImage = numpy.zeros(leftImage.shape, numpy.int32) for i in range(i1, i2): tmpData = self._stack.data[i] tmpData.shape = leftImage.shape if i == i1: minImageData = tmpData * 1.0 maxImageData = tmpData * 1.0 minImage[:,:] = i1 maxImage[:,:] = i1 else: tmpIndex = numpy.where(tmpData < minImageData) minImage[tmpIndex] = i minImageData[tmpIndex] = tmpData[tmpIndex] tmpIndex = numpy.where(tmpData > maxImageData) maxImage[tmpIndex] = i maxImageData[tmpIndex] = tmpData[tmpIndex] minImage = energy[minImage] maxImage = energy[maxImage] isUsingSuppliedEnergyAxis = True background = 0.5 * (i2 - i1) * (leftImage + rightImage) roiImage = numpy.sum(dataImage, axis=0, dtype=numpy.float64) logger.debug("Case 3 ROI image calculation elapsed = %f ", time.time() - t0) else: shape = self._stack.data.shape roiImage = numpy.zeros(self._stackImageData.shape, numpy.float64) background = roiImage * 1 leftImage = roiImage * 1 middleImage = roiImage * 1 rightImage = roiImage * 1 maxImage = numpy.zeros(roiImage.shape, numpy.int32) minImage = numpy.zeros(roiImage.shape, numpy.int32) istep = 1 for i in range(i1, i2): tmpData = self._stack.data[i:i + istep] tmpData.shape = roiImage.shape if i == i1: minImageData = tmpData * 1.0 maxImageData = tmpData * 1.0 minImage[:,:] = i1 maxImage[:,:] = i1 else: tmpIndex = numpy.where(tmpData < minImageData) minImage[tmpIndex] = i minImageData[tmpIndex] = tmpData[tmpIndex] tmpIndex = numpy.where(tmpData > maxImageData) maxImage[tmpIndex] = i maxImageData[tmpIndex] = tmpData[tmpIndex] numpy.add(roiImage, tmpData, roiImage) if (i == i1): leftImage = tmpData elif (i == imiddle): middleImage = tmpData elif i == (i2 - 1): rightImage = tmpData # the used approach is twice slower than argmax, but it # requires much less memory isUsingSuppliedEnergyAxis = True minImage = energy[minImage] maxImage = energy[maxImage] if i2 > i1: background = (leftImage + rightImage) * 0.5 * (i2 - i1) logger.debug("Case 4 Dynamic ROI elapsed = %f", time.time() - t0) else: t0 = time.time() if self._tryNumpy and\ isinstance(self._stack.data, numpy.ndarray): leftImage = self._stack.data[:, :, i1] middleImage = self._stack.data[:, :, imiddle] rightImage = self._stack.data[:, :, i2 - 1] dataImage = self._stack.data[:, :, i1:i2] background = 0.5 * (i2 - i1) * (leftImage + rightImage) roiImage = numpy.sum(dataImage, axis=2, dtype=numpy.float64) maxImage = energy[numpy.argmax(dataImage, axis=2) + i1] minImage = energy[numpy.argmin(dataImage, axis=2) + i1] isUsingSuppliedEnergyAxis = True logger.debug("Case 5 ROI Image elapsed = %f", time.time() - t0) else: shape = self._stack.data.shape roiImage = numpy.zeros(self._stackImageData.shape, numpy.float64) background = roiImage * 1 leftImage = roiImage * 1 middleImage = roiImage * 1 rightImage = roiImage * 1 maxImage = roiImage * 1 minImage = roiImage * 1 step = 1 for i in range(shape[0]): tmpData = self._stack.data[i:i+step,:, i1:i2] * 1 numpy.add(roiImage[i:i+step,:], numpy.sum(tmpData, axis=2, dtype=numpy.float64), roiImage[i:i+step,:]) numpy.add(minImage[i:i+step,:], numpy.min(tmpData, 2), minImage[i:i+step,:]) numpy.add(maxImage[i:i+step,:], numpy.max(tmpData, 2), maxImage[i:i+step,:]) leftImage[i:i+step, :] += tmpData[:, :, 0] middleImage[i:i+step, :] += tmpData[:, :, imiddle-i1] rightImage[i:i+step, :] += tmpData[:, :,-1] background = 0.5*(i2-i1)*(leftImage+rightImage) logger.debug("Case 6 Dynamic ROI image calculation elapsed = %f", time.time() - t0) else: #self.fileIndex = 2 t0 = time.time() if self.mcaIndex == 0: leftImage = self._stack.data[i1] middleImage = self._stack.data[imiddle] rightImage = self._stack.data[i2 - 1] background = 0.5 * (i2 - i1) * (leftImage + rightImage) dataImage = self._stack.data[i1:i2] roiImage = numpy.sum(dataImage, axis=0, dtype=numpy.float64) minImage = energy[numpy.argmin(dataImage, axis=0) + i1] maxImage = energy[numpy.argmax(dataImage, axis=0) + i1] isUsingSuppliedEnergyAxis = True logger.debug("Case 7 Default ROI image calculation elapsed = %f", time.time() - t0) else: leftImage = self._stack.data[:, i1, :] middleImage = self._stack.data[:, imiddle, :] rightImage = self._stack.data[:, i2 - 1, :] background = 0.5 * (i2 - i1) * (leftImage + rightImage) dataImage = self._stack.data[:, i1:i2, :] roiImage = numpy.sum(dataImage, axis=1, dtype=numpy.float64) minImage = energy[numpy.argmin(dataImage, axis=1) + i1] maxImage = energy[numpy.argmax(dataImage, axis=1) + i1] isUsingSuppliedEnergyAxis = True logger.debug("Case 8 Default ROI image calculation elapsed = %f", time.time() - t0) imageDict = {'ROI': roiImage, 'Maximum': maxImage, 'Minimum': minImage, 'Left': leftImage, 'Middle': middleImage, 'Right': rightImage, 'Background': background} self.__ROIImageCalculationIsUsingSuppliedEnergyAxis = isUsingSuppliedEnergyAxis logger.debug("ROI images calculated") return imageDict def setSelectionMask(self, mask): logger.debug("setSelectionMask called") goodData = numpy.isfinite(self._mcaData0.y[0].sum()) if goodData: self._selectionMask = mask else: if (self._ROIImageDict["ROI"] is not None) and\ (self.mcaIndex != 0): self._selectionMask = mask * numpy.isfinite(self._ROIImageDict["ROI"]) else: self._selectionMask = mask * numpy.isfinite(self._stackImageData) for key in self.pluginInstanceDict.keys(): self.pluginInstanceDict[key].selectionMaskUpdated() def getSelectionMask(self): logger.debug("getSelectionMask called") return self._selectionMask def addImage(self, image, name, info=None, replace=False, replot=True): """ Add image data to the RGB correlator """ print("Add image data not implemented") def removeImage(self, name, replace=True): """ Remove image data from the RGB correlator """ print("Remove image data not implemented") def addCurve(self, x, y, legend=None, info=None, replace=False, replot=True): """ Add the 1D curve given by x an y to the graph. """ print("addCurve not implemented") return None def removeCurve(self, legend, replot=True): """ Remove the curve associated to the supplied legend from the graph. The graph will be updated if replot is true. """ print("removeCurve not implemented") return None def getGraphXLabel(self): print("getGraphXLabel not implemented") return None def getGraphYLabel(self): print("getGraphYLabel not implemented") return None def getActiveCurve(self): """ Function to access the currently active curve. It returns None in case of not having an active curve. Default output has the form: xvalues, yvalues, legend, dict where dict is a dictionary containing curve info. For the time being, only the plot labels associated to the curve are warranted to be present under the keys xlabel, ylabel. If just_legend is True: The legend of the active curve (or None) is returned. """ logger.debug("getActiveCurve default implementation") info = {} info['xlabel'] = 'Channel' info['ylabel'] = 'Counts' legend = 'ICR Spectrum' return self._mcaData0.x[0], self._mcaData0.y[0], legend, info def getGraphXLimits(self): logger.debug("getGraphXLimits default implementation") return self._mcaData0.x[0].min(), self._mcaData0.x[0].max() def getGraphYLimits(self): logger.debug("getGraphYLimits default implementation") return self._mcaData0.y[0].min(), self._mcaData0.y[0].max() def getStackDataObject(self): return self._stack def getStackData(self): return self._stack.data def getStackInfo(self): return self._stack.info def setPositioners(self, positioners, copy=True): """Updates the "positioners" field in the stack info, after checking that the provided positioners have the proper number of values. :param dict positioners: Dictionary of positioners. The keys are the motor names. The values should preferably be arrays with the same number of values as there are stack pixels. Scalars are accepted if the positioner has a constant value. :param bool copy: If *True* (default), store a copy of the data. If *False*, store the original data whenever it is possible. If the dictionary contains lists, they need to be converted to numpy arrays and a copy is mandatory. :raise: TypeError if positioners is not a dict or if any positioner is not a scalar, list or numpy array. :raise: RuntimeError if any positioner is a list and copy=False """ if not hasattr(positioners, "items"): raise TypeError("Dictionary expected for positioners") npixels = self.getStackOriginalImage().size stackPositioners = {} for motorName, motorValues in positioners.items(): if numpy.isscalar(motorValues) or \ (hasattr(motorValues, "ndim") and motorValues.ndim == 0): stackPositioners[motorName] = motorValues elif hasattr(motorValues, "size"): # numpy array numMotorValues = motorValues.size if numMotorValues == npixels: stackPositioners[motorName] = numpy.array(motorValues, copy=copy) elif hasattr(motorValues, "__len__") and numpy.isscalar(motorValues[0]): # list: convert to numpy array before storing in info numMotorValues = len(motorValues) if numMotorValues == npixels: stackPositioners[motorName] = numpy.array(motorValues) else: raise TypeError( "Wrong type for positioner %s. " % motorName + "Valid types are numpy arrays, scalars or list of scalars.") if len(stackPositioners) != len(positioners): ignored_motors = list(set(positioners.keys()) - set(stackPositioners.keys())) logger.debug("Ignored motors due to mismatch in number of values: %s", ', '.join(ignored_motors)) self._stack.info["positioners"] = stackPositioners def _clearPositioners(self): """Removes the "positioners" field in the stack info""" if "positioners" in self._stack.info: self._stack.info["positioners"] = {} def getPositionersFromIndex(self, index): """Return the value of all positioners for the spectrum identified by its 1D index. If positioners are stored as 1D arrays ``a``, the value returned is simply ``a[index]``. If positioners are stored as 2D arrays, the index is applied to the flattened array. :param int index: Index of spectrum for which the the positioners value is to be returned. :return: dictionary of positioners values whose keys are the motor name :rtype: dict """ positioners = self._stack.info.get("positioners", {}) positionersAtIdx = {} if index >= self.getStackOriginalImage().size or index < 0: raise IndexError("index out of bounds") for motorName, motorValues in positioners.items(): # assert numpy.isscalar(motorValues) or hasattr(motorValues, "shape") if numpy.isscalar(motorValues): positionersAtIdx[motorName] = motorValues elif len(motorValues.shape) == 1: positionersAtIdx[motorName] = motorValues[index] else: positionersAtIdx[motorName] = motorValues.reshape((-1,))[index] return positionersAtIdx def test(): #create a dummy stack nrows = 100 ncols = 200 nchannels = 1024 a = numpy.ones((nrows, ncols), numpy.float64) stackData = numpy.zeros((nrows, ncols, nchannels), numpy.float64) for i in range(nchannels): stackData[:, :, i] = a * i stack = StackBase() stack.setStack(stackData, mcaindex=2) print("This should be 0 = %f" % stack.calculateROIImages(0, 0)['ROI'].sum()) print("This should be 0 = %f" % stack.calculateROIImages(0, 1)['ROI'].sum()) print("%f should be = %f" %\ (stackData[:, :, 0:10].sum(), stack.calculateROIImages(0, 10)['ROI'].sum())) if __name__ == "__main__": test()