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:	1081
Committed:	Thu Oct 19 20:49:05 2006 UTC (18 years, 6 months ago) by gezelter
File size:	9733 byte(s)
Log Message:	updated OpenBabel to version 2.0.2
#	Content
1	/**********************************************************************
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	The vector3 class was designed to simplify operations with floating
34	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	random results; details may depend on your particular floating
146	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	floating point implementation. The use of this method is
196	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	/*! 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	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.