qgisSpectre.py

# -*- coding: utf-8 -*-
from os import sys
sys.path.append("/usr/lib/python3/dist-packages/")
"""
/***************************************************************************
 qgisSpectre
                                 A QGIS plugin


View spectra stored in a geodataset. The spectral data must be stored in an array, e.g. as an postgres array datafield or as a comma separated string. When features in the dataset are selected, the integrated spectra over these features will be displayed. 


 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2019-10-07
        git sha              : $Format:%H$
        copyright            : (C) 2019 by Morten Sickel
        email                : morten@sickel.net
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from qgis.PyQt.QtCore import QSettings, QTranslator, QCoreApplication, Qt, QSize, QRectF
import qgis.PyQt.QtCore
from qgis.PyQt.QtGui import QIcon, QImage, QPainter
from qgis.PyQt.QtWidgets import QAction,QGraphicsScene,QApplication,QGraphicsView,QCheckBox, QFileDialog, QTableWidgetItem, QHeaderView
from PyQt5.QtGui import QIcon
from PyQt5.QtGui import QColor
# Initialize Qt resources from file resources.py
from .resources import *
from operator import add # To add spectra

from PyQt5 import QtCore,QtGui
from qgis.core import QgsProject, Qgis, QgsMapLayerType, QgsMapLayer,QgsMapLayerProxyModel,QgsFieldProxyModel,QgsSettings
from qgis.PyQt.QtGui import QPen, QBrush
# Import the code for the DockWidget
from .qgisSpectre_dockwidget import qgisSpectreDockWidget
import os.path
import math 
import json
import yaml
import numpy as np

class qgisSpectre:
    """QGIS Plugin Implementation."""

    def __init__(self, iface):
        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # Save reference to the QGIS interface
        self.iface = iface
        self.dlg=qgisSpectreDockWidget(self.iface.mainWindow())
        # stops it from showing up when starting QGIS
        try:    
            self.dlg.close()
        except:
            print("Should not be open here")
            # It should not open when QGIS is started
        self.pluginIsActive = None
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'qgisSpectre_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)
            QCoreApplication.installTranslator(self.translator)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr(u'&Spectral data')
        self.toolbar = self.iface.addToolBar(u'Spectre viewer')
        self.toolbar.setObjectName(u'Spectre viewer')
        self.pluginname="mortensickel_Spectrumviewer"

        self.view = MouseReadGraphicsView(self.iface)
            
        

    # noinspection PyMethodMayBeStatic
    def tr(self, message):
        """Get the translation for a string using Qt translation API.

        We implement this ourselves since we do not inherit QObject.

        :param message: String for translation.
        :type message: str, QString

        :returns: Translated version of message.
        :rtype: QString
        """
        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('qgisSpectre', message)


    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=False,
        status_tip = None,
        whats_this = None,
        parent=None):
        """Add a toolbar icon to the toolbar.

        :param icon_path: Path to the icon for this action. Can be a resource
            path (e.g. ':/plugins/foo/bar.png') or a normal file system path.
        :type icon_path: str

        :param text: Text that should be shown in menu items for this action.
        :type text: str

        :param callback: Function to be called when the action is triggered.
        :type callback: function

        :param enabled_flag: A flag indicating if the action should be enabled
            by default. Defaults to True.
        :type enabled_flag: bool

        :param add_to_menu: Flag indicating whether the action should also
            be added to the menu. Defaults to True.
        :type add_to_menu: bool

        :param add_to_toolbar: Flag indicating whether the action should also
            be added to the toolbar. Defaults to True.
        :type add_to_toolbar: bool

        :param status_tip: Optional text to show in a popup when mouse pointer
            hovers over the action.
        :type status_tip: str

        :param parent: Parent widget for the new action. Defaults None.
        :type parent: QWidget

        :param whats_this: Optional text to show in the status bar when the
            mouse pointer hovers over the action.

        :returns: The action that was created. Note that the action is also
            added to self.actions list.
        :rtype: QAction
        """

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            self.toolbar.addAction(action)

        if add_to_menu:
            self.iface.addPluginToVectorMenu(
                self.menu,
                action)

        self.actions.append(action)

        return action


    def initGui(self):
        """Create the menu entries and toolbar icons inside the QGIS GUI."""
        self.pluginIsActive = False
        icon_path = ':/plugins/qgisSpectre/icon.png'
        self.add_action(
            icon_path,
            text=self.tr(u'View Spectra'),
            callback=self.run,
            parent=self.iface.mainWindow())
        # Moved to __init__
        # self.dlg=qgisSpectreDockWidget(self.iface.mainWindow())
        self.view.setParent(self.dlg) 
        self.dlg.hlMain.addWidget(self.view)
        #self.dlg.cbLayer.currentIndexChanged['QString'].connect(self.listfields)
        self.dlg.qgLayer.setFilters(QgsMapLayerProxyModel.VectorLayer)
        self.dlg.qgField.setLayer(self.dlg.qgLayer.currentLayer())
       # self.dlg.qgField.setFilters(QgsFieldProxyModel.Numeric)
        self.dlg.qgLayer.layerChanged.connect(lambda: self.dlg.qgField.setLayer(self.dlg.qgLayer.currentLayer()))   
    #--------------------------------------------------------------------------

    def onClosePlugin(self):
        """Cleanup necessary items here when plugin dockwidget is closed"""

        #print "** CLOSING qgisSpectre"

        # disconnects
        self.dlg.closingPlugin.disconnect(self.onClosePlugin)

        # remove this statement if dockwidget is to remain
        # for reuse if plugin is reopened
        # Commented next statement since it causes QGIS crashe
        # when closing the docked window:
        # self.dockwidget = None
        
        self.pluginIsActive = False


    def unload(self):
        """Removes the plugin menu item and icon from QGIS GUI."""

        #print "** UNLOAD qgisSpectre"

        for action in self.actions:
            self.iface.removePluginVectorMenu(
                self.tr(u'Spectral data'),
                action)
            self.iface.removeToolBarIcon(action)
        # remove the toolbar
        del self.toolbar
    
    def drawspectra(self,data=None):
        """ Drawing the spectra on the graphicsstage """
        if data is None:
            spectrepen=QPen(Qt.black)
            layername=self.dlg.qgLayer.currentText()
            fieldname=self.dlg.qgField.currentText()
            layer=self.dlg.qgLayer.currentLayer()
            if layer==None:
                return # Happens some times, just as well to return
            try:
                acalib = float(self.dlg.leA.text())
                bcalib = float(self.dlg.leB.text())
            except ValueError:
                self.iface.messageBar().pushMessage(
                       'Invalid calibration value',
                        level=Qgis.Warning, duration=15)
                return
            self.scene.acalib = acalib
            self.scene.bcalib = bcalib
            self.updateUnit()
            data=self.view.spectreval
            # Prepares background and axis
            self.scene.h=330 # Height of scene
            self.scene.clear()
            self.scene.crdtext=None
            self.scene.markerline=None
            backgroundbrush=QBrush(Qt.white)
            outlinepen=QPen(Qt.white)
            self.scene.bottom=20 # Bottom clearing (for x tick marks and labels)
            self.scene.left=self.scene.bottom # Left clearing (for y tick marks and labels)
            self.scene.top = 30
            left=self.scene.left
            h=self.scene.h
            bt=self.scene.bottom
            # Needs this when saving as png, or background from map will shine through
            self.scene.addRect(0,0,1200,h,outlinepen,backgroundbrush) 
            # Y-axis:
            self.scene.addLine(float(left-1),float(h-bt),float(left-1),10.0) 
            # X-axis:
            self.scene.addLine(float(left-1),float(h-bt-1),float(len(data)+10),float(h-bt-1)) 
            if self.debug:
                self.scene.addLine(1.0,float(h),1000.0,float(h))
                self.scene.addLine(1.0,0.0,1000.0,0.0   )
            acalib=self.scene.acalib
            bcalib=self.scene.bcalib
            maxval=acalib*len(data)+bcalib
            tickval=self.tickinterval
            tickval = round(maxval/10)
            # set up some not too bad tick values. This will for a typical spectre
            # going a bit beyond 3000 keV give a tick distance of 300. Should try
            # to get 250 in stead. - there must be some standard algorithm for this
            # TODO: Fix better tick values
            oom=10**math.floor(math.log10(tickval))
            tickval=round(tickval/oom)*oom
            tickdist=tickval
            tickTextY = float(h-bt+2)
            while tickval < maxval:
                tickch=(tickval-bcalib)/acalib+left
                self.scene.addLine(float(tickch),float(h-bt),float(tickch),float(h-bt+4)) 
                # Ticklines
                text=self.scene.addText(str(tickval))
                text.setPos(tickch+left-40, tickTextY)
                tickval+=tickdist
        # Setting unit for ticks further down, after the size of the spectra is found    
        else:
            spectrepen=QPen(Qt.red)
        logscale=self.dlg.cbLog.isChecked()
        if logscale:
            dataset=[]
            for ch in data:
                if ch < self.minvalue and ch > 0:
                    ch = self.minvalue
                if ch==0:
                    ch=self.logoffset
                dataset.append(math.log(ch)-math.log(self.logoffset))
                
        else:
            dataset=data
        
        #DONE: Add x and y axis
        #DONE: Add scale factors to scale x axis from channel number to keV
        #DONE: Add settings to have custom unit
        #TODO: Custom scales
        #DONE: Keep spectra to compare - i.e. paste spectra
        #DONE: Draw spectra as line, not "line-histogram"
        #TODO: Select different drawing styles
        #DONE: Save as file 
        #DONE: export data to clipboard
        #TODO: export image to clipboard
        #DONE: Paste in a spectre copied spectre (i.e. commaseparated list) to show a second spectre 
        #DONE: Peak detection
        #TODO: Save different set of calibration values
        
        # Scales the spectra to fit with the size of the graphicsview
        
        bt = self.scene.bottom
        top = self.scene.top
        h = self.scene.h-self.scene.bottom
        fact = 1.0
        try:
            fact=(h-bt-top)/max(dataset)
        except ZeroDivisionError:
            self.iface.messageBar().pushMessage("Data Loader", f"No valid data in spectre'", level=Qgis.Critical)
            return
        prevvalue=0
        ch=self.scene.left
        for chvalue in dataset:
            # TODO: User selectable type of plot
       #     self.scene.addLine(float(n),float(h-bt),float(n),(h-bt-fact*ch))
       #     self.scene.addLine(float(n),float(h-(bt+4)-fact*ch),float(n),(h-bt-fact*ch))
            self.scene.addLine(float(ch),float(h-bt-fact*prevvalue),float(ch+1),(h-bt-fact*chvalue),spectrepen)
            prevvalue=chvalue
            ch+=1
        self.scene.end=ch-1
        text=self.scene.addText(self.scene.unit)
        text.setPos(self.scene.end+1,tickTextY)
        s = QgsSettings()
        layername=self.dlg.qgLayer.currentText()
        fieldname=self.dlg.qgField.currentText()
        # This is already read in from the fields
        #acalib=s.value(self.pluginname+"/"+layername+"_"+fieldname+"_a",s.value(self.pluginname+"/defaulta", 1))
        #bcalib=s.value(self.pluginname+"/"+layername+"_"+fieldname+"_b",s.value(self.pluginname+"/defaultb", 0))
        #self.scene.unit=s.value(self.pluginname+"/"+layername+"_"+fieldname+"_unit",s.value(self.pluginname+"/defaultunit", 0))
        
        ntext=self.scene.addText(f"n = {self.view.n}")
        ntext.setPos(left + 2,1)
        if self.dlg.cBautodetect.isChecked():
            self.detectpeaks(data)
            
    def smoothsum(self,s,m):
        smoothed=[]
        lastch=len(s)
        for ch in range(len(s)):
            tot=0
            for i in range(-1*m,m+1):
            #    print(ch,i,ch+i,tot)
                try:
                    if ch+i < 0 or s[ch+i] is None:
                        tot=None
                        break
                    else:
                        tot += s[ch+i]
                except IndexError:
                    tot=None
                    break
            #print(f"------{tot}")
            smoothed.append(tot)
        return(smoothed)
    
    def peak_finder(self,x,y,window,treshold):
        # Using Mariscotti’s second difference method
        # Mariscotti, M. A method for automatic identification of peaks in the presence of background and its application 
        # to spectrum analysis. Nucl. Instrum. Methods 1967, 50, 309–320. 
        if window %2 == 0:
            window += 1
            # Need window to be an uneven number to calculate values for centroid
        s=[] # 2nd differential
        f=[] # sd of 2nd differetial
        lastch=len(y)-1
        for ch in range(len(y)):
            if ch > 0 and ch < lastch:
                s.append(y[ch-1]-2*y[ch]+y[ch+1])
                f.append(math.sqrt(y[ch-1]+4*y[ch]+y[ch+1]))
            else:
                s.append(None)
                f.append(None)
        m = window // 2
        # TODO: Window should vary throughout the dataset as the peak width depends on the energy
        for i in range(5):
            # Does a five times sum smoothing
            s=self.smoothsum(s,m)
        # s holds a Mariscotti data set. Where s < 0 there is a peak
        peak=[]
        peakstart=0
        peaks=[]
        peakranges=[]
        # Going through the Mariscotti data set to look for peaks
        for ch in range(len(s)):
            if peak != []:
                # Within a peak
                # A peak may end into the Nulls at the end
                if s[ch] is None or s[ch] >= 0:
                    # Peak has finished
                    # Does peak refinement as in 
                    # Sam Fearn (2022): An Open-Source Iterative python Module for the Automated Identification of Photopeaks in Photon Spectra v2.0. 
                    # https://doi.org/10.5523/bris.n3cm8fnce5ri2k55dlipee3st
                    # Need to find the average value before the peak
                    pre  = sum(y[peakstart-10:peakstart])/10
                    # And after the peak, more complicated, since this may be too close to the end of the spectrum
                    try:
                        post = sum(y[ch:ch+10])/10
                    except TypeError:
                        # TODO: Needs to refine this. This may happen towards the end of the spectrum
                        # DONE: Needs to handle index error 
                        post = 0
                    except IndexError:
                        # Trying to read past end of spectre
                        if ch < len(y)-1:
                            post = sum(y[ch:])/(len(y)-1-ch)
                        else:
                            post = 0
                    # Adjusting for non-flat background through peak
                    a = (post-pre)/(ch-peakstart)
                    peakvalues=y[peakstart:ch]
                    peakadj = []
                    for i in range(len(peakvalues)):
                        val=peakvalues[i]-(i*a+pre)
                        peakadj.append(val)
                    coefs = [-2,3,6,7,6,3,-2]
                    peakadj = self.savgolsmooth(peakadj,coefs)
                    maxval = max(peakadj)
                    maxch = peakadj.index(maxval)+peakstart+1
                    # TODO: Fit a gaussian to peakadj to find peak channel
                    peaks.append(maxch)
                    peakranges.append([peakstart,ch-1])
                    peak = []
            if s[ch] is not None and s[ch] < 0:
                # This may start a new peak, if so peakstart must be stored,
                # if it is an ongoing peak, just append the channel value
                if peak == []:
                    peakstart=ch
                    
                peak.append(s[ch])
        return(list(zip(peaks, peakranges)))
    
    def savgolsmooth(self,data,coefs):
        # Savitzky Golay smoothing
        win = len(coefs)
        halfwin = int((win-1)/2)
        buff=[0] * halfwin
        data = buff + data + buff
        smoothed = []
        fact= sum(coefs)
        for i in range(halfwin,len(data)-halfwin):
            sgsum=0
            for j in range(-1*halfwin,halfwin+1):
                sgsum += data[i+j] * coefs[j+halfwin]
            smoothed.append(sgsum/fact)
        return(smoothed)
    
        
    def detectpeaks(self,spectre=None):
        # DONE: Another color for marker
        # DONE: Recalculate peak with correct baseline
        # DONE: Calculate peaks on smoothed spectrum
        # DONE: Find nuclides with correct energy
        # DONE: Print channel# or energy
        # DONE: Remove exsisting markers
        if spectre is None or not spectre:
            spectre=self.view.spectreval
        x=list(range(len(spectre)))
        window = int(self.dlg.leWindow.text())
        treshold = int(self.dlg.leTreshold.text())
        self.peaks=self.peak_finder(x,spectre,window,treshold)
        # DONE: Combine this two removeitem calls - no reason to have two differet lists
        if hasattr(self.scene,'peakdescriptions'):
            try:
                for pl in self.scene.peakdescriptions:
                    if pl.scene==self.scene:
                        self.scene.removeItem(pl)
            except:
                print("peakdescriptions problem")
                # Need some better handling?
        self.scene.peakdescriptions=[]
        bt=self.scene.bottom
        h=self.scene.h
        top = self.scene.top
        maxval=max(spectre)
        if self.dlg.cbLog.isChecked():
            maxval=math.log(maxval)-math.log(self.logoffset)
        fact=(h-bt-top)/maxval
        bluepen = QPen(QBrush(QColor(0,0,255,100)), 2, Qt.DashLine)
        peaktablewidget = self.dlg.tWpeaktable 
        peaktablewidget.setRowCount(len(self.peaks))
        peaktablewidget.setColumnCount(3)
        header = peaktablewidget.horizontalHeader()
        # Making correct widths for table columns
        header.setSectionResizeMode(0, QHeaderView.ResizeToContents)
        header.setSectionResizeMode(1, QHeaderView.ResizeToContents)
        header.setSectionResizeMode(2, QHeaderView.Stretch)
        peaktablewidget.setHorizontalHeaderLabels(["Channel","Energy","Target"])
        line = 0
        for(x,ends) in self.peaks:
            n=spectre[int(x)]
            y=n
            if self.dlg.cbLog.isChecked():
                if n==0:
                    n=self.logoffset
                y=math.log(n)-math.log(self.logoffset)
            ycoord = h-bt-fact*y
            xcoord = float(self.scene.left+x)
            pl=self.scene.addLine(xcoord,ycoord-10,xcoord,ycoord+10,bluepen)
            self.scene.peakdescriptions.append(pl)
            try:
                xval = x*self.scene.acalib+self.scene.bcalib
            except:
                xval = x
            pt = self.scene.addText(str(round(xval,1)))
            pt.setPos(xcoord+1,ycoord-30)
            self.scene.peakdescriptions.append(pt)
            peaktablewidget.setItem(line,0,QTableWidgetItem(str(x)))
            peaktablewidget.setItem(line,1,QTableWidgetItem('{:.1f}'.format(xval)))
            # Markers for debugging, to be turned off in production
            if self.debug:
                for end in ends:
                    xcoord = float(self.scene.left+end)
                    pl=self.scene.addLine(xcoord,ycoord-10,xcoord,ycoord+10,bluepen)
                    self.scene.peakdescriptions.append(pl)    
            line += 1
        self.drawnuclidelines()
        
    def drawnuclidelines(self):
        if hasattr(self.scene,'nuklines'):
            # been here before, may have drawn lines
            try:
                for pt in self.scene.nuklines:
                    if pt.scene == self.scene:
                        self.scene.removeItem(pt)
            except:
                # Need some better handling!
                print('Some problems here, removing nukline')
        self.scene.nuklines=[]
        #yellowpen = QPen(QBrush(QColor(255,255,0,100)), 2, Qt.DashLine)
        linepen = QPen(QBrush(QColor(255,0,0,100)), 1, Qt.DashLine) # red! 
        chaccuracy = float(self.dlg.leAccuracy.text())/100
        for nuc in self.gammas:
            for e in self.gammas[nuc]:
                try:
                    x = round((float(e) - self.scene.bcalib)/self.scene.acalib)
                except:
                    self.iface.messageBar().pushMessage(
                       f'Invalid energy {e} for {nuc}',
                        level=Qgis.Warning, duration=15)
                draw = False
                try:
                    for peak in self.peaks:
                        channel=peak[0]
                        draw = draw or (channel-channel*chaccuracy <= x and channel+channel*chaccuracy >=x)
                    # print(x,peak[0],draw)
                except e as TypeError:
                    # Fails out for some reason
                    self.iface.messageBar().pushMessage(
                       f'Problem plotting {e} for {nuc}',
                        level=Qgis.Warning, duration=15)
                if not draw:
                    continue
                try:
                    xcoord = float(self.scene.left+x)
                    ycoord = 2
                    pl=self.scene.addLine(xcoord,ycoord,xcoord,300 - self.scene.bottom+5 ,linepen)
                    self.scene.nuklines.append(pl)
                    # To keep the overview and be able to remove later o
                    pt = self.scene.addText(nuc)
                    pt.setPos(xcoord+1,1)
                    pt.setTextWidth(50)
                    self.scene.nuklines.append(pt)
                    # As above
                except NameError as e:
                    # Gets occasional 'name scene not defined', 
                    self.iface.messageBar().pushMessage(
                       f'Problem plotting {e} for {nuc}',
                        level=Qgis.Warning, duration=15)

            
    def updatecalib(self):
        """ Prepares newly typed in  calibration values for use """
        self.dlg.cbDefault.setChecked(False)
        layername=self.dlg.qgLayer.currentText()
        fieldname=self.dlg.qgField.currentText()
        try:
            self.scene.acalib=float(self.dlg.leA.text())
            self.scene.bcalib=float(self.dlg.leB.text())
        except ValueError:
            if self.dlg.leA.text()=='-' or self.dlg.leB.text()=='-' or self.dlg.leA.text()=='' or self.dlg.leB.text()=='' :
                pass
            else:
              self.iface.messageBar().pushMessage(
                   "Calibrating", "Invalid value(s)",
                    level=Qgis.Warning, duration=3)
    
    
    def setdefault(self):
        """ Stores actual values as defaults """
        s=QgsSettings()
        if self.dlg.cbDefault.isChecked():
            s.setValue(self.pluginname+"/defaulta",self.scene.acalib)
            s.setValue(self.pluginname+"/defaultb",self.scene.bcalib)
            s.setValue(self.pluginname+"/defaultunit",self.scene.unit)
        
    def usedefault(self):
        """Fetches default values for the plot """
        s=QgsSettings()
        self.scene.bcalib=s.value(self.pluginname+"/defaultb", 0)      
        self.scene.acalib=s.value(self.pluginname+"/defaulta", 1)
        self.scene.unit=s.value(self.pluginname+"/defaultunit", 'Ch')
        self.dlg.leA.setText(str(self.scene.acalib))
        self.dlg.leB.setText(str(self.scene.bcalib))
        self.dlg.leUnit.setText(str(self.scene.unit))

    def prepareplot(self):
        """ Reads in default values for selected layer before plotting"""
        layername=self.dlg.qgLayer.currentText()
        fieldname=self.dlg.qgField.currentText()
        s=QgsSettings()
        self.scene.acalib=float(s.value(self.pluginname+"/"+layername+"_"+fieldname+"_a",s.value(self.pluginname+"/defaulta", 1)))
        self.scene.bcalib=float(s.value(self.pluginname+"/"+layername+"_"+fieldname+"_b",s.value(self.pluginname+"/defaultb", 0)))
        self.scene.unit=s.value(self.pluginname+"/"+layername+"_"+fieldname+"_unit",s.value(self.pluginname+"/defaultunit", 'Ch'))
        
        self.dlg.leA.setText(str(self.scene.acalib))
        self.dlg.leB.setText(str(self.scene.bcalib))
        self.dlg.leUnit.setText(str(self.scene.unit))
        self.findselected()
    
    def findselected(self):
        """ Is being run when points have been selected. Makes a sum spectra from selected points"""
        layer=self.dlg.qgLayer.currentLayer()
        if layer is None:
          return
        sels=layer.selectedFeatures() # The selected features in the active (from this plugin's point of view) layer
        n=len(sels)
        if n>0:
            #self.iface.messageBar().pushMessage(
            #        "Drawing spectra", "Integrated over {} measurements".format(str(n)),
            #        level=Qgis.Success, duration=3)
            fieldname=self.dlg.qgField.currentText()
            # DONE: Rewrite to make it possible to read in a spectra as a string of comma-separated numbers
            if fieldname == '' or fieldname is None:
                return # Invalid fieldname, probably not selected yet
            try:
                stringspec = isinstance(sels[0][fieldname],str)
            except KeyError:
                return # Invalid fieldname, probably set from a dataset with other fields
            stringspec = stringspec and (sels[0][fieldname].find(',') != -1)
            if isinstance(sels[0][fieldname],list) or stringspec:
                # Only draw if a list field is selected
                sumspectre = None
                for sel in sels:
                    spectre=sel[fieldname]
                    if stringspec:
                        vals=spectre.split(',')
                        spectre = list(map(float, vals))
                    del spectre[-1] # To get rid of last channel i.e. cosmic from RSI-spectra
                                    # TODO: customable removal of channels at top and/or bottom
                    if sumspectre is None:
                        sumspectre = spectre
                    else:
                        sumspectre = list( map(add, spectre, sumspectre))
                self.view.spectreval=sumspectre
                self.view.n=n
                self.drawspectra()
            else:
                # This is coming up too often
                self.iface.messageBar().pushMessage(
                    "Warning", "Use an array field or a comma separated string",
                    level=Qgis.Warning, duration=3)
                    
    def spectreToClipboard(self):
        """ Copies the channel values to the clipboard as a comma separated string"""
        clipboard = QApplication.clipboard()
        text=",".join(str(x) for x in self.view.spectreval)
        clipboard.setText(text)
        
        # TODO: Make a graphical copy
    
    def spectreFromClipboard(self):
        clipboard = QApplication.clipboard()
        try:
            text = clipboard.text()
            textspec= text.split(",")
            if len(textspec) > 10:
                spectre=[float(x) for x in textspec]
                self.drawspectra(spectre)
            # Check if the input makes sense,
            # If so, call self.drawspectre(dataset)
        except:
            self.iface.messageBar().pushMessage(
                    "Warning", "Invalid data pasted",
                    level=Qgis.Warning, duration=3)
    def saveCalibration(self):
        """ Saves the calibration data """
        layername=self.dlg.qgLayer.currentText()
        fieldname=self.dlg.qgField.currentText()
        s=QgsSettings()
        s.setValue(self.pluginname+"/"+layername+"_"+fieldname+"_a", self.scene.acalib)
        s.setValue(self.pluginname+"/"+layername+"_"+fieldname+"_b", self.scene.bcalib)
        s.setValue(self.pluginname+"/"+layername+"_"+fieldname+"_unit", self.scene.unit)
    
    def updateUnit(self):
        self.scene.unit=self.dlg.leUnit.text()

    def resolve(self,name, basepath=None):
        if not basepath:
            basepath = os.path.dirname(os.path.realpath(__file__))
        return os.path.join(basepath, name)
    
    def estimateCoef(self,x, y):
        x = np.array(x)
        y = np.array(y)
        # number of observations/points
        n = np.size(x)
     
        # mean of x and y vector
        m_x = np.mean(x)
        m_y = np.mean(y)
     
        # calculating cross-deviation and deviation about x
        SS_xy = np.sum(y*x) - n*m_y*m_x
        SS_xx = np.sum(x*x) - n*m_x*m_x
     
        # calculating regression coefficients
        b_1 = SS_xy / SS_xx
        b_0 = m_y - b_1*m_x
     
        return (b_0, b_1)
    
    
    def calccalibrate(self):
        #Reads in values 
        self.dlg.cbUseCalibration.setChecked(False)
        tablewidget = self.dlg.tWpeaktable
        data = []
        chs = []
        targets = []
        maxtarget = -1
        #print(tablewidget.rowCount())
        for row in range(tablewidget.rowCount()):
            target = tablewidget.item(row,2)
            if target is None:
                continue
            ch = tablewidget.item(row,0).text()
            target = target.text()
            try:
                target = float(target)
                ch = int(ch)
                data.append([ch,target])
                chs.append(ch)
                targets.append(target)
            except:
                if target =='':
                    continue
                self.iface.messageBar().pushMessage(
                    f"'{ch}' or '{target}' is not numeric",
                    level=Qgis.Warning, duration=15)
                continue
            if target < maxtarget:
                self.iface.messageBar().pushMessage(
                    f'Channel {ch}: {target} &lt; {maxtarget}: targetenergy must be increasing',
                    level=Qgis.Warning, duration=15)
                continue
            maxtarget = target
        if len(data) >= 2:    
            (b,a) = self.estimateCoef(chs,targets)
            self.dlg.labA.setText('{:.4f}'.format(a))
            self.dlg.labB.setText('{:.4f}'.format(b))
        else:
            print(data)
            self.iface.messageBar().pushMessage(
                'Too few valid points, cannot calculate calibration',
                level=Qgis.Warning, duration=15)    
        
    
    def run(self):
        print("starting spectreviewer")
        """Run method that loads and starts the plugin"""
        if not self.pluginIsActive:
            self.pluginIsActive = True
            self.scene = QGraphicsScene()
            self.debug = False
            # Storing the spectra to be able to read out values later on
            # Setting the values storing line and text shown when the mouse button is clicked
            self.scene.crdtext = None
            self.scene.markerline = None
            self.scene.left = None
            # TODO: The four next settings to be user-settable
            self.tickinterval = 100
            s = QgsSettings()
            self.scene.acalib = s.value(self.pluginname+"/defaulta", 1)
            self.scene.bcalib = s.value(self.pluginname+"/defaultb", 0)
            self.scene.unit = s.value(self.pluginname+"/defaultunit","Ch")
            self.dlg.leA.setText(str(self.scene.acalib))
            self.dlg.leB.setText(str(self.scene.bcalib))
            self.dlg.leUnit.setText(str(self.scene.unit))
            showch = False # Set to True to show channel values
            if showch:
                self.scene.unit = 'Ch'
                self.scene.acalib = 1
                self.scene.bcalib = 0
            self.view.setScene(self.scene)
            self.scene.setSceneRect(0,0,1200,350)
            self.scene.top = 40
            # Replotting spectre when a new selection is made
            self.iface.mapCanvas().selectionChanged.connect(self.findselected)        
            # Listing layers
            # DONE: Only list vector layers
            # DONE: Repopulate when layers are added or removed
            # DONE both by using qgisWidget
            self.dlg.pBCopy.clicked.connect(self.spectreToClipboard)
            self.dlg.pBPaste.clicked.connect(self.spectreFromClipboard)
            self.dlg.pBUseDefault.clicked.connect(self.usedefault)
            self.dlg.pBSaveCalib.clicked.connect(self.saveCalibration)
            self.dlg.pBSave.clicked.connect(self.view.saveImage)
            self.dlg.pBPeakDetection.clicked.connect(self.detectandfind)
            self.dlg.leUnit.textChanged.connect(self.updateUnit)
            self.dlg.btRefresh.clicked.connect(self.findselected)
            self.dlg.btRefresh_2.clicked.connect(self.findselected)
            self.dlg.pBCalibrate.clicked.connect(self.calccalibrate)
            
            # connect to provide cleanup on closing of dockwidget
            self.dlg.closingPlugin.connect(self.onClosePlugin)
            # show the dockwidget
            self.iface.mainWindow().addDockWidget(Qt.BottomDockWidgetArea, self.dlg)
            self.dlg.show()
            self.dlg.cbLog.stateChanged.connect(self.findselected)
            self.dlg.cbUseCalibration.stateChanged.connect(self.copycalibdata)
            self.dlg.cbDefault.stateChanged.connect(self.setdefault)
            self.dlg.qgField.currentIndexChanged['QString'].connect(self.prepareplot)
            self.dlg.qgLayer.currentIndexChanged['QString'].connect(self.prepareplot)
            self.dlg.leA.textChanged['QString'].connect(self.updatecalib)
            self.dlg.leB.textChanged['QString'].connect(self.updatecalib)
            self.findselected()
            try:
                with open(self.resolve('gammas.yml')) as stream:
                    self.gammas=yaml.safe_load(stream)
            except FileNotFoundError:
                self.gammas={}
            #self.gammaenergies={}
            #for nuk in self.gammas:
            #    for energy in self.gammas[nuk]
            #        self.
            # Need to be able to print 0-values on a log scale.
            # Defines a offset where the 0 values are to be plotted.
            # A minimum value is set, all values below are adjusted up to this
            self.logoffset = 0.45
            self.minvalue = 0.5
            
            
    def copycalibdata(self):        
        if self.dlg.cbUseCalibration.isChecked():
            self.dlg.leA.setText(self.dlg.labA.text())
            self.dlg.leB.setText(self.dlg.labB.text())
            self.findselected()
            
    def detectandfind(self):
        self.detectpeaks()
        self.findselected()
          
class MouseReadGraphicsView(QGraphicsView):
    """ A class based on QGraphicsView to enable capture of mouse events"""
    
    def __init__(self, iface):
        self.iface = iface
        QGraphicsView.__init__(self)
        self.linex=0
    #DONE: Use arrowkeys to move marker up and down in spectra
        
    
    def drawline(self):
        """ Prints a marker line and reads out energy and number of counts"""
        #TODO: Show list of nuclides with peak at actual energy
        #      Maybe as a further extention as this is radionuclide specific.
        scene=self.scene()
        x=self.linex
        ch=x-scene.left
        unit=scene.unit
        energy=ch*scene.acalib+scene.bcalib
        # DONE: draw a vertical line where clicked. Mark energy
        if unit == 'Ch':
            message="{} {} (n={})".format(unit,int(energy),self.spectreval[int(ch)])
        else:
            message="Ch {}, {} {} (n={})".format(int(ch),int(energy),unit,self.spectreval[int(ch)])
        if self.scene().crdtext is not None:
            self.scene().removeItem(self.scene().crdtext)
        if self.scene().markerline is not None:
            self.scene().removeItem(self.scene().markerline)
        self.scene().crdtext=self.scene().addText(message)
        linetop=scene.top-20
        self.scene().crdtext.setPos(x,linetop)
        self.scene().markerline=self.scene().addLine(x,linetop,x,350-(scene.bottom+20))
    
    def keyPressEvent(self,event):
        ### Reads key presses to move marker line """
        #TODO: Use proper key constants
        if event.key()==Qt.Key_Space:
            self.saveImage()
            return
        if event.key()==Qt.Key_Right: #16777236: #right arrowkey
            self.linex+=1
        if event.key()==Qt.Key_Left: #16777234: # left arrowkey
            self.linex-=1
        if event.key()==Qt.Key_Up: #16777235: # up arrow
            self.linex+=10
        if event.key()==Qt.Key_Down: #16777237: # down arrow
            self.linex-=10
        self.linex=max(self.scene().left,self.linex)
        self.linex=min(self.scene().end,self.linex)
        self.drawline()
        if event.key()==Qt.Key_Escape: # To  be set to Esc 
            if self.scene().crdtext is not None:
                self.scene().removeItem(self.scene().crdtext)
            if self.scene().markerline is not None:
                self.scene().removeItem(self.scene().markerline)
        
    def saveImage(self):
        options = QFileDialog.Options()
        options |= QFileDialog.DontUseNativeDialog
        fileName, _ = QFileDialog.getSaveFileName(self,"QFileDialog.getSaveFileName()","","Image files (*.png);;All Files (*)", options=options)
        if not fileName:
            return
        if not fileName.endswith('.png'):
            fileName += '.png'
        # Get region of scene to capture from somewhere.
        area = self.scene().sceneRect()
        # Create a QImage to render to and fix up a QPainter for it.
        image = QImage(area.toRect().size(), QImage.Format_ARGB32_Premultiplied)
        painter = QPainter(image)
        # Render the region of interest to the QImage.
        self.scene().render(painter, QRectF(image.rect()),  # target
                    area)
        painter.end()
        # Save the image to a file.
        image.save(fileName)
        
    def mousePressEvent(self, event):
        """ Press the left mouse button to draw a line and print the energy at the point"""
        
        if event.button() == 1:
            if self.scene() is None or self.scene().left is None: # Not yet initialized
                return
            coords=self.mapToScene(event.pos())    
            x = coords.x()
            self.linex=x
            # Make sure the data not is read out when being outside the spectra
            if x is not None and x > self.scene().left and x < self.scene().end:
                self.drawline()