Skip to content

L
libcifpp
Commit 0d503056 by maarten
Options

double support 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@322 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent 6999b7a1
Showing with 124 additions and 94 deletions
include/cif++/Point.h
...@@ -19,42 +19,60 @@ const long double ...@@ -19,42 +19,60 @@ const long double
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// Point, a location with x, y and z coordinates as float. // Point, a location with x, y and z coordinates as floating point.
// This one is derived from a tuple<float,float,float> so // This one is derived from a tuple<float,float,float> so
// you can do things like: // you can do things like:
// //
// float x, y, z; // float x, y, z;
// tie(x, y, z) = atom.loc(); // tie(x, y, z) = atom.loc();
struct Point template<typename F>
struct PointF
{ {
float mX, mY, mZ; typedef F FType;
FType mX, mY, mZ;
PointF() : mX(0), mY(0), mZ(0) {}
PointF(FType x, FType y, FType z) : mX(x), mY(y), mZ(z) {}
PointF(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {}
Point() : mX(0), mY(0), mZ(0) {} template<typename PF>
Point(float x, float y, float z) : mX(x), mY(y), mZ(z) {} PointF(const PointF<PF>& pt)
Point(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {} : mX(static_cast<F>(pt.mX))
, mY(static_cast<F>(pt.mY))
, mZ(static_cast<F>(pt.mZ)) {}
Point& 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];
mZ = rhs[2]; mZ = rhs[2];
return *this; return *this;
} }
template<typename PF>
PointF& operator=(const PointF<PF>& rhs)
{
mX = static_cast<F>(rhs.mX);
mY = static_cast<F>(rhs.mY);
mZ = static_cast<F>(rhs.mZ);
return *this;
}
float& getX() { return mX; } FType& getX() { return mX; }
float getX() const { return mX; } FType getX() const { return mX; }
void setX(float x) { mX = x; } void setX(FType x) { mX = x; }
float& getY() { return mY; } FType& getY() { return mY; }
float getY() const { return mY; } FType getY() const { return mY; }
void setY(float y) { mY = y; } void setY(FType y) { mY = y; }
float& getZ() { return mZ; } FType& getZ() { return mZ; }
float getZ() const { return mZ; } FType getZ() const { return mZ; }
void setZ(float z) { mZ = z; } void setZ(FType z) { mZ = z; }
Point& operator+=(const Point& rhs) PointF& operator+=(const PointF& rhs)
{ {
mX += rhs.mX; mX += rhs.mX;
mY += rhs.mY; mY += rhs.mY;
...@@ -63,7 +81,7 @@ struct Point ...@@ -63,7 +81,7 @@ struct Point
return *this; return *this;
} }
Point& operator+=(float d) PointF& operator+=(FType d)
{ {
mX += d; mX += d;
mY += d; mY += d;
...@@ -72,7 +90,7 @@ struct Point ...@@ -72,7 +90,7 @@ struct Point
return *this; return *this;
} }
Point& operator-=(const Point& rhs) PointF& operator-=(const PointF& rhs)
{ {
mX -= rhs.mX; mX -= rhs.mX;
mY -= rhs.mY; mY -= rhs.mY;
...@@ -81,7 +99,7 @@ struct Point ...@@ -81,7 +99,7 @@ struct Point
return *this; return *this;
} }
Point& operator-=(float d) PointF& operator-=(FType d)
{ {
mX -= d; mX -= d;
mY -= d; mY -= d;
...@@ -90,7 +108,7 @@ struct Point ...@@ -90,7 +108,7 @@ struct Point
return *this; return *this;
} }
Point& operator*=(float rhs) PointF& operator*=(FType rhs)
{ {
mX *= rhs; mX *= rhs;
mY *= rhs; mY *= rhs;
...@@ -98,7 +116,7 @@ struct Point ...@@ -98,7 +116,7 @@ struct Point
return *this; return *this;
} }
Point& operator/=(float rhs) PointF& operator/=(FType rhs)
{ {
mX /= rhs; mX /= rhs;
mY /= rhs; mY /= rhs;
...@@ -106,7 +124,7 @@ struct Point ...@@ -106,7 +124,7 @@ struct Point
return *this; return *this;
} }
float normalize() FType normalize()
{ {
auto length = mX * mX + mY * mY + mZ * mZ; auto length = mX * mX + mY * mY + mZ * mZ;
if (length > 0) if (length > 0)
...@@ -117,9 +135,9 @@ struct Point ...@@ -117,9 +135,9 @@ struct Point
return length; return length;
} }
void rotate(const boost::math::quaternion<float>& q) void rotate(const boost::math::quaternion<FType>& q)
{ {
boost::math::quaternion<float> 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);
...@@ -133,73 +151,84 @@ struct Point ...@@ -133,73 +151,84 @@ struct Point
return clipper::Coord_orth(mX, mY, mZ); return clipper::Coord_orth(mX, mY, mZ);
} }
operator std::tuple<const float&, const float&, const float&>() 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<float&,float&,float&>() 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 Point& 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?
float lengthsq() const FType lengthsq() const
{ {
return mX * mX + mY * mY + mZ * mZ; return mX * mX + mY * mY + mZ * mZ;
} }
float length() const FType length() const
{ {
return sqrt(mX * mX + mY * mY + mZ * mZ); return sqrt(mX * mX + mY * mY + mZ * mZ);
} }
}; };
inline std::ostream& operator<<(std::ostream& os, const Point& pt) typedef PointF<float> Point;
typedef PointF<double> DPoint;
template<typename F>
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;
} }
inline Point operator+(const Point& lhs, const Point& rhs) template<typename F>
inline PointF<F> operator+(const PointF<F>& lhs, const PointF<F>& rhs)
{ {
return Point(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);
} }
inline Point operator-(const Point& lhs, const Point& rhs) template<typename F>
inline PointF<F> operator-(const PointF<F>& lhs, const PointF<F>& rhs)
{ {
return Point(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);
} }
inline Point operator-(const Point& pt) template<typename F>
inline PointF<F> operator-(const PointF<F>& pt)
{ {
return Point(-pt.mX, -pt.mY, -pt.mZ); return PointF<F>(-pt.mX, -pt.mY, -pt.mZ);
} }
inline Point operator*(const Point& pt, float f) template<typename F>
inline PointF<F> operator*(const PointF<F>& pt, F f)
{ {
return Point(pt.mX * f, pt.mY * f, pt.mZ * f); return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
} }
inline Point operator*(float f, const Point& pt) template<typename F>
inline PointF<F> operator*(F f, const PointF<F>& pt)
{ {
return Point(pt.mX * f, pt.mY * f, pt.mZ * f); return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
} }
inline Point operator/(const Point& pt, float f) template<typename F>
inline PointF<F> operator/(const PointF<F>& pt, F f)
{ {
return Point(pt.mX / f, pt.mY / f, pt.mZ / f); return PointF<F>(pt.mX / f, pt.mY / f, pt.mZ / f);
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// several standard 3d operations // several standard 3d operations
inline double DistanceSquared(const Point& a, const Point& b) template<typename F>
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) +
...@@ -207,7 +236,8 @@ inline double DistanceSquared(const Point& a, const Point& b) ...@@ -207,7 +236,8 @@ inline double DistanceSquared(const Point& a, const Point& b)
(a.mZ - b.mZ) * (a.mZ - b.mZ); (a.mZ - b.mZ) * (a.mZ - b.mZ);
} }
inline double Distance(const Point& a, const Point& b) template<typename F>
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) +
...@@ -215,27 +245,69 @@ inline double Distance(const Point& a, const Point& b) ...@@ -215,27 +245,69 @@ inline double Distance(const Point& a, const Point& b)
(a.mZ - b.mZ) * (a.mZ - b.mZ)); (a.mZ - b.mZ) * (a.mZ - b.mZ));
} }
inline float DotProduct(const Point& a, const Point& b) template<typename F>
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;
} }
inline Point CrossProduct(const Point& a, const Point& b) template<typename F>
inline PointF<F> CrossProduct(const PointF<F>& a, const PointF<F>& b)
{ {
return Point(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);
} }
float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4); template<typename F>
float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& 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> v43 = p4 - p3; // vector from p3 to p4
PointF<F> z = p2 - p3; // vector from p3 to p2
PointF<F> p = CrossProduct(z, v12);
PointF<F> x = CrossProduct(z, v43);
PointF<F> y = CrossProduct(z, x);
double u = DotProduct(x, x);
double v = DotProduct(y, y);
double result = 360;
if (u > 0 and v > 0)
{
u = DotProduct(p, x) / sqrt(u);
v = DotProduct(p, y) / sqrt(v);
if (u != 0 or v != 0)
result = atan2(v, u) * 180 / kPI;
}
return result;
}
template<typename F>
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> v34 = p3 - p4;
double result = 0;
double x = DotProduct(v12, v12) * DotProduct(v34, v34);
if (x > 0)
result = DotProduct(v12, v34) / sqrt(x);
return result;
}
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// For e.g. simulated annealing, returns a new point that is moved in // For e.g. simulated annealing, returns a new point that is moved in
// 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.
Point Nudge(Point p, float offset); template<typename F>
PointF<F> Nudge(PointF<F> p, F offset);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// We use quaternions to do rotations in 3d space // We use quaternions to do rotations in 3d space
......
src/Edia.cpp
...@@ -195,7 +195,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, ...@@ -195,7 +195,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
for (auto iv = iu; iv.coord().v() <= gridMax[1]; iv.next_v()) for (auto iv = iu; iv.coord().v() <= gridMax[1]; iv.next_v())
for (auto iw = iv; iw.coord().w() <= gridMax[2]; iw.next_w()) for (auto iw = iv; iw.coord().w() <= gridMax[2]; iw.next_w())
{ {
auto p = iw.coord_orth(); Point p = iw.coord_orth();
float z = (xmap[iw] - meanDensity) / rmsDensity; float z = (xmap[iw] - meanDensity) / rmsDensity;
......
src/Point.cpp
...@@ -34,48 +34,6 @@ quaternion Normalize(quaternion q) ...@@ -34,48 +34,6 @@ quaternion Normalize(quaternion q)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
{
Point v12 = p1 - p2; // vector from p2 to p1
Point v43 = p4 - p3; // vector from p3 to p4
Point z = p2 - p3; // vector from p3 to p2
Point p = CrossProduct(z, v12);
Point x = CrossProduct(z, v43);
Point y = CrossProduct(z, x);
double u = DotProduct(x, x);
double v = DotProduct(y, y);
double result = 360;
if (u > 0 and v > 0)
{
u = DotProduct(p, x) / sqrt(u);
v = DotProduct(p, y) / sqrt(v);
if (u != 0 or v != 0)
result = atan2(v, u) * 180 / kPI;
}
return result;
}
float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
{
Point v12 = p1 - p2;
Point v34 = p3 - p4;
double result = 0;
double x = DotProduct(v12, v12) * DotProduct(v34, v34);
if (x > 0)
result = DotProduct(v12, v34) / sqrt(x);
return result;
}
// --------------------------------------------------------------------
tuple<double,Point> QuaternionToAngleAxis(quaternion q) tuple<double,Point> QuaternionToAngleAxis(quaternion q)
{ {
if (q.R_component_1() > 1) if (q.R_component_1() > 1)
......
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