text: # -*- coding: utf-8 -*-
  r"""
  readmccd module is made for reading data contained in binary image file
  fully or partially.
  It can process a peak or blob search by various methods
  and refine the peak by a gaussian or lorentzian 2D model
  
  More tools can be found in LaueTools package at sourceforge.net and gitlab.esrf.fr
  March 2020
  """
  __author__ = "Jean-Sebastien Micha, CRG-IF BM32 @ ESRF"
  
  # built-in modules
  import sys
  import os
  import copy
  import time as ttt
  import numpy as np
  import configparser as CONF
  
  try:
      from lauetools import fit2Dintensity as fit2d
      from lauetools import fit2Dintensity_Lorentz as fit2d_l
      from lauetools import IOLaueTools as IOLT
      from lauetools import generaltools as GT
      from lauetools import IOimagefile as IOimage
      from lauetools import dict_LaueTools as DictLT
      from lauetools import imageprocessing as ImProc
  except:
      import lauetoolsnn.lauetools.fit2Dintensity as fit2d
      import lauetoolsnn.lauetools.fit2Dintensity_Lorentz as fit2d_l
      import lauetoolsnn.lauetools.IOLaueTools as IOLT
      import lauetoolsnn.lauetools.generaltools as GT
      import lauetoolsnn.lauetools.IOimagefile as IOimage
      import lauetoolsnn.lauetools.dict_LaueTools as DictLT
      import lauetoolsnn.lauetools.imageprocessing as ImProc
      
  listfile = os.listdir(os.curdir)
  
  # Default dictionary peak search parameters:
  
  PEAKSEARCHDICT_Convolve = {"PixelNearRadius": 10,
                          "thresholdConvolve": 500,
                          "IntensityThreshold": 10,
                          "removeedge": 2,
                          "local_maxima_search_method": 2,
                          "boxsize": 15,
                          "position_definition": 1,
                          "verbose": 0,
                          "fit_peaks_gaussian": 1,
                          "xtol": 0.001,
                          "FitPixelDev": 2.0,
                          "return_histo": 0,
                          "write_execution_time": 0,
                          "Data_for_localMaxima": "auto_background",
                          "NumberMaxofFits": 5000}
  # 'MaxPeakSize':3.0,
  # 'MinPeakSize':0.01
  # }
  
  IOimagefuncs = ['Add_Images', 'Add_Images2', 'SumImages', 'getIndex_fromfilename', 'get_imagesize',
  'getpixelValue', 'getwildcardstring',# 'libtiff_ctypes'
   'listfile', 'readCCDimage', 'read_header_marccd',
   'read_header_marccd2', 'read_header_scmos', 'read_motorsposition_fromheader', 'readheader', 'readoneimage', 'readoneimage_band', 'readoneimage_crop', 'readoneimage_crop_fast', 'readoneimage_full', 'readoneimage_manycrops',
   'readrectangle_in_image', 'setfilename', 'stringint', 'write_rawbinary', 'writeimage','getfilename']
  
  imageprocfuncs = ['ConvolvebyKernel', 'LoG', 'LoGArr', 'LocalMaxima_KernelConvolution', 'LocalMaxima_ShiftArrays',
   'LocalMaxima_from_thresholdarray', 'LocalMaxima_ndimage', 'applyformula_on_images', 'blurCCD', 'blurCCD_with_binning',
   'blur_image', 'check_array_indices', 'circularMask', 'compute_autobackground_image', 'computefilteredimage',
   'diff_pix', 'filter_minimum', 'filterimage', 'gauss_kern', 'getExtrema',
   'getIntegratedIntensity', 'getMinMax', 'getindices2cropArray', 'listfile', 'localmaxima',
   'minmax', 'minmax_fast', 'myfromfunction', 'normalize_shape',
   'plot_image_markers', 'pp', 'radialArr', 'rebin2Darray', 'remove_minimum_background',
   'shiftarrays_accum', 'to8bits']
  
  
  thismodule = sys.modules[__name__]
  
  def warn_func_is_now_in(func, modulename):
      def inner(*args, **kwargs):
          print("\n Warning!  %s is no longer in readmccd.py but in %s\n" % (func.__name__, modulename))
          return func(*args, **kwargs)
      return inner
  
  for iofunc in IOimagefuncs:
      setattr(thismodule, iofunc, warn_func_is_now_in(getattr(IOimage, iofunc), 'IOimagefile.py'))
  
  for ipfunc in imageprocfuncs:
      setattr(thismodule, ipfunc, warn_func_is_now_in(getattr(ImProc, ipfunc), 'imageprocessing.py'))
  
  def readoneimage_multiROIfit(filename, centers, boxsize, stackimageindex=-1, CCDLabel="PRINCETON",
                                                                          baseline="auto",
                                                                          startangles=0.0,
                                                                          start_sigma1=1.0,
                                                                          start_sigma2=1.0,
                                                                          position_start="max",
                                                                          fitfunc="gaussian",
                                                                          showfitresults=1,
                                                                          offsetposition=0,
                                                                          verbose=0,
                                                                          xtol=0.00000001,
                                                                          addImax=False,
                                                                          use_data_corrected=None):
      r"""
      Fit several peaks in one image
  
      :param filename: string, full path to image file
      :param centers: list or array like with shape=(n,2)
                      list of centers of selected ROI
      :param boxsize: (Truly HALF boxsize: fuill boxsize= 2(halfboxsize) +1), iterable 2 elements or integer
                      boxsizes [in x, in y] direction or integer to set a square ROI
      :param baseline: string, 'auto' (ie minimum intensity in ROI) or array of floats
      :param startangles: float or iterable of 2 floats, elliptic gaussian angle (major axis with respect to X direction),
                          one value or array of values
      :param start_sigma1, start_sigma2: floats, gaussian standard deviation (major and minor axis) in pixel,
      :param position_start:  string, starting gaussian center:'max' (position of maximum intensity in ROI),
                              "centers" (centre of each ROI)
      :param offsetposition: integer, 0 for no offset, 1  XMAS compatible, since XMAS consider first pixel as index 1 (in array, index starts with 0), 2  fit2d, since fit2d for peaksearch put pixel labelled n at the position n+0.5 (between n and n+1)
      :param use_data_corrected: tuple of 3 elements, Enter data instead of reading data from file:
                           fulldata, framedim, fliprot
                           where fulldata is a 2D ndarray
      :return: list of results:   bkg,  amp  (gaussian height-bkg), X , Y, major axis standard deviation, minor axis standard deviation,
                                  major axis tilt angle / Ox
      .. todo:: setting list of initial guesses can be improve with
          scipy.ndimages of a concatenate array of multiple slices?
      """
      # if 1:#verbose > 0:
      #     print("addImax", addImax)
      
      # read data (several arrays)
      ResData = IOimage.readoneimage_manycrops(filename,
                                      centers,
                                      boxsize,
                                      stackimageindex,
                                      CCDLabel=CCDLabel,
                                      addImax=addImax,
                                      use_data_corrected=use_data_corrected)
      if addImax:
          Datalist, Imax = ResData
      else:
          Datalist = ResData
  
      framedim = DictLT.dict_CCD[CCDLabel][0]
      saturation_value = DictLT.dict_CCD[CCDLabel][2]
     
      Data = np.array(Datalist)
  
      # setting initial guessed values for each center
      nb_Images = len(Data)
      #     print "nb of images to fitdata ... in  readoneimage_multiROIfit()", nb_Images
      if baseline in ("auto", None):  # background height or baseline level
          list_min = []
          for _, dd in enumerate(Data):
              #            print "k, dd.shape", k, dd.shape
              list_min.append(np.amin(dd))
          start_baseline = list_min
      else:  # input numerical value array
          start_baseline = baseline
  
      if isinstance(startangles, (int, float)):
          start_anglerot = startangles * np.ones(nb_Images)
      else:
          start_anglerot = startangles
  
      if isinstance(start_sigma1, (int, float)):
          start_sigma1 = start_sigma1 * np.ones(nb_Images)
      if isinstance(start_sigma2, (int, float)):
          start_sigma2 = start_sigma2 * np.ones(nb_Images)
  
      if isinstance(boxsize, (int, float)):
          xboxsize, yboxsize = int(boxsize), int(boxsize)
      else:
          xboxsize, yboxsize = boxsize
  
      Xboxsize = xboxsize * np.ones(nb_Images)
      Yboxsize = yboxsize * np.ones(nb_Images)
  
      start_j = []
      start_i = []
      start_amplitude = []
  
      if position_start in ("centers", "center"):  # starting position  from input center
          start_j, start_i = yboxsize, xboxsize
          start_amplitude = []
          d = 0
          for dd in Data:
              start_amplitude.append(dd[int(start_j[d]), int(start_i[d])] - start_baseline[d])
              d += 1
      elif position_start == "max":  # starting position  from maximum intensity in dat
  
          d = 0
          for dd in Data:
              start_j.append(np.argmax(dd) // dd.shape[1])
              start_i.append(np.argmax(dd) % dd.shape[1])
              start_amplitude.append(np.amax(dd) - start_baseline[d])
              d += 1
  
      startingparams_zip = np.array([start_baseline, start_amplitude,
                                      start_j, start_i,
                                      start_sigma1, start_sigma2,
                                      start_anglerot])
  
      RES_params = []
      RES_cov = []
      RES_infodict = []
      RES_errmsg = []
  
      if verbose:
          print("startingparams_zip", startingparams_zip.T)
  
      # consider that ROi shape will be constant over all ROIs,
      # so no need to recompute np.indices in gaussfit
      ROIshape = Data[0].shape
      ijindices_array = np.indices(ROIshape)
  
      k_image = 0
      for startingparams in startingparams_zip.T:
          # if (k_image%25) == 0: print "%d/%d"%(k_image,nb_Images)
          if verbose:
              print("startingparams", startingparams)
          if fitfunc == "gaussian":
              if Data[k_image].shape != ROIshape:
                  ijindices_array = None
  
              # ROIdata = Data[k_image]
              #             print 'ROIdata', ROIdata
              #             print 'np.amax(ROIdata)', np.amax(ROIdata)
              #             print 'np.amin (ROIdata)', np.amin(ROIdata)
              #             print 'np.argmax (ROIdata)', np.argmax(ROIdata)
  
              params, cov, infodict, errmsg = fit2d.gaussfit(Data[k_image],
                                                      err=None,
                                                      params=startingparams,
                                                      autoderiv=1,
                                                      return_all=1,
                                                      circle=0,
                                                      rotate=1,
                                                      vheight=1,
                                                      xtol=xtol,
                                                      Acceptable_HighestValue=saturation_value,
                                                      Acceptable_LowestValue=0,
                                                      ijindices_array=ijindices_array)
  
              if cov is not None and verbose:
                  print('params  solution ', params)
                  print('\n\n *****\n covariance matrix    --- ', cov.tolist())
  
          elif fitfunc == "lorentzian":
              params, cov, infodict, errmsg = fit2d_l.lorentzfit(Data[k_image],
                                                                  err=None,
                                                                  params=startingparams,
                                                                  autoderiv=1,
                                                                  return_all=1,
                                                                  circle=0,
                                                                  rotate=1,
                                                                  vheight=1,
                                                                  xtol=xtol)
  
          if showfitresults:
              # print "startingparams"
              # print startingparams
              print("\n *****fitting results ************\n")
              print(params)
              print("background intensity:                        {:.2f}".format(params[0]))
              print("Peak amplitude above background              {:.2f}".format(params[1]))
              print("pixel position (X)                   {:.2f}".format(
                      params[3] - Xboxsize[k_image] + centers[k_image][0])
              )  # WARNING Y and X are exchanged in params !
              print("pixel position (Y)                   {:.2f}".format(
                      params[2] - Yboxsize[k_image] + centers[k_image][1]))
              print("std 1,std 2 (pix)                    ( {:.2f} , {:.2f} )".format(
                      params[4], params[5]))
              print("e=min(std1,std2)/max(std1,std2)              {:.3f}".format(
                      min(params[4], params[5]) / max(params[4], params[5])))
              print("Rotation angle (deg)                 {:.2f}".format(params[6] % 360))
              print("************************************\n")
          bkg_sol, amp_sol, Y_sol, X_sol, std1_sol, std2_sol, ang_sol = params
  
          RES_cov.append(cov)
          RES_infodict.append(infodict)
          RES_errmsg.append(errmsg)
  
          params_sol = np.array([bkg_sol, amp_sol,
                  X_sol - Xboxsize[k_image] + centers[k_image][0],
                  Y_sol - Yboxsize[k_image] + centers[k_image][1],
                  std1_sol,
                  std2_sol,
                  ang_sol])  # now X,Y in safest order
  
          if offsetposition == 1:
              # PATCH: To match peak positions given by XMAS
              # (confusion coming from counting array indices from 0 or 1...)
              # array fitted by python module see pixels at lower position
              params_sol[3] = params_sol[3] + 1.0
              params_sol[2] = params_sol[2] + 1.0
              # End of PATCH
  
          elif offsetposition == 2:  # see Compute_data2thetachi() in find2thetachi.py
              # PATCH: To match peak positions given by peaksearch of fit2D
              # in fit2D graphics window first pixel labelled 1 is for the
              # peaksearch located at position in  between 0 and 1 (ie 0.5)
              params_sol[3] = params_sol[3] + 0.5
              params_sol[2] = (framedim[0] - params_sol[2]) + 0.5  # TODO: tocheck dim[0] or dim[1]
              # End of PATCH
  
          RES_params.append(params_sol)
  
          k_image += 1
  
      if addImax:
          return RES_params, RES_cov, RES_infodict, RES_errmsg, start_baseline, Imax
      else:
          return RES_params, RES_cov, RES_infodict, RES_errmsg, start_baseline
  
  
  def fitPeakMultiROIs(Data, centers, FittingParametersDict, showfitresults=True, verbose=False):
      r""" refine all peaks  guessed to be at center of several ROIs
  
      :param Data: list of Data array centered on peaks
      :param centers: list of pixels (x,y) positions of ROI centers
      :param FittingParametersDict: dict of fitting parameters
  
      :return: RES_params, RES_cov, RES_infodict, RES_errmsg, start_baseline
              which are all list of refinement results
      """
      boxsize = FittingParametersDict["boxsize"]
      if verbose: print("bosize", boxsize)
      if (boxsize[0] % 2 == 0) or (boxsize[1] % 2 == 0):
          print("boxsize", boxsize)
          raise ValueError("boxsizes are not odd !!")
  
      framedim = FittingParametersDict["framedim"]
      saturation_value = FittingParametersDict["saturation_value"]
      baseline = FittingParametersDict["baseline"]
  
      startangles = FittingParametersDict["startangles"]
      position_start = FittingParametersDict["position_start"]
      start_sigma1 = FittingParametersDict["start_sigma1"]
      start_sigma2 = FittingParametersDict["start_sigma2"]
  
      fitfunc = FittingParametersDict["fitfunction"]
      xtol = FittingParametersDict["xtol"]
      offsetposition = FittingParametersDict["offsetposition"]
  
      # setting initial guessed values for each center
      nb_Images = len(Data)
      if verbose:
          print("nb of images fitPeakMultiROIs", nb_Images)
          print("shape of Data", Data.shape)
  
      if baseline in ("auto", None):  # background height or baseline level
          list_min = []
          for _, dd in enumerate(Data):
              #            print "k, dd.shape", k, dd.shape
              list_min.append(np.amin(dd))
          start_baseline = list_min
      else:  # input numerical value array
          start_baseline = baseline
  
      if isinstance(startangles, (int, float)):
          start_anglerot = startangles * np.ones(nb_Images)
      else:
          start_anglerot = startangles
  
      if isinstance(start_sigma1, (int, float)):
          start_sigma1 = start_sigma1 * np.ones(nb_Images)
      if isinstance(start_sigma2, (int, float)):
          start_sigma2 = start_sigma2 * np.ones(nb_Images)
  
      if isinstance(boxsize, (int, float)):
          xboxsize, yboxsize = int(boxsize), int(boxsize)
      else:
          xboxsize, yboxsize = boxsize
  
      halfxboxsize = int(xboxsize / 2.0)
      halfyboxsize = int(yboxsize / 2.0)
  
      Xhalfboxsize = halfxboxsize * np.ones(nb_Images)
      Yhalfboxsize = halfyboxsize * np.ones(nb_Images)
  
      start_j = []
      start_i = []
      start_amplitude = []
  
      if position_start in ("centers", "center"):  # starting position  from input center
          start_j, start_i = halfyboxsize, halfxboxsize
          start_amplitude = []
          d = 0
          for dd in Data:
              start_amplitude.append(dd[int(start_j[d]), int(start_i[d])] - start_baseline[d])
              d += 1
      elif position_start == "max":  # starting position  from maximum intensity in dat
  
          d = 0
          for dd in Data:
              start_j.append(np.argmax(dd) // dd.shape[1])
              start_i.append(np.argmax(dd) % dd.shape[1])
              start_amplitude.append(np.amax(dd) - start_baseline[d])
              d += 1
  
      startingparams_zip = np.array([start_baseline, start_amplitude, start_j, start_i,
                              start_sigma1, start_sigma2, start_anglerot])
  
      RES_params = []
      RES_cov = []
      RES_infodict = []
      RES_errmsg = []
  
      if verbose:
          print("startingparams_zip", startingparams_zip.T)
  
      # consider that ROi shape will be constant over all ROIs,
      # so no need to recompute np.indices in gaussfit
      ROIshape = Data[0].shape
      ijindices_array = np.indices(ROIshape)
  
      k_image = 0
      for startingparams in startingparams_zip.T:
          # if (k_image%25) == 0: print "%d/%d"%(k_image,nb_Images)
          if verbose:
              print("startingparams", startingparams)
          if fitfunc == "gaussian":
              if Data[k_image].shape != ROIshape:
                  ijindices_array = None
  
              params, cov, infodict, errmsg = fit2d.gaussfit(Data[k_image],
                                                              err=None,
                                                              params=startingparams,
                                                              autoderiv=1,
                                                              return_all=1,
                                                              circle=0,
                                                              rotate=1,
                                                              vheight=1,
                                                              xtol=xtol,
                                                              Acceptable_HighestValue=saturation_value,
                                                              Acceptable_LowestValue=0,
                                                              ijindices_array=ijindices_array)
  
          elif fitfunc == "lorentzian":
              params, cov, infodict, errmsg = fit2d_l.lorentzfit(Data[k_image],
                                                                  err=None,
                                                                  params=startingparams,
                                                                  autoderiv=1,
                                                                  return_all=1,
                                                                  circle=0,
                                                                  rotate=1,
                                                                  vheight=1,
                                                                  xtol=xtol)
  
          if showfitresults:
              print("\n *****fitting results ************\n")
              print("  for k_image = ", k_image)
              print(params)
              print("background intensity:                        {:.2f}".format(params[0]))
              print("Peak amplitude above background              {:.2f}".format(params[1]))
              print("pixel position (X)                   {:.2f}".format(
                      params[3] - Xhalfboxsize[k_image] + centers[k_image][0]))  # WARNING Y and X are exchanged
              print("pixel position (Y)                   {:.2f}".format(
                      params[2] - Yhalfboxsize[k_image] + centers[k_image][1]))
              print("std 1,std 2 (pix)                    ( {:.2f} , {:.2f} )".format(
                      params[4], params[5]))
              print("e=min(std1,std2)/max(std1,std2)              {:.3f}".format(
                      min(params[4], params[5]) / max(params[4], params[5])))
              print("Rotation angle (deg)                 {:.2f}".format(params[6] % 360))
              print("- Xboxsize[k_image]", -Xhalfboxsize[k_image])
              print("centers[k_image][0]", centers[k_image][0])
              print("************************************\n")
          bkg_sol, amp_sol, Y_sol, X_sol, std1_sol, std2_sol, ang_sol = params
  
          RES_cov.append(cov)
          RES_infodict.append(infodict)
          RES_errmsg.append(errmsg)
  
          params_sol = np.array([bkg_sol, amp_sol,
                  X_sol - Xhalfboxsize[k_image] + centers[k_image][0],
                  Y_sol - Yhalfboxsize[k_image] + centers[k_image][1],
                  std1_sol,
                  std2_sol,
                  ang_sol])  # now X,Y in safest order
  
          if offsetposition == 1:
              # PATCH: To match peak positions given by XMAS
              # (confusion coming from counting array indices from 0 or 1...)
              # array fitted by python module see pixels at lower position
              params_sol[3] = params_sol[3] + 1.0
              params_sol[2] = params_sol[2] + 1.0
              # End of PATCH
  
          elif offsetposition == 2:  # see Compute_data2thetachi() in find2thetachi.py
              # PATCH: To match peak positions given by peaksearch of fit2D
              # in fit2D graphics window first pixel labelled 1 is for the
              # peaksearch located at position in  between 0 and 1 (ie 0.5)
              params_sol[3] = params_sol[3] + 0.5
              params_sol[2] = (framedim[0] - params_sol[2]) + 0.5  # TODO: tocheck dim[0] or dim[1]
              # End of PATCH
  
          RES_params.append(params_sol)
  
          k_image += 1
  
      return RES_params, RES_cov, RES_infodict, RES_errmsg, start_baseline
  
  
  def getIntegratedIntensities(fullpathimagefile,
                              list_centers,
                              boxsize,
                              CCDLabel="MARCCD165",
                              thresholdlevel=0.2,
                              flipxycenter=True):
      r"""
      read binary image file and compute integrated intensities of peaks
      whose center is given in list_centers
  
      :return: array whose columns are:
          - integrated intensity
          - absolute minimum intensity threshold
          - nb of pixels composing the peak
      """
      dataimage, framedim, _ = IOimage.readCCDimage(fullpathimagefile, CCDLabel, None, 0)
      res = []
      for center in list_centers:
          res.append(ImProc.getIntegratedIntensity(dataimage, center, boxsize, framedim,
                                                          thresholdlevel, flipxycenter))
      return np.array(res)
  
  
  
  def Find_optimal_thresholdconvolveValue(filename, IntensityThreshold, CCDLabel="PRINCETON", verbose=0):
      r"""
      give the lowest value for thresholdconvolve according to IntensityThreshold
      in order not to miss any blob detection in the first filtering
      (by convolution)
      that would be accepted in the second filtering
      (thresholding by IntensityThreshold)
      """
      # image = scipy.misc.lena().astype(float32)
      #    d = readoneimage(filename).reshape((2048,2048))
  
      res = []
      for tc in (0, 100, 200, 500, 750, 1000, 1500, 2000, 2500, 3000, 4000, 5000, 10000):
          # tstart = ttt.time()
          # Isorted, fitpeak, localpeak, nbrawblobs
          _, fitpeak, _, nbrawblobs = PeakSearch(filename,
                                                  CCDLabel=CCDLabel,
                                                  PixelNearRadius=20,
                                                  removeedge=2,
                                                  IntensityThreshold=IntensityThreshold,
                                                  thresholdConvolve=tc,
                                                  boxsize=15,
                                                  position_definition=1,
                                                  verbose=0,
                                                  fit_peaks_gaussian=0,
                                                  xtol=0.001,
                                                  return_histo=2)
  
          # res.append([tc,IntensityThreshold, nbrawblobs,len(fitpeak), ttt.time()-tstart])
          res.append([tc, IntensityThreshold, nbrawblobs, len(fitpeak)])
  
      Res = np.array(res)
      cond = (Res[:, 2] > Res[:, 3])  # if False  ==> threshold after convolution is too high, some blob may have been rejected before intensity thresholding
      indices_False = np.where(cond == False)[0]
      optim_value = 0
      if len(indices_False) > 1:
  
          if indices_False[0] != 0:
              optim_value = Res[indices_False[0] - 1, 0]
  
      if verbose: print("optim value for thresholdConvolve", optim_value)
  
      return optim_value, Res
  
  
  def writepeaklist(tabpeaks, output_filename,
                                          outputfolder=None, comments=None, initialfilename=None):
      r"""
      write peaks properties and comments in file with extension .dat added
      """
      outputfilefullpath = IOLT.writefile_Peaklist(output_filename,
                                                  tabpeaks,
                                                  dirname=outputfolder,
                                                  overwrite=1,
                                                  initialfilename=initialfilename,
                                                  comments=comments)
  
      return outputfilefullpath
  
  
  def fitoneimage_manypeaks(filename, peaklist, boxsize, stackimageindex=-1,
                                                      CCDLabel="PRINCETON",
                                                      dirname=None,
                                                      position_start="max",
                                                      type_of_function="gaussian",
                                                      guessed_peaksize=(1.0, 1.0),
                                                      xtol=0.001,
                                                      FitPixelDev=2.0,
                                                      Ipixmax=None,
                                                      MaxIntensity=100000000000,
                                                      MinIntensity=0,
                                                      PeakSizeRange=(0, 200),
                                                      verbose=0,
                                                      position_definition=1,
                                                      NumberMaxofFits=500,
                                                      ComputeIpixmax=False,
                                                      use_data_corrected=None,
                                                      reject_negative_baseline=True,
                                                      purgeDuplicates=True):
  
      r"""
      fit multiple ROI data to get peaks position in a single image
  
      Ipixmax  :  highest intensity above background in every ROI centered on element of peaklist
  
      use_data_corrected   :  enter data instead of reading data from file
                          must be a tuple of 3 elements:
                          fulldata, framedim, fliprot
                          where fulldata is an ndarray
  
      purgeDuplicates    : True   remove duplicates that are close within pixel distance of 'boxsize' and keep the most intense peak
  
      use_data_corrected   :  enter data instead of reading data from file
                          must be a tuple of 3 elements:
                          fulldata, framedim, fliprot
                          where fulldata  ndarray
  
      .. note:: used in PeakSearchGUI
      """
      #     print 'Ipixmax in fitoneimage_manypeaks', Ipixmax
      if len(peaklist) >= NumberMaxofFits:
          print("TOO MUCH peaks to fitdata.")
          print("(in fitoneimage_manypeaks) It may stuck the computer.")
          print("Try to reduce the number of Local Maxima or reduce NumberMaxofFits in fitoneimage_manypeaks()")
          return
  
      if dirname is not None:
          filename = os.path.join(dirname, filename)
  
      start_sigma1, start_sigma2 = guessed_peaksize
  
      tstart = ttt.time()
  
      ResFit = readoneimage_multiROIfit(filename,
                                          peaklist,
                                          boxsize,
                                          stackimageindex,
                                          CCDLabel=CCDLabel,
                                          baseline="auto",  # min in ROI box
                                          startangles=0.0,
                                          start_sigma1=start_sigma1,
                                          start_sigma2=start_sigma2,
                                          position_start=position_start,  # 'centers' or 'max'
                                          showfitresults=0,
                                          offsetposition=0,  # offset are applied after fit
                                          fitfunc=type_of_function,
                                          xtol=xtol,
                                          addImax=ComputeIpixmax,
                                          use_data_corrected=use_data_corrected)
  
      if verbose:
          print("fitting time for {} peaks is : {:.4f}".format(len(peaklist), ttt.time() - tstart))
          print("nb of results: ", len(ResFit[0]))
  
      if ComputeIpixmax:
          params, _, info, _, baseline, Ipixmax = ResFit
      else:
          params, _, info, _, baseline = ResFit
  
      par = np.array(params)
  
      #    print "par in fitoneimage_manypeaks", par
  
      if par == []:
          print("no fitted peaks")
          return
  
      peak_bkg = par[:, 0]
      peak_I = par[:, 1]
      peak_X = par[:, 2]
      peak_Y = par[:, 3]
      peak_fwaxmaj = par[:, 4]
      peak_fwaxmin = par[:, 5]
      peak_inclination = par[:, 6] % 360
  
      # pixel deviations from guessed initial position before fitting
      Xdev = peak_X - peaklist[:, 0]
      Ydev = peak_Y - peaklist[:, 1]
  
      #     print "peak_X",peak_X
      #     print "peaklist[:, 0]",peaklist[:, 0]
      #     print 'Xdev', Xdev
      #     print "Ydev", Ydev
  
      # --- --- PEAKS REJECTION -------------------------------
      # print "peaklist[:20]",peaklist[:]
      # number of iteration screening
      to_reject = []
      k = 0
      for inf in info:
          if inf["nfev"] > 1550:
              if verbose:
                  print("k= {}   too much iteration".format(k))
              to_reject.append(k)
          k += 1
  
      if CCDLabel == "FRELONID15_corrected":
          reject_negative_baseline = False
      # negative intensity rejection
      if reject_negative_baseline:
          to_reject2 = np.where((peak_bkg - baseline) < 0)[0]
      else:
          to_reject2 = []
  
      # too far found peak rejection
      to_reject3 = np.where(np.sqrt(Xdev ** 2 + Ydev ** 2) > FitPixelDev)[0]
      #     print 'to_reject3', to_reject3
  
      #     print "peak_I",peak_I
  
      # too intense compared to given threshold (saturation)
      to_reject4 = np.where(peak_I >= MaxIntensity)[0]
  
      # too weak compared to given threshold
      to_reject5 = np.where(peak_I <= MinIntensity)[0]
  
      maxpeaksize = np.maximum(peak_fwaxmaj, peak_fwaxmin)
  
      # too small peak compared to given threshold
      to_reject6 = np.where(maxpeaksize <= PeakSizeRange[0])[0]
  
      # too large peak compared to given threshold
      to_reject7 = np.where(maxpeaksize >= PeakSizeRange[1])[0]
  
      if verbose:
          print("to_reject", type(to_reject))
          print("to_reject ...(len)", len(to_reject))
          print(np.take(peaklist, to_reject, axis=0))
          print("to_reject2 ...(len)", len(to_reject2))
          print(np.take(peaklist, to_reject2, axis=0))
          print(np.take(peaklist, to_reject3, axis=0))
  
          print("After fitting, {}/{} peaks have been rejected\n due to (final - initial position)> FitPixelDev = {}".format(
                  len(to_reject3), len(peaklist), FitPixelDev))
          print("{} spots have been rejected\n due to negative baseline".format(len(to_reject2)))
          print("{} spots have been rejected\n due to much intensity ".format(len(to_reject4)))
          print("{} spots have been rejected\n due to weak intensity ".format(len(to_reject5)))
          print("{} spots have been rejected\n due to small peak size".format(len(to_reject6)))
          print("{} spots have been rejected\n due to large peak size".format(len(to_reject7)))
  
      # spots indices to reject
      ToR = (set(to_reject)
          | set(to_reject2)
          | set(to_reject3)
          | set(to_reject4)
          | set(to_reject5)
          | set(to_reject6)
          | set(to_reject7))  # to reject
  
      # spot indices to take
      ToTake = set(np.arange(len(peaklist))) - ToR
  
      if verbose:
          print("index ToTake", ToTake)
          print("nb indices in ToTake", len(ToTake))
      if len(ToTake) < 1:
          return None, par, peaklist
  
      #     print "Ipixmax",Ipixmax
      if Ipixmax is None:
          Ipixmax = peak_I
      else:
          # ask for maximum intensity in ROI, see
          pass
  
      # all peaks list building
      tabpeak = np.array([peak_X, peak_Y, peak_I, peak_fwaxmaj, peak_fwaxmin, peak_inclination,
                          Xdev, Ydev, peak_bkg, Ipixmax]).T
  
      # print("Results of all fits in tabpeak", tabpeak)
  
      tabpeak = np.take(tabpeak, list(ToTake), axis=0)
  
      #     print "tabpeak.shape",tabpeak.shape
      if len(tabpeak.shape) > 1:  # several peaks
          intense_rank = np.argsort(tabpeak[:, 2])[::-1]  # sort by decreasing intensity-bkg
          #    print "intense_rank", intense_rank
  
          tabIsorted = tabpeak[intense_rank]
      #         print "tabIsorted.shape case 1",tabIsorted.shape
  
      else:  # single peak
          #         tabIsorted = np.array(tabpeak)[:,0]
          print("tabIsorted.shape case 2", tabIsorted.shape)
  
      if position_definition == 1:  # XMAS offset
          tabIsorted[:, :2] = tabIsorted[:, :2] + np.array([1, 1])
  
      if verbose>1:
          print("\n\nIntensity sorted\n\n")
          print(tabIsorted[:10])
          print("X,Y", tabIsorted[:10, :2])
      if verbose:
          print("\n{} fitted peak(s)\n".format(len(tabIsorted)))
  
      if purgeDuplicates and len(tabIsorted) > 2:
          if verbose: print("Removing duplicates from fit")
  
          # remove duplicates (close points), the most intense pixel is kept
          # minimum distance fit solutions
          pixeldistance = boxsize
  
          # tabXY, index_todelete
          _, index_todelete = GT.purgeClosePoints2(tabIsorted[:, :2], pixeldistance)
  
          #         print tabXY
          #         print index_todelete
  
          tabIsorted = np.delete(tabIsorted, tuple(index_todelete), axis=0)
          if verbose:
              print(
                  "\n{} peaks found after removing duplicates minimum intermaxima distance = {})".format(
                      len(tabIsorted), pixeldistance))
  
      return tabIsorted, par, peaklist
  
  
  def PeakSearch(filename, stackimageindex=-1, CCDLabel="PRINCETON", center=None,
                                                  boxsizeROI=(200, 200),  # use only if center != None
                                                  PixelNearRadius=5,
                                                  removeedge=2,
                                                  IntensityThreshold=400,
                                                  thresholdConvolve=200,
                                                  paramsHat=(4, 5, 2),
                                                  boxsize=15,
                                                  verbose=0,
                                                  position_definition=1,
                                                  local_maxima_search_method=1,
                                                  peakposition_definition="max",
                                                  fit_peaks_gaussian=1,
                                                  xtol=0.00001,
                                                  return_histo=1,
                                                  FitPixelDev=25,  # to_reject3 parameter
                                                  write_execution_time=1,
                                                  Saturation_value=65535,  # to be merged in CCDLabel
                                                  Saturation_value_flatpeak=65535,
                                                  MinIntensity=0,
                                                  PeakSizeRange=(0, 200),
                                                  Data_for_localMaxima=None,
                                                  Fit_with_Data_for_localMaxima=False,
                                                  Remove_BlackListedPeaks_fromfile=None,
                                                  maxPixelDistanceRejection=15.0,
                                                  NumberMaxofFits=5000,
                                                  reject_negative_baseline=True,
                                                  formulaexpression="A-1.1*B",
                                                  listrois=None,
                                                  outputIpixmax=True):
      r"""
      Find local intensity maxima as starting position for fittinng and return peaklist.
  
      :param filename: string, full path to image data file
  
      :param stackimageindex: integer, index corresponding to the position of image data on a stacked images file
                  if -1  means single image data w/o stacking
  
      :param CCDLabel: string, label for CCD 2D detector used to read the image data file see dict_LaueTools.py
  
      :param center: position
  
      .. todo:: to be removed: position of the ROI center in CCD frame
  
      :param boxsizeROI: dimensions of the ROI to crop the data array
                      only used if center != None
  
      :param boxsize: half length of the selected ROI array centered on each peak, used for:
                  - fitting a peak
                  - estimating the background around a peak
                  - shifting array in second method of local maxima search (shifted arrays)
  
      :param IntensityThreshold: integer, pixel intensity level above which potential peaks are kept for fitting position procedure. For local maxima method 0 and 1, this level is relative to zero intensity. For local maxima method 2, this level is relative to lowest intensity in the ROI (local background).
      
      .. note:: Start with high value, because if too high, few peaks are found (only the most important), and if too low, too many local maxima are found leading to time consuming fitting procedure.
  
      :param thresholdConvolve: integer, pixel intensity level in convolved image above which potential peaks are kept for fitting position procedure. This threshold step on convolved image is applied prior to the local threshold step with IntensityThreshold on initial image (with respect to the local background)
  
      :param paramsHat: mexican hat kernel parameters (see :func:`LocalMaxima_ndimage`)
  
      :param PixelNearRadius: integer, pixel distance between two regions considered as peaks.
      
      .. note:: Start rather with a large value. If too low, there are very much peaks duplicates and this is very time consuming.
  
      :param local_maxima_search_method: integer, Select method for find the local maxima, each of them will fitted
                              - 0   extract all pixel above intensity threshold
                              - 1   find pixels are highest than their neighbours in horizontal, vertica and diagonal direction (up to a given pixel distance)
                              - 2   find local hot pixels which after numerical convolution give high intensity above threshold (thresholdConvolve) then threshold (IntensityThreshold) on raw data
  
      :param peakposition_definition: 'max' or 'center'  for local_maxima_search_method == 2 to assign to the blob position its hottest pixel position or its center (no weight)
  
      :param Saturation_value_flatpeak: saturation value of detector for local maxima search method 1
  
      :param Remove_BlackListedPeaks_fromfile:
          - None
          - file fullpath, str,  to a peaklist file containing peaks that will be deleted in peak list resulting from
          the local maxima search procedure (prior to peak refinement)
          - ndarray of nx2 X Y pixels cooordinates (avoid reading file in peaksearch series)
  
      :param maxPixelDistanceRejection: maximum distance between black listed peaks and current peaks
                                      (found by peak search) to be rejected
  
      :param NumberMaxofFits: highest acceptable number of local maxima peak to be refined with a 2D modelPeakSearch
  
      :param fit_peaks_gaussian:
          - 0  no position and shape refinement procedure performed from local maxima (or blob) result
          - 1  2D gaussian peak refinement
          - 2  2D lorentzian peak refinement
  
      :param xtol: relative error on solution (x vector)  see args for leastsq in scipy.optimize
      :param FitPixelDev: largest pixel distance between initial (from local maxima search)
                              and refined peak position
  
      :param position_definition: due to various conventional habits when reading array, add some offset to fitdata XMAS or fit2d peak search values:
          - 0   no offset (python numpy convention)
          - 1   XMAS offset (first pixel is counted as located at 1 instead of 0)
          - 2   fit2d offset (obsolete)
  
      :param return_histo: - 0   3 output elements
                           - 1   4 elemts, last one is histogram of data
                           - 2   4 elemts, last one is the nb of raw blob found after convolution and threshold
  
      :param Data_for_localMaxima:  object to be used only for initial step of finding local maxima (blobs) search
                                  (and not necessarly for peaks fitting procedure):
                                -  ndarray     = array data
                                - 'auto_background'  = calculate and remove background computed from image data itself (read in file 'filename')
                                - path to image file (string)  = B image to be used in a mathematical operation with Ato current image
  
      :param Fit_with_Data_for_localMaxima: use 'Data_for_localMaxima' object as image when refining peaks position and shape
                                         with initial peak position guess from local maxima search
  
      :param formulaexpression: string containing A (raw data array image) and B (other data array image)
                              expressing mathematical operation,e.g:
                              'A-3.2*B+10000'
                              for simple background substraction (with B as background data):
                              'A-B' or 'A-alpha*B' with alpha > 1.
  
      :param reject_negative_baseline:  True  reject refined peak result if intensity baseline (local background) is negative
                                          (2D model is maybe not suitable)
  
      :param outputIpixmax: compute maximal pixel intensity for all peaks found
  
      :return:  -  peak list sorted by decreasing (integrated intensity - fitted bkg)
                  -peak_X,peak_Y,peak_I,peak_fwaxmaj,peak_fwaxmin,peak_inclination,Xdev,Ydev,peak_bkg
  
      for fit_peaks_gaussian == 0 (no fitdata) and local_maxima_search_method==2 (convolution)
          if peakposition_definition ='max' then X,Y,I are from the hottest pixels
          if peakposition_definition ='center' then X,Y are blob center and I the hottest blob pixel
  
      .. warning:: nb of output elements depends on 'return_histo' argument
      """
  
      if return_histo in (0, 1):
          return_nb_raw_blobs = 0
      if return_histo in (2,):
          return_nb_raw_blobs = 1
      if write_execution_time:
          t0 = ttt.time()
  
      # user input its own shaped Data array
      if isinstance(Data_for_localMaxima, np.ndarray):
          # if verbose:
          #     print("Using 'Data_for_localMaxima' ndarray for finding local maxima")
          Data = Data_for_localMaxima
  
          #         print "min, max intensity", np.amin(Data), np.amax(Data)
          # TODO to test with VHR
          framedim = Data.shape
          ttread = ttt.time()
  
      # Data are read from image file
      elif isinstance(Data_for_localMaxima, str) or Data_for_localMaxima is None:
  
          Data, framedim, fliprot = IOimage.readCCDimage(filename,
                                                  stackimageindex=stackimageindex,
                                                  CCDLabel=CCDLabel,
                                                  dirname=None,
                                                  verbose=verbose)
          
          # if verbose: print("image from filename {} read!".format(filename))
  
          # peak search in a single and particular region of image
          if center is not None:
  
              framediminv = (framedim[1], framedim[0])
              imin, imax, jmin, jmax = ImProc.getindices2cropArray(center, boxsizeROI, framediminv)
              Data = Data[jmin:jmax, imin:imax]
  
          if write_execution_time:
              dtread = ttt.time() - t0
              ttread = ttt.time()
              # if verbose: print("Read Image. Execution time : {:.3f} seconds".format(dtread))
      if isinstance(Data_for_localMaxima, str):
          # if verbose: print("Using Data_for_localMaxima for local maxima search: --->", Data_for_localMaxima)
          # compute and remove background from this image
          if Data_for_localMaxima == "auto_background":
              # if verbose:
              #     print("computing background from current image ", filename)
              backgroundimage = ImProc.compute_autobackground_image(Data, boxsizefilter=10)
              # basic substraction
              usemask = True
          # path to a background image file
          else:
              if stackimageindex == -1:
                  raise ValueError("Use stacked images as background is not implement")
              path_to_bkgfile = Data_for_localMaxima
              if verbose: print("Using image file {} as background".format(path_to_bkgfile))
              try:
                  backgroundimage, _, _ = IOimage.readCCDimage(path_to_bkgfile,
                                                                          CCDLabel=CCDLabel)
              except IOError:
                  raise ValueError("{} does not seem to be a path file ".format(path_to_bkgfile))
  
              usemask = False
  
          # if verbose: print("Removing background for local maxima search")
          Data = ImProc.computefilteredimage(Data, backgroundimage, CCDLabel, usemask=usemask,
                                                              formulaexpression=formulaexpression)
  
      # if verbose > 1: print("Data.shape for local maxima", Data.shape)
  
      # --- PRE SELECTION OF HOT PIXELS as STARTING POINTS FOR FITTING ---------
      # first method ---------- "Basic Intensity Threshold"
      if local_maxima_search_method in (0, "0"):
  
          # if verbose: print("Using simple intensity thresholding to detect local maxima (method 1/3)")
          res = ImProc.LocalMaxima_from_thresholdarray(Data, IntensityThreshold=IntensityThreshold,
                                                      rois=listrois,
                                                      framedim=framedim,
                                                      outputIpixmax=outputIpixmax)
          if res is not None:
              if outputIpixmax:
                  peaklist, Ipixmax = res
                  ComputeIpixmax = True
              else:
                  peaklist = res
                  Ipixmax = np.ones(len(peaklist)) * IntensityThreshold
                  ComputeIpixmax = False
  
      # second method ----------- "Local Maxima in a box by shift array method"
      if local_maxima_search_method in (1, "1"):
          # flat top peaks (e.g. saturation) are NOT well detected
          # if verbose: print("Using shift arrays to detect local maxima (method 2/3)")
          peaklist, Ipixmax = ImProc.LocalMaxima_ShiftArrays(Data,
                                          framedim=framedim,
                                          IntensityThreshold=IntensityThreshold,
                                          Saturation_value=Saturation_value_flatpeak,
                                          boxsize_for_probing_minimal_value_background=boxsize,  # 30
                                          pixeldistance_remove_duplicates=PixelNearRadius,  # 25
                                          nb_of_shift=boxsize)  # 25
  
          ComputeIpixmax = True
  
      # third method: ------------ "Convolution by a gaussian kernel"
      if local_maxima_search_method in (2, "2"):
  
          # if verbose: print("Using mexican hat convolution to detect local maxima (method 3/3)")
  
          peakValConvolve, boxsizeConvolve, central_radiusConvolve = paramsHat
  
          Candidates = ImProc.LocalMaxima_KernelConvolution(Data,
                                          framedim=framedim,
                                          peakValConvolve=peakValConvolve,
                                          boxsizeConvolve=boxsizeConvolve,
                                          central_radiusConvolve=central_radiusConvolve,
                                          thresholdConvolve=thresholdConvolve,  # 600 for CdTe
                                          connectivity=1,
                                          IntensityThreshold=IntensityThreshold,
                                          boxsize_for_probing_minimal_value_background=PixelNearRadius,
                                          return_nb_raw_blobs=return_nb_raw_blobs,
                                          peakposition_definition=peakposition_definition)
  
          if Candidates is None:
              # if verbose: print("No local maxima found, change peak search parameters !!!")
              return None
  
          if return_nb_raw_blobs == 1:
              peaklist, Ipixmax, nbrawblobs = Candidates
          else:
              peaklist, Ipixmax = Candidates
  
          #         print "len(peaklist)", peaklist.shape
          #         print "Ipixmax", Ipixmax.shape
          ComputeIpixmax = True
      #         print "Ipixmax after convolution method", Ipixmax
      # -------------------------------------------------------------
      # --- END of blobs search methods calls
  
      if (peaklist is None
          or peaklist is []
          or peaklist is np.array([])
          or (len(peaklist) == 0)):
          # if verbose: print("No local maxima found, change peak search parameters !!!")
          return None
      # pixel origin correction due to ROI croping
      if center is not None:
          x1, y1 = center  # TODO: to ne checked !!
          peaklist = peaklist + np.array([x1, y1])
  
      if write_execution_time:
          dtsearch = ttt.time() - float(ttread)
  
          # if verbose: print("Local maxima search. Execution time : {:.3f} seconds".format(dtsearch))
  
      # removing some duplicates ------------
      if len(peaklist) >= 2:
          nb_peaks_before = len(peaklist)
          #         print "%d peaks in peaklist before purge" % nb_peaks_before
          #         print 'peaklist',in peaklist before purge
  
          if len(peaklist) >= NumberMaxofFits:
              if verbose:
                  print("TOO MUCH peaks to handle.")
                  print("(in PeakSearch) It may stuck the computer.")
                  print("Try to reduce the number of Local Maxima or\n reduce "
                          "NumberMaxofFits in PeakSearch()")
              return None
  
          Xpeaklist, Ypeaklist, tokeep = GT.removeClosePoints(peaklist[:, 0], peaklist[:, 1],
                                                                                  dist_tolerance=2)
  
          peaklist = np.array([Xpeaklist, Ypeaklist]).T
          Ipixmax = np.take(Ipixmax, tokeep)
  
          if verbose: print("Keep {} from {} initial peaks (ready for peak positions and shape fitting)".format(
                  len(peaklist), nb_peaks_before))
      # -----------------------------------------------
  
      #-------------------------------------------------
      # remove black listed peaks option
      # and update peaklist
      if Remove_BlackListedPeaks_fromfile is not None and len(peaklist)>1:
          if not isinstance(Remove_BlackListedPeaks_fromfile, str):
              # array of XY  shape = (n,2)
              XY_blacklisted = Remove_BlackListedPeaks_fromfile
  
          elif Remove_BlackListedPeaks_fromfile.endswith(('.dat', '.fit')):
              XY_blacklisted = Get_blacklisted_spots(Remove_BlackListedPeaks_fromfile)
          
          if XY_blacklisted is None: print('No or only 1 Blacklisted spots found...')
          else:  #
              X, Y = peaklist[:, :2].T
              (peakX, _, tokeep) = GT.removeClosePoints_two_sets([X, Y], XY_blacklisted,
                                                          dist_tolerance=maxPixelDistanceRejection,
                                                          verbose=0)
  
              npeak_before = len(X)
              npeak_after = len(peakX)
  
              if verbose: print("\n Removed {} (over {}) peaks belonging to the blacklist {}\n".format(
                      npeak_before - npeak_after,
                      npeak_before,
                      Remove_BlackListedPeaks_fromfile))
  
              peaklist = np.take(peaklist, tokeep, axis=0)
              Ipixmax = Ipixmax[tokeep]
      #-------------------------------------------------
  
      # ---- ----------- no FITTING ----------------------------
      # NO FIT  and return raw list of local maxima
      if fit_peaks_gaussian == 0:
  
          if position_definition == 1:  # XMAS like offset
              peaklist[:, :2] = peaklist[:, :2] + np.array([1, 1])
  
          if position_definition == 2:  # fit2D offset
              peaklist[:, 0] = peaklist[:, 0] + 0.5
              peaklist[:, 1] = framedim[0] - peaklist[:, 1] + 0.5
  
          if verbose:
              print("{} local maxima found".format(len(peaklist)))
              print("20 first peaks", peaklist[:20])
  
          # tabpeak mimics the array built after fitting procedures
          tabpeak = np.zeros((len(peaklist[:, 0]), 10))
          tabpeak[:, 0] = peaklist[:, 0]
          tabpeak[:, 1] = peaklist[:, 1]
          tabpeak[:, 2] = Ipixmax
          # return tabpeak, peaklist, peaklist, peaklist  # no fitdata return raw list of local maxima
          lastelem = peaklist
          if return_nb_raw_blobs == 1:
              lastelem = nbrawblobs
  
          return tabpeak, peaklist, peaklist, lastelem
  
      if (peaklist is None
          or peaklist is []
          or peaklist is np.array([])
          or (len(peaklist) == 0)):
          # print("No local maxima found, no peaks to fit !!!")
          return None
  
      # ----  ---------------FITTING ----------------------------
      # gaussian fitdata
      elif fit_peaks_gaussian == 1:
          type_of_function = "gaussian"
  
      # lorentzian fitdata
      elif fit_peaks_gaussian == 2:
          type_of_function = "lorentzian"
  
      else:
          raise ValueError("optional fit_peaks_gaussian value is not understood! Must be 0,1 or 2")
  
      if verbose:
          print("\n*****************")
          print("{} local maxima found".format(len(peaklist)))
          print("\n Fitting of each local maxima\n")
  
      if center is not None:
          position_start = "centers"
      else:
          position_start = "max"
  
      # if Data_for_localMaxima will be used for refining peak positions
      if Fit_with_Data_for_localMaxima:
          Data_to_Fit = (Data, framedim, fliprot)
      else:
          Data_to_Fit = None
  
      return fitoneimage_manypeaks(filename,
                                  peaklist,
                                  boxsize,
                                  stackimageindex,
                                  CCDLabel=CCDLabel,
                                  dirname=None,
                                  position_start=position_start,
                                  type_of_function=type_of_function,
                                  xtol=xtol,
                                  FitPixelDev=FitPixelDev,
                                  Ipixmax=Ipixmax,
                                  MaxIntensity=Saturation_value,
                                  MinIntensity=MinIntensity,
                                  PeakSizeRange=PeakSizeRange,
                                  verbose=verbose,
                                  position_definition=position_definition,
                                  NumberMaxofFits=NumberMaxofFits,
                                  ComputeIpixmax=ComputeIpixmax,
                                  use_data_corrected=Data_to_Fit,
                                  reject_negative_baseline=reject_negative_baseline)
  
  def Get_blacklisted_spots(filename):
      XY_blacklisted = None
      if filename.endswith('.dat'):
  
          data_peak_blacklisted = IOLT.read_Peaklist(filename, dirname=None)
          if len(data_peak_blacklisted) > 1:
  
              XY_blacklisted = data_peak_blacklisted[:, :2].T
  
      elif filename.endswith('.fit'):
          resdata = IOLT.readfile_fit(filename)
          if resdata is not None:
              allgrainsspotsdata = resdata[4]
              nspots = len(allgrainsspotsdata)
          else:
              nspots = 0
  
          if nspots > 1:
  
              XY_blacklisted = allgrainsspotsdata[:, 7:9].T
              
      return XY_blacklisted
  
  def peaksearch_on_Image(filename_in, pspfile, background_flag="no", blacklistpeaklist=None,
                                                          dictPeakSearch={},
                                                          CCDLabel="MARCCD165",
                                                          outputfilename=None,
                                                          psdict_Convolve=PEAKSEARCHDICT_Convolve,
                                                          verbose=0):
      r"""
      Perform a peaksearch by using .psp file
  
      # still not very used and checked?
      # missing dictPeakSearch   as function argument for formulaexpression  or dict_param??
      """
  
      dict_param = readPeakSearchConfigFile(pspfile)
  
      Data_for_localMaxima, formulaexpression = read_background_flag(background_flag)
  
      blacklistedpeaks_file = set_blacklist_filepath(blacklistpeaklist)
  
      dict_param["Data_for_localMaxima"] = Data_for_localMaxima
      dict_param["formulaexpression"] = formulaexpression
      dict_param["Remove_BlackListedPeaks_fromfile"] = blacklistedpeaks_file
  
      # create a data considered as background from an imagefile
      BackgroundImageCreated = False
      flag_for_backgroundremoval = dict_param["Data_for_localMaxima"]
  
      # flag_for_backgroundremoval is a file path to an imagefile
      # create background data: dataimage_bkg
      if flag_for_backgroundremoval not in ("auto_background", None) and not isinstance(
          flag_for_backgroundremoval, np.ndarray):
  
          fullpath_backgroundimage = psdict_Convolve["Data_for_localMaxima"]
  
          #         print "fullpath_backgroundimage ", fullpath_backgroundimage
  
          # dirname_bkg, imagefilename_bkg = os.path.split(fullpath_backgroundimage)
  
          # CCDlabel_bkg = CCDLabel
  
          BackgroundImageCreated = True
  
          if verbose: print("consider dataimagefile {} as background".format(fullpath_backgroundimage))
  
          if "formulaexpression" in dictPeakSearch:
              formulaexpression = dictPeakSearch["formulaexpression"]
          else:
              raise ValueError('Missing "formulaexpression" to operate on images before peaksearch in '                                    'peaksearch_fileseries()')
  
          # saturationlevel = DictLT.dict_CCD[CCDLabel][2]
  
          # dataimage_corrected = applyformula_on_images(dataimage_raw,
          #                                             dataimage_bkg,
          #                                             formulaexpression=formulaexpression,
          #                                             SaturationLevel=saturationlevel,
          #                                             clipintensities=True)
  
          if verbose: print("using {} in peaksearch_fileseries".format(formulaexpression))
  
          # for finding local maxima in image from formula
          psdict_Convolve["Data_for_localMaxima"] = fullpath_backgroundimage
  
          # for fitting peaks in image from formula
          psdict_Convolve["reject_negative_baseline"] = False
          psdict_Convolve["formulaexpression"] = formulaexpression
          psdict_Convolve["Fit_with_Data_for_localMaxima"] = True
  
      Res = PeakSearch(filename_in,
                      CCDLabel=CCDLabel,
                      Saturation_value=DictLT.dict_CCD[CCDLabel][2],
                      Saturation_value_flatpeak=DictLT.dict_CCD[CCDLabel][2],
                      **psdict_Convolve)
  
      if Res in (False, None):
          print("No peak found for image file: ", filename_in)
          return None
      # write file with comments
      Isorted, _, _ = Res[:3]
  
      if outputfilename:
  
          params_comments = "Peak Search and Fit parameters\n"
  
          params_comments += "# {}: {}\n".format("CCDLabel", CCDLabel)
  
          for key, val in list(psdict_Convolve.items()):
              if not BackgroundImageCreated or key not in ("Data_for_localMaxima",):
                  params_comments += "# " + key + " : " + str(val) + "\n"
  
          if BackgroundImageCreated:
              params_comments += ("# "
                                  + "Data_for_localMaxima"
                                  + " : {} \n".format(fullpath_backgroundimage))
          # .dat file extension is done in writefile_Peaklist()
  
          IOLT.writefile_Peaklist("{}".format(outputfilename),
                                  Isorted,
                                  overwrite=1,
                                  initialfilename=filename_in,
                                  comments=params_comments)
  
      return Isorted
  
  
  # -------------------  CONFIG file functions (.psp)
  # import configparser as CONF
  
  # --- ---- Local maxima and fit parameters
  CONVERTKEY_dict = {"fit_peaks_gaussian": "fit_peaks_gaussian",
                      "position_definition": "position_definition",
                      "local_maxima_search_method": "local_maxima_search_method",
                      "intensitythreshold": "IntensityThreshold",
                      "thresholdconvolve": "thresholdConvolve",
                      "boxsize": "boxsize",
                      "pixelnearradius": "PixelNearRadius",
                      "xtol": "xtol",
                      "fitpixeldev": "FitPixelDev",
                      "maxpixeldistancerejection": "maxPixelDistanceRejection",
                      "maxpeaksize": "MaxPeakSize",
                      "minpeaksize": "MinPeakSize"}
  
  LIST_OPTIONS_PEAKSEARCH = ["local_maxima_search_method",
                              "IntensityThreshold",
                              "thresholdConvolve",
                              "boxsize",
                              "PixelNearRadius",
                              "fit_peaks_gaussian",
                              "xtol",
                              "FitPixelDev",
                              "position_definition",
                              "maxPixelDistanceRejection",
                              "MinPeakSize",
                              "MaxPeakSize"]
  
  LIST_OPTIONS_TYPE_PEAKSEARCH = ["integer flag",
                                  "integer count",
                                  "float count",
                                  "integer pixel",
                                  "integer pixel",
                                  "integer flag",
                                  "float",
                                  "float pixel",
                                  "integer flag",
                                  "float",
                                  "float",
                                  "float"]
  
  LIST_OPTIONS_VALUESPARAMS = [1, 1000, 5000, 15, 10, 1, 0.001, 2.0, 1, 15.0, 0.01, 3.0]
  
  if (len(CONVERTKEY_dict) != len(LIST_OPTIONS_PEAKSEARCH)
                              != LIST_OPTIONS_TYPE_PEAKSEARCH != LIST_OPTIONS_VALUESPARAMS):
      raise ValueError(
          "Lists of parameters for config .psp file do not have the same length (readmccd.py)")
  
  
  def savePeakSearchConfigFile(dict_param, outputfilename=None):
      r""" save peak search parameters in .psp file"""
      # save peaksearch parameter in config file
      config = CONF.RawConfigParser()
      config.add_section("PeakSearch")
  
      params_comments = "Peak Search and Fit parameters\n"
  
      for key, val in list(dict_param.items()):
          params_comments += "# " + key + " : " + str(val) + "\n"
          config.set("PeakSearch", key, str(val))
  
      if outputfilename is None:
          outputfilename = "PeakSearch.psp"
  
      if not outputfilename.endswith(".psp"):
          if outputfilename.count(".") > 0:
  
              outputfilename = "".join(outputfilename.split(".")[:-1] + ".psp")
          else:
              outputfilename += ".psp"
  
      # Writing configuration file to 'PeakSearch.cfg'
      with open(outputfilename, "w") as configfile:
          config.write(configfile)
  
      return outputfilename
  
  
  def readPeakSearchConfigFile(filename):
      r""" read peak search parameters in .psp file"""
      config = CONF.RawConfigParser()
      config.optionxform = str
      #    config = MyCasePreservingConfigParser()
  
      config.read(filename)
  
      section = config.sections()[0]
  
      if section not in ("PeakSearch",):
          raise ValueError(
              "wrong section name in config file {}. Must be in {}".format(
                  filename, "IndexRefine"))
  
      dict_param = {}
  
      list_options = config.options(section)
  
      for option in list_options:
  
          #         print "\n option\n", option
          for option_ref, option_type in zip(LIST_OPTIONS_PEAKSEARCH, LIST_OPTIONS_TYPE_PEAKSEARCH):
  
              #             print "option_ref, option_type", option_ref, option_type
  
              if option_ref == option or option_ref.lower() == option:
  
                  #                 print "BINGO! I m able to read %s" % option_ref
                  #                 print "data type should be: %s" % option_type
  
                  try:
                      optionkey = CONVERTKEY_dict[option_ref]
                  except KeyError:
                      optionkey = option_ref
  
                  option_lower = option_ref.lower()
  
                  try:
                      if option_type.startswith("int"):
                          dict_param[optionkey] = int(config.getint(section, option_lower))
                      elif option_type.startswith(("float",)):
                          dict_param[optionkey] = float(config.getfloat(section, option_lower))
                      else:
                          dict_param[optionkey] = config.get(section, option_lower)
  
                  except ValueError:
                      print("Value of option '{}' has not the correct type".format(option))
                      return None
  
                  break
  
      #     print "Finally, I ve read these parameters"
      #     print dict_param
      return dict_param
  
  
  def read_background_flag(background_flag, verbose=0):
      r"""
      interpret the background flag (field used in FileSeries/Peak_Search.py)
  
      return two values to put in dict_param of peaksearch_series
      """
      formulaexpression = "A-B"
  
      if background_flag in ("auto", "AUTO", "yes", "YES", "y"):
          Data_for_localMaxima = "auto_background"
      elif background_flag in ("n", "no", "NO", None, "None", "NONE"):
          Data_for_localMaxima = None
      else:
          # should be a path to a single imagefile
          # (plus optionally a mathematical formula expression)
  
          ressplit = background_flag.split(";")
  
          if len(ressplit) != 2:
              filepath = background_flag
          elif len(ressplit) == 2:
              filepath, formulaexpression = ressplit
  
          if not os.path.exists(filepath):
              raise ValueError('File %s for background is does not exist'%filepath)
  
          Data_for_localMaxima = filepath
  
          if verbose: print("Image file path used for background", Data_for_localMaxima)
  
      return Data_for_localMaxima, formulaexpression
  
  
  def set_blacklist_filepath(filepathstr):
      r""" return None or path to file containing black listed spots list"""
      if filepathstr == "None":
          Remove_BlackListedPeaks_fromfile = None
      # fullpath
      else:
          Remove_BlackListedPeaks_fromfile = filepathstr
      return Remove_BlackListedPeaks_fromfile
  
  
  def set_rois_file(filepathstr):
      """ return None or path to file containing rois list"""
      if filepathstr == "None":
          roisfilepath = None
      # fullpath
      else:
          roisfilepath = filepathstr
      return roisfilepath
  
  
  # --- -------------- multiple file peak search
  def peaksearch_fileseries(fileindexrange,
                              filenameprefix="",
                              suffix="",
                              nbdigits=4,
                              dirname_in="/home/micha/LaueProjects/AxelUO2",
                              outputname=None,
                              dirname_out=None,
                              CCDLABEL="MARCCD165",
                              KF_DIRECTION="Z>0",  # not used yet
                              dictPeakSearch=None,
                              verbose=0,
                              writeResultDicts=0,
                              computetime=0):
      r"""
      peaksearch function to be called for multi or single processing
      """
      if computetime:
          t0 = ttt.time()
      print('\n\n ***** Starting peaksearch_fileseries()  *****\n\n')
      # peak search Parameters update from .psp file
      if isinstance(dictPeakSearch, dict):
          for key, val in list(dictPeakSearch.items()):
              PEAKSEARCHDICT_Convolve[key] = val
  
      if "MinPeakSize" not in PEAKSEARCHDICT_Convolve:
          PEAKSEARCHDICT_Convolve["MinPeakSize"] = 0.65
          PEAKSEARCHDICT_Convolve["MaxPeakSize"] = 3.
          if verbose: print("Default values for minimal and maximal peaksize are used!. Resp. 0.65 and 3 pixels.")
  
      PEAKSEARCHDICT_Convolve["PeakSizeRange"] = (copy.copy(PEAKSEARCHDICT_Convolve["MinPeakSize"]),
                                                  copy.copy(PEAKSEARCHDICT_Convolve["MaxPeakSize"]))
      del PEAKSEARCHDICT_Convolve["MinPeakSize"]
      del PEAKSEARCHDICT_Convolve["MaxPeakSize"]
  
      # ----handle reading of filename
      # special case for _mar.tif files...
      if nbdigits in ("varying",):
          pass
      # normal case
      else:
          nbdigits = int(nbdigits)
  
      if suffix == "":
          suffix = ".mccd"
  
      if dirname_in != None:
          filenameprefix_in = os.path.join(dirname_in, filenameprefix)
      else:
          filenameprefix_in = filenameprefix
  
      if outputname != None:
          prefix_outputname = outputname
      else:
          prefix_outputname = filenameprefix  # filename_wo_path[:-len(file_extension) - 1]
  
      if dirname_out != None:
          prefix_outputname = os.path.join(dirname_out, prefix_outputname)
  
      if len(fileindexrange) == 2:
          fileindexrange = fileindexrange[0], fileindexrange[1], 1
  
      # create a data considered as background from an imagefile
      BackgroundImageCreated = False
      flag_for_backgroundremoval = PEAKSEARCHDICT_Convolve["Data_for_localMaxima"]
  
      # flag_for_backgroundremoval is a file path to an imagefile
      # create background data: dataimage_bkg
      if flag_for_backgroundremoval not in ("auto_background", None) and not isinstance(
          flag_for_backgroundremoval, np.ndarray):
  
          fullpath_backgroundimage = PEAKSEARCHDICT_Convolve["Data_for_localMaxima"]
  
          BackgroundImageCreated = True
  
      # create
      blspots = PEAKSEARCHDICT_Convolve["Remove_BlackListedPeaks_fromfile"]
      if blspots is not None:
          PEAKSEARCHDICT_Convolve["Remove_BlackListedPeaks_fromfile"] = Get_blacklisted_spots(blspots)
          if PEAKSEARCHDICT_Convolve["Remove_BlackListedPeaks_fromfile"] is None:
              print('Warning!! blacklisted spots file may not be well read! blacked spots list is empty', PEAKSEARCHDICT_Convolve["Remove_BlackListedPeaks_fromfile"])
  
      DictPeaksList = {}
      file_ix, nb_empty_files = 0, 0  # nb of probed file, nb of zero peaks file
      listimageindices = list(range(fileindexrange[0],
                          fileindexrange[1] + 1,
                          fileindexrange[2]))
      nbimages = len(listimageindices)
      nbstepsprogress = 5
      progressstep = 0
      for fileindex in listimageindices:
          # TODO to move this branching elsewhere (readmccd)
          if suffix.endswith("_mar.tif"):
              filename_in = IOimage.setfilename(filenameprefix_in + "{}".format(fileindex) + suffix,
                                                                                          fileindex)
          else:
              #             filename_in = filenameprefix_in + encodingdigits % fileindex + suffix
              filename_in = filenameprefix_in + str(fileindex).zfill(nbdigits) + suffix
  
          tirets = "-" * 15
          if verbose: print("\n\n {} PeakSearch on filename {}\n{}\n{}{}{}n\n".format(
                  tirets, tirets, filename_in, tirets, tirets, tirets))
  
          if not os.path.exists(filename_in):
              raise ValueError("\n\n*******\nSomething wrong with the filename: {}. Please check "
                                              "carefully the filename!".format(filename_in))
  
          # remove a single image (considered as background) to current image
          if BackgroundImageCreated:
  
              if verbose: print("consider dataimagefile {} as background".format(fullpath_backgroundimage))
              # (dataimage_raw, _, _) = IOimage.readCCDimage(filename_in, CCDLabel=CCDLABEL,
              #                                                                         dirname=None)
  
              if "formulaexpression" in dictPeakSearch:
                  formulaexpression = dictPeakSearch["formulaexpression"]
              else:
                  raise ValueError('Missing "formulaexpression" to operate on images before '
                                  'peaksearch in peaksearch_fileseries()')
  
              if verbose: print("using {} in peaksearch_fileseries".format(formulaexpression))
  
              # for finding local maxima in image from formula
              PEAKSEARCHDICT_Convolve["Data_for_localMaxima"] = fullpath_backgroundimage
  
              # for fitting peaks in image from formula
              PEAKSEARCHDICT_Convolve["reject_negative_baseline"] = False
              PEAKSEARCHDICT_Convolve["formulaexpression"] = formulaexpression
              PEAKSEARCHDICT_Convolve["Fit_with_Data_for_localMaxima"] = True
  
          # --------------------------
          # launch peaksearch
          # -----------------------
          Res = PeakSearch(filename_in, CCDLabel=CCDLABEL,
                              Saturation_value=DictLT.dict_CCD[CCDLABEL][2],
                              Saturation_value_flatpeak=DictLT.dict_CCD[CCDLABEL][2],
                              **PEAKSEARCHDICT_Convolve)
  
          if Res in (False, None):
              print("No peak found for image file: ", filename_in)
  
              nb_empty_files += 1
          #             Isorted, fitpeak, localpeak = None, None, None
          else:  # write file with comments
              Isorted, _, _ = Res[:3]
  
              params_comments = "Peak Search and Fit parameters\n"
  
              params_comments += "# {}: {}\n".format("CCDLabel", CCDLABEL)
  
              for key, val in list(PEAKSEARCHDICT_Convolve.items()):
                  if not BackgroundImageCreated or key not in ("Data_for_localMaxima",):
                      params_comments += "# " + key + " : " + str(val) + "\n"
  
              if BackgroundImageCreated:
                  params_comments += ("# " + "Data_for_localMaxima"
                                      + " : {} \n".format(fullpath_backgroundimage))
              # .dat file extension is done in writefile_Peaklist()
              # filename_out = prefix_outputname + encodingdigits % fileindex
              # TODO valid whatever
              filename_out = prefix_outputname + str(fileindex).zfill(nbdigits)
              IOLT.writefile_Peaklist("{}".format(filename_out),
                                          Isorted,
                                          overwrite=1,
                                          initialfilename=filename_in,
                                          comments=params_comments)
              if writeResultDicts:                            
                  DictPeaksList[fileindex] = Isorted
  
          progress = int(np.floor(file_ix/nbimages*nbstepsprogress))
          if progress > progressstep:
              print('Imageindex: %d, Task Progress : %.2f %%' % (fileindex, file_ix / nbimages * 100))
              progressstep += 1
          file_ix += 1
  
      print("\n\n\n*******************\n\n\n task of peaksearch COMPLETED!")
      if computetime:
          print('Execution time %.2f sec'%(ttt.time()-t0))
      return DictPeaksList, file_ix, nb_empty_files
  
  # def peaksearch_multiprocessing(fileindexrange, filenameprefix, suffix="", nbdigits=4,
  #                                                     dirname_in="/home/micha/LaueProjects/AxelUO2",
  #                                                     outputname=None,
  #                                                     dirname_out=None,
  #                                                     CCDLABEL="MARCCD165",
  #                                                     KF_DIRECTION="Z>0",
  #                                                     dictPeakSearch=None,
  #                                                     nb_of_cpu=2,
  #                                                     verbose=0,
  #                                                     writeResultDicts=0):
  #     r"""
  #     launch several processes in parallel
  #     """
  #     import multiprocessing
  
  #     try:
  #         if len(fileindexrange) > 2:
  #             print("\n\n STEP INDEX is SET to 1 \n\n")
  #         index_start, index_final = fileindexrange[:2]
  #     except:
  #         raise ValueError("Need 2 file indices integers in fileindexrange=(indexstart, indexfinal)")
  
  #     t00 = ttt.time()
  
  #     max_nb_cpus = multiprocessing.cpu_count()
  #     nb_cpus = min(nb_of_cpu, max_nb_cpus)
  
  #     fileindexdivision = GT.getlist_fileindexrange_multiprocessing(index_start, index_final, nb_of_cpu)
  
  #     if nb_cpus > 1:
  #         print('using %d cpu(s)'%nb_cpus)
  #         fileindexdivision = GT.getlist_fileindexrange_multiprocessing(index_start, index_final, nb_cpus)
  #         nbimages = index_final - index_start + 1
  #         print('dispatch of fileindex ', fileindexdivision)
          
  #         peaksearch_fileseries.__defaults__ = (filenameprefix,
  #                                             suffix,
  #                                             nbdigits,
  #                                             dirname_in,
  #                                             outputname,
  #                                             dirname_out,
  #                                             CCDLABEL,
  #                                             KF_DIRECTION,
  #                                             dictPeakSearch,
  #                                             verbose,
  #                                             writeResultDicts,
  #                                             0) # compute execution time / task
  
  #         pool = multiprocessing.Pool(nb_of_cpu)
  #         multiple_results = pool.map(peaksearch_fileseries, fileindexdivision)
  
  #         # DictPeaksList, file_ix, nb_empty_files = multiple_results
  
  #     t_mp = ttt.time() - t00
  #     print("Execution time : %.2f" % t_mp)
  
  #     if nb_cpus > 1:
  #         nbtreatedimages = 0
  #         nbzeropeaksimages = 0
  #         for mres in multiple_results:
  #             print('nb treated files, nb zero peaks file', mres[1], mres[2])
  #             nbtreatedimages += mres[1]
  #             nbzeropeaksimages += mres[2]
  #         if nbzeropeaksimages != 0:
  #             print('total nb of zero peaks file', nbzeropeaksimages)
  #         else:
  #             print('all %s images contain at least one peak'%nbtreatedimages)
  
  #         # TODO  see end of indexFilesSeries()  to write a log file or hdf5 file with peaks props and other stas, nb of peaks per file, average nb , min and max number
  
  #     return multiple_results, nbtreatedimages, nbzeropeaksimages
  
  
  def peaklist_dict(prefixfilename, startindex, finalindex, dirname=None):
      r""" create a dict with key=image index and value=list of peaks """
      dict_peaks = {}
      for k in list(range(startindex, finalindex + 1)):
          filename = prefixfilename + "{:04d}.dat".format(k)
  
          array_peaks = IOLT.read_Peaklist(filename, dirname=dirname)
          dict_peaks[k] = array_peaks
  
      return dict_peaks
  
  
  def purgePeaksListFile(filename1, blacklisted_XY, dist_tolerance=0.5, dirname=None):
      r"""
      remove in peaklist .dat file peaks that are in blacklist
  
      :param blacklisted_XY:         [X1,Y1],[X2,Y2]
      """
      data_peak = IOLT.read_Peaklist(filename1, dirname=dirname)
  
      XY = data_peak[:, 0:2].T
  
      blacklisted_XY = np.array(blacklisted_XY).T
  
      peakX, peakY, tokeep = GT.removeClosePoints_two_sets(XY, blacklisted_XY,
                                              dist_tolerance=dist_tolerance, verbose=0)
  
      return peakX, peakY, tokeep
  
  def write_PurgedPeakListFile(filename1, blacklisted_XY, outputfilename, dist_tolerance=0.5,
                                                                          dirname=None):
      r"""
      write a new .dat file where peaks in blacklist are omitted
      """
      #peakX, peakY, tokeep
      _, _, tokeep = purgePeaksListFile(filename1, blacklisted_XY, dist_tolerance=0.5, dirname=dirname)
  
      data_peak = IOLT.read_Peaklist(filename1, dirname=dirname)
  
      new_data_peak = np.take(data_peak, tokeep, axis=0)
  
      if dirname is not None:
          outputfilename = os.path.join(dirname, outputfilename)
  
      IOLT.writefile_Peaklist(outputfilename,
                              new_data_peak,
                              overwrite=1,
                              initialfilename=filename1,
                              comments="Some peaks have been removed by write_PurgedPeakListFile",
                              dirname=dirname)
  
      print("New peak list file {} has been written".format(outputfilename))
  
  
  def removePeaks_inPeakList(PeakListfilename,
                              BlackListed_PeakListfilename,
                              outputfilename,
                              dist_tolerance=0.5,
                              dirname=None):
      r"""
      read peaks PeakListfilename and remove those in BlackListed_PeakListfilename
      and write a new peak list file
  
      .. note:: Not used ??
      """
      data_peak_blacklisted = IOLT.read_Peaklist(BlackListed_PeakListfilename, dirname=dirname)
  
      XY_blacklisted = data_peak_blacklisted[:, 0:2].T
  
      write_PurgedPeakListFile(PeakListfilename,
                              XY_blacklisted,
                              outputfilename,
                              dist_tolerance=0.5,
                              dirname=dirname)
  
  
  def merge_2Peaklist(filename1, filename2, dist_tolerance=5, dirname1=None, dirname2=None, verbose=0):
      r"""
      return merge spots data from two peaklists and removed duplicates within dist_tolerance (pixel)
      """
      data_peak_1 = IOLT.read_Peaklist(filename1, dirname=dirname1)
      data_peak_2 = IOLT.read_Peaklist(filename2, dirname=dirname2)
  
      XY1 = data_peak_1[:, 0:2]
      XY2 = data_peak_2[:, 0:2]
  
      #XY, ind_delele_1, ind_delele_2
      _, ind_delele_1, ind_delele_2 = GT.mergelistofPoints(XY1, XY2, dist_tolerance=dist_tolerance,
                                                              verbose=verbose)
  
      data1 = np.delete(data_peak_1, ind_delele_1, axis=0)
      data2 = np.delete(data_peak_2, ind_delele_2, axis=0)
  
      return np.concatenate((data1, data2), axis=0)
  
  
  def writefile_mergedPeaklist(filename1, filename2, outputfilename, dist_tolerance=5,
                                                                      dirname1=None,
                                                                      dirname2=None,
                                                                      verbose=0):
      r"""
      write peaklist file from the merge of spots data from two peaklists
      (and removed duplicates within dist_tolerance (pixel))
      """
      merged_data = merge_2Peaklist(
          filename1, filename2, dist_tolerance, dirname1, dirname2, verbose)
      comments = "Peaks from merging {} and {} with pixel tolerance {:.3f}".format(
          filename1, filename2, dist_tolerance)
      IOLT.writefile_Peaklist(outputfilename, merged_data, 1, None, comments, None)
  
      print("Merged peak lists written in file: ", outputfilename)
  