src/openbabel/vector3.cpp

/**********************************************************************
vector3.cpp - Handle 3D coordinates.
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#include <math.h>

#include "mol.hpp"
#include "vector3.hpp"

using namespace std;

namespace OpenBabel
{

/*! \class vector3
   \brief Represents a vector in the 3-dimensional real space.

The vector3 class was designed to simplify operations with floating
point coordinates. To this end many of the common operations have been
overloaded for simplicity. Vector addition, subtraction, scalar
multiplication, dot product, cross product, magnitude and a number of
other utility functions are built in to the vector class. For a full
description of the class member functions please consult the header
file vector3.h. The following code demonstrates several of the
functions of the vector class:
\code
vector3 v1,v2,v3;
v1 = VX;
v2 = VY;
v3 = cross(v1,v2);
v3 *= 2.5;
v3.normalize();
\endcode
*/

/*! This (slow) method allows to access the elements of the
  vector as if it were an array of doubles. If the index is > 2,
  then a warning is printed, and the program is terminated via
  exit(-1). Otherwise, if i is 0, 1 or 2, then a reference to x,
  y or z is returned, respectively.
  
  \warning This method is primarily designed to facilitate the
  integration ('Open Babelization') of code that uses arrays of
  doubles rather than the vector class. Due to the error checks
  the method is of course very slow and should therefore be
  avoided in production code.
*/
double& vector3::operator[] ( unsigned int i)
{
    if (i > 2)
    {
        cerr << "ERROR in OpenBabel::vector3::operator[]" << endl
        << "The method has been called with an illegal index i=" << i << "." << endl
        << "Please contact the author of the offending program immediately." << endl;
        exit(-1);
    }
    if (i == 0)
        return _vx;
    if (i == 1)
        return _vy;
    return _vz;
}

/*! replaces *this with a random unit vector, which is (supposed
  to be) uniformly distributed over the unit sphere. Uses the
  random number generator obRand, or uses the system number
  generator with a time seed if obRand == NULL.
     
  @param obRandP random number generator to use, or 0L, if the
  system random number generator (with time seed) should be used
*/
void vector3::randomUnitVector(OBRandom *obRandP)
{
    OBRandom *ptr;
    if (!obRandP)
    {
        ptr = new OBRandom(true);
        ptr->TimeSeed();
    }
    else
        ptr = obRandP;

    // obtain a random vector with 0.001 <= length^2 <= 1.0, normalize
    // the vector to obtain a random vector of length 1.0.
    double l;
    do
    {
        this->Set(ptr->NextFloat()-0.5, ptr->NextFloat()-0.5, ptr->NextFloat()-0.5);
        l = length_2();
    }
    while ( (l > 1.0) || (l < 1e-4) );
    this->normalize();

    if (!obRandP)
        delete ptr;
}

OBAPI ostream& operator<< ( ostream& co, const vector3& v )
{
    co << "< " << v._vx << ", " << v._vy << ", " << v._vz << " >" ;
    return co ;
}

OBAPI int operator== ( const vector3& v1, const vector3& v2 )
{
    if ( ( v1._vx == v2._vx ) &&
            ( v1._vy == v2._vy ) &&
            ( v1._vz == v2._vz ) )
        return ( true ) ;
    else
        return ( false ) ;
}

OBAPI int operator!= ( const vector3& v1, const vector3& v2 )
{
    if ( ( v1._vx != v2._vx ) ||
            ( v1._vy != v2._vy ) ||
            ( v1._vz != v2._vz ) )
        return ( true ) ;
    else
        return ( false ) ;
}

/*! This method checks if the current vector has length() ==
  0.0.  If so, *this remains unchanged. Otherwise, *this is
  scaled by 1.0/length().

  \warning If length() is very close to zero, but not == 0.0,
  this method may behave in unexpected ways and return almost
  random results; details may depend on your particular floating
  point implementation. The use of this method is therefore
  highly discouraged, unless you are certain that length() is in
  a reasonable range, away from 0.0 (Stefan Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version.

  @returns a reference to *this
*/
vector3& vector3 :: normalize ()
{
    double l = length ();

    if (IsNearZero(l))
        return(*this);

    _vx = _vx / l ;
    _vy = _vy / l ;
    _vz = _vz / l ;

    return(*this);
}

OBAPI double dot ( const vector3& v1, const vector3& v2 )
{
    return v1._vx*v2._vx + v1._vy*v2._vy + v1._vz*v2._vz ;
}

OBAPI vector3 cross ( const vector3& v1, const vector3& v2 )
{
    vector3 vv ;

    vv._vx =   v1._vy*v2._vz - v1._vz*v2._vy ;
    vv._vy = - v1._vx*v2._vz + v1._vz*v2._vx ;
    vv._vz =   v1._vx*v2._vy - v1._vy*v2._vx ;

    return ( vv ) ;
}


/*! This method calculates the angle between two vectors
     
  \warning If length() of any of the two vectors is == 0.0,
  this method will divide by zero. If the product of the
  length() of the two vectors is very close to 0.0, but not ==
  0.0, this method may behave in unexpected ways and return
  almost random results; details may depend on your particular
  floating point implementation. The use of this method is
  therefore highly discouraged, unless you are certain that the
  length()es are in a reasonable range, away from 0.0 (Stefan
  Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version.

  @returns the angle in degrees (0-360)
*/
OBAPI double vectorAngle ( const vector3& v1, const vector3& v2 )
{
    double mag;
    double dp;

    mag = v1.length() * v2.length();
    dp = dot(v1,v2)/mag;

    if (dp < -0.999999)
        dp = -0.9999999;

    if (dp > 0.9999999)
        dp = 0.9999999;

    if (dp > 1.0)
        dp = 1.0;

    return((RAD_TO_DEG * acos(dp)));
}

/*!  This function calculates the torsion angle of three vectors, represented
  by four points A--B--C--D, i.e. B and C are vertexes, but none of A--B,
  B--C, and C--D are colinear.  A "torsion angle" is the amount of "twist"
  or torsion needed around the B--C axis to bring A--B into the same plane
  as B--C--D.  The torsion is measured by "looking down" the vector B--C so
  that B is superimposed on C, then noting how far you'd have to rotate
  A--B to superimpose A over D.  Angles are + in the anticlockwise
  direction.  The operation is symmetrical in that if you reverse the image
  (look from C to B and rotate D over A), you get the same answer.
*/

OBAPI double CalcTorsionAngle(const vector3 &a, const vector3 &b,
                        const vector3 &c, const vector3 &d)
{
    double torsion;
    vector3 b1,b2,b3,c1,c2,c3;

    b1 = a - b;
    b2 = b - c;
    b3 = c - d;

    c1 = cross(b1,b2);
    c2 = cross(b2,b3);
    c3 = cross(c1,c2);

    if (c1.length() * c2.length() < 0.001)
        torsion = 0.0;
    else
    {
        torsion = vectorAngle(c1,c2);
        if (dot(b2,c3) > 0.0)
            torsion *= -1.0;
    }

    return(torsion);
}

/*! This method checks if the current vector *this is zero
  (i.e. if all entries == 0.0). If so, a warning message is
  printed, and the whole program is aborted with exit(0).
  Otherwise, a vector of length one is generated, which is
  orthogonal to *this, and stored in v. The resulting vector is
  not random.

  \warning If the entries of the *this (in particular the
  z-component) are very close to zero, but not == 0.0, this
  method may behave in unexpected ways and return almost random
  results; details may depend on your particular floating point
  implementation. The use of this method is therefore highly
  discouraged, unless you are certain that all components of
  *this are in a reasonable range, away from 0.0 (Stefan
  Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version that
  does not call exit()

  @param res a reference to a vector where the result will be
  stored
*/
void vector3::createOrthoVector(vector3 &res) const
{
    vector3 cO;

    if ( ( IsNearZero(this->x())) && (IsNearZero(this->y())) )
    {
        if ( IsNearZero(this->z()) )
        {
            cerr << "makeorthovec zero vector" << endl;
            exit(0);
        }
        cO.SetX(1.0);
    }
    else
    {
        cO.SetZ(1.0);
    }
    res= cross(cO,*this);
    res.normalize();
}

const vector3 VZero ( 0.0, 0.0, 0.0 ) ;
const vector3 VX    ( 1.0, 0.0, 0.0 ) ;
const vector3 VY    ( 0.0, 1.0, 0.0 ) ;
const vector3 VZ    ( 0.0, 0.0, 1.0 ) ;

/* Calculate the distance of point a to the plane determined by b,c,d */
double Point2Plane(vector3 a, vector3 b, vector3 c, vector3 d)
{
  double angle =0;
  double dist_ab =0;
  vector3 v_ba = a-b;
  vector3 v_normal = cross(c-b, d-b).normalize();
  angle = vectorAngle(v_normal, v_ba);
  dist_ab = v_ba.length();
  return fabs(dist_ab * cos(DEG_TO_RAD * angle));
}

} // namespace OpenBabel

//! \file vector3.cpp
//! \brief Handle 3D coordinates.
Revision:	3057
Committed:	Thu Oct 19 20:49:05 2006 UTC (17 years, 9 months ago) by gezelter
File size:	9733 byte(s)
Log Message:	updated OpenBabel to version 2.0.2
#	User	Rev	Content
1	tim	2440	/**********************************************************************
2			vector3.cpp - Handle 3D coordinates.
3
4			Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5			Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6
7			This file is part of the Open Babel project.
8			For more information, see <http://openbabel.sourceforge.net/>
9
10			This program is free software; you can redistribute it and/or modify
11			it under the terms of the GNU General Public License as published by
12			the Free Software Foundation version 2 of the License.
13
14			This program is distributed in the hope that it will be useful,
15			but WITHOUT ANY WARRANTY; without even the implied warranty of
16			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17			GNU General Public License for more details.
18			***********************************************************************/
19
20			#include <math.h>
21
22			#include "mol.hpp"
23			#include "vector3.hpp"
24
25			using namespace std;
26
27			namespace OpenBabel
28			{
29
30			/*! \class vector3
31			\brief Represents a vector in the 3-dimensional real space.
32
33	gezelter	3057	The vector3 class was designed to simplify operations with floating
34	tim	2440	point coordinates. To this end many of the common operations have been
35			overloaded for simplicity. Vector addition, subtraction, scalar
36			multiplication, dot product, cross product, magnitude and a number of
37			other utility functions are built in to the vector class. For a full
38			description of the class member functions please consult the header
39			file vector3.h. The following code demonstrates several of the
40			functions of the vector class:
41			\code
42			vector3 v1,v2,v3;
43			v1 = VX;
44			v2 = VY;
45			v3 = cross(v1,v2);
46			v3 *= 2.5;
47			v3.normalize();
48			\endcode
49			*/
50
51			/*! This (slow) method allows to access the elements of the
52			vector as if it were an array of doubles. If the index is > 2,
53			then a warning is printed, and the program is terminated via
54			exit(-1). Otherwise, if i is 0, 1 or 2, then a reference to x,
55			y or z is returned, respectively.
56
57			\warning This method is primarily designed to facilitate the
58			integration ('Open Babelization') of code that uses arrays of
59			doubles rather than the vector class. Due to the error checks
60			the method is of course very slow and should therefore be
61			avoided in production code.
62			*/
63			double& vector3::operator[] ( unsigned int i)
64			{
65			if (i > 2)
66			{
67			cerr << "ERROR in OpenBabel::vector3::operator[]" << endl
68			<< "The method has been called with an illegal index i=" << i << "." << endl
69			<< "Please contact the author of the offending program immediately." << endl;
70			exit(-1);
71			}
72			if (i == 0)
73			return _vx;
74			if (i == 1)
75			return _vy;
76			return _vz;
77			}
78
79			/! replaces this with a random unit vector, which is (supposed
80			to be) uniformly distributed over the unit sphere. Uses the
81			random number generator obRand, or uses the system number
82			generator with a time seed if obRand == NULL.
83
84			@param obRandP random number generator to use, or 0L, if the
85			system random number generator (with time seed) should be used
86			*/
87			void vector3::randomUnitVector(OBRandom *obRandP)
88			{
89			OBRandom *ptr;
90			if (!obRandP)
91			{
92			ptr = new OBRandom(true);
93			ptr->TimeSeed();
94			}
95			else
96			ptr = obRandP;
97
98			// obtain a random vector with 0.001 <= length^2 <= 1.0, normalize
99			// the vector to obtain a random vector of length 1.0.
100			double l;
101			do
102			{
103			this->Set(ptr->NextFloat()-0.5, ptr->NextFloat()-0.5, ptr->NextFloat()-0.5);
104			l = length_2();
105			}
106			while ( (l > 1.0) \|\| (l < 1e-4) );
107			this->normalize();
108
109			if (!obRandP)
110			delete ptr;
111			}
112
113			OBAPI ostream& operator<< ( ostream& co, const vector3& v )
114			{
115			co << "< " << v._vx << ", " << v._vy << ", " << v._vz << " >" ;
116			return co ;
117			}
118
119			OBAPI int operator== ( const vector3& v1, const vector3& v2 )
120			{
121			if ( ( v1._vx == v2._vx ) &&
122			( v1._vy == v2._vy ) &&
123			( v1._vz == v2._vz ) )
124			return ( true ) ;
125			else
126			return ( false ) ;
127			}
128
129			OBAPI int operator!= ( const vector3& v1, const vector3& v2 )
130			{
131			if ( ( v1._vx != v2._vx ) \|\|
132			( v1._vy != v2._vy ) \|\|
133			( v1._vz != v2._vz ) )
134			return ( true ) ;
135			else
136			return ( false ) ;
137			}
138
139			/*! This method checks if the current vector has length() ==
140			0.0. If so, this remains unchanged. Otherwise, this is
141			scaled by 1.0/length().
142
143			\warning If length() is very close to zero, but not == 0.0,
144			this method may behave in unexpected ways and return almost
145	gezelter	3057	random results; details may depend on your particular floating
146	tim	2440	point implementation. The use of this method is therefore
147			highly discouraged, unless you are certain that length() is in
148			a reasonable range, away from 0.0 (Stefan Kebekus)
149
150			\deprecated This method will probably replaced by a safer
151			algorithm in the future.
152
153			\todo Replace this method with a more fool-proof version.
154
155			@returns a reference to *this
156			*/
157			vector3& vector3 :: normalize ()
158			{
159			double l = length ();
160
161			if (IsNearZero(l))
162			return(*this);
163
164			_vx = _vx / l ;
165			_vy = _vy / l ;
166			_vz = _vz / l ;
167
168			return(*this);
169			}
170
171			OBAPI double dot ( const vector3& v1, const vector3& v2 )
172			{
173			return v1._vxv2._vx + v1._vyv2._vy + v1._vz*v2._vz ;
174			}
175
176			OBAPI vector3 cross ( const vector3& v1, const vector3& v2 )
177			{
178			vector3 vv ;
179
180			vv._vx = v1._vyv2._vz - v1._vzv2._vy ;
181			vv._vy = - v1._vxv2._vz + v1._vzv2._vx ;
182			vv._vz = v1._vxv2._vy - v1._vyv2._vx ;
183
184			return ( vv ) ;
185			}
186
187
188			/*! This method calculates the angle between two vectors
189
190			\warning If length() of any of the two vectors is == 0.0,
191			this method will divide by zero. If the product of the
192			length() of the two vectors is very close to 0.0, but not ==
193			0.0, this method may behave in unexpected ways and return
194			almost random results; details may depend on your particular
195	gezelter	3057	floating point implementation. The use of this method is
196	tim	2440	therefore highly discouraged, unless you are certain that the
197			length()es are in a reasonable range, away from 0.0 (Stefan
198			Kebekus)
199
200			\deprecated This method will probably replaced by a safer
201			algorithm in the future.
202
203			\todo Replace this method with a more fool-proof version.
204
205			@returns the angle in degrees (0-360)
206			*/
207			OBAPI double vectorAngle ( const vector3& v1, const vector3& v2 )
208			{
209			double mag;
210			double dp;
211
212			mag = v1.length() * v2.length();
213			dp = dot(v1,v2)/mag;
214
215			if (dp < -0.999999)
216			dp = -0.9999999;
217
218			if (dp > 0.9999999)
219			dp = 0.9999999;
220
221			if (dp > 1.0)
222			dp = 1.0;
223
224			return((RAD_TO_DEG * acos(dp)));
225			}
226
227	gezelter	3057	/*! This function calculates the torsion angle of three vectors, represented
228			by four points A--B--C--D, i.e. B and C are vertexes, but none of A--B,
229			B--C, and C--D are colinear. A "torsion angle" is the amount of "twist"
230			or torsion needed around the B--C axis to bring A--B into the same plane
231			as B--C--D. The torsion is measured by "looking down" the vector B--C so
232			that B is superimposed on C, then noting how far you'd have to rotate
233			A--B to superimpose A over D. Angles are + in the anticlockwise
234			direction. The operation is symmetrical in that if you reverse the image
235			(look from C to B and rotate D over A), you get the same answer.
236			*/
237
238	tim	2440	OBAPI double CalcTorsionAngle(const vector3 &a, const vector3 &b,
239			const vector3 &c, const vector3 &d)
240			{
241			double torsion;
242			vector3 b1,b2,b3,c1,c2,c3;
243
244			b1 = a - b;
245			b2 = b - c;
246			b3 = c - d;
247
248			c1 = cross(b1,b2);
249			c2 = cross(b2,b3);
250			c3 = cross(c1,c2);
251
252			if (c1.length() * c2.length() < 0.001)
253			torsion = 0.0;
254			else
255			{
256			torsion = vectorAngle(c1,c2);
257			if (dot(b2,c3) > 0.0)
258			torsion *= -1.0;
259			}
260
261			return(torsion);
262			}
263
264			/! This method checks if the current vector this is zero
265			(i.e. if all entries == 0.0). If so, a warning message is
266			printed, and the whole program is aborted with exit(0).
267			Otherwise, a vector of length one is generated, which is
268			orthogonal to *this, and stored in v. The resulting vector is
269			not random.
270
271			\warning If the entries of the *this (in particular the
272			z-component) are very close to zero, but not == 0.0, this
273			method may behave in unexpected ways and return almost random
274			results; details may depend on your particular floating point
275			implementation. The use of this method is therefore highly
276			discouraged, unless you are certain that all components of
277			*this are in a reasonable range, away from 0.0 (Stefan
278			Kebekus)
279
280			\deprecated This method will probably replaced by a safer
281			algorithm in the future.
282
283			\todo Replace this method with a more fool-proof version that
284			does not call exit()
285
286			@param res a reference to a vector where the result will be
287			stored
288			*/
289			void vector3::createOrthoVector(vector3 &res) const
290			{
291			vector3 cO;
292
293			if ( ( IsNearZero(this->x())) && (IsNearZero(this->y())) )
294			{
295			if ( IsNearZero(this->z()) )
296			{
297			cerr << "makeorthovec zero vector" << endl;
298			exit(0);
299			}
300			cO.SetX(1.0);
301			}
302			else
303			{
304			cO.SetZ(1.0);
305			}
306			res= cross(cO,*this);
307			res.normalize();
308			}
309
310			const vector3 VZero ( 0.0, 0.0, 0.0 ) ;
311			const vector3 VX ( 1.0, 0.0, 0.0 ) ;
312			const vector3 VY ( 0.0, 1.0, 0.0 ) ;
313			const vector3 VZ ( 0.0, 0.0, 1.0 ) ;
314
315			/* Calculate the distance of point a to the plane determined by b,c,d */
316			double Point2Plane(vector3 a, vector3 b, vector3 c, vector3 d)
317			{
318			double angle =0;
319			double dist_ab =0;
320			vector3 v_ba = a-b;
321			vector3 v_normal = cross(c-b, d-b).normalize();
322			angle = vectorAngle(v_normal, v_ba);
323			dist_ab = v_ba.length();
324			return fabs(dist_ab * cos(DEG_TO_RAD * angle));
325			}
326
327			} // namespace OpenBabel
328
329			//! \file vector3.cpp
330			//! \brief Handle 3D coordinates.