Python SpecfitFuns.interpol方法代码示例

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
def estimateXANESEdge(spectrum, energy=None, full=False):
    if energy is None:
        x = numpy.arange(len(spectrum)).astype(numpy.float)
    else:
        x = numpy.array(energy, dtype=numpy.float, copy=True)
    y = numpy.array(spectrum, dtype=numpy.float, copy=True)

    # make sure data are sorted
    idx = energy.argsort(kind='mergesort')
    x = numpy.take(energy, idx)
    y = numpy.take(spectrum, idx)

    # make sure data are strictly increasing
    delta = x[1:] - x[:-1]
    dmin = delta.min()
    dmax = delta.max()
    if delta.min() <= 1.0e-10:
        # force data are strictly increasing
        # although we do not consider last point
        idx = numpy.nonzero(delta>0)[0]
        x = numpy.take(x, idx)
        y = numpy.take(y, idx)
        delta = None

    sortedX = x
    sortedY = y

    # use a regularly spaced spectrum
    if dmax != dmin:
        # choose the number of points or deduce it from
        # the input data length?
        nchannels = 2 * len(spectrum)
        xi = numpy.linspace(x[1], x[-2], nchannels).reshape(-1, 1)
        x.shape = -1
        y.shape = -1
        y = SpecfitFuns.interpol([x], y, xi, y.min())
        x = xi
        x.shape = -1
        y.shape = -1

    # take the first derivative
    npoints = 7
    xPrime = x[npoints:-npoints]
    yPrime = SGModule.getSavitzkyGolay(y, npoints=npoints, degree=2, order=1)

    # get the index at maximum value
    iMax = numpy.argmax(yPrime)

    # get the center of mass
    w = 2 * npoints
    selection = yPrime[iMax-w:iMax+w+1]
    edge = (selection * xPrime[iMax-w:iMax+w+1]).sum(dtype=numpy.float)/\
           selection.sum(dtype=numpy.float)

    if full:
        # return intermediate information
        return edge, sortedX, sortedY, xPrime, yPrime
    else:
        # return the corresponding x value
        return edge

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
def estimateXANESEdge(spectrum, energy=None, npoints=5, full=False,
                      sanitize=True):
    if sanitize:
        if energy is None:
            x = numpy.arange(len(spectrum)).astype(numpy.float)
        else:
            x = numpy.array(energy, dtype=numpy.float, copy=False)
        y = numpy.array(spectrum, dtype=numpy.float, copy=False)
        # make sure data are sorted
        idx = energy.argsort(kind='mergesort')
        x = numpy.take(energy, idx)
        y = numpy.take(spectrum, idx)

        # make sure data are strictly increasing
        delta = x[1:] - x[:-1]
        dmin = delta.min()
        dmax = delta.max()
        if delta.min() <= 1.0e-10:
            # force data are strictly increasing
            # although we do not consider last point
            idx = numpy.nonzero(delta>0)[0]
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)
            delta = None

        # use a regularly spaced spectrum
        if dmax != dmin:
            # choose the number of points or deduce it from
            # the input data length?
            nchannels = 10 * x.size
            xi = numpy.linspace(x[1], x[-2], nchannels).reshape(-1, 1)
            x.shape = -1
            y.shape = -1
            y = SpecfitFuns.interpol([x], y, xi, y.min())
            x = xi
    else:
        # take views
        x = energy[:]
        y = spectrum[:]

    x.shape = -1
    y.shape = -1

    # Sorted and regularly spaced values
    sortedX = x
    sortedY = y    

    ddict = getE0SavitzkyGolay(sortedX,
                               sortedY,
                               points=npoints,
                               full=full)
    if full:
        # return intermediate information
        return ddict["edge"], sortedX, sortedY, ddict["xPrime"], ddict["yPrime"]
    else:
        # return the corresponding x value
        return ddict

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
    def interpolate(self, factor=1.):
        """
        Input
        -----
        factor : float
            factor used to oversample existing data, use
            with caution.

        Interpolates all existing curves to an equidistant
        x-range using the either the active or the first
        curve do determine the number of data points.
        Use this method instead of self.getAllCurves() when
        performin FFT related tasks.

        Returns
        -------
        interpCurves : ndarray
            Array containing the interpolated curves shown
            in the plot window.
            Format: [(x0, y0, legend0, info0), ...]
        """
        curves = self.getAllCurves()
        if len(curves) < 1:
            if DEBUG == 1:
                print('interpolate -- no curves present')
            raise ValueError("At least 1 curve needed")
            return

        activeCurve = self.getActiveCurve()
        if not activeCurve:
            activeCurve = curves[0]
        else:
            activeLegend = activeCurve[2]
            idx = list.index([curve[2] for curve in curves],
                             activeLegend)
            activeCurve = curves[idx]
        activeX, activeY, activeLegend, activeInfo = activeCurve[0:4]

        # Determine average spaceing between Datapoints
        step = numpy.average(numpy.diff(activeX))
        xmin, xmax = self.getXLimits([x for (x,y,leg,info) in curves],
                                     overlap=False)
        num  = factor * numpy.ceil((xmax-xmin)/step)

        # Create equidistant x-range, exclude first and last point
        xeq = numpy.linspace(xmin, xmax, num, endpoint=False)[:-1]

        # Interpolate on sections of xeq
        interpCurves = []
        for (x,y,legend,info) in curves:
            idx = numpy.nonzero((x.min()<xeq) & (xeq<x.max()))[0]
            xi = numpy.take(xeq, idx)
            yi = SpecfitFuns.interpol([x], y, xi.reshape(-1,1), y.min())
            yi.shape = -1
            interpCurves += [(xi, yi, legend, info)]
        return interpCurves

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
    def fftAlignment(self):
        curves = self.getMonotonicCurves()
        nCurves = len(curves)
        if nCurves < 2:
            raise ValueError("At least 2 curves needed")
            return

        # get legend of active curve
        try:
            activeCurveLegend = self.getActiveCurve()[2]
            if activeCurveLegend is None:
                activeCurveLegend = curves[0][2]
            for curve in curves:
                if curve[2] == activeCurveLegend:
                    activeCurve = curve
                    break
        except:
            activeCurve = curves[0]
            activeCurveLegend = activeCurve[2]

        # apply between graph limits
        x0, y0 = activeCurve[0:2]
        xmin, xmax =self.getGraphXLimits()

        for curve in curves:
            xmax = float(min(xmax, curve[0][-1]))
            xmin = float(max(xmin, curve[0][0]))

        idx = numpy.nonzero((x0 >= xmin) & (x0 <= xmax))[0]
        x0 = numpy.take(x0, idx)
        y0 = numpy.take(y0, idx)

        #make sure values are regularly spaced
        xi = numpy.linspace(x0[0], x0[-1], len(idx)).reshape(-1, 1)
        yi = SpecfitFuns.interpol([x0], y0, xi, y0.min())
        x0 = xi * 1
        y0 = yi * 1

        y0.shape = -1
        fft0 = numpy.fft.fft(y0)
        y0.shape = -1, 1
        x0.shape = -1, 1
        nChannels = x0.shape[0]

        # built a couple of temporary array of spectra for handy access
        shiftList = []
        curveList = []
        i = 0
        activeCurveIndex = 0
        for idx in range(nCurves):
            x, y, legend, info = curves[idx][0:4]
            if legend == activeCurveLegend:
                needInterpolation = False
                activeCurveIndex = idx
            elif x.size != x0.size:
                needInterpolation = True
            elif numpy.allclose(x, x0):
                # no need for interpolation
                needInterpolation = False
            else:
                needInterpolation = True

            if needInterpolation:
                # we have to interpolate
                x.shape = -1
                y.shape = -1
                xi = x0[:] * 1
                y = SpecfitFuns.interpol([x], y, xi, y0.min())
                x = xi

            y.shape = -1
            i += 1

            # now calculate the shift
            ffty = numpy.fft.fft(y)

            if numpy.allclose(fft0, ffty):
                shiftList.append(0.0)
                x.shape = -1
            else:
                shift = numpy.fft.ifft(fft0 * ffty.conjugate()).real
                shift2 = numpy.zeros(shift.shape, dtype=shift.dtype)
                m = shift2.size//2
                shift2[m:] = shift[:-m]
                shift2[:m] = shift[-m:]
                threshold = 0.80*shift2.max()
                #threshold = shift2.mean()
                idx = numpy.nonzero(shift2 > threshold)[0]
                #print("max indices = %d" % (m - idx))
                shift = (shift2[idx] * idx/shift2[idx].sum()).sum()
                #print("shift = ", shift - m, "in x units = ", (shift - m) * (x[1]-x[0]))

                # shift the curve
                shift = (shift - m) * (x[1]-x[0])
                x.shape = -1
                y = numpy.fft.ifft(numpy.exp(-2.0*numpy.pi*numpy.sqrt(numpy.complex(-1))*\
                                numpy.fft.fftfreq(len(x), d=x[1]-x[0])*shift)*ffty)
                y = y.real
                y.shape = -1
            curveList.append([x, y, legend + "SHIFT", False, False])
#.........这里部分代码省略.........

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
    def divideByActiveCurve(self):
        #all curves
        curves = self.getAllCurves()
        nCurves = len(curves)
        if nCurves < 2:
            raise ValueError("At least two curves needed")
            return

        #get active curve
        activeCurve = self.getActiveCurve()
        if activeCurve is None:
            raise ValueError("Please select an active curve")
            return

        x, y, legend0, info = activeCurve
        xmin, xmax = self.getGraphXLimits()
        y = y.astype(numpy.float)

        #get the nonzero values
        idx = numpy.nonzero(abs(y) != 0.0)[0]
        if not len(idx):
            raise ValueError("All divisor values are zero!")
        x0 = numpy.take(x, idx)
        y0 = numpy.take(y, idx)

        #sort the values
        idx = numpy.argsort(x0, kind='mergesort')
        x0 = numpy.take(x0, idx)
        y0 = numpy.take(y0, idx)

        i = 0
        for curve in curves:
            x, y, legend, info = curve[0:4]
            if legend == legend0:
                continue
            #take the portion ox x between limits
            idx = numpy.nonzero((x>=xmin) & (x<=xmax))[0]
            if not len(idx):
                #no overlap
                continue
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)

            idx = numpy.nonzero((x0>=x.min()) & (x0<=x.max()))[0]
            if not len(idx):
                #no overlap
                continue
            xi = numpy.take(x0, idx)
            yi = numpy.take(y0, idx)

            #perform interpolation
            xi.shape = -1, 1
            yw = SpecfitFuns.interpol([x], y, xi, yi.min())
            y = yw / yi
            if i == 0:
                replace = True
                replot = True
                i = 1
            else:
                replot = False
                replace = False
            self.addCurve(x, y,
                          legend=legend,
                          info=info,
                          replot=replot,
                          replace=replace)
            lastCurve = [x, y, legend]
        self.addCurve(lastCurve[0],
                      lastCurve[1],
                      legend=lastCurve[2],
                      info=info,
                      replot=True,
                      replace=False)

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
def _getROILineProfileCurve(image, roiStart, roiEnd, roiWidth,
                            lineProjectionMode):
    """Returns the profile values and the polygon in the given ROI.

    Works in image coordinates.

    See :func:`getAlignedROIProfileCurve`.
    """
    row0, col0 = roiStart
    row1, col1 = roiEnd

    deltaRow = abs(row1 - row0)
    deltaCol = abs(col1 - col0)

    if (lineProjectionMode == 'X' or
            (lineProjectionMode == 'D' and deltaCol >= deltaRow)):
        nPoints = deltaCol + 1
        coordsRange = col0, col1
    else:  # 'Y' or ('D' and deltaCol < deltaRow)
        nPoints = deltaRow + 1
        coordsRange = row0, row1

    if nPoints == 1:  # all points are the same
        if DEBUG:
            print("START AND END POINT ARE THE SAME!!")
        return None

    # the coordinates of the reference points
    x0 = numpy.arange(image.shape[0], dtype=numpy.float)
    y0 = numpy.arange(image.shape[1], dtype=numpy.float)

    if roiWidth < 1:
        x = numpy.zeros((nPoints, 2), numpy.float)
        x[:, 0] = numpy.linspace(row0, row1, nPoints)
        x[:, 1] = numpy.linspace(col0, col1, nPoints)
        # perform the interpolation
        ydata = SpecfitFuns.interpol((x0, y0), image, x)

        roiPolygonCols = numpy.array((col0, col1), dtype=numpy.float)
        roiPolygonRows = numpy.array((row0, row1), dtype=numpy.float)

    else:
        # find m and b in the line y = mx + b
        m = (row1 - row0) / float((col1 - col0))
        # Not used: b = row0 - m * col0
        alpha = numpy.arctan(m)
        # imagine the following sequence
        # - change origin to the first point
        # - clock-wise rotation to bring the line on the x axis of a new system
        #   so that the points (col0, row0) and (col1, row1)
        #   become (x0, 0) (x1, 0).
        # - counter clock-wise rotation to get the points (x0, -0.5 width),
        #   (x0, 0.5 width), (x1, 0.5 * width) and (x1, -0.5 * width) back to
        #   the original system.
        # - restore the origin to (0, 0)
        # - if those extremes are inside the image the selection is acceptable
        cosalpha = numpy.cos(alpha)
        sinalpha = numpy.sin(alpha)
        newCol0 = 0.0
        newCol1 = (col1 - col0) * cosalpha + (row1 - row0) * sinalpha
        newRow0 = 0.0
        newRow1 = - (col1 - col0) * sinalpha + (row1 - row0) * cosalpha

        if DEBUG:
            print("new X0 Y0 = %f, %f  " % (newCol0, newRow0))
            print("new X1 Y1 = %f, %f  " % (newCol1, newRow1))

        tmpX = numpy.linspace(newCol0, newCol1,
                              nPoints).astype(numpy.float)
        rotMatrix = numpy.zeros((2, 2), dtype=numpy.float)
        rotMatrix[0, 0] = cosalpha
        rotMatrix[0, 1] = - sinalpha
        rotMatrix[1, 0] = sinalpha
        rotMatrix[1, 1] = cosalpha
        if DEBUG:
            # test if I recover the original points
            testX = numpy.zeros((2, 1), numpy.float)
            colRow = numpy.dot(rotMatrix, testX)
            print("Recovered X0 = %f" % (colRow[0, 0] + col0))
            print("Recovered Y0 = %f" % (colRow[1, 0] + row0))
            print("It should be = %f, %f" % (col0, row0))
            testX[0, 0] = newCol1
            testX[1, 0] = newRow1
            colRow = numpy.dot(rotMatrix, testX)
            print("Recovered X1 = %f" % (colRow[0, 0] + col0))
            print("Recovered Y1 = %f" % (colRow[1, 0] + row0))
            print("It should be = %f, %f" % (col1, row1))

        # find the drawing limits
        testX = numpy.zeros((2, 4), numpy.float)
        testX[0, 0] = newCol0
        testX[0, 1] = newCol0
        testX[0, 2] = newCol1
        testX[0, 3] = newCol1
        testX[1, 0] = newRow0 - 0.5 * roiWidth
        testX[1, 1] = newRow0 + 0.5 * roiWidth
        testX[1, 2] = newRow1 + 0.5 * roiWidth
        testX[1, 3] = newRow1 - 0.5 * roiWidth
        colRow = numpy.dot(rotMatrix, testX)
        colLimits0 = colRow[0, :] + col0
#.........这里部分代码省略.........

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
    def XASNormalize(self):
        #all curves
        curves = self.getAllCurves()
        nCurves = len(curves)
        if nCurves < 1:
            raise ValueError("At least one curve needed")
            return

        #get active curve
        activeCurve = self.getActiveCurve()
        if activeCurve is None:
            raise ValueError("Please select an active curve")
            return

        x, y, legend0, info = activeCurve

        #sort the values
        idx = numpy.argsort(x, kind='mergesort')
        x0 = numpy.take(x, idx)
        y0 = numpy.take(y, idx)
        xmin, xmax = self.getGraphXLimits()

        # get calculation parameters
        if self.widget is None:
            self._createWidget(y0, x0)

        parameters = self.parameters
        if parameters['auto_edge']:
            edge = None
        else:
            edge = parameters['edge_energy']
        energy = x
        pre_edge_regions = parameters['pre_edge']['regions']
        post_edge_regions = parameters['post_edge']['regions']
        algorithm ='polynomial'
        algorithm_parameters = {}
        algorithm_parameters['pre_edge_order'] = parameters['pre_edge']\
                                                         ['polynomial']
        algorithm_parameters['post_edge_order'] = parameters['post_edge']\
                                                         ['polynomial']
        i = 0
        lastCurve = None
        for curve in curves:
            x, y, legend, info = curve[0:4]
            #take the portion ox x between limits
            idx = numpy.nonzero((x>=xmin) & (x<=xmax))[0]
            if not len(idx):
                #no overlap
                continue
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)

            idx = numpy.nonzero((x0>=x.min()) & (x0<=x.max()))[0]
            if not len(idx):
                #no overlap
                continue
            xi = numpy.take(x0, idx)
            yi = numpy.take(y0, idx)

            #perform interpolation
            xi.shape = -1, 1
            yw = SpecfitFuns.interpol([x], y, xi, yi.min())

            # try: ... except: here?
            yw.shape = -1
            xi.shape = -1
            x, y = XASNormalization.XASNormalization(yw,
                                energy=xi,
                                edge=edge,
                                pre_edge_regions=pre_edge_regions,
                                post_edge_regions=post_edge_regions,
                                algorithm=algorithm,
                                algorithm_parameters=algorithm_parameters)[0:2]
            #
            if i == 0:
                replace = True
                replot = True
                i = 1
            else:
                replot = False
                replace = False
            newLegend = " ".join(legend.split(" ")[:-1])
            if not newLegend.startswith('Norm.'):
                newLegend = "Norm. " + newLegend
            self.addCurve(x, y,
                          legend=newLegend,
                          info=info,
                          replot=replot,
                          replace=replace)
            lastCurve = [x, y, newLegend]
        self.addCurve(lastCurve[0],
                      lastCurve[1],
                      legend=lastCurve[2],
                      info=info,
                      replot=True,
                      replace=False)

# 需要导入模块: from PyMca5.PyMcaMath.fitting import SpecfitFuns [as 别名]
# 或者: from PyMca5.PyMcaMath.fitting.SpecfitFuns import interpol [as 别名]
    def applyMedianFilter(self, width=3):
        curves = self.getAllCurves()
        nCurves = len(curves)
        if nCurves < width:
            raise ValueError("At least %d curves needed" % width)
            return

        if self.__randomization:
            indices = numpy.random.permutation(nCurves)
        else:
            indices = range(nCurves)

        # get active curve
        activeCurve = self.getActiveCurve()
        if activeCurve is None:
            activeCurve = curves[0]

        # apply between graph limits
        x0 = activeCurve[0][:]
        y0 = activeCurve[1][:]
        xmin, xmax =self.getGraphXLimits()
        idx = numpy.nonzero((x0 >= xmin) & (x0 <= xmax))[0]
        x0 = numpy.take(x0, idx)
        y0 = numpy.take(y0, idx)

        #sort the values
        idx = numpy.argsort(x0, kind='mergesort')
        x0 = numpy.take(x0, idx)
        y0 = numpy.take(y0, idx)

        #remove duplicates
        x0 = x0.ravel()
        idx = numpy.nonzero((x0[1:] > x0[:-1]))[0]
        x0 = numpy.take(x0, idx)
        y0 = numpy.take(y0, idx)

        x0.shape = -1, 1
        nChannels = x0.shape[0]

        # built a couple of temporary array of spectra for handy access
        tmpArray = numpy.zeros((nChannels, nCurves), numpy.float)
        medianSpectra = numpy.zeros((nChannels, nCurves), numpy.float)
        i = 0
        for idx in indices:
            x, y, legend, info = curves[idx][0:4]
            #sort the values
            x = x[:]
            idx = numpy.argsort(x, kind='mergesort')
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)

            #take the portion of x between limits
            idx = numpy.nonzero((x>=xmin) & (x<=xmax))[0]
            if not len(idx):
                # no overlap
                continue
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)

            #remove duplicates
            x = x.ravel()
            idx = numpy.nonzero((x[1:] > x[:-1]))[0]
            x = numpy.take(x, idx)
            y = numpy.take(y, idx)
            x.shape = -1, 1
            if numpy.allclose(x, x0):
                # no need for interpolation
                pass
            else:
                # we have to interpolate
                x.shape = -1
                y.shape = -1
                xi = x0[:]
                y = SpecfitFuns.interpol([x], y, xi, y0.min())
            y.shape = -1
            tmpArray[:, i] = y
            i += 1

        # now perform the median filter
        for i in range(nChannels):
            medianSpectra[i, :] = medfilt1d(tmpArray[i,:],
                                            kernel_size=width)
        tmpArray = None
        # now get the final spectrum
        y = medianSpectra.sum(axis=1) / nCurves
        x0.shape = -1
        y.shape = x0.shape
        legend = "%d Median from %s to %s" % (width,
                                              curves[0][2],
                                              curves[-1][2])
        self.addCurve(x0,
                      y,
                      legend=legend,
                      info=None,
                      replot=True,
                      replace=True)