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 doubleing
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 doubleing
  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
  doubleing 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)));
}

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:	2440
Committed:	Wed Nov 16 19:42:11 2005 UTC (18 years, 7 months ago) by tim
File size:	9072 byte(s)
Log Message:	adding openbabel
#	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 doubleing
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 doubleing
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	doubleing 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	OBAPI double CalcTorsionAngle(const vector3 &a, const vector3 &b,
228	const vector3 &c, const vector3 &d)
229	{
230	double torsion;
231	vector3 b1,b2,b3,c1,c2,c3;
232
233	b1 = a - b;
234	b2 = b - c;
235	b3 = c - d;
236
237	c1 = cross(b1,b2);
238	c2 = cross(b2,b3);
239	c3 = cross(c1,c2);
240
241	if (c1.length() * c2.length() < 0.001)
242	torsion = 0.0;
243	else
244	{
245	torsion = vectorAngle(c1,c2);
246	if (dot(b2,c3) > 0.0)
247	torsion *= -1.0;
248	}
249
250	return(torsion);
251	}
252
253	/! This method checks if the current vector this is zero
254	(i.e. if all entries == 0.0). If so, a warning message is
255	printed, and the whole program is aborted with exit(0).
256	Otherwise, a vector of length one is generated, which is
257	orthogonal to *this, and stored in v. The resulting vector is
258	not random.
259
260	\warning If the entries of the *this (in particular the
261	z-component) are very close to zero, but not == 0.0, this
262	method may behave in unexpected ways and return almost random
263	results; details may depend on your particular floating point
264	implementation. The use of this method is therefore highly
265	discouraged, unless you are certain that all components of
266	*this are in a reasonable range, away from 0.0 (Stefan
267	Kebekus)
268
269	\deprecated This method will probably replaced by a safer
270	algorithm in the future.
271
272	\todo Replace this method with a more fool-proof version that
273	does not call exit()
274
275	@param res a reference to a vector where the result will be
276	stored
277	*/
278	void vector3::createOrthoVector(vector3 &res) const
279	{
280	vector3 cO;
281
282	if ( ( IsNearZero(this->x())) && (IsNearZero(this->y())) )
283	{
284	if ( IsNearZero(this->z()) )
285	{
286	cerr << "makeorthovec zero vector" << endl;
287	exit(0);
288	}
289	cO.SetX(1.0);
290	}
291	else
292	{
293	cO.SetZ(1.0);
294	}
295	res= cross(cO,*this);
296	res.normalize();
297	}
298
299	const vector3 VZero ( 0.0, 0.0, 0.0 ) ;
300	const vector3 VX ( 1.0, 0.0, 0.0 ) ;
301	const vector3 VY ( 0.0, 1.0, 0.0 ) ;
302	const vector3 VZ ( 0.0, 0.0, 1.0 ) ;
303
304	/* Calculate the distance of point a to the plane determined by b,c,d */
305	double Point2Plane(vector3 a, vector3 b, vector3 c, vector3 d)
306	{
307	double angle =0;
308	double dist_ab =0;
309	vector3 v_ba = a-b;
310	vector3 v_normal = cross(c-b, d-b).normalize();
311	angle = vectorAngle(v_normal, v_ba);
312	dist_ab = v_ba.length();
313	return fabs(dist_ab * cos(DEG_TO_RAD * angle));
314	}
315
316	} // namespace OpenBabel
317
318	//! \file vector3.cpp
319	//! \brief Handle 3D coordinates.