Branch data Line data Source code
1 : : //-------------------------------------------------------------------------
2 : : // COPYRIGHT 1998 CATERPILLAR INC. ALL RIGHTS RESERVED
3 : : //
4 : : // Filename : AnalyticGeometryTool.hpp
5 : : //
6 : : // Purpose : This class performs calculations on analytic geometry
7 : : // (points, lines, arcs, planes, polygons). Capabilities
8 : : // include vector and point math, matrix operations,
9 : : // measurements, intersections and comparison/containment checks.
10 : : //
11 : : // Special Notes : As most of these functions were taken from Caterpillar's
12 : : // Cmath library with minimal modification for this class,
13 : : // the vectors, planes, etc. are represented in non-CUBIT
14 : : // format. These could be converted with some effort.
15 : : //
16 : : // Other Notes:
17 : : // ------------
18 : : // 1. Points and vectors are represented as double[3].
19 : : //
20 : : // 2. Lines (unless denoted as 'segments') are defined as
21 : : // a point and a vector (point-direction form). The parametric
22 : : // equations of a line are:
23 : : //
24 : : // x = xo + t*a
25 : : // y = yo + t*a
26 : : // z = zo + t*a
27 : : //
28 : : // 3. Planes are also defined as a point and a vector
29 : : // (point-normal form). The equation of a plane is:
30 : : //
31 : : // a*(x-xo) + b*(y-yo) + c*(z-zo) = 0
32 : : //
33 : : // For example, the equation of a plane passing through the
34 : : // point (3,-1,7) and perpendicular to the vector n = (4,2,-5)
35 : : // is: 4(x-3) + 2(y+1) -5(z-7) = 0
36 : : //
37 : : // This can be reduced to:
38 : : //
39 : : // A*x + B*y + C*z + D = 0
40 : : //
41 : : // For the example, the equation is: 4x + 2y -5z + 25 = 0
42 : : //
43 : : // 4. Most functions which can accept a 3x3 matrix have
44 : : // an analogous 4x4 function. This is for convenience
45 : : // only (the extra rows/columns are ignored within the function).
46 : : //
47 : : //
48 : : // Creator : Steve Storm
49 : : //
50 : : // Creation Date : 10/16/98
51 : : //-------------------------------------------------------------------------
52 : :
53 : : #ifndef ANALYTIC_GEOMETRY_TOOL_HPP
54 : : #define ANALYTIC_GEOMETRY_TOOL_HPP
55 : :
56 : : #include <math.h>
57 : : #include "CubitDefines.h"
58 : : #include "CGMGeomConfigure.h"
59 : :
60 : : #include <iostream>
61 : :
62 : : class CubitBox;
63 : : class CubitPlane;
64 : : class CubitVector;
65 : : template <class X> class DLIList;
66 : :
67 : : // Multiplication constants
68 : : const double AGT_E = 2.7182818284590452354; // e
69 : : const double AGT_LOG_2E = 1.4426950408889634074; // log2(e)
70 : : const double AGT_LOG_10E = 0.43429448190325182765; // log10(e)
71 : : const double AGT_LN_2 = 0.69314718055994530942; // ln(2)
72 : : const double AGT_LN_10 = 2.30258509299404568402; // ln(10)
73 : : const double AGT_PI = 3.14159265358979323846; // PI
74 : : const double AGT_PI_DIV_2 = 1.57079632679489661923; // PI/2
75 : : const double AGT_PI_DIV_4 = 0.78539816339744830962; // PI/4
76 : : const double AGT_1_DIV_PI = 0.31830988618379067154; // 1/PI
77 : : const double AGT_2_DIV_PI = 0.63661977236758134308; // 2/PI
78 : : const double AGT_2_DIV_SQRT_PI = 1.12837916709551257390; // 2/(sqrt(PI))
79 : : const double AGT_SQRT_2 = 1.41421356237309504880; // sqrt(2)
80 : : const double AGT_SQRT_1_DIV_2 = 0.70710678118654752440; // sqrt(1/2)
81 : : const double AGT_PI_DIV_180 = 0.017453292519943295; // PI/180
82 : : const double AGT_180_DIV_PI = 57.295779513082323; // 180/PI
83 : :
84 : : const double DEGtoRAD = AGT_PI_DIV_180;
85 : : const double RADtoDEG = AGT_180_DIV_PI;
86 : :
87 : : typedef struct {
88 : : double Vec1[3]; /* First vector that defines plane of arc */
89 : : double Vec2[3]; /* Second vector that defines plane of arc */
90 : : double Origin[3]; /* Center of arc (on plane) */
91 : : double StartAngle; /* Starting angle (in radians) of arc */
92 : : double EndAngle; /* End angle (in radians) of arc */
93 : : double Radius; /* Radius of arc */
94 : : } AGT_3D_Arc;
95 : :
96 : : // Return values */
97 : : enum AgtEquality {
98 : : AGT_EQUAL = 0, // define for comparisons
99 : : AGT_LESSTHAN = 1, // define for comparisons
100 : : AGT_GREATERTHAN = 2 // define for comparisons
101 : : };
102 : :
103 : : enum AgtSide {
104 : : AGT_ON_PLANE = 0, // define for which side of plane
105 : : AGT_POS_SIDE = 1, // define for which side of plane
106 : : AGT_NEG_SIDE = 2, // define for which side of plane
107 : : AGT_INT_PLANE = 3, // define for intersects plane
108 : : AGT_CROSS = 4 // define for line crossing plane
109 : : };
110 : :
111 : : enum AgtOrientation {
112 : : AGT_OPP_DIR = 0, // define for vector direction comparision
113 : : AGT_SAME_DIR = 1 // define for vector direction comparision
114 : : };
115 : :
116 : : enum AgtDistanceMethod {
117 : : AGT_FRACTION = 0, // define for distance along a curve
118 : : AGT_ABSOLUTE = 1 // define for distance along a curve
119 : : };
120 : :
121 : : // Macros
122 : : #define number_sign(number) (number<0 ? -1 : 1)
123 : : #define copy_vec(vec1,vec2) copy_pnt(vec1,vec2)
124 : : #define set_vec(pnt,x,y,z) set_pnt(pnt,x,y,z)
125 : :
126 : : ///////////////////////////////////////////////////////////////////////////////
127 : : // MAGIC SOFTWARE //
128 : : // See note later in file regarding MAGIC softare. //
129 : : ///////////////////////////////////////////////////////////////////////////////
130 : : // Minimal 2D bounding box - parameters
131 : : const int maxPartition = 32;
132 : : const int invInterp = 32;
133 : : const double angleMin = -0.25*AGT_PI;
134 : : const double angleMax = +0.25*AGT_PI;
135 : : const double angleRange = 0.5*AGT_PI;
136 : :
137 : : typedef struct
138 : : {
139 : : double x, y;
140 : : }
141 : : Point2;
142 : :
143 : : typedef struct
144 : : {
145 : : double x, y, z;
146 : : }
147 : : Point3;
148 : :
149 : : typedef struct
150 : : {
151 : : Point2 center;
152 : : Point2 axis[2];
153 : : double extent[2];
154 : : }
155 : : OBBox2;
156 : :
157 : : typedef struct
158 : : {
159 : : Point3 center;
160 : : Point3 axis[3];
161 : : double extent[3];
162 : : }
163 : : OBBox3;
164 : :
165 : : // threshold on determinant to conclude lines/rays/segments are parallel
166 : : const double par_tolerance = 1e-06;
167 : : #define DIST(x) (Sqrt(Abs(x)))
168 : :
169 : : // Line in 2D
170 : : typedef struct
171 : : {
172 : : // line is Dot(N,X) = c, N not necessarily unit length
173 : : Point2 N;
174 : : double c;
175 : : }
176 : : AgtLine;
177 : :
178 : : // Triangle in 2D
179 : : typedef struct
180 : : {
181 : : Point2 v[3];
182 : : }
183 : : Triangle;
184 : :
185 : : // Linked list of triangle
186 : : typedef struct AgtTriList
187 : : {
188 : : Triangle* tri;
189 : : AgtTriList* next;
190 : : }
191 : : AgtTriList;
192 : :
193 : : //---------------------------------------------------------------------------
194 : : // Lines in 3D
195 : : //---------------------------------------------------------------------------
196 : : typedef struct
197 : : {
198 : : // Line is L(t) = b+t*m for any real-valued t
199 : : // Ray has constraint t >= 0, b is the origin of the ray
200 : : // Line segment has constraint 0 <= t <= 1, b and b+m are end points
201 : :
202 : : Point3 b, m;
203 : : }
204 : : Line3;
205 : : //---------------------------------------------------------------------------
206 : :
207 : : //---------------------------------------------------------------------------
208 : : // Circles in 3D
209 : : //---------------------------------------------------------------------------
210 : : typedef struct
211 : : {
212 : : // Plane containing circle is Dot(N,X-C) = 0 where X is any point in the
213 : : // plane. Vectors U, V, and N form an orthonormal right-handed set
214 : : // (matrix [U V N] is orthonormal and has determinant 1). Circle within
215 : : // the plane is parameterized by X = C + R*(cos(A)*U + sin(A)*V) where
216 : : // A is an angle in [0,2*pi).
217 : :
218 : : Point3 U, V, N; // coordinate system of plane containing circle
219 : : Point3 C; // center
220 : : double R; // radius > 0
221 : : }
222 : : Circle3;
223 : : //---------------------------------------------------------------------------
224 : :
225 : : //---------------------------------------------------------------------------
226 : : // Triangles in 3D
227 : : //---------------------------------------------------------------------------
228 : : typedef struct
229 : : {
230 : : // Triangle points are tri(s,t) = b+s*e0+t*e1 where 0 <= s <= 1,
231 : : // 0 <= t <= 1, and 0 <= s+t <= 1.
232 : :
233 : : // If the vertices of the triangle are v0, v1, and v2, then b = v0,
234 : : // e0 = v1 - v0, and e1 = v2 - v0.
235 : :
236 : : Point3 b, e0, e1;
237 : : }
238 : : Triangle3;
239 : : //---------------------------------------------------------------------------
240 : :
241 : : //---------------------------------------------------------------------------
242 : : // Parallelograms in 3D
243 : : //---------------------------------------------------------------------------
244 : : typedef struct
245 : : {
246 : : // Rectoid points are rect(s,t) = b+s*e0+t*e1 where 0 <= s <= 1
247 : : // and 0 <= t <= 1. Could have called this a parallelogram, but
248 : : // I do not like typing long words...
249 : :
250 : : Point3 b, e0, e1;
251 : : }
252 : : Rectoid3;
253 : : //---------------------------------------------------------------------------
254 : :
255 : : //---------------------------------------------------------------------------
256 : : // Planes in 3D
257 : : //---------------------------------------------------------------------------
258 : : typedef struct
259 : : {
260 : : // plane is Dot(N,X) = d
261 : : Point3 N;
262 : : double d;
263 : : }
264 : : Plane3;
265 : : //---------------------------------------------------------------------------
266 : :
267 : : // END OF MAGIC SOFTWARE
268 : : ///////////////////////////////////////////////////////////////////////////////
269 : :
270 : : //! Performs calculations on analytic geometry.
271 : : class CUBIT_GEOM_EXPORT AnalyticGeometryTool
272 : : {
273 : : public:
274 : :
275 : : ~AnalyticGeometryTool();
276 : : static AnalyticGeometryTool* instance();
277 : :
278 : 322 : static void delete_instance()
279 : : {
280 [ + + ]: 322 : if( instance_ )
281 [ + - ]: 271 : delete instance_;
282 : 322 : instance_ = NULL;
283 : 322 : };
284 : :
285 : : //*********************************************************
286 : : // Double numbers
287 : : //*********************************************************
288 : :
289 : : /*! Check to see if double numbers are equal within epsilon.
290 : : Values are equal if fabs(val_1 - val_2) < epsilon.
291 : : Epsilon is determined by next two functions. */
292 : : CubitBoolean dbl_eq( double val_1, double val_2 );
293 : :
294 : : //! Get epsilon used for double number comparisons.
295 : : //! Default value is 1e-6.
296 : : double get_epsilon();
297 : : //! Set epsilon used for double number comparisons.
298 : : //! Default value is 1e-6.
299 : : double set_epsilon( double new_epsilon );
300 : :
301 : : /*! Round value to either -1, 0, or 1 if within epsilon to
302 : : to one of these. Epsilon determined by get_epsilon/set_epsilon
303 : : functions. */
304 : : void round_near_val( double& val );
305 : :
306 : : //**************************************************************************
307 : : // Matrices & Transforms
308 : : //**************************************************************************
309 : : // Note: For these functions the matrix format is defined as follows:
310 : : //
311 : : // Consider the transformation from [x,y,z] to [x',y',z']:
312 : : // _ _
313 : : // | x1 y1 z1 0 |
314 : : // [x',y',z',1] = [x,y,z,1]| x2 y2 z2 0 |
315 : : // | x3 y3 z3 0 |
316 : : // | ox oy oz 1 |
317 : : // - -
318 : : /*! This functions applies the transformation matrix to the specified
319 : : point and returns the new point coordinates through the call list.
320 : : Point in and point out can be the same memory address or different. */
321 : : void transform_pnt( double m[4][4], double pin[3], double pout[3] );
322 : :
323 : : /*! This routine applies a transformation matrix to a specified vector
324 : : and returns the new vector orientation. Vector in and vector out can
325 : : be the same memory address or different. */
326 : : void transform_vec( double m3[3][3], double vin[3], double vout[3] );
327 : :
328 : : /*! This routine applies a transformation matrix to a specified vector
329 : : and returns the new vector orientation. Vector in and vector out can
330 : : be the same memory address or different. */
331 : : void transform_vec( double m4[4][4], double vin[3], double vout[3] );
332 : :
333 : : /*! Transforms a line using the given transformation mtx. The mtx is typically
334 : : built using add_origin_to_rotation_mtx. */
335 : : void transform_line( double rot_mtx[4][4], double origin[3], double axis[3] );
336 : : /*! Transforms a line using the given transformation mtx. The mtx is typically
337 : : built using add_origin_to_rotation_mtx. */
338 : : void transform_line( double rot_mtx[4][4], CubitVector &origin, CubitVector &axis );
339 : :
340 : :
341 : : //@{
342 : : /*! This routine simply copies one matrix to another. If a NULL is passed in
343 : : for the from matrix, then the identity matrix is copied into the out matrix.
344 : : Last function allows you to specify 1st 3 doubles of row 4 (origin) of the
345 : : 4x4 matrix.*/
346 : : void copy_mtx( double from[3][3],double to[3][3] );
347 : : void copy_mtx( double from[4][4], double to[4][4] );
348 : : void copy_mtx( double from[4][4], double to[3][3] );
349 : : void copy_mtx(double from[3][3], double to[4][4], double *origin = NULL );
350 : : //@}
351 : :
352 : : //@{
353 : : //! This routine determines the tranformation matrix given the theta and
354 : : //! the vector to cross through.
355 : : void create_rotation_mtx( double theta, double v[3], double mtx3x3[3][3] );
356 : : void create_rotation_mtx( double theta, double v[3], double mtx4x4[4][4] );
357 : : //@}
358 : :
359 : : /*! Adds origin to rotation matrix, so you can rotate points about a line.
360 : : Line is defined as original vector used in create_rotation_mtx and
361 : : the origin. */
362 : : void add_origin_to_rotation_mtx( double rot_mtx[4][4], double origin[3] );
363 : :
364 : : //@{
365 : : //! \brief Simply sets the given matrix to the identity matrix.
366 : : void identity_mtx( double mtx3x3[3][3] );
367 : : void identity_mtx( double mtx4x4[4][4] );
368 : : //@}
369 : :
370 : : //@{
371 : : //! Gets rotation angles to rotate one system to another - returned rotation
372 : : //! angles are about global system.
373 : : //! To use the result from this function, align the unoriented object's origin
374 : : //! with the global origin, then apply the rotation angles returned from this
375 : : //! routine to the object in the order of ax,ay,az about the global origin.
376 : : //! The object will be oriented such that its xyz axes will point in the same
377 : : //! direction as the object whose transformation matrix was given. The object
378 : : //! can then be translated.
379 : : void mtx_to_angs( double mtx3x3[3][3], double &ax, double &ay, double &az );
380 : : void mtx_to_angs( double mtx4x4[4][4], double &ax, double &ay, double &az );
381 : : //@}
382 : :
383 : :
384 : : //@{
385 : : // This routine rotates a 3x3 system given a transformation matrix. The
386 : : // matrix can be rotated in place by sending same variable in & out, or a
387 : : // new matrix can be created. In the 4x4 case, the translation portion
388 : : // of the matrix is unchanged.
389 : : void rotate_system( double mtx[3][3], double sys_in[3][3],
390 : : double sys_out[3][3] );
391 : : void rotate_system( double mtx[4][4], double sys_in[4][4],
392 : : double sys_out[4][4] );
393 : : //@}
394 : :
395 : : //! Find determinant of matrix.
396 : : double det_mtx( double m[3][3] );
397 : :
398 : : //@{
399 : : // Multiply matrices together. If any input is NULL, the identity
400 : : // matrix is used in its place.
401 : : void mult_mtx( double a[3][3],double b[3][3], double d[3][3] );
402 : : void mult_mtx( double a[4][4], double b[4][4], double d[4][4] );
403 : : //@}
404 : :
405 : : //!This routine inverts a 3x3 matrix using the adjoint method. If NULL
406 : : //! is sent in for first matrix, the second matrix is set to the identity
407 : : //! matrix. If the same memory address is passed in for both incoming and
408 : : //! outgoing matrices, the matrix is inverted in place. Returns CUBIT_FAILURE
409 : : //! if no inverse exists.
410 : : CubitStatus inv_mtx_adj( double mtx[3][3], double inv_mtx[3][3] );
411 : :
412 : : //! This routine inverts a 4x4 transformation matrix. If NULL is sent in
413 : : //! for first matrix, the second matrix is set to the identity matrix. If
414 : : //! the same memory address is passed in for both incoming and outgoing matrices,
415 : : //! the matrix is inverted in place. Uses LU decomposition.
416 : : CubitStatus inv_trans_mtx( double transf[4][4], double inv_transf[4][4] );
417 : :
418 : : //@{
419 : : //! Creates a transformation matrix given three vectors (and origin
420 : : //! for 4x4 case).
421 : : void vecs_to_mtx( double xvec[3], double yvec[3], double zvec[3],
422 : : double matrix[3][3] );
423 : : void vecs_to_mtx( double xvec[3], double yvec[3], double zvec[3],
424 : : double origin[3], double matrix[4][4] );
425 : : //@}
426 : :
427 : : //@{
428 : : //! Prints matrix values, for debugging.
429 : : void print_mtx( double mtx[3][3] );
430 : : void print_mtx( double mtx[4][4] );
431 : : //@}
432 : :
433 : : //**************************************************************************
434 : : // 3D Points
435 : : //**************************************************************************
436 : : //! Copy one double[3] point to another double[3] point (uses memcpy).
437 : : //! If first point in NULL then second point is set to 0,0,0.
438 : : void copy_pnt( double pnt_in[3], double pnt_out[3] );
439 : :
440 : : //@{
441 : : //! For going back and forth from CubitVectors
442 : : void copy_pnt( double pnt_in[3], CubitVector &cubit_vec );
443 : : void copy_pnt( CubitVector &cubit_vec, double pnt_out[3] );
444 : : //@}
445 : :
446 : : //! Sets the value of pnt to x,y,z (inline function).
447 : : void set_pnt( double pnt[3], double x, double y, double z );
448 : :
449 : : //! Compares two points to determine if they are equivalent. The
450 : : //! equivalence is based on epsilon (get_epsilon, set_epsilon).
451 : : //! If the distance between them is less than epsilon, the points
452 : : //! are considered equal.
453 : : CubitBoolean pnt_eq( double pnt1[3],double pnt2[3] );
454 : :
455 : : //! Function to mirror a point about a plane. The mirror point
456 : : //! is the same distance from the plane as the original, but on
457 : : //! the opposite side of the plane. Function returns CUBIT_FAILURE
458 : : //! if the point is on the plane - in which case the point is
459 : : //! just copied.
460 : : CubitStatus mirror_pnt( double pnt[3], double pln_orig[3],
461 : : double pln_norm[3], double m_pnt[3]);
462 : :
463 : : //! Given start pnt, vec dir and length find next point in space. The
464 : : //! vector direction is first unitized then the following formula is used:
465 : : //! new_point[] = start_point[] + (length * unit_vector[])
466 : : //! new_pnt can be the same as input point (overwrites old location).
467 : : CubitStatus next_pnt( double str_pnt[3], double vec_dir[3], double len,
468 : : double new_pnt[3]);
469 : :
470 : : //***************************************************************************
471 : : // 3D Vectors
472 : : //***************************************************************************
473 : : //! Finds unit vector of input vector (vout can equal vin - converts in place).
474 : : CubitStatus unit_vec( double vin[3], double vout[3] );
475 : :
476 : : //! Dots two vectors. Property of dot product is:
477 : : //! angle between vectors acute if dot product > 0
478 : : //! angle between vectors obtuse if dot product < 0
479 : : //! angle between vectors 90 deg if dot product = 0
480 : : double dot_vec( double uval[3], double vval[3] );
481 : :
482 : : //! Finds cross product of two vectors:
483 : : //! cross[0] = u[1] * v[2] - u[2] * v[1];
484 : : //! cross[1] = u[2] * v[0] - u[0] * v[2];
485 : : //! cross[2] = u[0] * v[1] - u[1] * v[0];
486 : : void cross_vec( double uval[3], double vval[3], double cross[3] );
487 : :
488 : : //! Finds unit vector that is cross product of two vectors.
489 : : void cross_vec_unit( double uval[3], double vval[3], double cross[3] );
490 : :
491 : : //! Finds 2 arbitrary orthoganal vectors to the first.
492 : : void orth_vecs( double uvect[3], double vvect[3], double wvect[3] );
493 : :
494 : : //! Finds the magnitude of a vector:
495 : : //! magnitude = sqrt (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
496 : : double mag_vec( double vec[3] );
497 : :
498 : : //! Finds a vector from str_pnt to stp_pnt.
499 : : //! Returns CUBIT_FAILURE if points are coincident.
500 : : CubitStatus get_vec ( double str_pnt[3], double stp_pnt[3],
501 : : double vector_out[3] );
502 : :
503 : : //! Finds a unit vector from str_pnt to stp_pnt.
504 : : //! Returns CUBIT_FAILURE if points are coincident.
505 : : CubitStatus get_vec_unit( double str_pnt[3], double stp_pnt[3],
506 : : double uv_out[3] );
507 : :
508 : : //! Multiples a vector by a constant (vec_out can equal vec).
509 : : void mult_vecxconst( double constant, double vec[3], double vec_out[3] );
510 : :
511 : : //! Multiples -1.0 by a vector (vout can equal vin).
512 : : void reverse_vec( double vin[3],double vout[3] );
513 : :
514 : : //! Finds angle in radians between two vectors. Angle will always be
515 : : //! between 0.0 and PI. For accuracy, the routine uses the cosine for large
516 : : //! angles and the sine for small angles.
517 : : double angle_vec_vec( double v1[3],double v2[3] );
518 : :
519 : : //***************************************************************************
520 : : // Distances
521 : : //***************************************************************************
522 : : //! Determines distance between two points in space.
523 : : double dist_pnt_pnt( double pnt1[3], double pnt2[3] );
524 : :
525 : : //! This routine determines the shortest distance between two planes. This
526 : : //! is the perpendicular distance between the two planes. Note that if the
527 : : //! two planes are not parallel, the returned distance is zero. The routine
528 : : //! can also be used to determine which side of the first plane the second
529 : : //! plane is on, and the orientation of the two planes to each other.
530 : : double dist_pln_pln( double pln_1_orig[3], double pln_1_norm[3],
531 : : double pln_2_orig[3], double pln_2_norm[3],
532 : : AgtSide *side = NULL, AgtOrientation *orien = NULL,
533 : : unsigned short *status = NULL );
534 : :
535 : : //***************************************************************************
536 : : // Intersections
537 : : //***************************************************************************
538 : : //! This routine calculates the intersection point of a line and a plane.
539 : : //! Returns CUBIT_FAILURE if no intersection exists (line is parallel to the
540 : : //! plane).
541 : : CubitStatus int_ln_pln( double ln_orig[3], double ln_vec[3], double pln_orig[3],
542 : : double pln_norm[3], double int_pnt[3] );
543 : : //! This function finds the intersection of two lines. If the lines cross,
544 : : //! it finds the one unique point. If the lines do not cross (but are non-
545 : : //! parallel), it finds the nearest intersection points (a line drawn from
546 : : //! each of these intersections would be perpendicular to both lines).
547 : : //! Returns number of intersections found:
548 : : //! 0 intersections found, lines are parallel
549 : : //! 1 intersection found, lines physically intersect
550 : : //! 2 intersections found, lines do not intersect but nearest
551 : : //! intersections found
552 : : int int_ln_ln( double p1[3], double v1[3], double p2[3], double v2[3],
553 : : double int_pnt1[3], double int_pnt2[3] );
554 : :
555 : : //! This routine finds the perpendicular intersection of a point & line. This
556 : : //! point will lie on the line.
557 : : //! Returns 0 if point is not on the line, otherwise 1.
558 : : int int_pnt_ln( double pnt[3], double ln_orig[3], double ln_vec[3],
559 : : double int_pnt[3] );
560 : :
561 : : //! This routine determines the perpendicular intersection of a point and a
562 : : //! plane. This is the perpendicular intersection point on the plane.
563 : : //! Returns 0 if point is not on the plane, otherwise 1.
564 : : int int_pnt_pln( double pnt[3], double pln_orig[3], double pln_norm[3],
565 : : double pln_int[3] );
566 : :
567 : : //! This routine finds intersection of two planes. This intersection results
568 : : //! in a line. Returns CUBIT_FAILURE if planes are parallel.
569 : : CubitStatus int_pln_pln( double pln_1_orig[3], double pln_1_norm[3],
570 : : double pln_2_orig[3], double pln_2_norm[3],
571 : : double ln_orig[3], double ln_vec[3] );
572 : :
573 : : //**************************************************************************
574 : : // Comparison/Containment Tests
575 : : //**************************************************************************
576 : : //! This routine checks to see if two vectors are parallel. Two vectors
577 : : //! are considered parallel if they are pointing in the same direction or
578 : : //! opposite directions.
579 : : CubitBoolean is_vec_par( double vec_1[3], double vec_2[3] );
580 : :
581 : : //! This routine checks to see if two vectors are perpendicular. Two vectors
582 : : //! are considered perpendicular if they are PI/2 radians to each other.
583 : : CubitBoolean is_vec_perp( double vec_1[3],double vec_2[3]);
584 : :
585 : : //! This routine checks to see if two vectors point in the same direction.
586 : : //! The vector magnitudes do not have to be equal for a successful return.
587 : : CubitBoolean is_vecs_same_dir( double vec_1[3], double vec_2[3] );
588 : :
589 : : //! This routined determines if a specified point is on a given infinite line
590 : : //! defined by a point and a vector.
591 : : CubitBoolean is_pnt_on_ln( double pnt[3], double ln_orig[3], double ln_vec[3] );
592 : :
593 : : //! This routine determines if a specified point is on a given line segment
594 : : //! defined by two endpoints. The point is on the line only if it lies on
595 : : //! the line between or on the two endpoints.
596 : : CubitBoolean is_pnt_on_ln_seg( double pnt1[3], double end1[3], double end2[3] );
597 : :
598 : : //! This routined determines if a specified point is on a given plane
599 : : //! which is defined by an origin and three vectors.
600 : : CubitBoolean is_pnt_on_pln( double pnt[3], double pln_orig[3], double pln_norm[3] );
601 : :
602 : : //! This routine determines if a specified line (defined by a point and
603 : : //! a vector) is in a given plane (defined by the two vectors in the plane
604 : : //! and the normal to that plane).
605 : : CubitBoolean is_ln_on_pln( double ln_orig[3], double ln_vec[3],
606 : : double pln_orig[3], double pln_norm[3] );
607 : :
608 : :
609 : : //! This routine checks to see if two infinite planes are coincident.
610 : : CubitBoolean is_pln_on_pln( double pln_orig1[3], double pln_norm1[3],
611 : : double pln_orig2[3], double pln_norm2[3] );
612 : :
613 : : //**************************************************************************
614 : : // Arcs/Circles
615 : : //**************************************************************************
616 : :
617 : : //@{
618 : : //! Functions to populate an arc structure. The arc is defined with two
619 : : //! orthogonal vectors in a plane, the arc origin, the beginning angle
620 : : //! in radians to rotate from the start vector (towards second vector), the
621 : : //! ending angle in radians to rotate from the start vector, and the radius
622 : : //! of the arc. The arc is parameterized from 0 to 1.
623 : : void setup_arc( double vec_1[3], double vec_2[3], double origin[3],
624 : : double start_angle, double end_angle, // Angles in radians
625 : : double radius, AGT_3D_Arc &arc );
626 : : void setup_arc( CubitVector& vec_1, CubitVector& vec_2,
627 : : CubitVector origin, double start_angle, // Angles in radians
628 : : double end_angle, double radius,
629 : : AGT_3D_Arc &arc );
630 : : //@}
631 : :
632 : : //@{
633 : : //! Given a parameter value from 0 to 1 on the arc, return the xyz location.
634 : : //! Arc is assumed to be parameterized from 0 to 1.
635 : : void get_arc_xyz( AGT_3D_Arc &arc, double param, double pnt[3] );
636 : : void get_arc_xyz( AGT_3D_Arc &arc, double param, CubitVector& pnt );
637 : : //@}
638 : :
639 : : //! Get number of tessellation points on the circle required to
640 : : //! approximate the length of the circle within len_tolerance. Can
641 : : //! be used to find the number of tessellations points to display a
642 : : //! circle - smaller circles will have fewer tessellation points, larger
643 : : //! circles will have more tessellation points with the same tolerance.
644 : : //! Minimum number of tessellations points returned is 8.
645 : : int get_num_circle_tess_pnts( double radius, double len_tolerance = 1e-4 );
646 : :
647 : : //**************************************************************************
648 : : // Miscellaneous
649 : : //**************************************************************************
650 : :
651 : : //! Finds an origin-normal format plane from reduced form
652 : : //! A*x + B*y + C*z + D = 0
653 : : void get_pln_orig_norm( double A, double B, double C, double D,
654 : : double pln_orig[3], double pln_norm[3] = NULL );
655 : : //! Find 8 corners of a box given minimum and maximum points. Box is
656 : : //! defined starting from left-bottom-front clockwise (at front of box),
657 : : //! then to the rear in same order.
658 : : void get_box_corners( double box_min[3],double box_max[3],double c[8][3] );
659 : :
660 : : //@{
661 : : //! Finds the 4 corner points of the input infinite plane that just
662 : : //! intersects the input box (defined by bottom-left and top-right corners,
663 : : //! axis aligned with the cubit coordinate system) - plane should have
664 : : //! minimal area. The resultant plane's corner points can be extended out
665 : : //! by a percentage of diagonal or absolute (making it bigger than the minimal
666 : : //! area plane). extension_type - 0=none, 1=percentage, 2=absolute
667 : : //! If silent is CUBIT_TRUE, no error messages are given. Note this returns
668 : : //! points even if the plane doesn't physically intersect the given box... it
669 : : //! does this by moving the plane to the center of the box, doing the
670 : : //! intersection, then projecting the points back to the original plane.
671 : : CubitStatus min_pln_box_int_corners( const CubitPlane& plane,
672 : : const CubitBox& box,
673 : : int extension_type, double extension,
674 : : CubitVector& p1, CubitVector& p2,
675 : : CubitVector& p3, CubitVector& p4,
676 : : CubitBoolean silent = CUBIT_FALSE );
677 : : CubitStatus min_pln_box_int_corners( CubitVector& vec1,
678 : : CubitVector& vec2,
679 : : CubitVector& vec3,
680 : : CubitVector& box_min, CubitVector& box_max,
681 : : int extension_type, double extension,
682 : : CubitVector& p1, CubitVector& p2,
683 : : CubitVector& p3, CubitVector& p4,
684 : : CubitBoolean silent = CUBIT_FALSE );
685 : : CubitStatus min_pln_box_int_corners( double plane_norm[3], double plane_coeff,
686 : : double box_min[3], double box_max[3],
687 : : int extension_type, double extension,
688 : : double p1[3], double p2[3], // OUT
689 : : double p3[3], double p4[3], // OUT
690 : : CubitBoolean silent = CUBIT_FALSE);
691 : : //@}
692 : :
693 : : //@{
694 : : //! Finds the minimum size bounding box that fits around the points. Box
695 : : //! will not necessarily be oriented in xyz directions.
696 : : CubitStatus get_tight_bounding_box( DLIList<CubitVector*> &point_list,
697 : : CubitVector& p1, CubitVector& p2,
698 : : CubitVector& p3, CubitVector& p4,
699 : : CubitVector& p5, CubitVector& p6,
700 : : CubitVector& p7, CubitVector& p8);
701 : : CubitStatus get_tight_bounding_box( DLIList<CubitVector*> &point_list,
702 : : CubitVector ¢er,
703 : : CubitVector axes[3],
704 : : CubitVector &extension );
705 : : //@}
706 : :
707 : :
708 : : //@{
709 : : //! Finds the start and end of a cylinder that just intersects the input
710 : : //! box (defined by bottom-left and top-right corners, axis aligned with
711 : : //! the cubit coordinate system). The resultant cylinder's start and end
712 : : //! points can be extended out by a percentage of cylinder's length or
713 : : //! absolute (making it longer in both directions).
714 : : //! extension_type - 0=none, 1=percentage, 2=absolute
715 : : //! The num_tess_pnts is the number of line points along the circle the
716 : : //! function projects to the box to calculate the minimum intersection
717 : : //! (more points could give a more accurate result for non-axis aligned
718 : : //! cylinders).
719 : : CubitStatus min_cyl_box_int( double radius, CubitVector& axis, CubitVector& center,
720 : : CubitBox& box,
721 : : int extension_type, double extension,
722 : : CubitVector &start, CubitVector &end,
723 : : int num_tess_pnts = 2048 );
724 : : CubitStatus min_cyl_box_int( double radius, double axis[3], double center[3],
725 : : double box_min[3], double box_max[3],
726 : : int extension_type, double extension,
727 : : double start[3], double end[3],
728 : : int num_tess_pnts = 2048 );
729 : : //@}
730 : :
731 : : //! Finds minimum distance between two triangles in 3D (MAGIC)
732 : : double MinTriangleTriangle (Triangle3& tri0,
733 : : Triangle3& tri1,
734 : : double& s, double& t, double& u, double& v);
735 : :
736 : : //! Finds area of intersection between two triangles (MAGIC)
737 : : double AreaIntersection (Triangle &tri1, Triangle &tri2 );
738 : :
739 : : //! Finds angle between normals of the two triangles (radians)
740 : : double Angle( Triangle3& tri1, Triangle3& tri2 );
741 : :
742 : : //! Finds minimum distance beween a triange and point in 3D (MAGIC)
743 : : double MinPointTriangle (const Point3& p, const Triangle3& tri, double& s, double& t);
744 : :
745 : : //! Finds normal vector of given triangle
746 : : void Normal( Triangle3& tri, double normal[3] );
747 : :
748 : : //! Projects tri2 to the plane of tri1 and finds the overlap area
749 : : double ProjectedOverlap( Triangle3& tri1, Triangle3& tri2, bool draw_overlap = false );
750 : :
751 : : protected:
752 : : AnalyticGeometryTool();
753 : : //- Class Constructor. (Not callable by user code. Class is constructed
754 : : //- by the {instance()} member function.
755 : :
756 : : private:
757 : :
758 : : static AnalyticGeometryTool* instance_;
759 : :
760 : : double agtEpsilon; // default = 1e-6
761 : :
762 : :
763 : : ///////////////////////////////////////////////////////////////////////////////
764 : : // MAGIC SOFTWARE
765 : : //This code was obtained from Dave Eberly, at www.magic-software.com. It
766 : : //has been modified only to be placed into AnalyticGeometryTool. This was
767 : : //done for convenience only. An alternate solution would be to create a
768 : : //separate library for these functions.
769 : : //
770 : : //Steve Storm, [email protected], 3-27-99
771 : : ///////////////////////////////////////////////////////////////////////////////
772 : : //MAGIC is an acronym for My Alternate Graphics and Image Code. The
773 : : //initial code base originated in 1991 as an attempt to answer questions that
774 : : //arise in my favorite news group, comp.graphics.algorithms, and has been
775 : : //growing ever since. Magic is intended to provide free source code for solving
776 : : //problems that commonly arise in computer graphics and image analysis. While
777 : : //the code at this web site is free, additional conditions are:
778 : : //
779 : : // * You may distribute the original source code to others at no charge. You
780 : : // got it for free, so do not charge others for it.
781 : : // * You may modify the original source code and distribute it to others at no
782 : : // charge. The modified code must be documented to indicate that it is not
783 : : // part of the original package. I do not want folks to blame me for bugs
784 : : // introduced by modifications. I do accept blame for bugs that are my
785 : : // doing and will do my best to fix them.
786 : : // * You may use this code for non-commercial purposes. You may also
787 : : // incorporate this code into commercial packages. However, you may not
788 : : // sell any of your source code which contains my original and/or modified
789 : : // source code (see items 1 and 2 in this list of conditions). In such a case,
790 : : // you would need to factor out my code and freely distribute it. Send me
791 : : // email if you use it. I am always interested in knowing where and how
792 : : // my code is being used.
793 : : // * The original code comes with absolutely no warranty and no guarantee is
794 : : // made that the code is bug-free. Caveat emptor.
795 : : //
796 : : // Dave Eberly, [email protected], www.magic-software.com
797 : : ///////////////////////////////////////////////////////////////////////////////
798 : : //FILE: minbox2.h
799 : : double Area (int N, Point2* pt, double angle);
800 : : void MinimalBoxForAngle (int N, Point2* pt, double angle, OBBox2& box);
801 : : OBBox2 MinimalBox2 (int N, Point2* pt);
802 : : // Functions to find minimal area rectangle that surrounds a set of points.
803 : :
804 : : #if 1
805 : : //FILE: minbox3.h
806 : : void MatrixToAngleAxis( double** R, double& angle, double axis[3] );
807 : : void AngleAxisToMatrix( double angle, double axis[3], double R[3][3] );
808 : : double Volume( int N, Point3* pt, double angle[3] );
809 : : void MinimalBoxForAngles( int N, Point3* pt, double angle[3], OBBox3& box );
810 : : void GetInterval( double A[3], double D[3], double& tmin, double& tmax );
811 : : void Combine( double result[3], double A[3], double t, double D[3] );
812 : : double MinimizeOnInterval( int N, Point3* pt, double A[3], double D[3] );
813 : : double MinimizeOnLattice( int N, Point3* pt, double A[3], int layers,
814 : : double thickness );
815 : : void InitialGuess( int N, Point3* pt, double angle[3] );
816 : : OBBox3 MinimalBox3( int N, Point3* pt );
817 : : // Functions to find minimal volume box that surrounds a set of points.
818 : : #endif
819 : :
820 : : double Abs (double x);
821 : : double ACos (double x);
822 : : double ATan2 (double y, double x);
823 : : double Cos (double x);
824 : : double Sign (double x);
825 : : double Sin (double x);
826 : : double Sqrt (double x);
827 : : double UnitRandom ();
828 : : double Dot (const Point2& p, const Point2& q);
829 : : double Length (const Point2& p);
830 : : CubitBoolean Unitize (Point2& p );
831 : : double Dot (const Point3& p, const Point3& q);
832 : : Point3 Cross (const Point3& p, const Point3& q);
833 : : double Length (const Point3& p);
834 : : CubitBoolean Unitize (Point3& p );
835 : : // Supporting code for triangle calculations
836 : :
837 : : double MinLineSegmentLineSegment (const Line3& seg0, const Line3& seg1, double& s, double& t);
838 : :
839 : : //get intersection between bounding box and plane
840 : : int get_plane_bbox_intersections( double box_min[3],
841 : : double box_max[3],
842 : : double pln_orig[3],
843 : : double pln_norm[3],
844 : : double int_array[6][3] );
845 : : // FILE: pt3tri3.cpp
846 : :
847 : :
848 : : // FILE: lin3tri3.cpp
849 : : double MinLineSegmentTriangle (const Line3& seg, const Triangle3& tri, double& r, double& s, double& t);
850 : : // Code to support distance between triangles
851 : :
852 : : //FILE: triasect.cpp
853 : : AgtLine EdgeToLine (Point2* v0, Point2* v1);
854 : : void TriangleLines (Triangle* T, AgtLine line[3]);
855 : : AgtTriList* SplitAndDecompose (Triangle* T, AgtLine* line);
856 : : AgtTriList* Intersection (Triangle* T0, Triangle* T1);
857 : : double AreaTriangle (Triangle* T);
858 : : double Area (AgtTriList* list);
859 : : double Area (Triangle &tri1, Triangle &tri2 );
860 : : // Code to support area of overlap of two triangles
861 : : };
862 : :
863 : : #if 1
864 : : // FILE: eigen.h
865 : : class mgcEigen
866 : : {
867 : : public:
868 : : mgcEigen (int _size);
869 : : ~mgcEigen ();
870 : :
871 : : mgcEigen& Matrix (float** inmat);
872 : :
873 : : private:
874 : : int size;
875 : : float** mat;
876 : : float* diag;
877 : : float* subd;
878 : :
879 : : // Householder reduction to tridiagonal form
880 : : void Tridiagonal2 (float** pmat, float* pdiag, float* psubd);
881 : : void Tridiagonal3 (float** pmat, float* pdiag, float* psubd);
882 : : void Tridiagonal4 (float** pmat, float* pdiag, float* psubd);
883 : : void TridiagonalN (int n, float** mat, float* diag, float* subd);
884 : :
885 : : // QL algorithm with implicit shifting, applies to tridiagonal matrices
886 : : void QLAlgorithm (int n, float* pdiag, float* psubd, float** pmat);
887 : :
888 : : // sort eigenvalues from largest to smallest
889 : : void DecreasingSort (int n, float* eigval, float** eigvec);
890 : :
891 : : // sort eigenvalues from smallest to largest
892 : : void IncreasingSort (int n, float* eigval, float** eigvec);
893 : :
894 : : // error handling
895 : : public:
896 : : static int verbose1;
897 : : static unsigned error;
898 : : static void Report (std::ostream& ostr);
899 : : private:
900 : : static const unsigned invalid_size;
901 : : static const unsigned allocation_failed;
902 : : static const unsigned ql_exceeded;
903 : : static const char* message[3];
904 : : static int Number (unsigned single_error);
905 : : static void Report (unsigned single_error);
906 : : };
907 : :
908 : : // FILE: eigen.h
909 : : class mgcEigenD
910 : : {
911 : : public:
912 : : mgcEigenD (int _size);
913 : : ~mgcEigenD ();
914 : :
915 : : // set the matrix for eigensolving
916 : 0 : double& Matrix (int row, int col) { return mat[row][col]; }
917 : : mgcEigenD& Matrix (double** inmat);
918 : :
919 : : // get the results of eigensolving
920 : : double Eigenvector (int row, int col) { return mat[row][col]; }
921 : 0 : const double** Eigenvector () { return (const double**) mat; }
922 : :
923 : : void EigenStuff3 (); // uses TriDiagonal3
924 : :
925 : : private:
926 : : int size;
927 : : double** mat;
928 : : double* diag;
929 : : double* subd;
930 : :
931 : : // Householder reduction to tridiagonal form
932 : : void Tridiagonal2 (double** mat, double* diag, double* subd);
933 : : void Tridiagonal3 (double** mat, double* diag, double* subd);
934 : : void Tridiagonal4 (double** mat, double* diag, double* subd);
935 : : void TridiagonalN (int n, double** mat, double* diag, double* subd);
936 : :
937 : : // QL algorithm with implicit shifting, applies to tridiagonal matrices
938 : : void QLAlgorithm (int n, double* diag, double* subd, double** mat);
939 : :
940 : : // sort eigenvalues from largest to smallest
941 : : void DecreasingSort (int n, double* eigval, double** eigvec);
942 : :
943 : : // sort eigenvalues from smallest to largest
944 : : void IncreasingSort (int n, double* eigval, double** eigvec);
945 : :
946 : : // error handling
947 : : public:
948 : : static int verbose1;
949 : : static unsigned error;
950 : : static void Report (std::ostream& ostr);
951 : : private:
952 : : static const unsigned invalid_size;
953 : : static const unsigned allocation_failed;
954 : : static const unsigned ql_exceeded;
955 : : static const char* message[3];
956 : : static int Number (unsigned single_error);
957 : : static void Report (unsigned single_error);
958 : : };
959 : : #endif
960 : :
961 : : ///// MAGIC SOFTWARE - INLINE FUNCTIONS
962 : : //---------------------------------------------------------------------------
963 : : // Wrapped math functions
964 : : //---------------------------------------------------------------------------
965 : 0 : inline double AnalyticGeometryTool::Abs (double x)
966 : : {
967 : 0 : return double(fabs(x));
968 : : }
969 : :
970 : : inline double AnalyticGeometryTool::ATan2 (double y, double x)
971 : : {
972 : : return double(atan2(y,x));
973 : : }
974 : :
975 : 0 : inline double AnalyticGeometryTool::Sqrt (double x)
976 : : {
977 : 0 : return double(sqrt(x));
978 : : }
979 : :
980 : : inline double AnalyticGeometryTool::UnitRandom ()
981 : : {
982 : : return double(rand())/double(RAND_MAX);
983 : : }
984 : :
985 : : //---------------------------------------------------------------------------
986 : : // Points in 2D
987 : : //---------------------------------------------------------------------------
988 : :
989 : 0 : inline double AnalyticGeometryTool::Dot (const Point3& p, const Point3& q)
990 : : {
991 : 0 : return double(p.x*q.x + p.y*q.y + p.z*q.z);
992 : : }
993 : :
994 : :
995 : 0 : inline double AnalyticGeometryTool::set_epsilon( double new_epsilon )
996 : : {
997 : 0 : double old_epsilon = agtEpsilon;
998 : 0 : agtEpsilon = new_epsilon;
999 : 0 : return old_epsilon;
1000 : : }
1001 : :
1002 : 0 : inline CubitBoolean AnalyticGeometryTool::dbl_eq(double val_1,double val_2)
1003 : : {
1004 : : CubitBoolean result;
1005 [ # # ]: 0 : if (fabs(val_1 - val_2) < agtEpsilon)
1006 : 0 : result = CUBIT_TRUE;
1007 : : else
1008 : 0 : result = CUBIT_FALSE;
1009 : 0 : return result;
1010 : : }
1011 : :
1012 : : #endif
1013 : :
|
