#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#    Project: Azimuthal integration
#             https://github.com/silx-kit/pyFAI
#
#    Copyright (C) 2013-2018 European Synchrotron Radiation Facility, Grenoble, France
#
#    Principal author:       Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

"""
pyFAI-calib

A tool for determining the geometry of a detector using a reference sample.

"""

__author__ = "Jerome Kieffer"
__contact__ = "Jerome.Kieffer@ESRF.eu"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
__date__ = "21/01/2021"
__status__ = "production"

import os
import sys
import time
import logging
logger = logging.getLogger(__name__)

import math
import numpy
from silx.image import marchingsquares
from scipy.stats import linregress
import fabio

from argparse import ArgumentParser
from urllib.parse import urlparse

from .matplotlib import pylab, matplotlib
from .utils import update_fig
from . import utils as gui_utils
from ..detectors import detector_factory, Detector
from ..geometryRefinement import GeometryRefinement
from .peak_picker import PeakPicker
from .. import units
from .. import average
from ..utils import measure_offset, expand_args, \
            readFloatFromKeyboard, FixedParameters, round_fft, \
            win32
from ..azimuthalIntegrator import AzimuthalIntegrator
from ..units import hc
from .. import version as PyFAI_VERSION
from .. import date as PyFAI_DATE
from ..calibrant import Calibrant, CALIBRANT_FACTORY
try:
    from ..ext._convolution import gaussian_filter
except ImportError:
    logger.debug("Backtrace", exc_info=True)
    from scipy.ndimage.filters import gaussian_filter

try:
    from ..ext import morphology
    pyFAI_morphology = True
except ImportError:
    logger.debug("Backtrace", exc_info=True)
    from scipy.ndimage import morphology
    pyFAI_morphology = False


def get_detector(detector, datafiles=None):
    """
    Detector factory taking into account the binning knowing the datafiles

    :param detector: string or detector or other junk
    :param datafiles: can be a list of images to be opened and their shape used
    :return: pyFAI.detector.Detector instance
    :raise RuntimeError: If no detector found
    """
    res = None
    if isinstance(detector, (str,)):
        try:
            res = detector_factory(detector)
        except RuntimeError:
            raise RuntimeError("Not a valid detector: %s" % detector)
    elif isinstance(detector, Detector):
        res = detector
    else:
        res = Detector()
    if datafiles and os.path.exists(datafiles[0]):
        shape = fabio.open(datafiles[0]).data.shape
        res.guess_binning(shape)
    return res


class AbstractCalibration(object):

    """
    Everything that is common to Calibration and Recalibration
    """

    win_error = "We are under windows with a 32 bit version of python,"\
                " matplotlib is not able to"\
                " display too many images without crashing, this"\
                " is why the window showing the diffraction image"\
                " is closed"
    _HELP = {"help": "Try to get the help of a given action, like 'refine?'. Use done when finished. "
             "Most command are composed of 'action parameter value' like 'set wavelength 1 A'.",
             "get": "print he value of a parameter",
             "set": "set the value of a parameter to the given value, i.e 'set wavelength 0.1 nm', units are optional",
             'fix': "fixes the value of a parameter so that its value will not be optimized, i.e. 'fix wavelength'",
             'free': "frees the parameter so that the value can be optimized, i.e. 'free wavelength'",
             'bound': "sets the upper and lower bound of a parameter: 'bound dist 0.1 0.2'",
             'bounds': "sets the upper and lower bound of all parameters",
             'refine': "performs a new cycle of refinement",
             'recalib': "extract a new set of rings and re-perform the calibration. One can specify how many rings to extract and the algorithm to use (blob, massif, watershed) and the nb_pts_per_deg in azimuth",
             'done': "finishes the processing, performs an integration and quits",
             'validate': "plot the offset between the calibrated image and the back-projected image",
             'validate2': "measures the offset of the center as function of azimuthal angle by cross-correlation of 2 plots, 180 deg appart. Option: number of azimuthal sliced, default: 36",
             'integrate': "perform the azimuthal integration and display results",
             'abort': "quit immediately, discarding any unsaved changes",
             'show': "Just print out the current parameter set. Optional parameters are units for length, rotation and wavelength, i.e. 'show mm deg A'",
             'reset': "Reset the geometry to the initial guess (rotation to zero, distance to 0.1m, poni at the center of the image)",
             'assign': "Change the assignment of a group of points to a rings",
             "weight": "toggle from weighted to unweighted mode...",
             "define": "Re-define the value for a constant internal parameter of the program like max_iter, nPt_1D, nPt_2D_azim, nPt_2D_rad, integrator_method, error_model. Warning: attribute change may be harmful !",
             "chiplot": "plot control point radial error as function of azimuthal angle, optional parameters: the rings for which this need to be plotted",
             "delete": "delete a group of points, provide the letter."
             }
    PARAMETERS = ["dist", "poni1", "poni2", "rot1", "rot2", "rot3", "wavelength"]
    UNITS = {"dist": "meter", "poni1": "meter", "poni2": "meter", "rot1": "radian",
             "rot2": "radian", "rot3": "radian", "wavelength": "meter"}
    VALID_URL = ["", 'file', 'hdf5', "nxs", "h5"]
    PTS_PER_DEG = 0.3

    def __init__(self, dataFiles=None, darkFiles=None, flatFiles=None, pixelSize=None,
                 splineFile=None, detector=None, wavelength=None, calibrant=None):
        """Constructor of AbstractCalibration

        :param dataFiles: list of filenames containing data images
        :param darkFiles: list of filenames containing dark current images
        :param flatFiles: list of filenames containing flat images
        :param pixelSize: size of the pixel in meter as 2 tuple
        :param splineFile: file containing the distortion of the taper
        :param detector: Detector name or instance
        :param wavelength: radiation wavelength in meter
        :param calibrant: pyFAI.calibrant.Calibrant instance
        """
        self.dataFiles = dataFiles
        self.darkFiles = darkFiles
        self.flatFiles = flatFiles
        self.pointfile = None

        self.detector = get_detector(detector, dataFiles)

        if splineFile and os.path.isfile(splineFile):
            self.detector.splineFile = os.path.abspath(splineFile)
        if pixelSize:
            if "__len__" in dir(pixelSize) and len(pixelSize) >= 2:
                self.detector.pixel1 = float(pixelSize[0])
                self.detector.pixel2 = float(pixelSize[1])
            else:
                self.detector.pixel1 = self.detector.pixel2 = float(pixelSize)

        self.cutBackground = None
        self.outfile = "merged.edf"
        self.peakPicker = None
        self.img = None
        self.ai = AzimuthalIntegrator(dist=0.1, detector=self.detector)
        self.wavelength = wavelength
        if wavelength:
            self.ai.wavelength = wavelength

        self.data = None
        self.basename = None
        self.geoRef = None
        self.reconstruct = False
        if calibrant:
            if isinstance(calibrant, Calibrant):
                self.calibrant = calibrant
            elif calibrant in CALIBRANT_FACTORY:
                self.calibrant = CALIBRANT_FACTORY(calibrant)
            elif os.path.isfile(calibrant) and os.path.isfile(calibrant):
                self.calibrant = Calibrant(calibrant)
            else:
                logger.error("Unable to handle such calibrant %s", calibrant)
                self.calibrant = None
        else:
            self.calibrant = None
        self.mask = None
        self.saturation = 0
        self.fixed = ["wavelength"]  # parameter fixed during optimization
        self.max_rings = None
        self.max_iter = 1000
        self.gui = True
        self.interactive = True
        self.filter = "mean"
        self.weighted = False
        self.polarization_factor = None
        self.parser = None
        self.nPt_1D = 1024
        self.nPt_2D_azim = 360
        self.nPt_2D_rad = 400
        self.unit = units.to_unit("2th_deg")
        self.keep = True
        self.check_calib = None
        self.fig_integrate = self.ax_xrpd_1d = self.ax_xrpd_2d = None
        self.fig_chiplot = self.ax_chiplot = None
        self.fig_center = self.ax_center = None
        self.integrator_method = "splitbbox"
        self.error_model = ""

    def __repr__(self):
        lst = ["Calibration object:"]
        if self.dataFiles:
            lst.append("data= " + ", ".join(self.dataFiles))
        else:
            lst.append("data= None")
        if self.darkFiles:
            lst.append("dark= " + ", ".join(self.darkFiles))
        else:
            lst.append("dark= None")
        if self.flatFiles:
            lst.append("flat= " + ", ".join(self.flatFiles))
        else:
            lst.append("flat= None")
        if self.fixed:
            lst.append("fixed=" + ", ".join(self.fixed))
        else:
            lst.append("fixed= None")
        lst.append(self.detector.__repr__())
        return os.linesep.join(lst)

    def configure_parser(self, version="calibration from pyFAI  version %s: %s" % (PyFAI_VERSION, PyFAI_DATE),
                         usage="pyFAI-calib [options] input_image.edf",
                         description=None, epilog=None):
        """Common configuration for parsers
        """
        self.parser = ArgumentParser(usage=usage, description=description, epilog=epilog)
        self.parser.add_argument("-V", "--version", action='version', version=version)
        self.parser.add_argument("args", metavar="FILE", help="List of files to calibrate", nargs='+')
        self.parser.add_argument("-o", "--out", dest="outfile",
                                 help="Filename where processed image is saved", metavar="FILE",
                                 default="merged.edf")
        self.parser.add_argument("-v", "--verbose",
                                 action="store_true", dest="debug", default=False,
                                 help="switch to debug/verbose mode")
        self.parser.add_argument("-c", "--calibrant", dest="spacing", metavar="FILE",
                                 help="Calibrant name or file containing d-spacing of the reference sample (MANDATORY, case sensitive !)",
                                 default=None)
        self.parser.add_argument("-w", "--wavelength", dest="wavelength", type=float,
                                 help="wavelength of the X-Ray beam in Angstrom. Mandatory ", default=None)
        self.parser.add_argument("-e", "--energy", dest="energy", type=float,
                                 help="energy of the X-Ray beam in keV (hc=%skeV.A)." % hc, default=None)
        self.parser.add_argument("-P", "--polarization", dest="polarization_factor",
                                 type=float, default=None,
                                 help="polarization factor, from -1 (vertical) to +1 (horizontal),"
                                 " default is None (no correction), synchrotrons are around 0.95")
        self.parser.add_argument("-i", "--poni", dest="poni", metavar="FILE",
                                 help="file containing the diffraction parameter (poni-file). MANDATORY for pyFAI-recalib!",
                                 default=None)
        self.parser.add_argument("-b", "--background", dest="background",
                                 help="Automatic background subtraction if no value are provided",
                                 default=None)
        self.parser.add_argument("-d", "--dark", dest="dark",
                                 help="list of comma separated dark images to average and subtract", default=None)
        self.parser.add_argument("-f", "--flat", dest="flat",
                                 help="list of comma separated flat images to average and divide", default=None)
        self.parser.add_argument("-s", "--spline", dest="spline",
                                 help="spline file describing the detector distortion", default=None)
        self.parser.add_argument("-D", "--detector", dest="detector_name",
                                 help="Detector name (instead of pixel size+spline)", default=None)
        self.parser.add_argument("-m", "--mask", dest="mask",
                                 help="file containing the mask (for image reconstruction)", default=None)
        self.parser.add_argument("-n", "--pt", dest="npt",
                                 help="file with datapoints saved. Default: basename.npt", default=None)
        self.parser.add_argument("--filter", dest="filter",
                                 help="select the filter, either mean(default), max or median",
                                 default="mean")
        self.parser.add_argument("-l", "--distance", dest="distance", type=float,
                                 help="sample-detector distance in millimeter. Default: 100mm", default=None)
        self.parser.add_argument("--dist", dest="dist", type=float,
                                 help="sample-detector distance in meter. Default: 0.1m", default=None)
        self.parser.add_argument("--poni1", dest="poni1", type=float,
                                 help="poni1 coordinate in meter. Default: center of detector", default=None)
        self.parser.add_argument("--poni2", dest="poni2", type=float,
                                 help="poni2 coordinate in meter. Default: center of detector", default=None)
        self.parser.add_argument("--rot1", dest="rot1", type=float,
                                 help="rot1 in radians. default: 0", default=None)
        self.parser.add_argument("--rot2", dest="rot2", type=float,
                                 help="rot2 in radians. default: 0", default=None)
        self.parser.add_argument("--rot3", dest="rot3", type=float,
                                 help="rot3 in radians. default: 0", default=None)

        self.parser.add_argument("--fix-dist", dest="fix_dist",
                                 help="fix the distance parameter", default=None, action="store_true")
        self.parser.add_argument("--free-dist", dest="fix_dist",
                                 help="free the distance parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-poni1", dest="fix_poni1",
                                 help="fix the poni1 parameter", default=None, action="store_true")
        self.parser.add_argument("--free-poni1", dest="fix_poni1",
                                 help="free the poni1 parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-poni2", dest="fix_poni2",
                                 help="fix the poni2 parameter", default=None, action="store_true")
        self.parser.add_argument("--free-poni2", dest="fix_poni2",
                                 help="free the poni2 parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-rot1", dest="fix_rot1",
                                 help="fix the rot1 parameter", default=None, action="store_true")
        self.parser.add_argument("--free-rot1", dest="fix_rot1",
                                 help="free the rot1 parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-rot2", dest="fix_rot2",
                                 help="fix the rot2 parameter", default=None, action="store_true")
        self.parser.add_argument("--free-rot2", dest="fix_rot2",
                                 help="free the rot2 parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-rot3", dest="fix_rot3",
                                 help="fix the rot3 parameter", default=None, action="store_true")
        self.parser.add_argument("--free-rot3", dest="fix_rot3",
                                 help="free the rot3 parameter. Default: Activated", default=None, action="store_false")

        self.parser.add_argument("--fix-wavelength", dest="fix_wavelength",
                                 help="fix the wavelength parameter. Default: Activated", default=True, action="store_true")
        self.parser.add_argument("--free-wavelength", dest="fix_wavelength",
                                 help="free the wavelength parameter. Default: Deactivated ", default=True, action="store_false")

        self.parser.add_argument("--tilt", dest="tilt",
                                 help="Allow initially detector tilt to be refined (rot1, rot2, rot3). Default: Activated", default=None, action="store_true")
        self.parser.add_argument("--no-tilt", dest="tilt",
                                 help="Deactivated tilt refinement and set all rotation to 0", default=None, action="store_false")

        self.parser.add_argument("--saturation", dest="saturation",
                                 help="consider all pixel>max*(1-saturation) as saturated and "
                                 "reconstruct them, default: 0 (deactivated)",
                                 default=0, type=float)
        self.parser.add_argument("--weighted", dest="weighted",
                                 help="weight fit by intensity, by default not.",
                                 default=False, action="store_true")
        self.parser.add_argument("--npt", dest="nPt_1D",
                                 help="Number of point in 1D integrated pattern, Default: 1024", type=int,
                                 default=1024)
        self.parser.add_argument("--npt-azim", dest="nPt_2D_azim",
                                 help="Number of azimuthal sectors in 2D integrated images. Default: 360", type=int,
                                 default=360)
        self.parser.add_argument("--npt-rad", dest="nPt_2D_rad",
                                 help="Number of radial bins in 2D integrated images. Default: 400", type=int,
                                 default=400)
        self.parser.add_argument("--unit", dest="unit",
                                 help="Valid units for radial range: 2th_deg, 2th_rad, q_nm^-1,"
                                 " q_A^-1, r_mm. Default: 2th_deg", type=str, default="2th_deg")
        self.parser.add_argument("--no-gui", dest="gui",
                                 help="force the program to run without a Graphical interface",
                                 default=True, action="store_false")
        self.parser.add_argument("--no-interactive", dest="interactive",
                                 help="force the program to run and exit without prompting"
                                 " for refinements", default=True, action="store_false")

    def analyse_options(self, options=None, args=None):
        """Analyzes options and arguments

        :return: option,arguments
        """
        if (options is None) and (args is None):
            options = self.parser.parse_args()
            args = options.args
        if options.debug:
            logger.setLevel(logging.DEBUG)
        self.outfile = options.outfile
        if options.dark:
            self.darkFiles = [f for f in options.dark.split(",") if os.path.isfile(f)]
            if not self.darkFiles:  # empty container !!!
                logger.error("No dark file exists !!!")
                self.darkFiles = None
        if options.flat:
            self.flatFiles = [f for f in options.flat.split(",") if os.path.isfile(f)]
            if not self.flatFiles:  # empty container !!!
                logger.error("No flat file exists !!!")
                self.flatFiles = None

        if options.detector_name:
            self.detector = get_detector(options.detector_name, args)
            self.ai.detector = self.detector
        if options.spline:
            if "Pilatus" in self.detector.name:
                self.detector.set_splineFile(options.spline)  # is as 2-tuple of path
            elif os.path.isfile(options.spline):
                self.detector.set_splineFile(os.path.abspath(options.spline))
            else:
                logger.error("Unknown spline file %s", options.spline)

        if options.mask and os.path.isfile(options.mask):
            self.mask = (fabio.open(options.mask).data != 0)
        else:  # Use default mask provided by detector
            self.mask = self.detector.mask

        self.pointfile = options.npt
        if options.spacing:
            if options.spacing in CALIBRANT_FACTORY:
                self.calibrant = CALIBRANT_FACTORY(options.spacing)
            elif os.path.isfile(options.spacing):
                self.calibrant = Calibrant(options.spacing)
            else:
                logger.error("No such Calibrant / d-Spacing file: %s", options.spacing)

        if self.calibrant is None:
            self.read_dSpacingFile(True)

        if options.poni:
            self.ai.load(options.poni)

        if options.wavelength:
            self.ai.wavelength = self.wavelength = 1e-10 * options.wavelength
        elif options.energy:
            self.ai.wavelength = self.wavelength = 1e-10 * hc / options.energy
#        else:
            # This should be read from the poni. It it is missing; it is called in preprocess.
#            self.read_wavelength()
#            pass
        if options.distance:
            self.ai.dist = 1e-3 * options.distance
        if options.dist:
            self.ai.dist = options.dist

        if options.tilt is False:
            self.ai.rot1 = 0.0
            self.ai.rot2 = 0.0
            self.ai.rot3 = 0.0
        if options.poni1 is not None:
            self.ai.poni1 = options.poni1
        if options.poni2 is not None:
            self.ai.poni2 = options.poni2
        if options.rot1 is not None:
            self.ai.rot1 = options.rot1
        if options.rot2 is not None:
            self.ai.rot2 = options.rot2
        if options.rot3 is not None:
            self.ai.rot3 = options.rot3
        self.dataFiles = expand_args(args)
        if not self.dataFiles:
            raise RuntimeError("Please provide some calibration images ... "
                               "if you want to analyze them. Try also the "
                               "--help option to see all options!")
        self.fixed = FixedParameters()
        if options.tilt is not None:
            for key in ["rot1", "rot2", "rot3"]:
                self.fixed.add_or_discard(key, not(options.tilt))

        self.fixed.add_or_discard("dist", options.fix_dist)
        self.fixed.add_or_discard("poni1", options.fix_poni1)
        self.fixed.add_or_discard("poni2", options.fix_poni2)
        self.fixed.add_or_discard("rot1", options.fix_rot1)
        self.fixed.add_or_discard("rot2", options.fix_rot2)
        self.fixed.add_or_discard("rot3", options.fix_rot3)
        self.fixed.add_or_discard("wavelength", options.fix_wavelength)
        print(self.fixed)
        self.saturation = options.saturation

        self.gui = options.gui
        self.interactive = options.interactive
        self.filter = options.filter
        self.weighted = options.weighted
        self.polarization_factor = options.polarization_factor
        self.detector = self.ai.detector
        self.nPt_1D = options.nPt_1D
        self.nPt_2D_azim = options.nPt_2D_azim
        self.nPt_2D_rad = options.nPt_2D_rad
        self.unit = units.to_unit(options.unit)
        if options.background is not None:
            try:
                self.cutBackground = float(options.background)
            except Exception:
                self.cutBackground = True

        return options, args

    def get_pixelSize(self, ans):
        """convert a comma separated sting into pixel size"""
        sp = ans.split(",")
        if len(sp) >= 2:
            try:
                pixelSizeXY = [float(i) * 1e-6 for i in sp[:2]]
            except Exception:
                logger.error("error in reading pixel size_2")
                return
        elif len(sp) == 1:
            px = sp[0]
            try:
                pixelSizeXY = [float(px) * 1e-6, float(px) * 1e-6]
            except Exception:
                logger.error("error in reading pixel size_1")
                return
        else:
            logger.error("error in reading pixel size_0")
            return
        self.detector.pixel1 = pixelSizeXY[1]
        self.detector.pixel2 = pixelSizeXY[0]

    def read_pixelsSize(self):
        """Read the pixel size from prompt if not available"""
        if (self.detector.pixel1 is None) and (self.detector.splineFile is None):
            pixelSize = [15, 15]
            ans = input("Please enter the pixel size (in micron, comma separated X,Y "
                        " i.e. %.2e,%.2e) or a spline file: " % tuple(pixelSize)).strip()
            if os.path.isfile(ans):
                self.detector.splineFile = ans
            else:
                self.get_pixelSize(ans)

    def read_dSpacingFile(self, verbose=True):
        """Read the name of the calibrant / file with d-spacing"""
        if (self.calibrant is None):
            comments = ["pyFAI calib has changed !!!",
                        "Instead of entering the 2theta value, which was tedious,"
                        "the program takes a calibrant name or a d-spacing file in input "
                        "(just a serie of number representing the inter-planar "
                        "distance in Angstrom)",
                        "and an associated wavelength",
                        "You will be asked to enter the ring number,"
                        " which is usually a simpler than the 2theta value."]
            if verbose:
                print(os.linesep.join(comments))
            valid = False
            while valid:
                ans = input("Please enter the calibrant name or the file"
                            " containing the d-spacing:\t").strip()
                if ans in CALIBRANT_FACTORY:
                    self.calibrant = CALIBRANT_FACTORY(ans)
                    valid = True
                elif os.path.isfile(ans):
                    self.calibrant = Calibrant(ans)
                    valid = True

    def read_wavelength(self):
        """Read the wavelength"""
        while not self.wavelength:
            ans = input("Please enter wavelength in Angstrom:\t").strip()
            try:
                self.wavelength = self.ai.wavelength = 1e-10 * float(ans)
            except Exception:
                self.wavelength = None

    def preprocess(self):
        """
        Common part:
        do dark, flat correction thresholding, ...
        and read missing data from keyboard if needed
        """
        # GF: self.saturation ignored if none of the other options active...
        if len(self.dataFiles) > 1 or self.cutBackground or self.darkFiles or self.flatFiles:
            self.outfile = average.average_images(self.dataFiles, self.outfile,
                                                  threshold=self.saturation,
                                                  minimum=self.cutBackground,
                                                  darks=self.darkFiles,
                                                  flats=self.flatFiles,
                                                  filter_=self.filter)
        else:
            self.outfile = self.dataFiles[0]

        url = urlparse(self.outfile)
        if (sys.platform == "win32") and (len(url.scheme) == 1):  # "c:" like path
            path = self.outfile
        else:
            if url.scheme not in self.VALID_URL:
                logger.warning("unexpected URL: %s", self.outfile)
            path = url.path
        self.basename, ext = os.path.splitext(path)
        if ext in [".gz", ".bz2"]:
            self.basename = os.path.splitext(self.basename)[0]

        if isinstance(self, Recalibration):
            self.keep = False
            self.pointfile = None
        else:
            self.pointfile = self.basename + ".npt"
        if self.wavelength is None:
            self.wavelength = self.ai.wavelength

        data = fabio.open(self.outfile).data
        self.peakPicker = PeakPicker(data, reconst=self.reconstruct, mask=self.mask,
                                     pointfile=self.pointfile, calibrant=self.calibrant,
                                     wavelength=self.ai.wavelength, detector=self.detector)
        if not self.keep:
            self.peakPicker.points.reset()
            if not self.peakPicker.points.calibrant.wavelength:
                self.peakPicker.points.calibrant.wavelength = self.ai.wavelength
            elif self.ai.wavelength != self.peakPicker.points.calibrant.wavelength:
                self.peakPicker.points.calibrant.setWavelength_change2th(self.ai.wavelength)
        if not self.peakPicker.points.calibrant.dSpacing:
            wl = self.peakPicker.points.calibrant.wavelength
            self.read_dSpacingFile()
            if wl:
                self.peakPicker.points.calibrant.wavelength = wl
        if not self.peakPicker.points.calibrant.wavelength:
            self.read_wavelength()
            self.peakPicker.points.calibrant.wavelength = self.wavelength

    def extract_cpt(self, method="massif", pts_per_deg=1.0):
        """
        Performs an automatic keypoint extraction:
        Can be used in recalib or in calib after a first calibration has been performed.

        :param method: method for keypoint extraction
        :param pts_per_deg: number of control points per azimuthal degree (increase for better precision)
        """
        logger.info("in extract_cpt with method %s", method)
        assert self.ai
        assert self.calibrant
        assert self.peakPicker
        self.peakPicker.reset()
        self.peakPicker.init(method, False)
        if self.geoRef:
            self.ai.setPyFAI(**self.geoRef.getPyFAI())
        tth = numpy.array([i for i in self.calibrant.get_2th() if i is not None])
        tth = numpy.unique(tth)
        tth_min = numpy.zeros_like(tth)
        tth_max = numpy.zeros_like(tth)
        delta = (tth[1:] - tth[:-1]) / 4.0
        tth_max[:-1] = delta
        tth_max[-1] = delta[-1]
        tth_min[1:] = -delta
        tth_min[0] = -delta[0]
        tth_max += tth
        tth_min += tth

        if self.geoRef:
            ttha = self.geoRef.get_ttha()
            chia = self.geoRef.get_chia()
            if (ttha is None) or (ttha.shape != self.peakPicker.data.shape):
                ttha = self.geoRef.twoThetaArray(self.peakPicker.data.shape)
            if (chia is None) or (chia.shape != self.peakPicker.data.shape):
                chia = self.geoRef.chiArray(self.peakPicker.data.shape)
        else:
            ttha = self.ai.twoThetaArray(self.peakPicker.data.shape)
            chia = self.ai.chiArray(self.peakPicker.data.shape)
        rings = 0
        self.peakPicker.sync_init()
        if self.max_rings is None:
            self.max_rings = tth.size

        ms = marchingsquares.MarchingSquaresMergeImpl(ttha, self.mask, use_minmax_cache=True)

        for i in range(tth.size):
            if rings >= self.max_rings:
                break
            mask = numpy.logical_and(ttha >= tth_min[i], ttha < tth_max[i])
            if self.mask is not None:
                mask = numpy.logical_and(mask, numpy.logical_not(self.mask))
            size = mask.sum(dtype=int)
            if (size > 0):
                rings += 1
                self.peakPicker.massif_contour(mask)
                if self.gui:
                    update_fig(self.peakPicker.fig)
                sub_data = self.peakPicker.data.ravel()[numpy.where(mask.ravel())]
                mean = sub_data.mean(dtype=numpy.float64)
                std = sub_data.std(dtype=numpy.float64)
                upper_limit = mean + std
                mask2 = numpy.logical_and(self.peakPicker.data > upper_limit, mask)
                size2 = mask2.sum(dtype=int)
                if size2 < 1000:
                    upper_limit = mean
                    mask2 = numpy.logical_and(self.peakPicker.data > upper_limit, mask)
                    size2 = mask2.sum()
                # length of the arc:
                points = ms.find_pixels(tth[i])
                seeds = set((i[0], i[1]) for i in points if mask2[i[0], i[1]])
                # max number of points: 360 points for a full circle
                azimuthal = chia[points[:, 0].clip(0, self.peakPicker.data.shape[0]), points[:, 1].clip(0, self.peakPicker.data.shape[1])]
                nb_deg_azim = numpy.unique(numpy.rad2deg(azimuthal).round()).size
                keep = int(nb_deg_azim * pts_per_deg)
                if keep == 0:
                    continue
                dist_min = len(seeds) / 2.0 / keep
                # why 3.0, why not ?

                logger.info("Extracting datapoint for ring %s (2theta = %.2f deg); "
                            "searching for %i pts out of %i with I>%.1f, dmin=%.1f" %
                            (i, numpy.degrees(tth[i]), keep, size2, upper_limit, dist_min))
                _res = self.peakPicker.peaks_from_area(mask=mask2, Imin=upper_limit, keep=keep, method=method, ring=i, dmin=dist_min, seed=seeds)

        self.peakPicker.points.save(self.basename + ".npt")
        if self.weighted:
            self.data = self.peakPicker.points.getWeightedList(self.peakPicker.data)
        else:
            self.data = self.peakPicker.points.getList()

    def refine(self):
        """
        Contains the common geometry refinement part
        """
        if win32 and self.peakPicker is not None:
            logging.info(self.win_error)
            self.peakPicker.closeGUI()
        print("Before refinement, the geometry is:")
        print(self.geoRef)
        previous = sys.maxsize
        finished = False
        fig2 = None
        while not finished:
            count = 0
            if "wavelength" in self.fixed:
                while (previous > self.geoRef.chi2()) and (count < self.max_iter):
                    if (count == 0):
                        previous = sys.maxsize
                    else:
                        previous = self.geoRef.chi2()
                    self.geoRef.refine2(1000000, fix=self.fixed)
                    print(self.geoRef)
                    count += 1
            else:
                while previous > self.geoRef.chi2_wavelength() and (count < self.max_iter):
                    if (count == 0):
                        previous = sys.maxsize
                    else:
                        previous = self.geoRef.chi2_wavelength()
                    self.geoRef.refine2_wavelength(1000000, fix=self.fixed)
                    print(self.geoRef)
                    count += 1
                self.peakPicker.points.setWavelength_change2th(self.geoRef.wavelength)
            self.geoRef.save(self.basename + ".poni")
            self.geoRef.del_ttha()
            self.geoRef.del_dssa()
            self.geoRef.del_chia()
            tth = self.geoRef.twoThetaArray(self.peakPicker.shape)
            dsa = self.geoRef.solidAngleArray(self.peakPicker.shape)
            # self.geoRef.chiArray(self.peakPicker.shape)
            # self.geoRef.corner_array(self.peakPicker.shape, unit=units.TTH_RAD, scale=False)
            if win32:
                logger.info(self.win_error)
            else:
                if self.gui:
                    self.peakPicker.contour(tth)
                    if self.interactive:
                        if fig2 is None:
                            fig2 = pylab.plt.figure()
                            sp = fig2.add_subplot(111)
                            im = sp.imshow(dsa, origin="lower")
                            _cbar = fig2.colorbar(im)  # Add color bar
                            sp.set_title("Pixels solid-angle (relative to PONI)")
                        else:
                            im.set_array(dsa)
                            im.autoscale()

                        fig2.show()
                        update_fig(fig2)

            if self.interactive:
                finished = self.prompt()
            else:
                finished = True
            if not finished:
                previous = sys.maxsize

    def prompt(self):
        """
        prompt for commands to guide the calibration process

        :return: True when the user is happy with what he has, False to request another refinement
        """

        while True:
            req_help = False
            print("Fixed: " + ", ".join(self.fixed))
            ans = input("Modify parameters (or ? for help)?\t ").strip()
            if "?" in ans:
                req_help = True
            if not ans:
                print("'done' to continue")
                continue
            words = ans.lower().split()
            action = words[0]
            if action in ["help", "?"]:
                req_help = True
            if req_help:
                for what in self._HELP.keys():
                    if action.startswith(what):
                        print("Help on %s" % what)
                        print(self._HELP[what])
                        break
                else:
                    print("Help on commands")
                    print(self._HELP["help"])
                    print("Valid actions: " + ", ".join(self._HELP.keys()))
                    print("Valid parameters: " + ", ".join(self.PARAMETERS))
            elif action == "get":  # get wavelength
                if (len(words) >= 2):
                    for param in words[1:]:
                        if param in self.PARAMETERS:
                            print("Value of parameter %s: %s  %s" % (param, self.geoRef.__getattribute__(param), self.UNITS[param]))
                        else:
                            print("No a parameter: %s" % param)
                else:
                    print(self._HELP[action])

            elif action == "set":  # set wavelength 1e-10
                if (len(words) in (3, 4)) and words[1] in self.PARAMETERS:
                    param = words[1]
                    try:
                        value = float(words[2])
                    except ValueError:
                        logger.warning("invalid value")
                    else:
                        scale = 1.0
                        if len(words) == 4:
                            unit = units.to_unit(words[3], units.LENGTH_UNITS + units.ANGLE_UNITS)
                            if unit:
                                scale = unit.scale
                        setattr(self.geoRef, param, value / scale)
                else:
                    print(self._HELP[action])
            elif action == "fix":  # fix wavelength
                if (len(words) >= 2):
                    for param in words[1:]:
                        if param in self.PARAMETERS:
                            print("Value of parameter %s: %s %s" % (param, self.geoRef.__getattribute__(param), self.UNITS[param]))
                            self.fixed.add(param)
                        else:
                            print("No a parameter: %s" % param)
                else:
                    print(self._HELP[action])
            elif action == "free":  # free wavelength
                if (len(words) >= 2):
                    for param in words[1:]:
                        if param in self.PARAMETERS:
                            print("Value of parameter %s: %s %s" % (param, self.geoRef.__getattribute__(param), self.UNITS[param]))
                            self.fixed.discard(param)
                        else:
                            print("No a parameter: %s" % param)
                else:
                    print(self._HELP[action])

            elif action == "recalib":
                max_rings = None
                pts_per_deg = self.PTS_PER_DEG
                if len(words) >= 2:
                    try:
                        max_rings = int(words[1])
                    except Exception:
                        logger.warning("specify the number of rings to extract")
                        max_rings = None
                    else:
                        self.max_rings = max_rings
                else:
                    self.max_rings = None
                if len(words) >= 3 and words[2] in PeakPicker.VALID_METHODS:
                    method = words[2]
                else:
                    method = "blob"
                if len(words) >= 4:
                    try:
                        pts_per_deg = float(words[3])
                    except ValueError:
                        pts_per_deg = self.PTS_PER_DEG
                self.extract_cpt(method, pts_per_deg)
                self.geoRef.data = numpy.array(self.data, dtype=numpy.float64)
                return False
            elif action == "bound":  # bound dist
                if len(words) >= 2 and words[1] in self.PARAMETERS:
                    param = words[1]
                    if len(words) == 2:
                        text = ("Enter %s in %s " % (param, self.UNITS[param]) +
                                "(or %s_min[%.3f] %s[%.3f] %s_max[%.3f]):\t " % (
                                    param, self.geoRef.__getattribute__("get_%s_min" % param)(),
                                    param, self.geoRef.__getattribute__("get_%s" % param)(),
                                    param, self.geoRef.__getattribute__("get_%s_max" % param)()))
                        values = {
                            1: [self.geoRef.__getattribute__("set_%s" % param)],
                            2: [self.geoRef.__getattribute__("set_%s_min" % param),
                                self.geoRef.__getattribute__("set_%s_max" % param)],
                            3: [self.geoRef.__getattribute__("set_%s_min" % param),
                                self.geoRef.__getattribute__("set_%s" % param),
                                self.geoRef.__getattribute__("set_%s_max" % param)]}
                        readFloatFromKeyboard(text, values)
                    elif len(words) == 3:
                        try:
                            value = float(words[2])
                        except ValueError:
                            logger.warning("invalid value")
                        else:
                            self.geoRef.__getattribute__("set_%s" % param)(value)
                    elif len(words) == 4:
                        try:
                            value_min = float(words[2])
                            value_max = float(words[3])
                        except ValueError:
                            logger.warning("invalid value")
                        else:
                            self.geoRef.__getattribute__("set_%s_min" % param)(value_min)
                            self.geoRef.__getattribute__("set_%s_max" % param)(value_max)
                    elif len(words) == 5:
                        try:
                            value_min = float(words[2])
                            value = float(words[3])
                            value_max = float(words[4])
                        except ValueError:
                            logger.warning("invalid value")
                        else:
                            self.geoRef.__getattribute__("set_%s_min" % param)(value_min)
                            self.geoRef.__getattribute__("set_%s" % param)(value)
                            self.geoRef.__getattribute__("set_%s_max" % param)(value_max)
                    else:
                        print(self._HELP[action])
                else:
                    print(self._HELP[action])
            elif action == "bounds":
                readFloatFromKeyboard("Enter Distance in meter "
                                      "(or dist_min[%.3f] dist[%.3f] dist_max[%.3f]):\t " %
                                      (self.geoRef.dist_min, self.geoRef.dist, self.geoRef.dist_max),
                                      {1: [self.geoRef.set_dist], 2: [self.geoRef.set_dist_min, self.geoRef.set_dist_max],
                                       3: [self.geoRef.set_dist_min, self.geoRef.set_dist, self.geoRef.set_dist_max]})
                readFloatFromKeyboard("Enter Poni1 in meter "
                                      "(or poni1_min[%.3f] poni1[%.3f] poni1_max[%.3f]):\t " %
                                      (self.geoRef.poni1_min, self.geoRef.poni1, self.geoRef.poni1_max),
                                      {1: [self.geoRef.set_poni1], 2: [self.geoRef.set_poni1_min, self.geoRef.set_poni1_max],
                                       3: [self.geoRef.set_poni1_min, self.geoRef.set_poni1, self.geoRef.set_poni1_max]})
                readFloatFromKeyboard("Enter Poni2 in meter "
                                      "(or poni2_min[%.3f] poni2[%.3f] poni2_max[%.3f]):\t " %
                                      (self.geoRef.poni2_min, self.geoRef.poni2, self.geoRef.poni2_max),
                                      {1: [self.geoRef.set_poni2], 2: [self.geoRef.set_poni2_min, self.geoRef.set_poni2_max],
                                       3: [self.geoRef.set_poni2_min, self.geoRef.set_poni2, self.geoRef.set_poni2_max]})
                readFloatFromKeyboard("Enter Rot1 in rad "
                                      "(or rot1_min[%.3f] rot1[%.3f] rot1_max[%.3f]):\t " %
                                      (self.geoRef.rot1_min, self.geoRef.rot1, self.geoRef.rot1_max),
                                      {1: [self.geoRef.set_rot1], 2: [self.geoRef.set_rot1_min, self.geoRef.set_rot1_max],
                                       3: [self.geoRef.set_rot1_min, self.geoRef.set_rot1, self.geoRef.set_rot1_max]})
                readFloatFromKeyboard("Enter Rot2 in rad "
                                      "(or rot2_min[%.3f] rot2[%.3f] rot2_max[%.3f]):\t " %
                                      (self.geoRef.rot2_min, self.geoRef.rot2, self.geoRef.rot2_max),
                                      {1: [self.geoRef.set_rot2], 2: [self.geoRef.set_rot2_min, self.geoRef.set_rot2_max],
                                       3: [self.geoRef.set_rot2_min, self.geoRef.set_rot2, self.geoRef.set_rot2_max]})
                readFloatFromKeyboard("Enter Rot3 in rad "
                                      "(or rot3_min[%.3f] rot3[%.3f] rot3_max[%.3f]):\t " %
                                      (self.geoRef.rot3_min, self.geoRef.rot3, self.geoRef.rot3_max),
                                      {1: [self.geoRef.set_rot3], 2: [self.geoRef.set_rot3_min, self.geoRef.set_rot3_max],
                                       3: [self.geoRef.set_rot3_min, self.geoRef.set_rot3, self.geoRef.set_rot3_max]})
            elif action == "done":
                self.postProcess()
                return True
            elif action == "quit":
                return True
            elif action == "refine":
                return False
            elif action == "fit":
                return False
            elif action == "validate":
                self.validate_calibration()
            elif action == "validate2":
                if len(words) > 1:
                    nb = int(words[1])
                else:
                    nb = 36
                self.validate_center(nb)
            elif action == "integrate":
                self.postProcess()
            elif action == "abort":
                sys.exit()
            elif action == "show":
                args = []
                print("The current parameter set is:")
                if len(words) > 1:
                    args = ans.split()[1:]
                print(self.geoRef.__repr__(*args))
            elif action == "reset":
                if len(words) > 1:
                    how = words[1]
                else:
                    how = "center"
                self.reset_geometry(how)
            elif action == "assign":
                # Re assign a group of point to a ring ...
                if self.peakPicker and self.peakPicker.points:
                    control_points = self.peakPicker.points
                    control_points.readRingNrFromKeyboard()
                    control_points.save(self.basename + ".npt")
                    if self.weighted:
                        self.data = self.peakPicker.points.getWeightedList(self.peakPicker.data)
                    else:
                        self.data = self.peakPicker.points.getList()
                    self.geoRef.data = numpy.array(self.data, dtype=numpy.float64)
            elif action == "weight":
                old = self.weighted
                if len(words) == 2:
                    value = words[1].lower()
                    if value in ("0", "off", "no", "none", "false"):
                        self.weighted = False
                    elif value in ("1", "on", "yes", "true"):
                        self.weighted = True
                    else:
                        logger.warning("Unrecognized argument for weight: %s", value)
                        continue
                print("Weights: %s" % self.weighted)
                if (old != self.weighted):
                    if self.weighted:
                        self.data = self.peakPicker.points.getWeightedList(self.peakPicker.data)
                    else:
                        self.data = self.peakPicker.points.getList()
                    self.geoRef.data = numpy.array(self.data, dtype=numpy.float64)
            elif action == "define":
                if len(words) == 3:
                    param = words[1]
                    sval = words[2]
                    for cs_param in dir(self):
                        if cs_param.lower() == param:
                            oldval = self.__getattribute__(cs_param)
                            t = type(oldval)
                            print("constant %s was %s of type %s, setting to %s" % (cs_param, oldval, t, sval))
                            try:
                                newval = t(sval)
                            except Exception as err:
                                print("Unable to convert type")
                                logger.warning(err)
                            self.__setattr__(cs_param, newval)
                            break
                    else:
                        print("No such parameter %s" % param)
                else:
                    print(self._HELP[action])
            elif action == "chiplot":
                    print(self._HELP[action])
                    rings = None
                    if len(words) > 1:
                        try:
                            rings = [int(i) for i in words[1:]]
                        except ValueError:
                            print("Please provide ring numbers ... ")
                    self.chiplot(rings)
            elif action == "delete":
                if len(words) < 2:
                    print(self._HELP[action])
                else:
                    for code in words[1:]:
                        self.peakPicker.remove_grp(code)
                    self.data = self.peakPicker.points.getList()
                    self.geoRef.data = numpy.array(self.data, dtype=numpy.float64)
            else:
                logger.warning("Unrecognized action: %s, type 'quit' to leave ", action)

    def chiplot(self, rings=None):
        """
        plot delta_2theta/2theta = f(chi) and fit the curve.

        :param rings: list of rings to consider
        """
        from scipy.optimize import leastsq
        model = lambda x, mean, amp, phase: mean + amp * numpy.sin(x + phase)
        error = lambda param, xdata, ydata: model(xdata, *param) - ydata

        def jacob(param, xdata, ydata):
            j = numpy.ones((param.size, xdata.size))
            j[1,:] = numpy.sin(xdata + param[2])
            j[2,:] = param[1] * numpy.cos(xdata + param[2])
            return j

        sqrt2 = math.sqrt(2.)
        ttha = self.geoRef.twoThetaArray(self.detector.shape)
        resolution = numpy.rad2deg(max(abs(ttha[1:] - ttha[:-1]).max(),
                                       abs(ttha[:, 1:] - ttha[:,:-1]).max()))
        if self.gui:
            if self.fig_chiplot:
                self.fig_chiplot.clf()
            else:
                self.fig_chiplot = pylab.plt.figure()
            self.ax_chiplot = self.fig_chiplot.add_subplot(1, 1, 1)
            self.ax_chiplot.set_xlim(-180, 180)
            self.ax_chiplot.set_xticks(numpy.linspace(-180, 180, 9))
            self.ax_chiplot.set_xlabel("Azimuthal angle $\chi$ ($^o$)")
            self.ax_chiplot.set_ylabel(r"Error in Radial angle $\Delta$ 2$\theta$/2$\theta$*10$^4$")
            self.ax_chiplot.set_title("Chi plot")

        else:
            print("chiplot display only possible with GUI")
        if rings is None:
            rings = list(set(int(i[2]) for i in self.data))
            rings.sort()
        for ring in rings:
            ref_2th = numpy.rad2deg(self.calibrant.get_2th()[ring])
            print("Fitting ring #%s (2th=%.3fdeg)" % (ring, ref_2th))
            d1 = []
            d2 = []
            for i in self.data:
                if i[2] == ring:
                    d1.append(i[0])
                    d2.append(i[1])
            if len(d1) < 5:
                print(" Skip group of length %i" % len(d1))
                continue
            d1 = numpy.array(d1)
            d2 = numpy.array(d2)
            tth = numpy.rad2deg(self.geoRef.tth(d1, d2))
            err4 = (tth - ref_2th) / ref_2th * 10000
            chi = self.geoRef.chi(d1, d2)
            mean = err4.mean()
            amp = err4.std() * sqrt2
            phase = 0.0
            param = numpy.array([mean, amp, phase])
            print(" guessed err4 = %.3f + %.3f *sin($\chi$+ %.3f )" % (mean, amp, phase))
            res = leastsq(error, param, (chi, err4), jacob, col_deriv=True)
            popt = res[0]
            str_res = "%.3f + %.3f *sin($\chi$+ %.3f )" % tuple(popt)
            print(" fitted err4 = " + str_res)
            chi = numpy.rad2deg(chi)
            if self.ax_chiplot:
                color = list(matplotlib.colors.cnames.keys())[ring]
                self.ax_chiplot.plot(chi, err4, "o", color=color, label="ring #%i (%.3f$^o$)" % (ring, ref_2th))
                chi2 = numpy.linspace(-180, 180, 360)
                self.ax_chiplot.plot(chi2, model(numpy.deg2rad(chi2), *popt), color=color, label=str_res)

        self.ax_chiplot.legend()
        if not gui_utils.main_loop:
            self.fig_chiplot.show()
        update_fig(self.fig_chiplot)
        logger.info("One pixel = %.3e deg", resolution)

    def postProcess(self):
        """
        Common part: shows the result of the azimuthal integration in 1D and 2D
        """
        if self.geoRef is None:
            self.refine()
        if "wavelength" not in self.fixed:
            self.peakPicker.points.setWavelength_change2th(self.geoRef.wavelength)
        self.peakPicker.points.save(self.basename + ".npt")
        self.geoRef.save(self.basename + ".poni")
        self.geoRef.mask = self.mask
        self.geoRef.del_ttha()
        self.geoRef.del_dssa()
        self.geoRef.del_chia()
        t0 = time.perf_counter()
        _tth = self.geoRef.twoThetaArray(self.peakPicker.shape)
        t1 = time.perf_counter()
        _dsa = self.geoRef.solidAngleArray(self.peakPicker.shape)
        t2 = time.perf_counter()
        self.geoRef.chiArray(self.peakPicker.shape)
        t2a = time.perf_counter()
        self.geoRef.corner_array(self.peakPicker.shape, units.TTH_DEG,
                                 scale=False)
        t2b = time.perf_counter()
        if self.gui:
            if self.fig_integrate is None:
                self.fig_integrate = pylab.plt.figure()
                self.ax_xrpd_1d = self.fig_integrate.add_subplot(1, 2, 1)
                self.ax_xrpd_2d = self.fig_integrate.add_subplot(1, 2, 2)
            else:
                self.ax_xrpd_1d.cla()
                self.ax_xrpd_2d.cla()
                update_fig(self.fig_integrate)

        t3 = time.perf_counter()
        res1 = self.geoRef.integrate1d_ng(self.peakPicker.data, self.nPt_1D,
                                          filename=self.basename + ".xy",
                                          unit=self.unit,
                                          polarization_factor=self.polarization_factor,
                                          method=self.integrator_method,
                                          error_model=self.error_model)
        t4 = time.perf_counter()
        res2 = self.geoRef.integrate2d(self.peakPicker.data,
                                          self.nPt_2D_rad, self.nPt_2D_azim,
                                          filename=self.basename + ".azim",
                                          unit=self.unit,
                                          polarization_factor=self.polarization_factor,
                                          method=self.integrator_method,
                                          error_model=self.error_model)
        t5 = time.perf_counter()
        logger.info(os.linesep.join(["Timings (%s):" % self.integrator_method,
                                     " * two theta array generation %.3fs" % (t1 - t0),
                                     " * diff Solid Angle           %.3fs" % (t2 - t1),
                                     " * chi array generation       %.3fs" % (t2a - t2),
                                     " * corner coordinate array    %.3fs" % (t2b - t2a),
                                     " * 1D Azimuthal integration   %.3fs" % (t4 - t3),
                                     " * 2D Azimuthal integration   %.3fs" % (t5 - t4)]))

        if self.gui:
            self.ax_xrpd_1d.plot(res1.radial, res1.intensity)
            # GF: Add vertical line for each used calibration ring:
            xValues = None
            twoTheta = numpy.array([i for i in self.peakPicker.points.calibrant.get_2th() if i])  # in radian
            if self.unit == units.TTH_DEG:
                xValues = numpy.rad2deg(twoTheta)
            elif self.unit == units.TTH_RAD:
                xValues = twoTheta
            elif self.unit == units.Q_A:
                xValues = (4.e-10 * numpy.pi / self.wavelength) * numpy.sin(.5 * twoTheta)
            elif self.unit == units.Q_NM:
                xValues = (4.e-9 * numpy.pi / self.wavelength) * numpy.sin(.5 * twoTheta)
            elif self.unit == units.R_MM:
                # GF: correct formula?
                dBeamCentre = self.geoRef.getFit2D()["directDist"]  # in mm!!
                xValues = dBeamCentre * numpy.tan(twoTheta)
            else:
                logger.warning("Unknown unit %s, do not plot calibration rings", self.unit)
            if xValues is not None:
                for x in xValues:
                    line = matplotlib.lines.Line2D([x, x], self.ax_xrpd_1d.axis()[2:4],
                                                   color='red', linestyle='--')
                    self.ax_xrpd_1d.add_line(line)
            self.ax_xrpd_1d.set_title("1D integration")
            self.ax_xrpd_1d.set_xlabel(self.unit.label)
            self.ax_xrpd_1d.set_ylabel("Intensity")
            img = res2.intensity
            pos_rad = res2.radial
            pos_azim = res2.azimuthal
            self.ax_xrpd_2d.imshow(numpy.log(img - img.min() + 1e-3), origin="lower",
                                   extent=[pos_rad.min(), pos_rad.max(), pos_azim.min(), pos_azim.max()],
                                   aspect="auto")
            self.ax_xrpd_2d.set_title("2D regrouping")
            self.ax_xrpd_2d.set_xlabel(self.unit.label)
            self.ax_xrpd_2d.set_ylabel(r"Azimuthal angle $\chi$ ($^{o}$)")
            if not gui_utils.main_loop:
                self.fig_integrate.show()
            update_fig(self.fig_integrate)

    def validate_calibration(self):
        """
        Validate the calibration and calculate the offset in the diffraction image
        """
        if not self.check_calib:
            self.check_calib = CheckCalib()
        if self.geoRef:
            self.ai.setPyFAI(**self.geoRef.getPyFAI())
            self.ai.wavelength = self.geoRef.wavelength
        self.check_calib.ai = self.ai
        self.check_calib.img = self.peakPicker.data
        self.check_calib.mask = self.peakPicker.mask
        self.check_calib.wavelength = self.check_calib.wavelength
        self.check_calib.integrate()
        self.check_calib.rebuild()
        self.check_calib.show()

    def validate_center(self, slices=36):
        """
        Validate the position of the center by cross-correlating two spectra 180 deg appart.
        Output values are in micron.

        Designed for orthogonal setup with centered beam...

        :param slices: number of slices on which perform
        """
        if slices <= 0:
            logger.warning("The number of slices should be strictly positive")
            slices = 1
        if slices % 2 == 1:
            logger.warning("Validate assumes the number of slices is even. adding one")
            slices += 1
        half_slices = slices // 2
        npt = round_fft(int(math.sqrt(self.peakPicker.data.shape[0] ** 2 + self.peakPicker.data.shape[1] ** 2) + 1))

        if self.geoRef:
            self.ai.setPyFAI(**self.geoRef.getPyFAI())
            self.ai.wavelength = self.geoRef.wavelength
        logger.info("Performing autocorrelation on %sx%s, Fourier analysis may take some time", slices, npt)
        img, tth, chi = self.ai.integrate2d(self.peakPicker.data, npt, slices, azimuth_range=(-180, 180), unit="r_mm", method="splitpixel")
        ft = numpy.fft.fft(img, npt * 2, axis=-1)
        crosscor = numpy.fft.ifft(ft[:half_slices,:] * (ft[half_slices:,:].conj()), axis=-1).real
        centered = numpy.empty_like(crosscor)
        centered[:,:npt] = crosscor[:, npt:]
        centered[:, npt:] = crosscor[:,:npt]

        center = numpy.zeros(slices)  # in micron
        dr = (tth[1] - tth[0]) * 1000.0  # ouput in r(mm) -> micron

        # sub-bin precision obtained by second order expantion of peak
        range_half_slices = range(half_slices)
        x0 = centered.argmax(axis=-1)
        f_x = centered[(range_half_slices, x0)]
        f_xm1 = centered[(range_half_slices, x0 - 1)]
        f_xp1 = centered[(range_half_slices, x0 + 1)]
        f_prime = (f_xp1 - f_xm1) / 2.0
        f_second = (f_xp1 + f_xm1 - 2.0 * f_x)
        dx = -f_prime / f_second
        if (abs(dx) >= 0.5).any():
            msk = abs(dx) > 1
            logger.info("Correction is important ! %s", msk)
            dx[msk] = 0.0
        center[half_slices:] = (x0 + dx - npt) * dr
        center[:half_slices] = -center[half_slices:]
        if self.gui:
            if self.fig_center:
                self.fig_center.clf()
            else:
                self.fig_center = pylab.plt.figure()
            self.ax_center = self.fig_center.add_subplot(1, 1, 1)
            self.ax_center.set_xlim(-180, 180)
            self.ax_center.set_xticks(numpy.linspace(-180, 180, 9))
            self.ax_center.set_xlabel("Azimuthal angle $\chi$ ($^o$)")
            self.ax_center.set_ylabel(r"Error of the center position along radius ($\mu$m)")
            self.ax_center.set_title("Center plot")
            self.ax_center.plot(chi, center, label="From pattern cross-correlation")
            self.fig_center.show()
            update_fig(self.fig_center)

    def set_data(self, data):
        """
        call-back function for the peak-picker
        """
        self.data = data
        if not self.weighted:
            self.data = numpy.array(self.data)[:,:-1]
        self.refine()

    def reset_geometry(self, how="center", refine=False):
        """
        Reset the geometry: no tilt in all cases

        :param how: multiple options
            * center: set the PONI at the center of the detector
            * ring: center the poni at the middle of the inner-most ring
            * best: try both option and keeps the best (this option is not available)
        :param refine: launch the refinement (argument not used)
        """
        if how not in ["center", "ring"]:  # ,"best"]:
            logger.warning("unknow geometry reset method: %s, fall back on detector center", how)
            how = "center"
        if self.data is None:
            logger.warning("No datapoint: fall back on detector center")
            how = "center"
        # this is true for all:
        self.ai.rot1 = 0.0
        self.ai.rot2 = 0.0
        self.ai.rot3 = 0.0

        if how == "ring":
            inner_ring = min(set(i[2] for i in self.data))
            print("inner ring: %s" % inner_ring)
            data = numpy.array([[i[0], i[1]] for i in self.data if i[2] == inner_ring])
            center = data.mean(axis=0)
            self.ai.poni1, self.ai.poni2 = data.mean(axis=0)
            tth = self.calibrant.get_2th()[int(inner_ring)]
            dist = (data - center)
            d = numpy.sqrt(dist[:, 0] ** 2 + dist[:, 1] ** 2).mean()
            self.ai.dist = d / numpy.tan(tth)
        elif how == "center":
            self.ai.dist = 0.1
            try:
                p1, p2, _p3 = self.detector.calc_cartesian_positions()
                self.ai.poni1 = p1.max() / 2.0
                self.ai.poni2 = p2.max() / 2.0
            except Exception as err:
                logger.warning(err)
                self.ai.poni1 = self.detector.pixel1 * (self.peakPicker.shape[0] / 2.)
                self.ai.poni2 = self.detector.pixel2 * (self.peakPicker.shape[1] / 2.)

        if self.geoRef:
            # reset geoRef object
            self.geoRef.set_dist_min(0)
            self.geoRef.set_dist_max(100)
            self.geoRef.set_dist(self.ai.dist)

            self.geoRef.set_poni1_min(-10.0 * self.ai.poni1)
            self.geoRef.set_poni1_max(10.0 * self.ai.poni1)
            self.geoRef.set_poni1(self.ai.poni1)

            self.geoRef.set_poni2_min(-10.0 * self.ai.poni2)
            self.geoRef.set_poni2_max(10.0 * self.ai.poni2)
            self.geoRef.set_poni2(self.ai.poni2)

            self.geoRef.set_rot1_min(-math.pi)
            self.geoRef.set_rot1_max(math.pi)
            self.geoRef.set_rot1(self.ai.rot1)

            self.geoRef.set_rot2_min(-math.pi)
            self.geoRef.set_rot2_max(math.pi)
            self.geoRef.set_rot2(self.ai.rot2)

            self.geoRef.set_rot3_min(-math.pi)
            self.geoRef.set_rot3_max(math.pi)
            self.geoRef.set_rot3(self.ai.rot3)

################################################################################
# Calibration
################################################################################


class Calibration(AbstractCalibration):
    """
    class doing the calibration of frames
    """

    def __init__(self, dataFiles=None, darkFiles=None, flatFiles=None, pixelSize=None,
                 splineFile=None, detector=None, gaussianWidth=None,
                 wavelength=None, calibrant=None):
        """
        Constructor for calibration:

        :param dataFiles: list of filenames containing data images
        :param darkFiles: list of filenames containing dark current images
        :param flatFiles: list of filenames containing flat images
        :param pixelSize: size of the pixel in meter as 2 tuple
        :param splineFile: file containing the distortion of the taper
        :param detector: Detector name or instance
        :param wavelength: radiation wavelength in meter
        :param calibrant: pyFAI.calibrant.Calibrant instance

        """
        AbstractCalibration.__init__(self, dataFiles=dataFiles,
                                     darkFiles=darkFiles,
                                     flatFiles=flatFiles,
                                     pixelSize=pixelSize,
                                     splineFile=splineFile,
                                     detector=detector,
                                     calibrant=calibrant,
                                     wavelength=wavelength)
        self.gaussianWidth = gaussianWidth
        self.labelPattern = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]

    def __repr__(self):
        return AbstractCalibration.__repr__(self) + \
            "%sgaussian= %s" % (os.linesep, self.gaussianWidth)

    def parse(self):
        """
        parse options from command line
        """
        description = """Calibrate the diffraction setup geometry based on Debye-Sherrer rings images
without a priori knowledge of your setup.
You will need to provide a calibrant or a "d-spacing" file containing the spacing of Miller plans in
Angstrom (in decreasing order).
%s
or search in the American Mineralogist database:
http://rruff.geo.arizona.edu/AMS/amcsd.php
The --calibrant option is mandatory !""" % str(CALIBRANT_FACTORY)

        epilog = """The output of this program is a "PONI" file containing the detector description
and the 6 refined parameters (distance, center, rotation) and wavelength.
An 1D and 2D diffraction patterns are also produced. (.dat and .azim files)
        """
        usage = "pyFAI-calib [options] -w 1 -D detector -c calibrant.D imagefile.edf"
        self.configure_parser(usage=usage, description=description, epilog=epilog)  # common
        self.parser.add_argument("-r", "--reconstruct", dest="reconstruct",
                                 help="Reconstruct image where data are masked or <0  (for Pilatus "
                                 "detectors or detectors with modules)",
                                 action="store_true", default=False)

        self.parser.add_argument("-g", "--gaussian", dest="gaussian",
                                 help="""Size of the gaussian kernel.
Size of the gap (in pixels) between two consecutive rings, by default 100
Increase the value if the arc is not complete;
decrease the value if arcs are mixed together.""", default=None)
        self.parser.add_argument("--square", dest="square", action="store_true",
                                 help="Use square kernel shape for neighbor search instead of diamond shape",
                                 default=False)
        self.parser.add_argument("-p", "--pixel", dest="pixel",
                                 help="size of the pixel in micron", default=None)

        (options, _) = self.analyse_options()
        # Analyse remaining aruments and options
        self.reconstruct = options.reconstruct
        self.gaussianWidth = options.gaussian
        if options.square:
            self.labelPattern = [[1] * 3] * 3
        else:
            self.labelPattern = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]

        if options.pixel is not None:
            self.get_pixelSize(options.pixel)

    def preprocess(self):
        """
        do dark, flat correction thresholding, ...
        """
        AbstractCalibration.preprocess(self)

        if self.gaussianWidth is not None:
            self.peakPicker.massif.valley_size = self.gaussianWidth
        else:
            self.peakPicker.massif.initValleySize()
        if self.gui:
            self.peakPicker.gui(log=True, maximize=True, pick=True)
            update_fig(self.peakPicker.fig)

    def gui_peakPicker(self):
        if self.peakPicker is None:
            self.preprocess()
#        self.peakPicker.gui(True)
        if os.path.isfile(self.pointfile):
            self.peakPicker.load(self.pointfile)
        if self.gui:
            update_fig(self.peakPicker.fig)
#        self.peakPicker.finish(self.pointfile, callback=self.set_data)
        self.set_data(self.peakPicker.finish(self.pointfile))
#        input("Please press enter when you are happy with your selection" + os.linesep)
#        while self.data is None:
#            update_fig(self.peakPicker.fig)
#            time.sleep(0.1)

    def initgeoRef(self):
        """
        Tries to initialise the GeometryRefinement (dist, poni, rot)
        Returns a dictionary of key value pairs
        """
        defaults = {"dist": 0.1, "poni1": 0.0, "poni2": 0.0,
                    "rot1": 0.0, "rot2": 0.0, "rot3": 0.0}
        if self.detector:
            try:
                p1, p2, _p3 = self.detector.calc_cartesian_positions()
                defaults["poni1"] = p1.max() / 2.
                defaults["poni2"] = p2.max() / 2.
            except Exception as err:
                logger.warning(err)
        if self.ai:
            for key in defaults.keys():  # not PARAMETERS which holds wavelength
                val = getattr(self.ai, key, None)
                if val is not None:
                    defaults[key] = val
        return defaults

    def refine(self):
        """
        Contains the geometry refinement part specific to Calibration
        Sets up the initial guess when starting pyFAI-calib
        """
        # First attempt
        defaults = self.initgeoRef()
        self.geoRef = GeometryRefinement(self.data,
                                         detector=self.detector,
                                         wavelength=self.wavelength,
                                         calibrant=self.calibrant,
                                         **defaults)
        self.geoRef.refine2(1000000, fix=self.fixed)
        scor = self.geoRef.chi2()
        pars = [getattr(self.geoRef, p) for p in self.PARAMETERS]

        scores = [(scor, pars), ]

        # Second attempt
        defaults = self.initgeoRef()
        self.geoRef = GeometryRefinement(self.data,
                                         detector=self.detector,
                                         wavelength=self.wavelength,
                                         calibrant=self.calibrant,
                                         **defaults)
        self.geoRef.guess_poni()
        self.geoRef.refine2(1000000, fix=self.fixed)
        scor = self.geoRef.chi2()
        pars = [getattr(self.geoRef, p) for p in self.PARAMETERS]

        scores.append((scor, pars))

        # Third attempt (can be from when a program bombed last time)
        paramfile = self.basename + ".poni"
        if os.path.isfile(paramfile):
            self.geoRef.load(paramfile)
            if self.wavelength:
                try:
                    old_wl = self.geoRef.wavelength
                except Exception as err:
                    logger.warning(err)
                else:
                    logger.warning("Overwriting wavelength from PONI file (%s) "
                                   "with the one from command line (%s)" %
                                   (old_wl, self.wavelength))
                self.geoRef.wavelength = self.wavelength
            if self.detector:
                gr_det = str(self.geoRef.detector)
                nw_det = str(self.detector)
                if gr_det != nw_det:
                    logger.warning("Overwriting detector from PONI file: %s%s "
                                   "with the one from command line %s%s" %
                                   (os.linesep, gr_det, os.linesep, nw_det))
                    self.geoRef.detector = self.detector

        # Third attempt
        self.geoRef.refine2(1000000, fix=self.fixed)
        scor = self.geoRef.chi2()
        pars = [getattr(self.geoRef, p) for p in self.PARAMETERS]

        scores.append((scor, pars))

        # Choose the best scoring method: At this point we might also ask
        # a user to just type the numbers in?
        scores.sort()
        scor, pars = scores[0]
        for parval, parname in zip(pars, self.PARAMETERS):
            setattr(self.geoRef, parname, parval)

        # Now continue as before
        AbstractCalibration.refine(self)

################################################################################
# Recalibration
################################################################################


class Recalibration(AbstractCalibration):
    """
    class doing the re-calibration of frames
    """

    def __init__(self, dataFiles=None, darkFiles=None, flatFiles=None, pixelSize=None,
                 splineFile=None, detector=None, wavelength=None, calibrant=None):
        """
        Constructor for Recalibration:

        :param dataFiles: list of filenames containing data images
        :param darkFiles: list of filenames containing dark current images
        :param flatFiles: list of filenames containing flat images
        :param pixelSize: size of the pixel in meter as 2 tuple
        :param splineFile: file containing the distortion of the taper
        :param detector: Detector name or instance
        :param wavelength: radiation wavelength in meter
        :param calibrant: pyFAI.calibrant.Calibrant instance
        """
        AbstractCalibration.__init__(self, dataFiles=dataFiles,
                                     darkFiles=darkFiles,
                                     flatFiles=flatFiles,
                                     pixelSize=pixelSize,
                                     splineFile=splineFile,
                                     detector=detector,
                                     wavelength=wavelength,
                                     calibrant=calibrant)

    def parse(self):
        """
        parse options from command line
        """
        description = """Calibrate the diffraction setup geometry based on Debye-Sherrer rings images
with a priori knowledge of your setup (an input PONI-file).
You will need to provide a calibrant or a "d-spacing" file containing the spacing of Miller plans in
Angstrom (in decreasing order).
%s
or search in the American Mineralogist database:
http://rruff.geo.arizona.edu/AMS/amcsd.php
The --calibrant option is mandatory !
""" % str(CALIBRANT_FACTORY)

        epilog = """The main difference with pyFAI-calib is the way control-point hence Debye-Sherrer
rings are extracted. While pyFAI-calib relies on the contiguity of a region of peaks
called massif; pyFAI-recalib knows approximatly the geometry and is able to select
the region where the ring should be. From this region it selects automatically
the various peaks; making pyFAI-recalib able to run without graphical interface and
without human intervention (--no-gui and --no-interactive options).


Note that `pyFAI-recalib` program is obsolete as the same functionality is
available from within pyFAI-calib, using the `recalib` command in the
refinement process.
Two option are available for recalib: the numbe of rings to extract (similar to the -r option of this program)
and a new option which lets you choose between the original `massif` algorithm and newer ones like `blob` and `watershed` detection.
        """
        usage = "pyFAI-recalib [options] -i ponifile -w 1 -c calibrant.D imagefile.edf"
        self.configure_parser(usage=usage, description=description, epilog=epilog)

        self.parser.add_argument("-r", "--ring", dest="max_rings", type=int,
                                 help="maximum number of rings to extract. Default: all accessible", default=None)
        self.parser.add_argument("-k", "--keep", dest="keep",
                                 help="Keep existing control point and append new",
                                 default=False, action="store_true")

        options = self.parser.parse_args()
        args = options.args
        # Analyse aruments and options
        if (not options.poni) or (not os.path.isfile(options.poni)):
            logger.error("You should provide a PONI file as starting point !!")
        else:
            self.ai = AzimuthalIntegrator.sload(options.poni)
        if self.wavelength:
            self.ai.wavelength = self.wavelength
        self.max_rings = options.max_rings
        self.detector = self.ai.detector
        self.keep = options.keep
        self.analyse_options(options, args)

    def read_dSpacingFile(self):
        """Read the name of the file with d-spacing"""
        AbstractCalibration.read_dSpacingFile(self, verbose=False)

    def preprocess(self):
        """
        do dark, flat correction thresholding, ...
        """
        AbstractCalibration.preprocess(self)

        if self.gui:
            self.peakPicker.gui(log=True, maximize=True, pick=False)
            update_fig(self.peakPicker.fig)

    def refine(self):
        """
        Contains the geometry refinement part specific to Recalibration
        """

        self.geoRef = GeometryRefinement(self.data, dist=self.ai.dist, poni1=self.ai.poni1,
                                         poni2=self.ai.poni2, rot1=self.ai.rot1,
                                         rot2=self.ai.rot2, rot3=self.ai.rot3,
                                         detector=self.ai.detector, calibrant=self.calibrant,
                                         wavelength=self.wavelength)
        self.ai = self.geoRef
        self.geoRef.set_tolerance(10)
        AbstractCalibration.refine(self)


class MultiCalib(object):

    def __init__(self, dataFiles=None, darkFiles=None, flatFiles=None, pixelSize=None, splineFile=None, detector=None):
        """
        """
        self.dataFiles = dataFiles or []
        self.darkFiles = darkFiles or []
        self.flatFiles = flatFiles or []
        self.data = {}

        self.detector = get_detector(detector, dataFiles)

        if splineFile and os.path.isfile(splineFile):
            self.detector.splineFile = os.path.abspath(splineFile)
        if pixelSize:
            if "__len__" in dir(pixelSize) and len(pixelSize) >= 2:
                self.detector.pixel1 = float(pixelSize[0])
                self.detector.pixel2 = float(pixelSize[1])
            else:
                self.detector.pixel1 = self.detector.pixel2 = float(pixelSize)
        self.cutBackground = None
        self.outfile = "merged.edf"
        self.peakPicker = "blob"
        self.basename = None
        self.geoRef = None
#        self.reconstruct = False
        self.mask = None
        self.max_iter = 1000
        self.filter = "mean"
        self.saturation = 0.1
        self.calibrant = None
        self.wavelength = None
        self.weighted = False
        self.polarization_factor = 0
        self.results = {}
        self.gui = True
        self.interactive = True
        self.poni1 = None
        self.poni2 = None
        self.dist = None
        self.fixed = FixedParameters()
        self.max_rings = None
        self.rot1 = 0.0
        self.rot2 = 0.0
        self.rot3 = 0.0

    def __repr__(self):
        lst = ["Multi-Calibration object:",
               "data= " + ", ".join(self.dataFiles),
               "dark= " + ", ".join(self.darkFiles),
               "flat= " + ", ".join(self.flatFiles)]
        lst.append(self.detector.__repr__())
#        lst.append("gaussian= %s" % self.gaussianWidth)
        return os.linesep.join(lst)

    def parse(self, exe=None, description=None, epilog=None):
        """
        parse options from command line
        :param exe: name of the program (MX-calibrate)
        :param description: Description of the program
        """
        if exe is None:
            exe = "MX-Calibrate"
            usage = "%s -w 1.54 -c CeO2 file1.cbf file2.cbf ..." % exe
            version = "%s from pyFAI version %s: %s" % (exe, PyFAI_VERSION, PyFAI_DATE)
            description = """
        Calibrate automatically a set of frames taken at various sample-detector distance.
        Return the linear regression of the fit in funtion of the sample-setector distance.
        """
            epilog = """This tool has been developed for ESRF MX-beamlines where an acceptable calibration is
        usually present is the header of the image. PyFAI reads it and does a "recalib" on
        each of them before exporting a linear regression of all parameters versus this distance.
        """
        else:
            description = description or ""
            epilog = epilog or ""
        parser = ArgumentParser(usage=usage, description=description, epilog=epilog)
        parser.add_argument("-V", "--version", action='version', version=version)
        parser.add_argument("args", metavar="FILE", help="List of files to calibrate", nargs='+')
#        parser.add_argument("-V", "--version", dest="version", action="store_true",
#                          help="print version of the program and quit", metavar="FILE", default=False)
#        parser.add_argument("-o", "--out", dest="outfile",
#                          help="Filename where processed image is saved", metavar="FILE", default="merged.edf")
        parser.add_argument("-v", "--verbose",
                            action="store_true", dest="debug", default=False,
                            help="switch to debug/verbose mode")
#        parser.add_argument("-g", "--gaussian", dest="gaussian", help="""Size of the gaussian kernel.
# Size of the gap (in pixels) between two consecutive rings, by default 100
# Increase the value if the arc is not complete;
# decrease the value if arcs are mixed together.""", default=None)
#        parser.add_argument("-c", "--square", dest="square", action="store_true",
#                      help="Use square kernel shape for neighbor search instead of diamond shape", default=False)
        parser.add_argument("-c", "--calibrant", dest="calibrant", metavar="FILE",
                            help="file containing d-spacing of the calibrant reference sample (MANDATORY)", default=None)
        parser.add_argument("-w", "--wavelength", dest="wavelength", type=float,
                            help="wavelength of the X-Ray beam in Angstrom", default=None)
        parser.add_argument("-e", "--energy", dest="energy", type=float,
                            help="energy of the X-Ray beam in keV (hc=%skeV.A)" % hc, default=None)
        parser.add_argument("-P", "--polarization", dest="polarization_factor",
                            type=float, default=0.0,
                            help="polarization factor, from -1 (vertical) to +1 (horizontal), default is 0, synchrotrons are around 0.95")
        parser.add_argument("-b", "--background", dest="background",
                            help="Automatic background subtraction if no value are provided", default=None)
        parser.add_argument("-d", "--dark", dest="dark",
                            help="list of dark images to average and subtract", default=None)
        parser.add_argument("-f", "--flat", dest="flat",
                            help="list of flat images to average and divide", default=None)
#        parser.add_argument("-r", "--reconstruct", dest="reconstruct",
#                      help="Reconstruct image where data are masked or <0  (for Pilatus detectors or detectors with modules)",
#                      action="store_true", default=False)
        parser.add_argument("-s", "--spline", dest="spline",
                            help="spline file describing the detector distortion", default=None)
        parser.add_argument("-p", "--pixel", dest="pixel",
                            help="size of the pixel in micron", default=None)
        parser.add_argument("-D", "--detector", dest="detector_name",
                            help="Detector name (instead of pixel size+spline)", default=None)
        parser.add_argument("-m", "--mask", dest="mask",
                            help="file containing the mask (for image reconstruction)", default=None)
#        parser.add_argument("-n", "--npt", dest="npt",
#                      help="file with datapoints saved", default=None)
        parser.add_argument("--filter", dest="filter",
                            help="select the filter, either mean(default), max or median",
                            default="mean")
        parser.add_argument("--saturation", dest="saturation",
                            help="consider all pixel>max*(1-saturation) as saturated and reconstruct them",
                            default=0.1, type=float)
        parser.add_argument("-r", "--ring", dest="max_rings", type=float,
                            help="maximum number of rings to extract", default=None)

        parser.add_argument("--weighted", dest="weighted",
                            help="weight fit by intensity",
                            default=False, action="store_true")
        parser.add_argument("-l", "--distance", dest="distance", type=float,
                            help="sample-detector distance in millimeter", default=None)

        parser.add_argument("--tilt", dest="tilt",
                            help="Allow initially detector tilt to be refined (rot1, rot2, rot3). Default: Activated",
                            default=None, action="store_true")
        parser.add_argument("--no-tilt", dest="tilt",
                            help="Deactivated tilt refinement and set all rotation to 0", default=None, action="store_false")

        parser.add_argument("--dist", dest="dist", type=float,
                            help="sample-detector distance in meter", default=None)
        parser.add_argument("--poni1", dest="poni1", type=float,
                            help="poni1 coordinate in meter", default=None)
        parser.add_argument("--poni2", dest="poni2", type=float,
                            help="poni2 coordinate in meter", default=None)
        parser.add_argument("--rot1", dest="rot1", type=float,
                            help="rot1 in radians", default=None)
        parser.add_argument("--rot2", dest="rot2", type=float,
                            help="rot2 in radians", default=None)
        parser.add_argument("--rot3", dest="rot3", type=float,
                            help="rot3 in radians", default=None)

        parser.add_argument("--fix-dist", dest="fix_dist",
                            help="fix the distance parameter", default=None, action="store_true")
        parser.add_argument("--free-dist", dest="fix_dist",
                            help="free the distance parameter", default=None, action="store_false")

        parser.add_argument("--fix-poni1", dest="fix_poni1",
                            help="fix the poni1 parameter", default=None, action="store_true")
        parser.add_argument("--free-poni1", dest="fix_poni1",
                            help="free the poni1 parameter", default=None, action="store_false")

        parser.add_argument("--fix-poni2", dest="fix_poni2",
                            help="fix the poni2 parameter", default=None, action="store_true")
        parser.add_argument("--free-poni2", dest="fix_poni2",
                            help="free the poni2 parameter", default=None, action="store_false")

        parser.add_argument("--fix-rot1", dest="fix_rot1",
                            help="fix the rot1 parameter", default=None, action="store_true")
        parser.add_argument("--free-rot1", dest="fix_rot1",
                            help="free the rot1 parameter", default=None, action="store_false")

        parser.add_argument("--fix-rot2", dest="fix_rot2",
                            help="fix the rot2 parameter", default=None, action="store_true")
        parser.add_argument("--free-rot2", dest="fix_rot2",
                            help="free the rot2 parameter", default=None, action="store_false")

        parser.add_argument("--fix-rot3", dest="fix_rot3",
                            help="fix the rot3 parameter", default=None, action="store_true")
        parser.add_argument("--free-rot3", dest="fix_rot3",
                            help="free the rot3 parameter", default=None, action="store_false")

        parser.add_argument("--fix-wavelength", dest="fix_wavelength",
                            help="fix the wavelength parameter", default=True, action="store_true")
        parser.add_argument("--free-wavelength", dest="fix_wavelength",
                            help="free the wavelength parameter", default=True, action="store_false")

        parser.add_argument("--no-gui", dest="gui",
                            help="force the program to run without a Graphical interface",
                            default=True, action="store_false")
        parser.add_argument("--gui", dest="gui",
                            help="force the program to run with a Graphical interface",
                            default=True, action="store_true")

        parser.add_argument("--no-interactive", dest="interactive",
                            help="force the program to run and exit without prompting for refinements",
                            default=True, action="store_false")
        parser.add_argument("--interactive", dest="interactive",
                            help="force the program to prompt for refinements",
                            default=True, action="store_true")
        parser.add_argument("--peak-picker", dest="peakPicker",
                            help="Uses the 'massif', 'blob' or 'watershed' peak-picker algorithm (default: blob)",
                            default="blob", type=str)
        options = parser.parse_args()

        # Analyse aruments and options
        if options.debug:
            logger.setLevel(logging.DEBUG)
        if options.background is not None:
            try:
                self.cutBackground = float(options.background)
            except Exception:
                self.cutBackground = True
        if options.dark:
            self.darkFiles = [f for f in options.dark.split(",") if os.path.isfile(f)]
        if options.flat:
            self.flatFiles = [f for f in options.flat.split(",") if os.path.isfile(f)]
        if options.mask and os.path.isfile(options.mask):
            self.mask = fabio.open(options.mask).data

        if options.detector_name:
            self.detector = get_detector(options.detector_name, options.args)
        if options.spline:
            if os.path.isfile(options.spline):
                self.detector.splineFile = os.path.abspath(options.spline)
            else:
                logger.error("Unknown spline file %s", options.spline)
        if options.pixel is not None:
            self.get_pixelSize(options.pixel)
        self.filter = options.filter
        self.saturation = options.saturation
        if options.wavelength:
            self.wavelength = 1e-10 * options.wavelength
        elif options.energy:
            self.wavelength = 1e-10 * hc / options.energy
        if not options.calibrant:
            logger.error("The calibrant is mandatory: please use the -c option")
        self.calibrant = options.calibrant
        self.polarization_factor = options.polarization_factor
        self.gui = options.gui
        self.interactive = options.interactive
        self.max_rings = options.max_rings

        self.fixed = FixedParameters()
        if options.tilt is not None:
            for key in ["rot1", "rot2", "rot3"]:
                self.fixed.add_or_discard(key, not(options.tilt))
        self.fixed.add_or_discard("dist", options.fix_dist)
        self.fixed.add_or_discard("poni1", options.fix_poni1)
        self.fixed.add_or_discard("poni2", options.fix_poni2)
        self.fixed.add_or_discard("rot1", options.fix_rot1)
        self.fixed.add_or_discard("rot2", options.fix_rot2)
        self.fixed.add_or_discard("rot3", options.fix_rot3)
        self.fixed.add_or_discard("wavelength", options.fix_wavelength)

        if options.distance:
            self.dist = 1e-3 * float(options.distance)
        if options.dist:
            self.dist = float(options.dist)
        if options.poni1:
            self.poni1 = float(options.poni1)
        if options.poni2:
            self.poni2 = float(options.poni2)
        if options.rot1:
            self.rot1 = float(options.rot1)
        if options.rot2:
            self.rot2 = float(options.rot2)
        if options.rot3:
            self.rot3 = float(options.rot3)

        self.dataFiles = [f for f in options.args if os.path.isfile(f)]
        if not self.dataFiles:
            raise RuntimeError("Please provide some calibration images ... "
                               "if you want to analyze them. Try also the --help option to see all options!")
        self.weighted = options.weighted
        if options.peakPicker.lower() in PeakPicker.VALID_METHODS:
            self.peakPicker = options.peakPicker.lower()

    def get_pixelSize(self, ans):
        """convert a comma separated sting into pixel size"""
        sp = ans.split(",")
        if len(sp) >= 2:
            try:
                pixelSizeXY = [float(i) * 1e-6 for i in sp[:2]]
            except Exception:
                logger.error("error in reading pixel size_2")
                return
        elif len(sp) == 1:
            px = sp[0]
            try:
                pixelSizeXY = [float(px) * 1e-6, float(px) * 1e-6]
            except Exception:
                logger.error("error in reading pixel size_1")
                return
        else:
            logger.error("error in reading pixel size_0")
            return
        self.detector.pixel1 = pixelSizeXY[1]
        self.detector.pixel2 = pixelSizeXY[0]

    def read_pixelsSize(self):
        """Read the pixel size from prompt if not available"""
        if (self.detector.pixel1 is None) and (self.detector.splineFile is None):
            pixelSize = [15, 15]
            ans = input("Please enter the pixel size (in micron, comma separated X, Y "
                        "i.e. %.2e,%.2e) or a spline file: " % tuple(pixelSize)).strip()
            if os.path.isfile(ans):
                self.detector.splineFile = ans
            else:
                self.get_pixelSize(ans)

    def read_dSpacingFile(self):
        """Read the name of the calibrant or the file with d-spacing"""
        if self.calibrant in CALIBRANT_FACTORY:
            self.calibrant = CALIBRANT_FACTORY(self.calibrant)
        elif os.path.isfile(self.calibrant):
            self.calibrant = Calibrant(filename=self.calibrant)
        else:
            comments = ["MX-calibrate has changed !!!",
                        "Instead of entering the 2theta value, which was tedious,"
                        "the program takes a calibrant as in input "
                        "(either a reference one like Ceo2, either a "
                        "d-spacing file with inter planar distance in Angstrom)",
                        "and an associated wavelength", ""
                        "You will be asked to enter the ring number, "
                        "which is usually a simpler than the 2theta value."]
            print(os.linesep.join(comments))
            ans = ""
            while not self.calibrant:
                ans = input("Please enter the name of the calibrant"
                            " or the file containing the d-spacing:\t").strip()
                if ans in CALIBRANT_FACTORY:
                    self.calibrant = CALIBRANT_FACTORY(ans)
                elif os.path.isfile(ans):
                    self.calibrant = Calibrant(filename=ans)

    def read_wavelength(self):
        """Read the wavelength"""
        while not self.wavelength:
            ans = input("Please enter wavelength in Angstrom:\t").strip()
            try:
                self.wavelength = 1e-10 * float(ans)
            except:
                self.wavelength = None

    def process(self):
        """

        """
        self.dataFiles.sort()
        for fn in self.dataFiles:
            fabimg = fabio.open(fn)
            wavelength = self.wavelength
            dist = self.dist
            if self.poni2:
                centerX = self.poni2 / self.detector.pixel2
            else:
                centerX = None
            if self.poni1:
                centerY = self.poni1 / self.detector.pixel1
            else:
                centerY = None
            if "_array_data.header_contents" in fabimg.header:
                headers = fabimg.header["_array_data.header_contents"].lower().split()
                if "detector_distance" in headers:
                    dist = float(headers[headers.index("detector_distance") + 1])
                if "wavelength" in headers:
                    wavelength = float(headers[headers.index("wavelength") + 1]) * 1e-10
                if "beam_xy" in headers:
                    centerX = float(headers[headers.index("beam_xy") + 1][1:-1])
                    centerY = float(headers[headers.index("beam_xy") + 2][:-1])
            if dist is None:
                digits = ""
                for i in os.path.basename(fn):
                    if i.isdigit() and not digits:
                        digits += i
                    elif i.isdigit():
                        digits += i
                    elif not i.isdigit() and digits:
                        break
                dist = int(digits) * 0.001
            if centerX is None:
                centerX = fabimg.data.shape[1] // 2
            if centerY is None:
                centerY = fabimg.data.shape[0] // 2
            self.results[fn] = {"wavelength": wavelength, "dist": dist}
            rec = Recalibration(dataFiles=[fn], darkFiles=self.darkFiles,
                                flatFiles=self.flatFiles, detector=self.detector,
                                calibrant=self.calibrant, wavelength=wavelength)
            rec.outfile = os.path.splitext(fn)[0] + ".proc.edf"
            rec.interactive = self.interactive
            rec.gui = self.gui
            rec.saturation = self.saturation
            rec.mask = self.mask
            rec.filter = self.filter
            rec.cutBackground = self.cutBackground
            rec.fixed = self.fixed
            rec.max_rings = self.max_rings
            rec.weighted = self.weighted
            if centerY:
                rec.ai.poni1 = centerY * self.detector.pixel1
            if centerX:
                rec.ai.poni2 = centerX * self.detector.pixel2
            if dist:
                rec.ai.dist = dist
            rec.preprocess()
            rec.extract_cpt(method=self.peakPicker)
            rec.refine()
            self.results[fn]["ai"] = rec.ai

    def regression(self):
        print(self.results)
        dist = numpy.zeros(len(self.results))
        x = dist.copy()
        poni1 = dist.copy()
        poni2 = dist.copy()
        rot1 = dist.copy()
        rot2 = dist.copy()
        rot3 = dist.copy()
        direct = dist.copy()
        tilt = dist.copy()
        trp = dist.copy()
        centerX = dist.copy()
        centerY = dist.copy()
        idx = 0
        print("")
        print("Results of linear regression for distance in mm")
        for key, dico in self.results.items():
            print(key, dico["dist"])
            print(dico["ai"])
            x[idx] = dico["dist"] * 1000
            dist[idx] = dico["ai"].dist
            poni1[idx] = dico["ai"].poni1
            poni2[idx] = dico["ai"].poni2
            rot1[idx] = dico["ai"].rot1
            rot2[idx] = dico["ai"].rot2
            rot3[idx] = dico["ai"].rot3
            f = dico["ai"].getFit2D()
            direct[idx] = f["directDist"]
            tilt[idx] = f["tilt"]
            trp[idx] = f["tiltPlanRotation"]
            centerX[idx] = f["centerX"]
            centerY[idx] = f["centerY"]
            idx += 1
        for name, elt in [("dist", dist),
                          ("poni1", poni1), ("poni2", poni2),
                          ("rot1", rot1), ("rot2", rot2), ("rot3", rot3),
                          ("direct", direct), ("tilt", tilt), ("trp", trp),
                          ("centerX", centerX), ("centerY", centerY)]:
            slope, intercept, r, _two, stderr = linregress(x, elt)

            print("%s = %s * dist_mm + %s \t R= %s\t stderr= %s" % (name, slope, intercept, r, stderr))


class CheckCalib(object):

    def __init__(self, poni=None, img=None, unit="2th_deg"):
        self.ponifile = poni
        if poni:
            self.ai = AzimuthalIntegrator.sload(poni)
        else:
            self.ai = None
        if img:
            self.img = fabio.open(img).data
        else:
            self.img = None
        self.mask = None
        self.r = None
        self.I = None
        self.wavelength = None
        self.resynth = None
        self.delta = None
        self.unit = unit
        self.masked_resynth = None
        self.masked_image = None
        self.offset = None
        self.data = None
        self.fig = None

    def __repr__(self, *args, **kwargs):
        if self.ai:
            return self.ai.__repr__()

    def parse(self):
        logger.debug("in parse")
        usage = "check_calib [options] -p param.poni image.edf"
        description = """Check_calib is a research tool aiming at validating both the geometric
calibration and everything else like flat-field correction, distortion
correction, at a sub-pixel level.

Note that `check_calib` program is obsolete as the same functionality is
available from within pyFAI-calib, using the `validate` command in the
refinement process.

        :returns: True if the parsing succeed, else False
        """
        version = "check_calib from pyFAI version %s: %s" % (PyFAI_VERSION, PyFAI_DATE)
        parser = ArgumentParser(usage=usage,
                                description=description)
        parser.add_argument("-V", "--version", action='version', version=version)
        parser.add_argument("args", metavar="FILE", help="Image file to check calibration for", nargs='+')
        parser.add_argument("-v", "--verbose",
                            action="store_true", dest="verbose", default=False,
                            help="switch to debug mode")
        parser.add_argument("-d", "--dark", dest="dark", metavar="FILE", type=str,
                            help="file containing the dark images to subtract", default=None)
        parser.add_argument("-f", "--flat", dest="flat", metavar="FILE", type=str,
                            help="file containing the flat images to divide", default=None)
        parser.add_argument("-m", "--mask", dest="mask", metavar="FILE", type=str,
                            help="file containing the mask", default=None)
        parser.add_argument("-p", "--poni", dest="poni", metavar="FILE", type=str,
                            help="file containing the diffraction parameter (poni-file)",
                            default=None)
        parser.add_argument("-e", "--energy", dest="energy", type=float,
                            help="energy of the X-Ray beam in keV (hc=%skeV.A)" % hc, default=None)
        parser.add_argument("-w", "--wavelength", dest="wavelength", type=float,
                            help="wavelength of the X-Ray beam in Angstrom", default=None)

        options = parser.parse_args()
        if options.verbose:
            logger.setLevel(logging.DEBUG)

        if options.mask is not None:
            self.mask = (fabio.open(options.mask).data != 0)
        args = expand_args(options.args)
        if len(args) > 0:
            f = args[0]
            if os.path.isfile(f):
                self.img = fabio.open(f).data.astype(numpy.float32)
            else:
                print("Please enter diffraction images as arguments")
                return False
            for f in args[1:]:
                self.img += fabio.open(f).data
        if options.dark and os.path.exists(options.dark):
            self.img -= fabio.open(options.dark).data
        if options.flat and os.path.exists(options.flat):
            self.img /= fabio.open(options.flat).data
        if options.poni:
            self.ai = AzimuthalIntegrator.sload(options.poni)
        self.data = [f for f in args if os.path.isfile(f)]
        if options.poni is None:
            logger.error("PONI parameter is mandatory")
            return False
        self.ai = AzimuthalIntegrator.sload(options.poni)
        if options.wavelength:
            self.ai.wavelength = 1e-10 * options.wavelength
        elif options.energy:
            self.ai.wavelength = 1e-10 * hc / options.energy
        # else:
        #     self.read_wavelength()
        return True

    def get_1dsize(self):
        logger.debug("in get_1dsize")
        return int(numpy.sqrt(self.img.shape[0] ** 2 + self.img.shape[1] ** 2))

    size1d = property(get_1dsize)

    def integrate(self):
        logger.debug("in integrate")
        self.r, self.I = self.ai.integrate1d_ng(self.img, self.size1d, mask=self.mask,
                                                unit=self.unit, method=("full", "histo", "cython"))

    def rebuild(self):
        """
        Rebuild the diffraction image and measures the offset with the reference
        :return: offset
        """
        logger.debug("in rebuild")
        if self.r is None:
            self.integrate()

        self.resynth = self.ai.calcfrom1d(self.r, self.I, shape=self.img.shape, mask=self.mask,
                                          dim1_unit=self.unit, correctSolidAngle=True)
        if self.mask is not None:
            self.img[numpy.where(self.mask)] = 0

        self.delta = self.resynth - self.img
        if self.mask is not None:
            smooth_mask = self.smooth_mask()
        else:
            smooth_mask = 1.0
        self.masked_resynth = self.resynth * smooth_mask
        self.masked_image = self.img * smooth_mask
        self.offset = measure_offset(self.masked_resynth, self.masked_image, withLog=0)
        print("Measured offset: %s" % str(self.offset))
        return self.offset

    def smooth_mask(self, hwhm=5):
        """
        smooth out around the mask to avoid aligning on the mask
        """
        logger.debug("in smooth_mask")
        fwhm = int(round(2.0 * hwhm))
        sigma = hwhm / math.sqrt(2 * math.log(2))

        if self.mask is not None:
            if not pyFAI_morphology:
                my, mx = numpy.ogrid[-fwhm: fwhm + 1, -fwhm:fwhm + 1]
                grow = (mx * mx + my * my) <= 4.0 * hwhm * hwhm
                big_mask = morphology.binary_dilation(self.mask, grow)
            else:
                big_mask = morphology.binary_dilation(self.mask.astype(numpy.int8), fwhm)
            smooth_mask = 1.0 - gaussian_filter(big_mask.astype(numpy.float32), sigma)
            return smooth_mask

    def show(self):
        """
        Show the image with the the errors
        """
        if self.fig is None:
            self.fig = pylab.figure()
            if not gui_utils.main_loop:
                self.fig.show()
        else:
            self.fig.clf()
        ax1 = self.fig.add_subplot(2, 2, 3)
        ax1.imshow(self.delta, aspect="auto", interpolation="nearest", origin="bottom")
        ax1.set_title("Difference image")
        ax2 = self.fig.add_subplot(2, 2, 1)
        ax2.imshow(self.masked_image, aspect="auto", interpolation="nearest", origin="bottom")
        ax2.set_title("Raw image")
        ax3 = self.fig.add_subplot(2, 2, 2)
        ax3.imshow(self.masked_resynth, aspect="auto", interpolation="nearest", origin="bottom")
        ax3.set_title("Rebuild image")
        ax4 = self.fig.add_subplot(2, 2, 4)
        ax4.plot(self.r, self.I)
        ax4.set_title("powder pattern")
        ax4.set_xlabel(r"2$\theta$ ($^o$)")
        ax4.set_ylabel("Intensity")
        update_fig(self.fig)


# Procedural version of calibration
def calib(img, calibrant, detector, basename="from_ipython", reconstruct=False, dist=0.1, gaussian=None, interactive=True):
    """
    Procedural interfact for calibration

    :param img: 2d array representing the calibration image
    :param calibrant: Instance of Calibrant, set-up with wavelength
    :param detector: Detector instance containing the mask
    :param basename: output file base
    :param recontruct: perform image reconstruction of masked pixel ?
    :param dist: initial distance
    :param gaussian: width of the gaussian used for difference of gaussian in the "massif" peak-picking algorithm
    :param interactive: set to False for testing
    :return: AzimuthalIntegrator instance
    """
    assert isinstance(detector, Detector)
    assert isinstance(calibrant, Calibrant)
    assert calibrant.wavelength

    if logging.root.level > logging.INFO:
        logging.warning("Lowering the log-level to INFO")
        logging.root.setLevel(logging.INFO)
    c = Calibration(wavelength=calibrant.wavelength,
                    detector=detector,
                    calibrant=calibrant,
                    gaussianWidth=gaussian)
    c.gui = interactive
    c.basename = basename
    c.pointfile = basename + ".npt"
    c.ai = AzimuthalIntegrator(dist=dist, detector=detector, wavelength=calibrant.wavelength)
    c.peakPicker = PeakPicker(img, reconst=reconstruct, mask=detector.mask,
                              pointfile=c.pointfile, calibrant=calibrant,
                              wavelength=calibrant.wavelength)
    if gaussian is not None:
        c.peakPicker.massif.setValleySize(gaussian)
    else:
        c.peakPicker.massif.initValleySize()

    if interactive:
        c.peakPicker.gui(log=True, maximize=True, pick=True)
        update_fig(c.peakPicker.fig)
    c.gui_peakPicker()
    c.ai.setPyFAI(**c.geoRef.getPyFAI())
    c.ai.wavelength = c.geoRef.wavelength

    return c.ai