Skip to content

L
libcifpp
Commit 31499b97 by Maarten L. Hekkelman
Options

Fix the 3d alignment code

parent 45f33e4b
Showing with 534 additions and 554 deletions
include/cif++/Matrix.hpp deleted 100644 → 0
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
* Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
* Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
// --------------------------------------------------------------------
// uBlas compatible matrix types
#pragma once
#include <iostream>
#include <vector>
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename T>
class MatrixBase
{
public:
using value_type = T;
virtual ~MatrixBase() {}
virtual uint32_t dim_m() const = 0;
virtual uint32_t dim_n() const = 0;
virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
MatrixBase &operator*=(const value_type &rhs);
MatrixBase &operator-=(const value_type &rhs);
};
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) *= rhs;
}
}
return *this;
}
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) -= rhs;
}
}
return *this;
}
template <typename T>
std::ostream &operator<<(std::ostream &lhs, const MatrixBase<T> &rhs)
{
lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
for (uint32_t i = 0; i < rhs.dim_m(); ++i)
{
lhs << '(';
for (uint32_t j = 0; j < rhs.dim_n(); ++j)
{
if (j > 0)
lhs << ',';
lhs << rhs(i, j);
}
lhs << ')';
}
lhs << ')';
return lhs;
}
template <typename T>
class Matrix : public MatrixBase<T>
{
public:
using value_type = T;
template <typename T2>
Matrix(const MatrixBase<T2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
{
m_data = new value_type[m_m * m_n];
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix()
: m_data(nullptr)
, m_m(0)
, m_n(0)
{
}
Matrix(const Matrix &m)
: m_m(m.m_m)
, m_n(m.m_n)
{
m_data = new value_type[m_m * m_n];
std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
}
Matrix &operator=(const Matrix &m)
{
value_type *t = new value_type[m.m_m * m.m_n];
std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
delete[] m_data;
m_data = t;
m_m = m.m_m;
m_n = m.m_n;
return *this;
}
Matrix(uint32_t m, uint32_t n, T v = T())
: m_m(m)
, m_n(n)
{
m_data = new value_type[m_m * m_n];
std::fill(m_data, m_data + (m_m * m_n), v);
}
virtual ~Matrix()
{
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_m; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
virtual value_type &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
template <typename Func>
void each(Func f)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
f(m_data[i]);
}
template <typename U>
Matrix &operator/=(U v)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
m_data[i] /= v;
return *this;
}
private:
value_type *m_data;
uint32_t m_m, m_n;
};
// --------------------------------------------------------------------
template <typename T>
class SymmetricMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
SymmetricMatrix(uint32_t n, T v = T())
: m_owner(true)
, m_n(n)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
m_data = new value_type[N];
std::fill(m_data, m_data + N, v);
}
SymmetricMatrix(const T *data, uint32_t n)
: m_owner(false)
, m_data(const_cast<T *>(data))
, m_n(n)
{
}
virtual ~SymmetricMatrix()
{
if (m_owner)
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
T operator()(uint32_t i, uint32_t j) const;
virtual T &operator()(uint32_t i, uint32_t j);
// erase two rows, add one at the end (for neighbour joining)
void erase_2(uint32_t i, uint32_t j);
template <typename Func>
void each(Func f)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
f(m_data[i]);
}
template <typename U>
SymmetricMatrix &operator/=(U v)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
m_data[i] /= v;
return *this;
}
private:
bool m_owner;
value_type *m_data;
uint32_t m_n;
};
template <typename T>
inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
template <typename T>
inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
template <typename T>
void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
{
uint32_t s = 0, d = 0;
for (uint32_t i = 0; i < m_n; ++i)
{
for (uint32_t j = 0; j < i; ++j)
{
if (i != di and j != dj and i != dj and j != di)
{
if (s != d)
m_data[d] = m_data[s];
++d;
}
++s;
}
}
--m_n;
}
template <typename T>
class IdentityMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
value_type result = 0;
if (i == j)
result = 1;
return result;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
result(i, j) += lhs(li, j) * rhs(i, rj);
}
}
return result;
}
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs);
result *= rhs;
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
result(i, j) = lhs(i, j) - rhs(i, j);
}
}
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs.dim_m(), lhs.dim_n());
result -= rhs;
return result;
}
// template <typename T>
// symmetric_matrix<T> hammingDistance(const MatrixBase<T> &lhs, T rhs);
// template <typename T>
// std::vector<T> sum(const MatrixBase<T> &m);
include/cif++/Point.hpp
/*- /*-
* SPDX-License-Identifier: BSD-2-Clause * SPDX-License-Identifier: BSD-2-Clause
* *
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
* *
* Redistribution and use in source and binary forms, with or without * Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met: * modification, are permitted provided that the following conditions are met:
* *
* 1. Redistributions of source code must retain the above copyright notice, this * 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer * list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice, * 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation * this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution. * and/or other materials provided with the distribution.
* *
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
...@@ -37,7 +37,7 @@ ...@@ -37,7 +37,7 @@
namespace mmcif namespace mmcif
{ {
typedef boost::math::quaternion<float> Quaternion; typedef boost::math::quaternion<float> Quaternion;
const double const double
kPI = 3.141592653589793238462643383279502884; kPI = 3.141592653589793238462643383279502884;
...@@ -51,26 +51,43 @@ const double ...@@ -51,26 +51,43 @@ const double
// float x, y, z; // float x, y, z;
// tie(x, y, z) = atom.loc(); // tie(x, y, z) = atom.loc();
template<typename F> template <typename F>
struct PointF struct PointF
{ {
typedef F FType; typedef F FType;
FType mX, mY, mZ; FType mX, mY, mZ;
PointF() : mX(0), mY(0), mZ(0) {}
PointF(FType x, FType y, FType z) : mX(x), mY(y), mZ(z) {}
template<typename PF> PointF()
PointF(const PointF<PF>& pt) : mX(0)
, mY(0)
, mZ(0)
{
}
PointF(FType x, FType y, FType z)
: mX(x)
, mY(y)
, mZ(z)
{
}
template <typename PF>
PointF(const PointF<PF> &pt)
: mX(static_cast<F>(pt.mX)) : mX(static_cast<F>(pt.mX))
, mY(static_cast<F>(pt.mY)) , mY(static_cast<F>(pt.mY))
, mZ(static_cast<F>(pt.mZ)) {} , mZ(static_cast<F>(pt.mZ))
{
}
#if HAVE_LIBCLIPPER #if HAVE_LIBCLIPPER
PointF(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {} PointF(const clipper::Coord_orth &pt)
: mX(pt[0])
, mY(pt[1])
, mZ(pt[2])
{
}
PointF& operator=(const clipper::Coord_orth& rhs) PointF &operator=(const clipper::Coord_orth &rhs)
{ {
mX = rhs[0]; mX = rhs[0];
mY = rhs[1]; mY = rhs[1];
...@@ -79,72 +96,72 @@ struct PointF ...@@ -79,72 +96,72 @@ struct PointF
} }
#endif #endif
template<typename PF> template <typename PF>
PointF& operator=(const PointF<PF>& rhs) PointF &operator=(const PointF<PF> &rhs)
{ {
mX = static_cast<F>(rhs.mX); mX = static_cast<F>(rhs.mX);
mY = static_cast<F>(rhs.mY); mY = static_cast<F>(rhs.mY);
mZ = static_cast<F>(rhs.mZ); mZ = static_cast<F>(rhs.mZ);
return *this; return *this;
} }
FType& getX() { return mX; } FType &getX() { return mX; }
FType getX() const { return mX; } FType getX() const { return mX; }
void setX(FType x) { mX = x; } void setX(FType x) { mX = x; }
FType& getY() { return mY; } FType &getY() { return mY; }
FType getY() const { return mY; } FType getY() const { return mY; }
void setY(FType y) { mY = y; } void setY(FType y) { mY = y; }
FType& getZ() { return mZ; } FType &getZ() { return mZ; }
FType getZ() const { return mZ; } FType getZ() const { return mZ; }
void setZ(FType z) { mZ = z; } void setZ(FType z) { mZ = z; }
PointF& operator+=(const PointF& rhs) PointF &operator+=(const PointF &rhs)
{ {
mX += rhs.mX; mX += rhs.mX;
mY += rhs.mY; mY += rhs.mY;
mZ += rhs.mZ; mZ += rhs.mZ;
return *this; return *this;
} }
PointF& operator+=(FType d) PointF &operator+=(FType d)
{ {
mX += d; mX += d;
mY += d; mY += d;
mZ += d; mZ += d;
return *this; return *this;
} }
PointF& operator-=(const PointF& rhs) PointF &operator-=(const PointF &rhs)
{ {
mX -= rhs.mX; mX -= rhs.mX;
mY -= rhs.mY; mY -= rhs.mY;
mZ -= rhs.mZ; mZ -= rhs.mZ;
return *this; return *this;
} }
PointF& operator-=(FType d) PointF &operator-=(FType d)
{ {
mX -= d; mX -= d;
mY -= d; mY -= d;
mZ -= d; mZ -= d;
return *this; return *this;
} }
PointF& operator*=(FType rhs) PointF &operator*=(FType rhs)
{ {
mX *= rhs; mX *= rhs;
mY *= rhs; mY *= rhs;
mZ *= rhs; mZ *= rhs;
return *this; return *this;
} }
PointF& operator/=(FType rhs) PointF &operator/=(FType rhs)
{ {
mX /= rhs; mX /= rhs;
mY /= rhs; mY /= rhs;
...@@ -162,18 +179,18 @@ struct PointF ...@@ -162,18 +179,18 @@ struct PointF
} }
return length; return length;
} }
void rotate(const boost::math::quaternion<FType>& q) void rotate(const boost::math::quaternion<FType> &q)
{ {
boost::math::quaternion<FType> p(0, mX, mY, mZ); boost::math::quaternion<FType> p(0, mX, mY, mZ);
p = q * p * boost::math::conj(q); p = q * p * boost::math::conj(q);
mX = p.R_component_2(); mX = p.R_component_2();
mY = p.R_component_3(); mY = p.R_component_3();
mZ = p.R_component_4(); mZ = p.R_component_4();
} }
#if HAVE_LIBCLIPPER #if HAVE_LIBCLIPPER
operator clipper::Coord_orth() const operator clipper::Coord_orth() const
{ {
...@@ -181,21 +198,21 @@ struct PointF ...@@ -181,21 +198,21 @@ struct PointF
} }
#endif #endif
operator std::tuple<const FType&, const FType&, const FType&>() const operator std::tuple<const FType &, const FType &, const FType &>() const
{ {
return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ)); return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
} }
operator std::tuple<FType&,FType&,FType&>() operator std::tuple<FType &, FType &, FType &>()
{ {
return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ)); return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
} }
bool operator==(const PointF& rhs) const bool operator==(const PointF &rhs) const
{ {
return mX == rhs.mX and mY == rhs.mY and mZ == rhs.mZ; return mX == rhs.mX and mY == rhs.mY and mZ == rhs.mZ;
} }
// consider point as a vector... perhaps I should rename Point? // consider point as a vector... perhaps I should rename Point?
FType lengthsq() const FType lengthsq() const
{ {
...@@ -211,45 +228,45 @@ struct PointF ...@@ -211,45 +228,45 @@ struct PointF
typedef PointF<float> Point; typedef PointF<float> Point;
typedef PointF<double> DPoint; typedef PointF<double> DPoint;
template<typename F> template <typename F>
inline std::ostream& operator<<(std::ostream& os, const PointF<F>& pt) inline std::ostream &operator<<(std::ostream &os, const PointF<F> &pt)
{ {
os << '(' << pt.mX << ',' << pt.mY << ',' << pt.mZ << ')'; os << '(' << pt.mX << ',' << pt.mY << ',' << pt.mZ << ')';
return os; return os;
} }
template<typename F> template <typename F>
inline PointF<F> operator+(const PointF<F>& lhs, const PointF<F>& rhs) inline PointF<F> operator+(const PointF<F> &lhs, const PointF<F> &rhs)
{ {
return PointF<F>(lhs.mX + rhs.mX, lhs.mY + rhs.mY, lhs.mZ + rhs.mZ); return PointF<F>(lhs.mX + rhs.mX, lhs.mY + rhs.mY, lhs.mZ + rhs.mZ);
} }
template<typename F> template <typename F>
inline PointF<F> operator-(const PointF<F>& lhs, const PointF<F>& rhs) inline PointF<F> operator-(const PointF<F> &lhs, const PointF<F> &rhs)
{ {
return PointF<F>(lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ); return PointF<F>(lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ);
} }
template<typename F> template <typename F>
inline PointF<F> operator-(const PointF<F>& pt) inline PointF<F> operator-(const PointF<F> &pt)
{ {
return PointF<F>(-pt.mX, -pt.mY, -pt.mZ); return PointF<F>(-pt.mX, -pt.mY, -pt.mZ);
} }
template<typename F> template <typename F>
inline PointF<F> operator*(const PointF<F>& pt, F f) inline PointF<F> operator*(const PointF<F> &pt, F f)
{ {
return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f); return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
} }
template<typename F> template <typename F>
inline PointF<F> operator*(F f, const PointF<F>& pt) inline PointF<F> operator*(F f, const PointF<F> &pt)
{ {
return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f); return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
} }
template<typename F> template <typename F>
inline PointF<F> operator/(const PointF<F>& pt, F f) inline PointF<F> operator/(const PointF<F> &pt, F f)
{ {
return PointF<F>(pt.mX / f, pt.mY / f, pt.mZ / f); return PointF<F>(pt.mX / f, pt.mY / f, pt.mZ / f);
} }
...@@ -257,17 +274,16 @@ inline PointF<F> operator/(const PointF<F>& pt, F f) ...@@ -257,17 +274,16 @@ inline PointF<F> operator/(const PointF<F>& pt, F f)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// several standard 3d operations // several standard 3d operations
template<typename F> template <typename F>
inline double DistanceSquared(const PointF<F>& a, const PointF<F>& b) inline double DistanceSquared(const PointF<F> &a, const PointF<F> &b)
{ {
return return (a.mX - b.mX) * (a.mX - b.mX) +
(a.mX - b.mX) * (a.mX - b.mX) + (a.mY - b.mY) * (a.mY - b.mY) +
(a.mY - b.mY) * (a.mY - b.mY) + (a.mZ - b.mZ) * (a.mZ - b.mZ);
(a.mZ - b.mZ) * (a.mZ - b.mZ);
} }
template<typename F> template <typename F>
inline double Distance(const PointF<F>& a, const PointF<F>& b) inline double Distance(const PointF<F> &a, const PointF<F> &b)
{ {
return sqrt( return sqrt(
(a.mX - b.mX) * (a.mX - b.mX) + (a.mX - b.mX) * (a.mX - b.mX) +
...@@ -275,44 +291,44 @@ inline double Distance(const PointF<F>& a, const PointF<F>& b) ...@@ -275,44 +291,44 @@ inline double Distance(const PointF<F>& a, const PointF<F>& b)
(a.mZ - b.mZ) * (a.mZ - b.mZ)); (a.mZ - b.mZ) * (a.mZ - b.mZ));
} }
template<typename F> template <typename F>
inline F DotProduct(const PointF<F>& a, const PointF<F>& b) inline F DotProduct(const PointF<F> &a, const PointF<F> &b)
{ {
return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ; return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ;
} }
template<typename F> template <typename F>
inline PointF<F> CrossProduct(const PointF<F>& a, const PointF<F>& b) inline PointF<F> CrossProduct(const PointF<F> &a, const PointF<F> &b)
{ {
return PointF<F>(a.mY * b.mZ - b.mY * a.mZ, return PointF<F>(a.mY * b.mZ - b.mY * a.mZ,
a.mZ * b.mX - b.mZ * a.mX, a.mZ * b.mX - b.mZ * a.mX,
a.mX * b.mY - b.mX * a.mY); a.mX * b.mY - b.mX * a.mY);
} }
template<typename F> template <typename F>
double Angle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3) double Angle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3)
{ {
PointF<F> v1 = p1 - p2; PointF<F> v1 = p1 - p2;
PointF<F> v2 = p3 - p2; PointF<F> v2 = p3 - p2;
return std::acos(DotProduct(v1, v2) / (v1.length() * v2.length())) * 180 / kPI; return std::acos(DotProduct(v1, v2) / (v1.length() * v2.length())) * 180 / kPI;
} }
template<typename F> template <typename F>
double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4) double DihedralAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
{ {
PointF<F> v12 = p1 - p2; // vector from p2 to p1 PointF<F> v12 = p1 - p2; // vector from p2 to p1
PointF<F> v43 = p4 - p3; // vector from p3 to p4 PointF<F> v43 = p4 - p3; // vector from p3 to p4
PointF<F> z = p2 - p3; // vector from p3 to p2 PointF<F> z = p2 - p3; // vector from p3 to p2
PointF<F> p = CrossProduct(z, v12); PointF<F> p = CrossProduct(z, v12);
PointF<F> x = CrossProduct(z, v43); PointF<F> x = CrossProduct(z, v43);
PointF<F> y = CrossProduct(z, x); PointF<F> y = CrossProduct(z, x);
double u = DotProduct(x, x); double u = DotProduct(x, x);
double v = DotProduct(y, y); double v = DotProduct(y, y);
double result = 360; double result = 360;
if (u > 0 and v > 0) if (u > 0 and v > 0)
{ {
...@@ -321,33 +337,33 @@ double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& ...@@ -321,33 +337,33 @@ double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>&
if (u != 0 or v != 0) if (u != 0 or v != 0)
result = atan2(v, u) * 180 / kPI; result = atan2(v, u) * 180 / kPI;
} }
return result; return result;
} }
template<typename F> template <typename F>
double CosinusAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4) double CosinusAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
{ {
PointF<F> v12 = p1 - p2; PointF<F> v12 = p1 - p2;
PointF<F> v34 = p3 - p4; PointF<F> v34 = p3 - p4;
double result = 0; double result = 0;
double x = DotProduct(v12, v12) * DotProduct(v34, v34); double x = DotProduct(v12, v12) * DotProduct(v34, v34);
if (x > 0) if (x > 0)
result = DotProduct(v12, v34) / sqrt(x); result = DotProduct(v12, v34) / sqrt(x);
return result; return result;
} }
template<typename F> template <typename F>
auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<F> &p) auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<F> &p)
{ {
auto line = l2 - l1; auto line = l2 - l1;
auto p_to_l1 = p - l1; auto p_to_l1 = p - l1;
auto p_to_l2 = p - l2; auto p_to_l2 = p - l2;
auto cross = CrossProduct(p_to_l1, p_to_l2); auto cross = CrossProduct(p_to_l1, p_to_l2);
return cross.length() / line.length(); return cross.length() / line.length();
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -355,7 +371,7 @@ auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF< ...@@ -355,7 +371,7 @@ auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<
// a random direction with a distance randomly chosen from a normal // a random direction with a distance randomly chosen from a normal
// distribution with a stddev of offset. // distribution with a stddev of offset.
template<typename F> template <typename F>
PointF<F> Nudge(PointF<F> p, F offset); PointF<F> Nudge(PointF<F> p, F offset);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -363,66 +379,77 @@ PointF<F> Nudge(PointF<F> p, F offset); ...@@ -363,66 +379,77 @@ PointF<F> Nudge(PointF<F> p, F offset);
Quaternion Normalize(Quaternion q); Quaternion Normalize(Quaternion q);
std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q); std::tuple<double, Point> QuaternionToAngleAxis(Quaternion q);
Point Centroid(std::vector<Point>& Points); Point Centroid(std::vector<Point> &Points);
Point CenterPoints(std::vector<Point>& Points); Point CenterPoints(std::vector<Point> &Points);
Quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
double RMSd(const std::vector<Point>& a, const std::vector<Point>& b); /// \brief Returns how the two sets of points \a a and \b b can be aligned
///
/// \param a The first set of points
/// \param b The second set of points
/// \result The quaternion which should be applied to the points in \a a to
/// obtain the best superposition.
Quaternion AlignPoints(const std::vector<Point> &a, const std::vector<Point> &b);
/// \brief The RMSd for the points in \a a and \a b
double RMSd(const std::vector<Point> &a, const std::vector<Point> &b);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// Helper class to generate evenly divided Points on a sphere // Helper class to generate evenly divided Points on a sphere
// we use a fibonacci sphere to calculate even distribution of the dots // we use a fibonacci sphere to calculate even distribution of the dots
template<int N> template <int N>
class SphericalDots class SphericalDots
{ {
public: public:
enum { P = 2 * N + 1 }; enum
typedef typename std::array<Point,P> array_type; {
typedef typename array_type::const_iterator iterator; P = 2 * N + 1
};
typedef typename std::array<Point, P> array_type;
typedef typename array_type::const_iterator iterator;
static SphericalDots& instance() static SphericalDots &instance()
{ {
static SphericalDots sInstance; static SphericalDots sInstance;
return sInstance; return sInstance;
} }
size_t size() const { return mPoints.size(); }
const Point operator[](uint32_t inIx) const { return mPoints[inIx]; }
iterator begin() const { return mPoints.begin(); }
iterator end() const { return mPoints.end(); }
double weight() const { return mWeight; } size_t size() const { return mPoints.size(); }
const Point operator[](uint32_t inIx) const { return mPoints[inIx]; }
iterator begin() const { return mPoints.begin(); }
iterator end() const { return mPoints.end(); }
double weight() const { return mWeight; }
SphericalDots() SphericalDots()
{ {
const double const double
kGoldenRatio = (1 + std::sqrt(5.0)) / 2; kGoldenRatio = (1 + std::sqrt(5.0)) / 2;
mWeight = (4 * kPI) / P; mWeight = (4 * kPI) / P;
auto p = mPoints.begin(); auto p = mPoints.begin();
for (int32_t i = -N; i <= N; ++i) for (int32_t i = -N; i <= N; ++i)
{ {
double lat = std::asin((2.0 * i) / P); double lat = std::asin((2.0 * i) / P);
double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio; double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio;
p->mX = sin(lon) * cos(lat); p->mX = sin(lon) * cos(lat);
p->mY = cos(lon) * cos(lat); p->mY = cos(lon) * cos(lat);
p->mZ = sin(lat); p->mZ = sin(lat);
++p; ++p;
} }
} }
private: private:
array_type mPoints;
array_type mPoints; double mWeight;
double mWeight;
}; };
typedef SphericalDots<50> SphericalDots_50; typedef SphericalDots<50> SphericalDots_50;
} } // namespace mmcif
src/Point.cpp
...@@ -28,12 +28,378 @@ ...@@ -28,12 +28,378 @@
#include <valarray> #include <valarray>
#include "cif++/Point.hpp" #include "cif++/Point.hpp"
#include "cif++/Matrix.hpp"
namespace mmcif namespace mmcif
{ {
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// uBlas compatible matrix types
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename T>
class MatrixBase
{
public:
using value_type = T;
virtual ~MatrixBase() {}
virtual uint32_t dim_m() const = 0;
virtual uint32_t dim_n() const = 0;
virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
MatrixBase &operator*=(const value_type &rhs);
MatrixBase &operator-=(const value_type &rhs);
};
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) *= rhs;
}
}
return *this;
}
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) -= rhs;
}
}
return *this;
}
template <typename T>
class Matrix : public MatrixBase<T>
{
public:
using value_type = T;
template <typename T2>
Matrix(const MatrixBase<T2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
{
m_data = new value_type[m_m * m_n];
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix()
: m_data(nullptr)
, m_m(0)
, m_n(0)
{
}
Matrix(const Matrix &m)
: m_m(m.m_m)
, m_n(m.m_n)
{
m_data = new value_type[m_m * m_n];
std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
}
Matrix &operator=(const Matrix &m)
{
value_type *t = new value_type[m.m_m * m.m_n];
std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
delete[] m_data;
m_data = t;
m_m = m.m_m;
m_n = m.m_n;
return *this;
}
Matrix(uint32_t m, uint32_t n, T v = T())
: m_m(m)
, m_n(n)
{
m_data = new value_type[m_m * m_n];
std::fill(m_data, m_data + (m_m * m_n), v);
}
virtual ~Matrix()
{
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_m; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
virtual value_type &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
template <typename Func>
void each(Func f)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
f(m_data[i]);
}
template <typename U>
Matrix &operator/=(U v)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
m_data[i] /= v;
return *this;
}
private:
value_type *m_data;
uint32_t m_m, m_n;
};
// --------------------------------------------------------------------
template <typename T>
class SymmetricMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
SymmetricMatrix(uint32_t n, T v = T())
: m_owner(true)
, m_n(n)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
m_data = new value_type[N];
std::fill(m_data, m_data + N, v);
}
SymmetricMatrix(const T *data, uint32_t n)
: m_owner(false)
, m_data(const_cast<T *>(data))
, m_n(n)
{
}
virtual ~SymmetricMatrix()
{
if (m_owner)
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
T operator()(uint32_t i, uint32_t j) const;
virtual T &operator()(uint32_t i, uint32_t j);
// erase two rows, add one at the end (for neighbour joining)
void erase_2(uint32_t i, uint32_t j);
template <typename Func>
void each(Func f)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
f(m_data[i]);
}
template <typename U>
SymmetricMatrix &operator/=(U v)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
m_data[i] /= v;
return *this;
}
private:
bool m_owner;
value_type *m_data;
uint32_t m_n;
};
template <typename T>
inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
template <typename T>
inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
template <typename T>
void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
{
uint32_t s = 0, d = 0;
for (uint32_t i = 0; i < m_n; ++i)
{
for (uint32_t j = 0; j < i; ++j)
{
if (i != di and j != dj and i != dj and j != di)
{
if (s != d)
m_data[d] = m_data[s];
++d;
}
++s;
}
}
--m_n;
}
template <typename T>
class IdentityMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
value_type result = 0;
if (i == j)
result = 1;
return result;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
result(i, j) += lhs(li, j) * rhs(i, rj);
}
}
return result;
}
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs);
result *= rhs;
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
result(i, j) = lhs(i, j) - rhs(i, j);
}
}
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs.dim_m(), lhs.dim_n());
result -= rhs;
return result;
}
template<class M1, typename T>
void cofactors(const M1& m, SymmetricMatrix<T>& cf)
{
const size_t ixs[4][3] =
{
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
for (size_t x = 0; x < 4; ++x)
{
const size_t* ix = ixs[x];
for (size_t y = x; y < 4; ++y)
{
const size_t* iy = ixs[y];
cf(x, y) =
m(ix[0], iy[0]) * m(ix[1], iy[1]) * m(ix[2], iy[2]) +
m(ix[0], iy[1]) * m(ix[1], iy[2]) * m(ix[2], iy[0]) +
m(ix[0], iy[2]) * m(ix[1], iy[0]) * m(ix[2], iy[1]) -
m(ix[0], iy[2]) * m(ix[1], iy[1]) * m(ix[2], iy[0]) -
m(ix[0], iy[1]) * m(ix[1], iy[0]) * m(ix[2], iy[2]) -
m(ix[0], iy[0]) * m(ix[1], iy[2]) * m(ix[2], iy[1]);
if ((x + y) % 2 == 1)
cf(x, y) *= -1;
}
}
}
// --------------------------------------------------------------------
Quaternion Normalize(Quaternion q) Quaternion Normalize(Quaternion q)
{ {
...@@ -225,6 +591,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p ...@@ -225,6 +591,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
M(1, 2) * M(2, 0) * M(0, 1) + M(1, 2) * M(2, 0) * M(0, 1) +
M(2, 1) * M(1, 0) * M(0, 2)); M(2, 1) * M(1, 0) * M(0, 2));
// E is the determinant of N:
double E = double E =
(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) + (N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) + (N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
...@@ -240,41 +607,18 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p ...@@ -240,41 +607,18 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
Matrix<double> li = IdentityMatrix<double>(4) * lm; Matrix<double> li = IdentityMatrix<double>(4) * lm;
Matrix<double> t = N - li; Matrix<double> t = N - li;
// calculate a Matrix of cofactors for t // calculate a Matrix of cofactors for t, since N is symmetric, t must be symmetric as well and so will be cf
Matrix<double> cf(4, 4); SymmetricMatrix<double> cf(4);
cofactors(t, cf);
const uint32_t ixs[4][3] = int maxR = 0;
{ for (int r = 1; r < 4; ++r)
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
uint32_t maxR = 0;
for (uint32_t r = 0; r < 4; ++r)
{ {
const uint32_t* ir = ixs[r]; if (cf(r, 0) > cf(maxR, 0))
for (uint32_t c = 0; c < 4; ++c)
{
const uint32_t* ic = ixs[c];
cf(r, c) =
t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
}
if (r > maxR and cf(r, 0) > cf(maxR, 0))
maxR = r; maxR = r;
} }
// NOTE the negation of the y here, why? Maybe I swapped r/c above? Quaternion q(cf(maxR, 0), cf(maxR, 1), cf(maxR, 2), cf(maxR, 3));
Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
q = Normalize(q); q = Normalize(q);
return q; return q;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment