"""
Operations on rotation matrices.
Copyright Schrodinger, LLC. All rights reserved.
"""
# This is a partial port of the cctbx/sgtbx/rot_mx.cpp (LEGAL-980)
# Please don't use this outside of the cctbx modules
import math
from fractions import Fraction
import numpy
import sympy
[docs]def basis_of_mirror_plane_with_normal(u):
""" Primitive setting assumed """
assert u != (0, 0, 0)
basis = ()
for t in ((1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 0), (1, 0, 1), (0, 1, 1),
(-1, 1, 0), (-1, 0, 1), (0, -1, 1), (1, 1, 1), (-1, 1, 1),
(1, -1, 1), (1, 1, -1)):
if len(basis) == 2:
break
if u[0] * t[0] + u[1] * t[1] + u[2] * t[2] == 0:
basis += (t,)
return basis
[docs]def back_substitution_int(re_mx, sol=None):
# From scitbx/matrix/row_echelon.h
nr, nc = re_mx.shape
flag_indep = [True] * nc
if sol is not None:
sol = numpy.array(sol, dtype=int)
d = 1
for ir in reversed(range(nr)):
for ic in range(nc):
if re_mx[ir, ic] != 0:
break
else:
continue
flag_indep[ic] = False
if sol is None:
continue
icp = ic + 1
nv = nc - icp
if nv != 0:
sol[ic:ic + nv] += numpy.matmul(re_mx[ir:ir + nv, icp:icp + nv],
sol[icp:icp + nv])
sol[ic] *= -1
else:
sol[ic] = 0
mrc = re_mx[ir, ic]
f = math.gcd(sol[ic], mrc)
if mrc < 0:
f *= -1
sol[ic] //= f
f = mrc // f
if f != 1:
for jc in range(nc):
if jc != ic:
sol[jc] *= f
d *= f
return d, flag_indep, sol
[docs]def homog_rank_2(re_mx):
# From cctbx/sgtbx/row_echelon_solve.cpp
assert re_mx.shape == (2, 3)
d, indep, tmp = back_substitution_int(re_mx)
assert d >= 1
assert len(list(x for x in indep if x)) == 1
ev = [0, 0, 0]
ev[indep.index(True)] = 1
d, flag_indep, ev = back_substitution_int(re_mx, ev)
assert d >= 1
return ev
[docs]class parse_string(object):
[docs] def __init__(self, operator):
self.op_ = operator
[docs] def string(self):
return self.op_
[docs] def parse_symm_op(self):
line = self.op_
rot_mx = numpy.zeros(shape=(3, 3), dtype=int)
tr_vec = numpy.zeros(3)
for row, oper in enumerate(line.split(',')):
orig_oper = oper
for col, lett in enumerate(['x', 'y', 'z']):
if '+' + lett in oper:
rot_mx[row, col] = 1
oper = oper.replace('+' + lett, '')
elif '-' + lett in oper:
rot_mx[row, col] = -1
oper = oper.replace('-' + lett, '')
elif lett in oper:
rot_mx[row, col] = 1
oper = oper.replace(lett, '')
if oper:
try:
tr_vec[row] = float(Fraction(oper))
except:
print(oper, orig_oper)
return rot_mx, tr_vec
[docs]class tr_vec(object):
# See cctbx/sgtbx/basic.h
sg_t_den = 12
[docs] def __init__(self, *args):
if len(args) == 0:
self.num_ = numpy.array([0, 0, 0], dtype=int)
self.den_ = self.sg_t_den
elif len(args) == 1:
if isinstance(args[0], int):
self.num_ = numpy.array([0, 0, 0], dtype=int)
self.den_ = args[0]
else:
self.num_ = numpy.array(args[0], dtype=int)
self.den_ = self.sg_t_den
elif len(args) == 2:
self.num_ = numpy.array(args[0], dtype=int)
self.den_ = args[1]
assert self.num_.shape == (3,)
def __eq__(self, other):
assert isinstance(other, self.__class__)
if (self.den_ != other.den_):
return False
return numpy.array_equal(self.num_, other.num_)
[docs] def den(self):
"""
Get denominator of the rotation matrix.
"""
return self.den_
[docs] def num(self):
"""
Get rotation matrix.
"""
return tuple(self.num_.flatten().tolist())
[docs] def is_valid(self):
"""
Check if the matrix is valid.
True only if den() != 0.
"""
return self.den() != 0
[docs] def is_zero(self):
"""
True only if all elements of the numertor are equal to zero.
"""
return not self.num_.all()
[docs] def cancel(self):
ret = self.copy()
tmp = numpy.append(ret.num_, ret.den_)
gcd = numpy.gcd.reduce(tmp)
if gcd == 1:
return ret
ret.num_ = ret.num_ / gcd
ret.den_ = ret.den_ // gcd
return ret
[docs] def copy(self):
ret = self.__class__()
ret.num_ = numpy.array(self.num_, dtype=int)
ret.den_ = self.den_
return ret
[docs] def new_denominator(self, new_den):
ret = self.__class__(new_den)
ret.num_ = self.num_ * new_den
if (ret.num_ % self.den()).any():
raise ValueError('Wrong new determinant %d' % new_den)
ret.num_ = ret.num_ / self.den()
return ret
[docs] def scale(self, fac):
assert fac >= 1
assert isinstance(fac, int)
ret = self.copy()
ret.num_ *= fac
ret.den_ *= fac
return ret
[docs] def divide(self, rhs):
assert not (self.num_ % rhs).any()
return tr_vec(self.num_ / rhs, self.den_)
[docs] def minus(self, rhs):
assert self.den_ == rhs.den_, '%s != %s' % (self.den_, rhs.den_)
return tr_vec(self.num_ - rhs.num_, self.den_)
[docs] def as_double(self):
return tuple((self.num_ / self.den_).tolist())
[docs]class rot_mx(object):
"""
3x3 rotation matrix.
The elements of the matrix are stored as integers and a common
denominator. The actual value of an element is obtained by
dividing the integer number by the denominator.
"""
[docs] def __init__(self, *args):
"""
Initialize matrix.
See cctbx/sgtbx/rot_mx.h for more documentation
"""
if len(args) == 0:
self.num_ = numpy.eye(3, dtype=int)
self.den_ = 1
elif len(args) == 1:
val = args[0]
if isinstance(val, int):
self.num_ = numpy.eye(3, dtype=int) * val
self.den_ = val
else:
assert len(val) == 9
self.num_ = numpy.array(val, dtype=int).reshape((3, 3))
self.den_ = 1
elif len(args) == 2:
val = args[0]
if isinstance(val, int):
self.den_ = val
self.num_ = numpy.eye(3, dtype=int) * self.den_ * args[1]
else:
assert len(val) == 9
self.num_ = numpy.array(val, dtype=int).reshape((3, 3))
self.den_ = args[1]
def __eq__(self, other):
"""
Check if two objects are equal.
True only if both the numerators and the denominators are equal.
"""
assert isinstance(other, self.__class__)
if (self.den_ != other.den()):
return False
return self.num() == other.num()
[docs] def den(self):
"""
Get denominator of the rotation matrix.
"""
return self.den_
[docs] def num(self):
"""
Get rotation matrix.
"""
return tuple(self.num_.flatten().tolist())
[docs] def is_valid(self):
"""
Check if the matrix is valid.
True only if den() != 0.
"""
return self.den() != 0
[docs] def is_unit_mx(self):
"""
"""
return numpy.array_equal(self.num_,
numpy.eye(3, dtype=int) * self.den())
[docs] def minus_unit_mx(self):
val = self.num_ - (numpy.eye(3, dtype=int) * self.den())
return self.__class__(tuple(val.flatten().tolist()))
[docs] def new_denominator(self, new_den):
ret = self.__class__(new_den)
ret.num_ = self.num_ * new_den
if (ret.num_ % self.den()).any():
raise ValueError('Wrong new determinant %d' % new_den)
ret.num_ = ret.num_ / self.den()
return ret
[docs] def matrix_cofactor(self):
m = self.num_.flatten().tolist()
ret = [
m[4] * m[8] - m[5] * m[7], -m[1] * m[8] + m[2] * m[7],
m[1] * m[5] - m[2] * m[4], -m[3] * m[8] + m[5] * m[6],
m[0] * m[8] - m[2] * m[6], -m[0] * m[5] + m[2] * m[3],
m[3] * m[7] - m[4] * m[6], -m[0] * m[7] + m[1] * m[6],
m[0] * m[4] - m[1] * m[3]
]
return numpy.array(ret, dtype=int).reshape((3, 3))
[docs] def copy(self, transpose=False, minus=False):
ret = self.__class__()
ret.num_ = numpy.array(self.num_, dtype=int)
if minus:
ret.num_ = -ret.num_
if transpose:
ret.num_ = ret.num_.T
ret.den_ = self.den_
return ret
[docs] def scale(self, fac):
assert fac >= 1
assert isinstance(fac, int)
ret = self.copy()
ret.num_ *= fac
ret.den_ *= fac
return ret
[docs] def determinant(self):
assert self.den_ > 0
return int(numpy.linalg.det(self.num_) / (self.den_**3))
[docs] def transpose(self):
return self.copy(transpose=True)
[docs] def inverse(self):
det = int(numpy.linalg.det(self.num_))
if det == 0:
raise ValueError('Cannot invert matrix')
ret = self.copy()
ret.num_ = self.matrix_cofactor() * (ret.den()**2)
return ret.opdiv(det)
[docs] def cancel(self):
ret = self.copy()
tmp = numpy.append(ret.num_.flatten(), ret.den())
gcd = numpy.gcd.reduce(tmp)
if gcd == 1:
return ret
ret.num_ = ret.num_ / gcd
ret.den_ = ret.den_ // gcd
return ret
[docs] def inverse_cancel(self):
det = int(numpy.linalg.det(self.num_))
if det == 0:
raise ValueError('Cannot invert matrix')
ret = self.copy()
gcd = numpy.gcd(det, self.den_)
ret.num_ = ret.matrix_cofactor() * self.den_ // gcd
ret.den_ = 1
return ret.divide(det // gcd)
[docs] def opdiv(self, rhs):
assert isinstance(rhs, int)
if (self.num_ % rhs).any():
raise ValueError('Cannot divide by %d' % rhs)
ret = self.copy()
ret.num_ = ret.num_ // rhs
return ret
[docs] def divide(self, rhs):
ret = self.copy()
if (rhs < 0):
ret.num *= -1
rhs = -rhs
ret.den_ *= rhs
return ret.cancel()
[docs] def multiply(self, rhs):
assert isinstance(rhs, (self.__class__, tr_vec))
if isinstance(rhs, tr_vec):
ret = tr_vec()
ret.num_ = numpy.dot(self.num_, rhs.num_)
ret.den_ = self.den_ * rhs.den_
return ret.cancel()
ret = self.copy()
ret.num_ = numpy.dot(ret.num_, rhs.num_)
ret.den_ *= rhs.den_
return ret
[docs] def oper_multiply(self, rhs):
assert isinstance(rhs, tr_vec)
return tr_vec(numpy.dot(self.num_, rhs.num_), self.den() * rhs.den())
[docs] def as_double(self):
return tuple((numpy.array(self.num_, dtype=float) /
self.den_).flatten().tolist())
[docs] def type(self):
det = numpy.linalg.det(self.num_)
if det == -1 or det == 1:
trace = self.num_.trace()
if trace == -3:
return -1
elif trace == -2:
return -6
elif trace == -1:
if det == -1:
return -4
else:
return 2
elif trace == 0:
if det == -1:
return -3
else:
return 3
elif trace == 1:
if det == -1:
return -2
else:
return 4
elif trace == 2:
return 6
elif trace == 3:
return 1
return 0
[docs] def order(self, type_in=0):
if type_in == 0:
type_in = self.type()
if type_in > 0:
return type_in
if type_in % 2:
return -type_in * 2
return -type_in
[docs] def accumulate(self, type_in=0):
assert self.den_ == 1
order = self.order(type_in)
if order == 1:
return self.copy()
assert order != 0
result = numpy.eye(3, dtype=int)
num_ = numpy.array(self.num_, dtype=int)
result += num_
for idx in range(2, order):
num_ = numpy.dot(num_, self.num_)
result += num_
ret = self.copy()
ret.num_ = result
return ret
[docs] def info(self):
return rot_mx_info(self)
[docs]class rot_mx_info(object):
[docs] def __init__(self, rot_mx_in):
assert rot_mx_in.den() == 1
assert rot_mx_in.type() != 0
self.type_ = rot_mx_in.type()
self.ev_ = (0, 0, 0)
self.sense_ = 0
proper_r = rot_mx_in.copy()
proper_order = self.type_
if proper_order < 0:
proper_order *= -1
proper_r = rot_mx_in.copy(minus=True)
if proper_order > 1:
rmi = proper_r.minus_unit_mx()
matrix, pivots = sympy.Matrix(rmi.num_).rref()
assert len(pivots) == 2
re_mx = numpy.array(matrix.tolist(), dtype=int)[:len(pivots)]
self.ev_ = tuple(homog_rank_2(re_mx).tolist())
self.sense_ = self.sense_of_rotation(rot_mx_in, self.type_,
self.ev_)
[docs] def type(self):
return self.type_
[docs] def ev(self):
return self.ev_
[docs] def basis_of_invariant(self):
if self.type() == 1:
basis = ((1, 0, 0), (0, 1, 0), (0, 0, 1))
elif self.type() == -2:
basis = basis_of_mirror_plane_with_normal(self.ev())
elif self.type() < 0:
basis = ()
else:
basis = (self.ev(),)
return basis
[docs] def sense_of_rotation(self, rot_mx, type_, ev):
# M.B. Boisen, Jr. & G.V. Gibbs
# Mathematical Crystallography, Revised Edition 1990
# pp. 348-349, 354-356
r = rot_mx.num()
f = -1 if type_ < 0 else 1
trace = f * (r[0] + r[4] + r[8])
if trace in [-1, 3]:
return 0
# 1-fold or 2-fold
if ev[1] == 0 and ev[2] == 0:
if ev[0] * f * r[7] > 0:
return 1
elif f * (r[3] * ev[2] - r[6] * ev[1]) > 0:
return 1
return -1
[docs] def sense(self):
return self.sense_
[docs]class rt_mx(object):
[docs] def __init__(self, *args):
if len(args) == 0:
self.r_ = rot_mx(1)
self.t_ = tr_vec()
elif len(args) == 1:
if isinstance(args[0], int):
self.r_ = rot_mx(args[0])
self.t_ = tr_vec()
elif isinstance(args[0], rot_mx):
self.r_ = args[0].copy()
self.t_ = tr_vec()
elif isinstance(args[0], tr_vec):
self.r_ = rot_mx(1)
self.t_ = args[0].copy()
elif isinstance(args[0], parse_string):
self.r_ = rot_mx()
self.t_ = tr_vec(1)
self.r_.num_, self.t_.num_ = args[0].parse_symm_op()
elif len(args) == 2:
if isinstance(args[0], int) and isinstance(args[1], int):
self.r_ = rot_mx(args[0])
self.t_ = tr_vec(args[1])
elif isinstance(args[0], rot_mx) and isinstance(args[1], int):
self.r_ = args[0].copy()
self.t_ = tr_vec(args[1])
elif isinstance(args[0], rot_mx) and isinstance(args[1], tr_vec):
self.r_ = args[0].copy()
self.t_ = args[1].copy()
elif isinstance(args[0], tr_vec) and isinstance(args[1], int):
self.r_ = rot_mx(args[1])
self.t_ = args[0].copy()
[docs] def r(self):
return self.r_
[docs] def t(self):
return self.t_
[docs] def is_unit_mx(self):
return self.r_.is_unit_mx() and self.t_.is_zero()
[docs] def t_intrinsic_part(self):
type_ = self.r().type()
mult = self.r().accumulate(type_).oper_multiply(self.t())
return mult.divide(self.r_.order(type_))
[docs] def t_location_part(self, wi):
return wi.minus(self.t())
[docs]class translation_part_info(object):
[docs] def __init__(self, rt_mx):
self.intrinsic_part_ = rt_mx.t_intrinsic_part()
self.location_part_ = rt_mx.t_location_part(self.intrinsic_part_)
[docs] def intrinsic_part(self):
return self.intrinsic_part_
[docs] def location_part(self):
return self.location_part_
[docs] def origin_shift(self):
raise NotImplementedError