Python azimuthalIntegrator.AzimuthalIntegrator类代码示例

    def scan_tilt(self, width=1.0, points=10):
        """
        ???

        @param width: ???
        @type width: float ???
        @param points: ???
        @type points: int ???

        @return: ???
        @rtype: ???
        """
        f = self.ai.getFit2D()
        out = []
        for x in numpy.linspace(f["tilt"] - width / 2.0, f["tilt"] + width / 2.0, points):
            ax = AzimuthalIntegrator()
            fx = f.copy()
            fx["tilt"] = x
            ax.setFit2D(**fx)
            #            print ax
            ref = Refinment2D(self.img, ax)
            res = ref.diff_tth_tilt()
            print "x= %.3f mean= %e" % (x, res)
            out.append(res)
        return numpy.linspace(f["tilt"] - width / 2.0, f["tilt"] + width / 2.0, points), out

class TestSetter(unittest.TestCase):
    def setUp(self):
        self.ai = AzimuthalIntegrator()
        shape = (10, 15)
        self.rnd1 = numpy.random.random(shape).astype(numpy.float32)
        self.rnd2 = numpy.random.random(shape).astype(numpy.float32)
        if not os.path.isdir(tmp_dir):
            os.mkdir(tmp_dir)

        fd, self.edf1 = tempfile.mkstemp(".edf", "testAI1", tmp_dir)
        os.close(fd)
        fd, self.edf2 = tempfile.mkstemp(".edf", "testAI2", tmp_dir)
        os.close(fd)
        fabio.edfimage.edfimage(data=self.rnd1).write(self.edf1)
        fabio.edfimage.edfimage(data=self.rnd2).write(self.edf2)

    def tearDown(self):
        recursive_delete(tmp_dir)

    def test_flat(self):
        self.ai.set_flatfiles((self.edf1, self.edf2), method="mean")
        self.assert_(self.ai.flatfiles == "%s(%s,%s)" % ("mean", self.edf1, self.edf2), "flatfiles string is OK")
        self.assert_(abs(self.ai.flatfield - 0.5 * (self.rnd1 + self.rnd2)).max() == 0, "Flat array is OK")

    def test_dark(self):
        self.ai.set_darkfiles((self.edf1, self.edf2), method="mean")
        self.assert_(self.ai.darkfiles == "%s(%s,%s)" % ("mean", self.edf1, self.edf2), "darkfiles string is OK")
        self.assert_(abs(self.ai.darkcurrent - 0.5 * (self.rnd1 + self.rnd2)).max() == 0, "Dark array is OK")

    def diff_Fit2D(self, axis="all", dx=0.1):
        """
        ???

        @param axis: ???
        @type axis: ???
        @param dx: ???
        @type dx: ???

        @return: ???
        @rtype: ???
        """
        tth, I = self.ai.xrpd(self.img, max(self.shape))
        dI = SGModule.getSavitzkyGolay(I, npoints=5, degree=2, order=1) / (tth[1] - tth[0])
        dImg = self.reconstruct(tth, dI)
        f = self.ai.getFit2D()
        tth2d_ref = self.ai.twoThetaArray(self.shape)  # useless variable ???

        keys = ["centerX", "centerY", "tilt", "tiltPlanRotation"]
        if axis != "all":
            keys = [i for i in keys if i == axis]
        grad = {}
        for key in keys:
            fp = f.copy()
            fp[key] += dx
            ap = AzimuthalIntegrator()
            ap.setFit2D(**fp)
            dtth = (ap.twoThetaArray(self.shape) - self.ai.twoThetaArray(self.shape)) / dx
            grad[key] = (dtth * dImg).sum()
        if axis == "all":
            return grad
        else:
            return grad[axis]

    def scan_Fit2D(self, width=1.0, points=10, axis="tilt", dx=0.1):
        """
        ???

        @param width: ???
        @type width: float ???
        @param points: ???
        @type points: int ???
        @param axis: ???
        @type axis: str ???
        @param dx: ???
        @type dx: float ???

        @return: ???
        @rtype: ???
        """
        logger.info("Scanning along axis %s" % axis)
        f = self.ai.getFit2D()
        out = []
        meas_pts = numpy.linspace(f[axis] - width / 2.0, f[axis] + width / 2.0, points)
        for x in meas_pts:
            ax = AzimuthalIntegrator()
            fx = f.copy()
            fx[axis] = x
            ax.setFit2D(**fx)
            ref = Refinment2D(self.img, ax)
            res = ref.diff_Fit2D(axis=axis, dx=dx)
            print "x= %.3f mean= %e" % (x, res)
            out.append(res)
        return meas_pts, out

class test_saxs(unittest.TestCase):
    saxsPilatus = "1492/bsa_013_01.edf"
    maskFile = "1491/Pcon_01Apr_msk.edf"
    maskRef = "1490/bioSaxsMaskOnly.edf"
    maskDummy = "1488/bioSaxsMaskDummy.edf"
    poniFile = "1489/bioSaxs.poni"
    ai = None
    tmp_dir = os.environ.get("PYFAI_TEMPDIR", os.path.join(os.path.dirname(os.path.abspath(__file__)), "tmp"))

    def setUp(self):
        self.edfPilatus = UtilsTest.getimage(self.__class__.saxsPilatus)
        self.maskFile = UtilsTest.getimage(self.__class__.maskFile)
        self.poniFile = UtilsTest.getimage(self.__class__.poniFile)
        self.maskRef = UtilsTest.getimage(self.__class__.maskRef)
        self.maskDummy = UtilsTest.getimage(self.__class__.maskDummy)
        self.ai = AzimuthalIntegrator()
        self.ai.load(self.poniFile)
        if not os.path.isdir(self.tmp_dir):
            os.mkdir(self.tmp_dir)

    def test_mask(self):
        """test the generation of mask"""
        print self.edfPilatus
        data = fabio.open(self.edfPilatus).data
        mask = fabio.open(self.maskFile).data
        assert abs(self.ai.makeMask(data, mask=mask).astype(int) - fabio.open(self.maskRef).data).max() == 0
        assert abs(self.ai.makeMask(data, mask=mask, dummy= -2, delta_dummy=1.1).astype(int) - fabio.open(self.maskDummy).data).max() == 0

    def __init__(self, data, dist=1, poni1=None, poni2=None,
                 rot1=0, rot2=0, rot3=0,
                 pixel1=None, pixel2=None, splineFile=None, detector=None,
                 wavelength=None, dSpacing=None):
        """
        @param data: ndarray float64 shape = n, 3
            col0: pos in dim0 (in pixels)
            col1: pos in dim1 (in pixels)
            col2: ring index in dSpacing file
        @param dist: guessed sample-detector distance (optional, in m)
        @param poni1: guessed PONI coordinate along the Y axis (optional, in m)
        @param poni2: guessed PONI coordinate along the X axis (optional, in m)
        @param rot1: guessed tilt of the detector around the Y axis (optional, in rad)
        @param rot2: guessed tilt of the detector around the X axis (optional, in rad)
        @param rot3: guessed tilt of the detector around the incoming beam axis (optional, in rad)
        @param pixel1: Pixel size along the vertical direction of the detector (in m), almost mandatory
        @param pixel2: Pixel size along the horizontal direction of the detector (in m), almost mandatory
        @param splineFile: file describing the detector as 2 cubic splines. Replaces pixel1 & pixel2
        @param detector: name of the detector or Detector instance. Replaces splineFile, pixel1 & pixel2
        @param wavelength: wavelength in m (1.54e-10)
        @param dSpacing: filename or list or array or vector containing the d-spacing (in Angstrom)

        """
        self.data = numpy.array(data, dtype="float64")
        assert self.data.ndim == 2
        assert self.data.shape[1] == 3
        assert self.data.shape[0]>0

        if (pixel1 is None) and (pixel2 is None) and (splineFile is None) and (detector is None):
            raise RuntimeError("Setting up the geometry refinement without knowing the detector makes little sense")
        AzimuthalIntegrator.__init__(self, dist, 0, 0,
                                     rot1, rot2, rot3,
                                     pixel1, pixel2, splineFile, detector, wavelength=wavelength)

        if (poni1 is None) or (poni2 is None):
            self.guess_poni()
        else:
            self.poni1 = float(poni1)
            self.poni2 = float(poni2)
        self._dist_min = 0
        self._dist_max = 10
        self._poni1_min = -10000 * self.pixel1
        self._poni1_max = 15000 * self.pixel1
        self._poni2_min = -10000 * self.pixel2
        self._poni2_max = 15000 * self.pixel2
        self._rot1_min = -pi
        self._rot1_max = pi
        self._rot2_min = -pi
        self._rot2_max = pi
        self._rot3_min = -pi
        self._rot3_max = pi
        self._wavelength_min = 1e-15
        self._wavelength_max = 100.e-10
        if dSpacing is not None:
            if type(dSpacing) in types.StringTypes:
                self.dSpacing = numpy.loadtxt(dSpacing)
            else:
                self.dSpacing = numpy.array(dSpacing, dtype=numpy.float64)
        else:
            self.dSpacing = numpy.array([])

class test_setter(unittest.TestCase):
    tmp_dir = os.environ.get("PYFAI_TEMPDIR", os.path.join(os.path.dirname(os.path.abspath(__file__)), "tmp"))

    def setUp(self):
        self.ai = AzimuthalIntegrator()
        shape = (10, 15)
        self.rnd1 = numpy.random.random(shape).astype(numpy.float32)
        self.rnd2 = numpy.random.random(shape).astype(numpy.float32)
        if not os.path.isdir(self.tmp_dir):
            os.mkdir(self.tmp_dir)

        fd, self.edf1 = tempfile.mkstemp(".edf", "testAI1", self.tmp_dir)
        os.close(fd)
        fd, self.edf2 = tempfile.mkstemp(".edf", "testAI2", self.tmp_dir)
        os.close(fd)
        fabio.edfimage.edfimage(data=self.rnd1).write(self.edf1)
        fabio.edfimage.edfimage(data=self.rnd2).write(self.edf2)
    def tearDown(self):
        if os.path.exists(self.edf1):
            os.unlink(self.edf1)
        if os.path.exists(self.edf2):
            os.unlink(self.edf2)
    def test_flat(self):
        self.ai.set_flatfiles((self.edf1,self.edf2), method="mean")
        self.assert_(self.ai.flatfiles == "%s(%s,%s)" % ("mean", self.edf1, self.edf2), "flatfiles string is OK")
        self.assert_(abs(self.ai.flatfield-0.5*(self.rnd1+self.rnd2)).max() == 0, "Flat array is OK")
    def test_dark(self):
        self.ai.set_darkfiles((self.edf1, self.edf2), method="mean")
        self.assert_(self.ai.darkfiles == "%s(%s,%s)" % ("mean", self.edf1, self.edf2), "darkfiles string is OK")
        self.assert_(abs(self.ai.darkcurrent-0.5*(self.rnd1+self.rnd2)).max() == 0, "Dark array is OK")

class TestAzimPilatus(unittest.TestCase):
    img = UtilsTest.getimage("1884/Pilatus6M.cbf")

    def setUp(self):
        """Download files"""
        self.data = fabio.open(self.img).data
        self.ai = AzimuthalIntegrator(detector="pilatus6m")
        self.ai.setFit2D(300, 1326, 1303)

    def test_separate(self):
        bragg, amorphous = self.ai.separate(self.data)
        self.assert_(amorphous.max() < bragg.max(), "bragg is more intense than amorphous")
        self.assert_(amorphous.std() < bragg.std(), "bragg is more variatic than amorphous")

    def diff_tth_X(self, dx=0.1):
        """
        Jerome peux-tu décrire de quoi il retourne ???

        @param dx: ???
        @type: float ???

        @return: ???
        @rtype: ???
        """
        f = self.ai.getFit2D()
        fp = f.copy()
        fm = f.copy()
        fm["centerX"] -= dx / 2.0
        fp["centerX"] += dx / 2.0
        ap = AzimuthalIntegrator()
        am = AzimuthalIntegrator()
        ap.setFit2D(**fp)
        am.setFit2D(**fm)
        dtthX = (ap.twoThetaArray(self.shape) - am.twoThetaArray(self.shape))\
            / dx
        tth, I = self.ai.xrpd(self.img, max(self.shape))
        dI = SGModule.getSavitzkyGolay(I, npoints=5, degree=2, order=1)\
            / (tth[1] - tth[0])
        dImg = self.reconstruct(tth, dI)
        return (dtthX * dImg).sum()

 def setUp(self):
     """Download files"""
     self.fit2dFile = UtilsTest.getimage(self.__class__.fit2dFile)
     self.halfFrelon = UtilsTest.getimage(self.__class__.halfFrelon)
     self.splineFile = UtilsTest.getimage(self.__class__.splineFile)
     poniFile = UtilsTest.getimage(self.__class__.poniFile)
     with open(poniFile) as f:
         data = []
         for line in f:
             if line.startswith("SplineFile:"):
                 data.append("SplineFile: " + self.splineFile)
             else:
                 data.append(line.strip())
     self.poniFile = os.path.join(self.tmp_dir, os.path.basename(poniFile))
     with open(self.poniFile, "w") as f:
         f.write(os.linesep.join(data))
     self.fit2d = numpy.loadtxt(self.fit2dFile)
     self.ai = AzimuthalIntegrator()
     self.ai.load(self.poniFile)
     self.data = fabio.open(self.halfFrelon).data
     if not os.path.isdir(self.tmp_dir):
         os.makedirs(self.tmp_dir)
     for tmpfile in self.tmpfiles.values():
         if os.path.isfile(tmpfile):
             os.unlink(tmpfile)

    def __init__(self, img_model=None):
        """
        :param img_model:
        :type img_model: ImgModel
        """
        self.img_model = img_model
        self.points = []
        self.points_index = []
        self.spectrum_geometry = AzimuthalIntegrator()
        self.cake_geometry = None
        self.calibrant = Calibrant()
        self.start_values = {'dist': 200e-3,
                             'wavelength': 0.3344e-10,
                             'pixel_width': 79e-6,
                             'pixel_height': 79e-6,
                             'polarization_factor': 0.99}
        self.orig_pixel1 = 79e-6
        self.orig_pixel2 = 79e-6
        self.fit_wavelength = False
        self.fit_distance = True
        self.is_calibrated = False
        self.use_mask = False
        self.filename = ''
        self.calibration_name = 'None'
        self.polarization_factor = 0.99
        self.supersampling_factor = 1
        self._calibrants_working_dir = os.path.dirname(calibrants.__file__)

        self.cake_img = np.zeros((2048, 2048))
        self.tth = np.linspace(0, 25)
        self.int = np.sin(self.tth)
        self.num_points = len(self.int)

        self.peak_search_algorithm = None

def image_test_rings():
    rings = 10
    mod = 50
    detector = detector_factory("Titan")
    sigma = detector.pixel1 * 4
    shape = detector.max_shape
    ai = AzimuthalIntegrator(detector=detector)
    ai.setFit2D(1000, 1000, 1000)
    r = ai.rArray(shape)
    r_max = r.max()
    chi = ai.chiArray(shape)
    img = numpy.zeros(shape)
    modulation = (1 + numpy.sin(5 * r + chi * mod))
    for radius in numpy.linspace(0, r_max, rings):
        img += numpy.exp(-(r - radius) ** 2 / (2 * (sigma * sigma)))
    return img * modulation

    def setUp(self):
        # defining geometry
        image_shape = [2048, 2048]  # pixel
        detector_distance = 200  # mm
        wavelength = 0.31  # angstrom
        center_x = 1024  # pixel
        center_y = 1024  # pixel
        self.tilt = 0  # degree
        self.rotation = 0  # degree
        pixel_size = 79  # um
        dummy_tth = np.linspace(0, 35, 2000)
        dummy_int = np.ones(dummy_tth.shape)
        self.geometry = AzimuthalIntegrator()
        self.geometry.setFit2D(directDist=detector_distance,
                               centerX=center_x,
                               centerY=center_y,
                               tilt=self.tilt,
                               tiltPlanRotation=self.rotation,
                               pixelX=pixel_size,
                               pixelY=pixel_size)
        self.geometry.wavelength = wavelength / 1e10
        self.dummy_img = self.geometry.calcfrom1d(dummy_tth, dummy_int, shape=image_shape, correctSolidAngle=True)

        self.tth_array = self.geometry.twoThetaArray(image_shape)
        self.azi_array = self.geometry.chiArray(image_shape)

class ObliqueAngleDetectorAbsorptionCorrectionTest(unittest.TestCase):
    def setUp(self):
        # defining geometry
        image_shape = [2048, 2048]  # pixel
        detector_distance = 200  # mm
        wavelength = 0.31  # angstrom
        center_x = 1024  # pixel
        center_y = 1024  # pixel
        self.tilt = 0  # degree
        self.rotation = 0  # degree
        pixel_size = 79  # um
        dummy_tth = np.linspace(0, 35, 2000)
        dummy_int = np.ones(dummy_tth.shape)
        self.geometry = AzimuthalIntegrator()
        self.geometry.setFit2D(directDist=detector_distance,
                               centerX=center_x,
                               centerY=center_y,
                               tilt=self.tilt,
                               tiltPlanRotation=self.rotation,
                               pixelX=pixel_size,
                               pixelY=pixel_size)
        self.geometry.wavelength = wavelength / 1e10
        self.dummy_img = self.geometry.calcfrom1d(dummy_tth, dummy_int, shape=image_shape, correctSolidAngle=True)

        self.tth_array = self.geometry.twoThetaArray(image_shape)
        self.azi_array = self.geometry.chiArray(image_shape)

    def tearDown(self):
        del self.azi_array
        del self.tth_array
        del self.dummy_img
        del self.geometry
        gc.collect()

    def test_that_it_is_correctly_calculating(self):
        oblique_correction = ObliqueAngleDetectorAbsorptionCorrection(
                tth_array=self.tth_array,
                azi_array=self.azi_array,
                detector_thickness=40,
                absorption_length=465.5,
                tilt=self.tilt,
                rotation=self.rotation
        )
        oblique_correction_data = oblique_correction.get_data()
        self.assertGreater(np.sum(oblique_correction_data), 0)
        self.assertEqual(oblique_correction_data.shape, self.dummy_img.shape)
        del oblique_correction

    def run(self):
        ai = AzimuthalIntegrator(
            dist=self.__distance,
            poni1=self.__poni1,
            poni2=self.__poni2,
            rot1=self.__rotation1,
            rot2=self.__rotation2,
            rot3=self.__rotation3,
            detector=self.__detector,
            wavelength=self.__wavelength)

        numberPoint1D = 1024
        numberPointRadial = 400
        numberPointAzimuthal = 360

        # FIXME error model, method

        self.__result1d = ai.integrate1d(
            data=self.__image,
            npt=numberPoint1D,
            unit=self.__radialUnit,
            mask=self.__mask,
            polarization_factor=self.__polarizationFactor)

        self.__result2d = ai.integrate2d(
            data=self.__image,
            npt_rad=numberPointRadial,
            npt_azim=numberPointAzimuthal,
            unit=self.__radialUnit,
            mask=self.__mask,
            polarization_factor=self.__polarizationFactor)

        if self.__calibrant:

            rings = self.__calibrant.get_2th()
            rings = filter(lambda x: x <= self.__result1d.radial[-1], rings)
            rings = list(rings)
            try:
                rings = utils.from2ThRad(rings, self.__radialUnit, self.__wavelength, ai)
            except ValueError:
                message = "Convertion to unit %s not supported. Ring marks ignored"
                _logger.warning(message, self.__radialUnit)
                rings = []
        else:
            rings = []
        self.__ringAngles = rings