mPlane.h

Engine/source/math/mPlane.h

More...

Classes:

class
class

A 3D plane defined by a normal and a distance along the normal.

Public Defines

define
define
PlaneSwitchCode(s, e) (s * 3 + e)

Detailed Description

Public Defines

PARALLEL_PLANE() 1e20f
PlaneSwitchCode(s, e) (s * 3 + e)
  1
  2//-----------------------------------------------------------------------------
  3// Copyright (c) 2012 GarageGames, LLC
  4//
  5// Permission is hereby granted, free of charge, to any person obtaining a copy
  6// of this software and associated documentation files (the "Software"), to
  7// deal in the Software without restriction, including without limitation the
  8// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
  9// sell copies of the Software, and to permit persons to whom the Software is
 10// furnished to do so, subject to the following conditions:
 11//
 12// The above copyright notice and this permission notice shall be included in
 13// all copies or substantial portions of the Software.
 14//
 15// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 16// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 17// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 18// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 19// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 20// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 21// IN THE SOFTWARE.
 22//-----------------------------------------------------------------------------
 23
 24#ifndef _MPLANE_H_
 25#define _MPLANE_H_
 26
 27#ifndef _MPOINT3_H_
 28#include "math/mPoint3.h"
 29#endif
 30
 31#ifndef _MBOX_H_
 32#include "math/mBox.h"
 33#endif
 34
 35#ifndef _MORIENTEDBOX_H_
 36#include "math/mOrientedBox.h"
 37#endif
 38
 39#ifndef _MSPHERE_H_
 40#include "math/mSphere.h"
 41#endif
 42
 43
 44/// A 3D plane defined by a normal and a distance along the normal.
 45///
 46/// @note The distance is stored negative.
 47class PlaneF : public Point3F
 48{
 49   public:
 50
 51      /// *NEGATIVE* distance along the xyz normal.
 52      F32 d;
 53
 54      /// Return the plane's normal.
 55      const Point3F& getNormal() const { return *this; }
 56
 57      /// Return the plane's position.
 58      Point3F getPosition() const { return Point3F( x, y, z ) * -d; }
 59
 60      bool isHorizontal() const;
 61      bool isVertical() const;
 62
 63      bool isParallelTo( const PlaneF& plane, F32 epsilon = POINT_EPSILON ) const;
 64
 65      bool isPerpendicularTo( const PlaneF& plane, F32 epsilon = POINT_EPSILON ) const;
 66
 67      /// @name Initialization
 68      /// @{
 69
 70      PlaneF() :d(0.0f) {}
 71      PlaneF( const Point3F& p, const Point3F& n );
 72      /// NOTE: d is the NEGATIVE distance along the xyz normal.
 73      PlaneF( F32 _x, F32 _y, F32 _z, F32 _d);
 74      PlaneF( const Point3F& j, const Point3F& k, const Point3F& l );
 75
 76      void set( F32 _x, F32 _y, F32 _z );
 77      /// NOTE: d is the NEGATIVE distance along the xyz normal.
 78      void set( F32 _x, F32 _y, F32 _z, F32 _d );
 79      void set( const Point3F& p, const Point3F& n);
 80      void set( const Point3F& k, const Point3F& j, const Point3F& l );
 81      void setPoint(const Point3F &p); // assumes the x,y,z fields are already set
 82                                           // creates an un-normalized plane
 83
 84      void setXY(F32 zz);
 85      void setYZ(F32 xx);
 86      void setXZ(F32 yy);
 87      void setXY(const Point3F& P, F32 dir);
 88      void setYZ(const Point3F& P, F32 dir);
 89      void setXZ(const Point3F& P, F32 dir);
 90
 91      /// @}
 92
 93      void shiftX(F32 xx);
 94      void shiftY(F32 yy);
 95      void shiftZ(F32 zz);
 96      void invert();
 97      void neg();
 98
 99      /// Return the distance of the given point from the plane.
100      F32 distToPlane( const Point3F& cp ) const;
101
102      /// Project the given point onto the plane.
103      Point3F project(const Point3F &pt) const;
104
105      /// @name Intersection
106      /// @{
107
108      enum Side
109      {
110         Front = 1,
111         On    = 0,
112         Back  = -1
113      };
114
115      /// Compute on which side of the plane the given point lies.
116      Side whichSide( const Point3F& cp ) const;
117
118      /// Compute which side the given sphere lies on.
119      Side whichSide( const SphereF& sphere ) const;
120
121      /// Compute which side the given AABB lies on.
122      Side whichSide( const Box3F& aabb ) const;
123
124      /// Compute which side the given OBB lies on.
125      Side whichSide( const OrientedBox3F& obb ) const;
126
127      /// Compute which side the given box lies on.
128      ///
129      /// @param center Center point.
130      /// @param axisx X-Axis vector.  Length is box half-extent along X.
131      /// @param axisy Y-Axis vector.  Length is box half-extent along Y.
132      /// @param axisz Z-Axis vector.  Length is box half-extent along Z.
133      ///
134      /// @return Side that the given box is on.
135      Side whichSideBox( const Point3F& center,
136                         const Point3F& axisx,
137                         const Point3F& axisy,
138                         const Point3F& axisz ) const;
139
140      /// @}
141
142      /// @name Intersection
143      /// @{
144
145      /// Compute the distance at which the ray traveling from @a start in the direction
146      /// of @end intersects the plane
147      /// @param start Starting point of the ray.
148      /// @param end Point in the direction of which the ray travels from @a start.
149      /// @return The distance (as a fraction/multiple of the length of the vector between
150      ///   @a start and @a end) at which the ray intersects the plane or PARALLEL_PLANE if the
151      ///   ray is parallel to the plane.
152      F32 intersect( const Point3F &start, const Point3F &end ) const;
153
154      /// Compute the intersection of the ray that emanates at @a start in the direction
155      /// @dir.
156      /// @param start Point where the ray emanates.
157      /// @param dir Direction in which the ray travels.  Must be normalized.
158      /// @param outHit Resulting intersection point.  Only set if there indeed is a hit.
159      /// @return True if the ray intersects the plane or false if not.
160      bool intersect( const Point3F &start, const Point3F &dir, Point3F *outHit ) const;
161
162      /// Compute the intersection between two planes.
163      ///
164      /// @param plane Plane to intersect with.
165      /// @param outLineOrigin Used to store the origin of the resulting intersection line.
166      /// @param outLineDirection Used to store the direction of the resulting intersection line.
167      ///
168      /// @return True if there is an intersection or false if the two planes are coplanar.
169      bool intersect( const PlaneF& plane, Point3F& outLineOrigin, Point3F& outLineDirection ) const;
170
171      /// @}
172
173      /// @name Clipping
174      /// @{
175
176      /// Clip a convex polygon by the plane.
177      ///
178      /// The resulting polygon will be the input polygon minus the part on the negative side
179      /// of the plane.  The input polygon must be convex and @a inVertices must be in CCW or CW order.
180      ///
181      /// @param inVertices Array holding the vertices of the input polygon.
182      /// @param inNumVertices Number of vertices in @a inVertices.
183      /// @param outVertices Array to hold the vertices of the clipped polygon.  Must have space for one additional
184      ///   vertex in case the polygon is split by the plane such that an additional vertex appears.  Must not
185      ///   be the same as @a inVertices.
186      /// @return Number of vertices in the clipped polygon, i.e. number of vertices in @a outVertices.
187      ///
188      /// @note Be aware that if the polygon fully lies on the negative side of the plane,
189      ///   the resulting @a outNumVertices will be zero, i.e. no polygon will result from the clip.
190      U32 clipPolygon( const Point3F* inVertices, U32 inNumVertices, Point3F* outVertices ) const;
191
192      /// Clip a line segment by the plane.
193      ///
194      /// @param start Start point of the line segment.
195      /// @param end End point of the line segment.
196      /// @param outNewEnd New end point if there is an intersection with the plane.
197      ///
198      /// @return True
199      bool clipSegment( const Point3F& start, const Point3F& end, Point3F& outNewEnd ) const;
200
201      /// @}
202};
203
204#define PARALLEL_PLANE  1e20f
205
206#define PlaneSwitchCode(s, e) (s * 3 + e)
207
208
209//---------------------------------------------------------------------------
210
211inline PlaneF::PlaneF( F32 _x, F32 _y, F32 _z, F32 _d )
212{
213   set( _x, _y, _z, _d );
214}
215
216//-----------------------------------------------------------------------------
217
218inline PlaneF::PlaneF( const Point3F& p, const Point3F& n )
219{
220   set(p,n);
221}
222
223//-----------------------------------------------------------------------------
224
225inline PlaneF::PlaneF( const Point3F& j, const Point3F& k, const Point3F& l )
226{
227   set(j,k,l);
228}
229
230//-----------------------------------------------------------------------------
231
232inline void PlaneF::setXY( F32 zz )
233{
234   x = y = 0.0f; z = 1.0f; d = -zz;
235}
236
237//-----------------------------------------------------------------------------
238
239inline void PlaneF::setYZ( F32 xx )
240{
241   x = 1.0f; z = y = 0.0f; d = -xx;
242}
243
244//-----------------------------------------------------------------------------
245
246inline void PlaneF::setXZ( F32 yy )
247{
248   x = z = 0.0f; y = 1.0f; d = -yy;
249}
250
251//-----------------------------------------------------------------------------
252
253inline void PlaneF::setXY(const Point3F& point, F32 dir)       // Normal = (0, 0, -1|1)
254{
255   x = y = 0.0f;
256   d = -((z = dir) * point.z);
257}
258
259//-----------------------------------------------------------------------------
260
261inline void PlaneF::setYZ(const Point3F& point, F32 dir)       // Normal = (-1|1, 0, 0)
262{
263   z = y = 0.0f;
264   d = -((x = dir) * point.x);
265}
266
267//-----------------------------------------------------------------------------
268
269inline void PlaneF::setXZ(const Point3F& point, F32 dir)       // Normal = (0, -1|1, 0)
270{
271   x = z = 0.0f;
272   d = -((y = dir) * point.y);
273}
274
275//-----------------------------------------------------------------------------
276
277inline void PlaneF::shiftX( F32 xx )
278{
279   d -= xx * x;
280}
281
282//-----------------------------------------------------------------------------
283
284inline void PlaneF::shiftY( F32 yy )
285{
286   d -= yy * y;
287}
288
289//-----------------------------------------------------------------------------
290
291inline void PlaneF::shiftZ( F32 zz )
292{
293   d -= zz * z;
294}
295
296//-----------------------------------------------------------------------------
297
298inline bool PlaneF::isHorizontal() const
299{
300   return (x == 0.0f && y == 0.0f) ? true : false;
301}
302
303//-----------------------------------------------------------------------------
304
305inline bool PlaneF::isVertical() const
306{
307    return ((x != 0.0f || y != 0.0f) && z == 0.0f) ? true : false;
308}
309
310//-----------------------------------------------------------------------------
311
312inline Point3F PlaneF::project(const Point3F &pt) const
313{
314   F32 dist = distToPlane(pt);
315   return Point3F(pt.x - x * dist, pt.y - y * dist, pt.z - z * dist);
316}
317
318//-----------------------------------------------------------------------------
319
320inline F32 PlaneF::distToPlane( const Point3F& cp ) const
321{
322   // return mDot(*this,cp) + d;
323   return (x * cp.x + y * cp.y + z * cp.z) + d;
324}
325
326//-----------------------------------------------------------------------------
327
328inline PlaneF::Side PlaneF::whichSide(const Point3F& cp) const
329{
330   F32 dist = distToPlane(cp);
331   if (dist >= 0.005f)                 // if (mFabs(dist) < 0.005f)
332      return Front;                    //    return On;
333   else if (dist <= -0.005f)           // else if (dist > 0.0f)
334      return Back;                     //    return Front;
335   else                                // else
336      return On;                       //    return Back;
337}
338
339//-----------------------------------------------------------------------------
340
341inline PlaneF::Side PlaneF::whichSide( const SphereF& sphere ) const
342{
343   const F32 dist = distToPlane( sphere.center );
344   if( dist > sphere.radius )
345      return Front;
346   else if( dist < - sphere.radius )
347      return Back;
348   else
349      return On;
350}
351
352//-----------------------------------------------------------------------------
353
354inline PlaneF::Side PlaneF::whichSide( const Box3F& aabb ) const
355{
356   // See Graphics Gems IV, 1.7 and "A Faster Overlap Test for a Plane and a Bounding Box"
357   // (http://replay.waybackmachine.org/19981203032829/http://www.cs.unc.edu/~hoff/research/vfculler/boxplane.html)
358   // for details.
359
360   Point3F pVertex;
361   Point3F nVertex;
362
363   pVertex.x = ( x > 0.0f ) ? aabb.maxExtents.x : aabb.minExtents.x;
364   pVertex.y = ( y > 0.0f ) ? aabb.maxExtents.y : aabb.minExtents.y;
365   pVertex.z = ( z > 0.0f ) ? aabb.maxExtents.z : aabb.minExtents.z;
366
367   if( whichSide( pVertex ) == Back )
368      return Back;
369
370   nVertex.x = ( x > 0.0f ) ? aabb.minExtents.x : aabb.maxExtents.x;
371   nVertex.y = ( y > 0.0f ) ? aabb.minExtents.y : aabb.maxExtents.y;
372   nVertex.z = ( z > 0.0f ) ? aabb.minExtents.z : aabb.maxExtents.z;
373
374   if( whichSide( nVertex ) == Front )
375      return Front;
376
377   return On;
378}
379
380//-----------------------------------------------------------------------------
381
382inline PlaneF::Side PlaneF::whichSide( const OrientedBox3F& obb ) const
383{
384   // Project the box onto the line defined by the plane center and normal.
385   // See "3D Game Engine Design" chapter 4.3.2.
386
387   Point3F mObbHalf = obb.getHalfExtents();
388   const F32 r = mObbHalf.x * mFabs( mDot( obb.getAxis( 0 ), *this ) ) +
389                 mObbHalf.y * mFabs( mDot( obb.getAxis( 1 ), *this ) ) +
390                 mObbHalf.z * mFabs( mDot( obb.getAxis( 2 ), *this ) );
391
392   const F32 dist = distToPlane( obb.getCenter() );
393   if( dist > r )
394      return Front;
395   else if( dist < - r )
396      return Back;
397   else
398      return On;
399}
400
401//-----------------------------------------------------------------------------
402
403inline PlaneF::Side PlaneF::whichSideBox(const Point3F& center,
404                                         const Point3F& axisx,
405                                         const Point3F& axisy,
406                                         const Point3F& axisz) const
407{
408   F32 baseDist = distToPlane(center);
409
410   F32 compDist = mFabs(mDot(axisx, *this)) +
411      mFabs(mDot(axisy, *this)) +
412      mFabs(mDot(axisz, *this));
413
414   if (baseDist >= compDist)
415      return Front;
416   else if (baseDist <= -compDist)
417      return Back;
418   else
419      return On;
420}
421
422inline void PlaneF::set( F32 _x, F32 _y, F32 _z )
423{
424    Point3F::set(_x,_y,_z);
425}
426
427inline void PlaneF::set( F32 _x, F32 _y, F32 _z, F32 _d )
428{
429    Point3F::set(_x,_y,_z);
430    d = _d;
431}
432
433//---------------------------------------------------------------------------
434/// Calculate the coefficients of the plane passing through
435/// a point with the given normal.
436inline void PlaneF::setPoint(const Point3F &p)
437{
438   d = -(p.x * x + p.y * y + p.z * z);
439}
440
441inline void PlaneF::set( const Point3F& p, const Point3F& n )
442{
443   x = n.x; y = n.y; z = n.z;
444   normalize();
445
446   // Calculate the last plane coefficient.
447
448   d = -(p.x * x + p.y * y + p.z * z);
449}
450
451//---------------------------------------------------------------------------
452// Calculate the coefficients of the plane passing through
453// three points.  Basically it calculates the normal to the three
454// points then calculates a plane through the middle point with that
455// normal.
456inline void PlaneF::set( const Point3F& k, const Point3F& j, const Point3F& l )
457{
458   // Point3F kj, lj, pv;
459   // kj = k - j;
460   // lj = l - j;
461   // mCross( kj, lj, &pv );
462   // set(j, pv);
463
464   // Inline for speed...
465   const F32 ax = k.x-j.x;
466   const F32 ay = k.y-j.y;
467   const F32 az = k.z-j.z;
468   const F32 bx = l.x-j.x;
469   const F32 by = l.y-j.y;
470   const F32 bz = l.z-j.z;
471   x = ay*bz - az*by;
472   y = az*bx - ax*bz;
473   z = ax*by - ay*bx;
474      
475   m_point3F_normalize( (F32 *)(&x) );
476   d = -(j.x * x + j.y * y + j.z * z);
477}
478
479inline void PlaneF::invert()
480{
481   x = -x;
482   y = -y;
483   z = -z;
484   d = -d;
485}
486
487inline void PlaneF::neg()
488{
489   invert();
490}
491
492inline F32 PlaneF::intersect(const Point3F &p1, const Point3F &p2) const
493{
494   const F32 den = mDot(p2 - p1, *this);
495   if( mIsZero( den ) )
496      return PARALLEL_PLANE;
497   return -distToPlane(p1) / den;
498}
499
500inline bool PlaneF::intersect( const Point3F &start, const Point3F &dir, Point3F *outHit ) const
501{
502   const F32 den = mDot( dir, *this );
503   if ( mIsZero( den ) )
504      return false;
505
506   F32 dist = -distToPlane( start ) / den;
507   *outHit = start + dir * dist;
508   return true;
509}
510
511class PlaneD: public Point3D
512{
513public:
514   /// NOTE: d is the NEGATIVE distance along the xyz normal.
515   F64 d;
516
517   PlaneD();
518   PlaneD( const PlaneF& copy);
519   PlaneD( const Point3D& p, const Point3D& n );
520   /// NOTE: d is the NEGATIVE distance along the xyz normal.
521   PlaneD( F64 _x, F64 _y, F64 _z, F64 _d);
522   PlaneD( const Point3D& j, const Point3D& k, const Point3D& l );
523
524   // Methods
525   //using Point3D::set;
526    void set(const F64 _x, const F64 _y, const F64 _z);
527
528   void     set( const Point3D& p, const Point3D& n);
529   void     set( const Point3D& k, const Point3D& j, const Point3D& l );
530   void     setPoint(const Point3D &p); // assumes the x,y,z fields are already set
531                                        // creates an un-normalized plane
532
533   void     setXY(F64 zz);
534   void     setYZ(F64 xx);
535   void     setXZ(F64 yy);
536   void     setXY(const Point3D& P, F64 dir);
537   void     setYZ(const Point3D& P, F64 dir);
538   void     setXZ(const Point3D& P, F64 dir);
539   void     shiftX(F64 xx);
540   void     shiftY(F64 yy);
541   void     shiftZ(F64 zz);
542   void     invert();
543   void     neg();
544   Point3D  project(const Point3D &pt) const; // projects the point onto the plane.
545
546   F64      distToPlane( const Point3D& cp ) const;
547
548   enum Side {
549      Front = 1,
550      On    = 0,
551      Back  = -1
552   };
553
554   Side whichSide(const Point3D& cp) const;
555   F64  intersect(const Point3D &start, const Point3D &end) const;
556   //DLLAPI bool     split( const Poly3F& poly, Poly3F* front, Poly3F* back );
557
558   bool     isHorizontal() const;
559   bool     isVertical() const;
560
561   Side whichSideBox(const Point3D& center,
562                     const Point3D& axisx,
563                     const Point3D& axisy,
564                     const Point3D& axisz) const;
565};
566//#define PARALLEL_PLANE  1e20f
567
568//#define PlaneSwitchCode(s, e) (s * 3 + e)
569
570
571//---------------------------------------------------------------------------
572
573inline PlaneD::PlaneD()
574{
575   d = 0.0;
576}
577
578inline PlaneD::
579PlaneD( F64 _x, F64 _y, F64 _z, F64 _d )
580{
581   x = _x; y = _y; z = _z; d = _d;
582}
583
584inline PlaneD::PlaneD( const PlaneF& copy)
585{
586   x = copy.x; y = copy.y; z = copy.z; d = copy.d;
587}
588
589inline PlaneD::PlaneD( const Point3D& p, const Point3D& n )
590{
591   set(p,n);
592}
593
594inline PlaneD::PlaneD( const Point3D& j, const Point3D& k, const Point3D& l )
595{
596   set(j,k,l);
597}
598
599inline void PlaneD::setXY( F64 zz )
600{
601   x = y = 0; z = 1; d = -zz;
602}
603
604inline void PlaneD::setYZ( F64 xx )
605{
606   x = 1; z = y = 0; d = -xx;
607}
608
609inline void PlaneD::setXZ( F64 yy )
610{
611   x = z = 0; y = 1; d = -yy;
612}
613
614inline void PlaneD::setXY(const Point3D& point, F64 dir)       // Normal = (0, 0, -1|1)
615{
616   x = y = 0; 
617   d = -((z = dir) * point.z);
618}
619
620inline void PlaneD::setYZ(const Point3D& point, F64 dir)       // Normal = (-1|1, 0, 0)
621{
622   z = y = 0; 
623   d = -((x = dir) * point.x);
624}
625
626inline void PlaneD::setXZ(const Point3D& point, F64 dir)       // Normal = (0, -1|1, 0)
627{
628   x = z = 0; 
629   d = -((y = dir) * point.y);
630}
631
632inline void PlaneD::shiftX( F64 xx )
633{
634   d -= xx * x;
635}
636
637inline void PlaneD::shiftY( F64 yy )
638{
639   d -= yy * y;
640}
641
642inline void PlaneD::shiftZ( F64 zz )
643{
644   d -= zz * z;
645}
646
647inline bool PlaneD::isHorizontal() const
648{
649   return (x == 0 && y == 0) ? true : false;
650}
651
652inline bool PlaneD::isVertical() const
653{
654    return ((x != 0 || y != 0) && z == 0) ? true : false;
655}
656
657inline Point3D PlaneD::project(const Point3D &pt) const
658{
659   F64 dist = distToPlane(pt);
660   return Point3D(pt.x - x * dist, pt.y - y * dist, pt.z - z * dist);
661}
662
663inline F64 PlaneD::distToPlane( const Point3D& cp ) const
664{
665   // return mDot(*this,cp) + d;
666   return (x * cp.x + y * cp.y + z * cp.z) + d;
667}
668
669inline PlaneD::Side PlaneD::whichSide(const Point3D& cp) const
670{
671   F64 dist = distToPlane(cp);
672   if (dist >= 0.005f)                 // if (mFabs(dist) < 0.005f)
673      return Front;                    //    return On;
674   else if (dist <= -0.005f)           // else if (dist > 0.0f)
675      return Back;                     //    return Front;
676   else                                // else
677      return On;                       //    return Back;
678}
679
680inline void PlaneD::set(const F64 _x, const F64 _y, const F64 _z)
681{
682    Point3D::set(_x,_y,_z); 
683}
684
685//---------------------------------------------------------------------------
686// Calculate the coefficients of the plane passing through 
687// a point with the given normal.
688
689////inline void PlaneD::set( const Point3D& p, const Point3D& n )
690inline void PlaneD::setPoint(const Point3D &p)
691{
692   d = -(p.x * x + p.y * y + p.z * z);
693}
694
695inline void PlaneD::set( const Point3D& p, const Point3D& n )
696{
697   x = n.x; y = n.y; z = n.z;
698   normalize();
699
700   // Calculate the last plane coefficient.
701
702   d = -(p.x * x + p.y * y + p.z * z);
703}
704
705//---------------------------------------------------------------------------
706// Calculate the coefficients of the plane passing through 
707// three points.  Basically it calculates the normal to the three
708// points then calculates a plane through the middle point with that
709// normal.
710
711inline void PlaneD::set( const Point3D& k, const Point3D& j, const Point3D& l )
712{
713//   Point3D kj,lj,pv;
714//   kj = k;
715//   kj -= j;
716//   lj = l;
717//   lj -= j;
718//   mCross( kj, lj, &pv );
719//   set(j,pv);
720
721// Above ends up making function calls up the %*#...
722// This is called often enough to be a little more direct
723// about it (sqrt should become intrinsic in the future)...
724   F64 ax = k.x-j.x;
725   F64 ay = k.y-j.y;
726   F64 az = k.z-j.z;
727   F64 bx = l.x-j.x;
728   F64 by = l.y-j.y;
729   F64 bz = l.z-j.z;
730   x = ay*bz - az*by;
731   y = az*bx - ax*bz;
732   z = ax*by - ay*bx;
733   F64 squared = x*x + y*y + z*z;
734   AssertFatal(squared != 0.0, "Error, no plane possible!");
735
736   // In non-debug mode
737   if (squared != 0) 
738   {
739      F64 invSqrt = 1.0 / (F64)mSqrt(x*x + y*y + z*z);
740      x *= invSqrt;
741      y *= invSqrt;
742      z *= invSqrt;
743      d = -(j.x * x + j.y * y + j.z * z);
744   }
745   else 
746   {
747      x = 0;
748      y = 0;
749      z = 1;
750      d = -(j.x * x + j.y * y + j.z * z);
751   }
752}
753
754inline void PlaneD::invert()
755{
756   x = -x;
757   y = -y;
758   z = -z;
759   d = -d;
760}
761
762inline void PlaneD::neg()
763{
764   invert();
765}
766
767inline F64 PlaneD::intersect(const Point3D &p1, const Point3D &p2) const
768{
769   F64 den = mDot(p2 - p1, *this);
770   if(den == 0)
771      return PARALLEL_PLANE;
772   return -distToPlane(p1) / den;
773}
774
775inline PlaneD::Side PlaneD::whichSideBox(const Point3D& center,
776                                         const Point3D& axisx,
777                                         const Point3D& axisy,
778                                         const Point3D& axisz) const
779{
780   F64 baseDist = distToPlane(center);
781
782   F64 compDist = mFabs(mDot(axisx, *this)) +
783                  mFabs(mDot(axisy, *this)) +
784                  mFabs(mDot(axisz, *this));
785
786   // Intersects
787   if (baseDist >= compDist)
788      return Front;
789   else if (baseDist <= -compDist)
790      return Back;
791   else
792      return On;
793
794//   if (compDist > mFabs(baseDist))
795//      return On;
796//   else
797//      return baseDist < 0.0 ? Back : Front;
798}
799
800#endif // _MPLANE_H_
801