Python pycgtypes.mat3函数代码示例

def mat3T(*args):
    if len(args) == 3:
        return mat3(args[0], args[1], args[2]).transpose()
    elif len(args) == 1:
        return mat3(args[0][0], args[0][1], args[0][2]).transpose()
    else:
        raise Exception, "Number of argument should be 1 or 3!"

    def UBxds_to_dnz(self):
        """ Convert the XDS direct space Orientation Matrix to a mosflm OM
        
        Denzo CAMERA coordinate frame has orthonormal axes with:

          y // to the rotation (spindel) axis
          z // to the beam
          x perpendicular to z and to the beam
        
        For more details see the denzo documentation:
        http://www.ccp4.ac.uk/dist/x-windows/Mosflm/doc/mosflm_user_guide.html#a3
        """
            
        if "UB" not in self.dict.keys():
            self.debut()
            
        BEAM = vec3(self.dict["beam"])
        ROT = vec3(self.dict["rot"])
        UBxds = self.dict["UB"]
                
        CAMERA_y = ROT.normalize()
        CAMERA_x = CAMERA_y.cross(BEAM).normalize()
        CAMERA_z = CAMERA_x.cross(CAMERA_y)
        CAMERA = mat3(CAMERA_x,CAMERA_y,CAMERA_z).transpose()
                    
        return  CAMERA * UBxds

def getPermutUB(PGoperators, UB, _epsilonCell=1e-2):
    "Apply Point group permutations to UB, return a list of permutated UB"

    cell = reciprocal(UB_to_cellParam(UB))
    permutedList = []

    for equiv in PGoperators:
        # Note that equiv_mat is the transpose of the normal equiv matrix!
        # because of the way mat3 matrices are contructed.
        equiv_mat = mat3(equiv[0],equiv[1],equiv[2])
        new_UB = UB * equiv_mat

        # One simple way to verify that the permutation is correct is
        # to extract the cell parameters from the permuted UB matrix
        if _debug:
            new_cell = reciprocal(UB_to_cellParam(new_UB))
            assert abs(new_cell[0] - cell[0]) < _epsilonCell and \
                   abs(new_cell[1] - cell[1]) < _epsilonCell and \
                   abs(new_cell[2] - cell[2]) < _epsilonCell and \
                   abs(new_cell[3] - cell[3]) < _epsilonCell and \
                   abs(new_cell[4] - cell[4]) < _epsilonCell and \
                   abs(new_cell[5] - cell[5]) < _epsilonCell 

        permutedList.append(new_UB)
    return permutedList

    def debut(self):
        "Do simple cristallographic calculations from XDS initial parameters"

        A = vec3(self.dict["A"])
        B = vec3(self.dict["B"])
        C = vec3(self.dict["C"])

        volum = A.cross(B)*C
        Ar = B.cross(C).__div__(volum)
        Br = C.cross(A).__div__(volum)
        Cr = A.cross(B).__div__(volum)

        """
        Ar = B.cross(C)/volum
        Br = C.cross(A)/volum
        Cr = A.cross(B)/volum
        """
        UBxds = mat3(Ar,Br,Cr)

        BEAM = vec3(self.dict["beam"])
        wavelength = 1/BEAM.length()

        self.dict["cell_volum"] = volum
        self.dict["wavelength"] = wavelength
        self.dict["Ar"] = Ar
        self.dict["Br"] = Br
        self.dict["Cr"] = Cr
        self.dict["UB"] = UBxds

    def __init__(self, init, rotationAxes=(ex, ey, ez), inversAxesOrder=0):
        self.inversAxesOrder = inversAxesOrder
        self.e1 = Vector(rotationAxes[0]).normalize()
        self.e2 = Vector(rotationAxes[1]).normalize()
        self.e3 = Vector(rotationAxes[2]).normalize()
        if inversAxesOrder:
            self.e1 = Vector(rotationAxes[2]).normalize()
            self.e3 = Vector(rotationAxes[0]).normalize()
        self.rotationAxes = self.e1, self.e2, self.e3
	if hasattr(init,'is_tensor') == 1:
	    self.tensor = init
	elif hasattr(init,'is_rotation') == 1:
	    self.tensor = init.tensor
	elif len(init) == 3:
            a1, a2, a3 = init
            if inversAxesOrder: a3, a2, a1 = init
            R1 = Rotation2(self.e1, a1)
            R2 = Rotation2(self.e2, a2)
            R3 = Rotation2(self.e3, a3)
	    self.tensor = R1*R2*R3
	elif len(init) == 9 and \
            (type(init) == type([]) or type(init) == type(())):
            "Used for compatibility with other tensor types, like cgtypes.mat3"
            "In that case, use mat3.mlist."
            self.tensor = mat3(list(init))
	else:
	    raise TypeError, 'no valid arguments'

 def UBxds_to_mos(self):
     """ Convert the XDS direct space Orientation Matrix to a mosflm OM
     
     Mosflm CAMERA coordinate frame has orthonormal axes with:
     
       z // rotation axis
       y perpendicular to z and to the beam
       x perpendicular to y and z (along the beam)
     
     For more details see the mosflm documentation:
     http://www.ccp4.ac.uk/dist/x-windows/Mosflm/doc/mosflm_user_guide.html#a3
     """
     
     if "UB" not in self.dict.keys():
         self.debut()
         
     BEAM = vec3(self.dict["beam"])
     ROT = vec3(self.dict["rot"])
     UBxds = self.dict["UB"]
     
     CAMERA_z = ROT.normalize()
     CAMERA_y = CAMERA_z.cross(BEAM).normalize()
     CAMERA_x = CAMERA_y.cross(CAMERA_z)
     CAMERA = mat3(CAMERA_x,CAMERA_y,CAMERA_z).transpose()
     
     return  CAMERA * UBxds * self.dict["wavelength"]

def test_2(axes_i):
    from ThreeAxisRotation import ThreeAxisRotation
    import Numeric, random
    
    r = random.uniform
    a1, a2, a3 = r(-pi,pi),r(-pi,pi),r(-pi,pi)
    angles_i = [a1, a2, a3]
    rot1 = ThreeAxisRotation (angles_i, rotationAxes=axes_i, inversAxesOrder=1)
    rot2 = ThreeAxisRotation2(angles_i, rotationAxes=axes_i, inversAxesOrder=1)
    rot1t = mat3(list(Numeric.ravel(rot1.tensor.array)))
    diffmat1 = math.sqrt(addReduceSq((rot1t - rot2.tensor).mlist))

    sol1A, sol1B = rot1.getAngles()
    sol2A, sol2B = rot2.getAngles()
    
    diff1I = math.sqrt(min(addReduceSq(map(diff, zip(sol1A,angles_i))), 
                           addReduceSq(map(diff, zip(sol1B,angles_i)))))
    diff2I = math.sqrt(min(addReduceSq(map(diff, zip(sol2A,angles_i))), 
                           addReduceSq(map(diff, zip(sol2B,angles_i)))))
    diffAA = math.sqrt(addReduceSq(map(diff, zip(sol2A,sol1A))))
    diffBB = math.sqrt(addReduceSq(map(diff, zip(sol2B,sol1B))))
    
    #if 1:
    ae = 1e-12
    if diffmat1 > 1e-12 or diffAA > ae or diffBB > ae or diff1I > ae or diff2I > ae:
        print "0 "+3*"%10.2f" % tuple(angles_i)
        print "1A"+3*"%10.2f" % tuple(sol1A)
        print "2A"+3*"%10.2f" % tuple(sol2A)
        print "1B"+3*"%10.2f" % tuple(sol1B)
        print "2B"+3*"%10.2f" % tuple(sol2B)
        print 'Matrix difference \t%.1e' % diffmat1
        print '1I Angle difference \t%.1e' % diff1I
        print '2I Angle difference \t%.1e' % diff1I
        print 'AA Angle difference \t%.1e' % diffAA
        print 'BB Angle difference \t%.1e' % diffBB

 def Adnz_to_rotxyz(self, Adnz, Bdnz, vertical, spindle):
     U2a = (Bdnz * spindle).normalize()
     U1a = Bdnz * vertical
     U1b = (U1a - (U1a * U2a) * U2a).normalize()
     Udnz0 = mat3(U1b, U2a, U1b.cross(U2a)).transpose()
     Adnz0 = Udnz0 * Bdnz
     RMAT = Adnz * Adnz0.inverse()
     dnzrmat = ThreeAxisRotation2(RMAT.mlist, self.DNZAxes)
     rotxyz = dnzrmat.getAngles()
     return rotxyz, Udnz0

def getPermutU(PGoperators, U, _epsilonCell=1e-4, debug=True):
    "Apply Point group permutations to U, return a list of permutated U"

    permutedList = []
    for equiv in PGoperators:
        # Note that equiv_mat is the transpose of the normal equiv matrix!
        # because of the way mat3 matrices are contructed.
        equiv_mat = mat3(equiv[0],equiv[1],equiv[2])
        new_U = equiv_mat.transpose() * U
        if is_orthogonal(new_U):
            permutedList.append(new_U)
        else:
            print "Internal Error: permuted U matrix is not orthogonal !"
            print new_U
    return permutedList

    def getOmega(self):
        """Calculate an Omega value (in radian) wich defines how the fast (X) and
        slow (Y) axis of detector files are  orientated toward the camera frame.
        The calculation of this omega value is supposed to reflect the Mosflm
        definition...

        But it seems that I get different values from the mosflm defaults... This
        may be due to: A) My missanderstanding of the mosflm documentation, B)
        Some tricks in the image reading routines.

        Nonetheless, this calculated value works for translating
        correctly the beam coordinates from XDS to mosflm [at least in the tested
        cases of MARCCD, MAR345 and ADSC detector images].

        Reference:
        http://www.ccp4.ac.uk/dist/x-windows/Mosflm/doc/mosflm_user_guide.html#a3
        """
        # Xd = CAMERA_y = beam
        # Yd = CAMERA_z = rot
        Xd =  vec3(self.dict["beam"]).normalize()
        Yd =  vec3(self.dict["rot"]).normalize()
        CAMERA_x = Xd.cross(Yd)
        CAMERA = mat3(CAMERA_x, Xd, Yd).transpose()
            
        # This is the definition of the fast:X and slow:Y axis for the detector files.
        XDSdetector_X = vec3(self.dict["detector_X"])
        XDSdetector_Y = vec3(self.dict["detector_Y"])
                
        # Now this axes are translated in the mosflm Camera frame
        Xs = XDSdetector_X*CAMERA
        Ys = XDSdetector_Y*CAMERA
        
        # Both angles should be identical.
        omegaX = Xd.angle(Xs)
        omegaY = Yd.angle(Ys)
        
        if _debug:
            print "DEBUG: X xds: fast =",XDSdetector_X
            print "DEBUG: Y xds: slow =",XDSdetector_Y
            print "DEBUG: Xs: fast =",XDSdetector_X,"->", Xs
            print "DEBUG: Ys: slow =",XDSdetector_Y,"->", Ys
            print "DEBUG: Xd: ", Xd
            print "DEBUG: Yd: ", Yd
            print "DEBUG: OmegaX:   %8.2f" % (omegaX*r2d)
            print "DEBUG: OmegaY:   %8.2f" % (omegaY*r2d)
            
        return omegaX

    def get_U0(self, rcell=None, vertical=None, spindle=None, clean=False):
        "Calculate denzo U0 from spindle, verctical"
        if not rcell: rcell = self.cell_r
        if not vertical: vertical = self.verticalAxis
        if not spindle: spindle = self.spindleAxis

        Bmat = self.get_B(rcell)
        vertical = vec3(vertical)
        spindle = vec3(spindle)

        U0y = (Bmat * spindle).normalize()
        U0xi = Bmat * vertical
        U0x = (U0xi - (U0xi * U0y) * U0y).normalize()

        U0 = mat3(U0x, U0y, U0x.cross(U0y)).transpose()

        # cleaning... Just cosmetic, not realy needed.
        if clean: U0 = cleanU0(U0)
        return U0

 def get_B(self, rcell=None):
     """ Denzo Orthogonalisation matrix
       b* is aligned with spindle axis (2-nd coordinate)
       a* is in the plane perpendicular to the beam (1-st,2-nd coords)
     """
     if not rcell:
         rcell = self.cell_r
     sr = map(sind, rcell[3:6])
     cr = map(cosd, rcell[3:6])
     B  = mat3()
     B13 = (cr[1]-cr[2]*cr[0])/sr[2]
     B[0,0] = rcell[0] * sr[2]
     B[1,0] = rcell[0] * cr[2]
     B[1,1] = rcell[1]
     B[0,2] = rcell[2] * B13
     B[1,2] = rcell[2] * cr[0]
     B[2,2] = rcell[2] * (sr[0]**2 - B13**2)**0.5
     
     return B

def axis_and_angle(mat_3):
    """From a rotation matrix return a corresponding rotation as an
       axis (a normalized vector) and angle (in radians).
       The angle is in the interval (-pi, pi]
    """
    asym = -asymmetrical_part(mat_3)
    axis = vec3(asym[1, 2], asym[2, 0], asym[0, 1])
    sine = axis.length()
    if abs(sine) > 1.e-10:
        axis = axis/sine
        projector = dyadic_product(axis, axis)
        cosine = trace((mat_3-projector))/(3.-axis*axis)
        angle = angle_from_sine_and_cosine(sine, cosine)
    else:
        tsr = 0.5*(mat_3+mat3(1))
        diag = tsr[0, 0], tsr[1, 1], tsr[3, 3] 
        i = tsr.index(max(diag))
        axis = vec3(tsr.getRow(i)/(tsr[i, i])**0.5)
        angle = 0.
        if trace(tsr) < 2.:
            angle = math.pi
    return axis, angle

        diag = tsr[0, 0], tsr[1, 1], tsr[3, 3] 
        i = tsr.index(max(diag))
        axis = vec3(tsr.getRow(i)/(tsr[i, i])**0.5)
        angle = 0.
        if trace(tsr) < 2.:
            angle = math.pi
    return axis, angle

# Test code
if __name__ == '__main__':

    from Scientific.Geometry import Vector ##.Transformation import *
    from Scientific.Geometry.Transformation import Rotation
    from random import random
    #
    Q = mat3(0.36, 0.48, -0.8, -0.8, 0.6, 0, 0.48, 0.64, 0.60)
    axis_q, angle_q = axis_and_angle(Q)
    print "Axis_q:  %9.6f%9.6f%9.6f" % tuple(axis_q),
    print "Angle_q: %10.5f" % (angle_q*R2D)
    #
    for iii in range(1e6):
        axis_i = list(vec3([random(), random(), random()]).normalize())
        angle_i = 3*random()
        rme = mat3().rotation(angle_i, vec3(axis_i))
        axis_1, angle_1 = axis_and_angle(rme)

        v = Vector(axis_i)
        r = Rotation(v, angle_i)
        axis_2, angle_2 = r.axisAndAngle()
        axis_d = (axis_1 - vec3(tuple(axis_2))).length()
        angle_d = abs(angle_1 - angle_2)

def dyadic_product(vector1, vector2):
    "Dyadic product of two vectors."
    matr1, matr2 = mat3(), mat3()
    matr1.setColumn(0, vector1)
    matr2.setRow(0, vector2)
    return matr1*matr2