source/elementarycoords.h

   1 /* This file is part of Shapes.
   2  *
   3  * Shapes is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * any later version.
   7  *
   8  * Shapes is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with Shapes.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright 2008, 2009 Henrik Tidefelt
  17  */
  18
  19 #pragma once
  20
  21 #include <cmath>
  22
  23 #include "Shapes_Concrete_decls.h"
  24 #include "Shapes_Lang_decls.h"
  25
  26 #include "elementarylength.h"
  27 #include "shapesexceptions.h"
  28
  29
  30 namespace Shapes
  31 {
  32
  33         namespace Concrete
  34         {
  35
  36                 class UnitFloatPair
  37                 {
  38                 public:
  39                         double x_;
  40                         double y_;
  41                         UnitFloatPair( double x, double y )
  42                         {
  43                                 if( x == 0 && y == 0 )
  44                                         {
  45                                                 throw Exceptions::InternalError( "A UnitFloatPair was initialized with a zero vector." );
  46                                         }
  47                                 double tmp = 1. / hypot( x, y );
  48                                 x_ = tmp * x;
  49                                 y_ = tmp * y;
  50                         }
  51                         UnitFloatPair( bool ) : x_( 1 ), y_( 0 ) { }
  52                         UnitFloatPair( bool, bool ) : x_( 0 ), y_( 1 ) { }
  53                         explicit UnitFloatPair( double x, double y, double precomputedNorm )
  54                         {
  55                                 double tmp = 1. / precomputedNorm;
  56                                 x_ = tmp * x;
  57                                 y_ = tmp * y;
  58                         }
  59                         explicit UnitFloatPair( double x, double y, bool )
  60                         {
  61                                 x_ = x;
  62                                 y_ = y;
  63                         }
  64                         double normSquaredThatOughtToBeOne( ) const
  65                         {
  66                                 return x_ * x_ + y_ * y_;
  67                         }
  68                         Concrete::UnitFloatPair unnormalizedScaling( double a ) const;
  69                         Concrete::UnitFloatPair reverse( ) const /* Does not normalize the result!  Result will have same normSquaredThatOughtToBeOne( ) as *this.*/
  70                         {
  71                                 return Concrete::UnitFloatPair( - x_, - y_, bool ( ) );
  72                         }
  73                         Concrete::UnitFloatPair quarterTurnCounterClockwise( ) const /* Does not normalize the result!  Result will have same normSquaredThatOughtToBeOne( ) as *this.*/
  74                         {
  75                                 return Concrete::UnitFloatPair( - y_, x_, bool ( ) );
  76                         }
  77                 };
  78
  79                 class Coords2D
  80                 {
  81                 public:
  82                         Length x_;
  83                         Length y_;
  84                         Coords2D( Concrete::Length x, Concrete::Length y ) : x_( x ), y_( y ) { }
  85                         Coords2D( double x, double y ) : x_( x ), y_( y ) { }
  86                         Coords2D( const Concrete::Coords2D & orig ) : x_( orig.x_ ), y_( orig.y_ ) { }
  87                         Coords2D( const Lang::Coords2D & orig );
  88                         Coords2D( const Concrete::Coords2D & base, const Lang::Coords2D & orig );
  89                         operator Lang::Coords2D ( ) const;
  90                         Coords2D * transformedPtr( const Lang::Transform2D & tf ) const;
  91                         Coords2D transformed( const Lang::Transform2D & tf ) const;
  92                         Concrete::Coords3D * transformedPtr( const Lang::Transform3D & tf ) const; // treat as z = 0
  93                         friend std::ostream & operator << ( std::ostream & os, const Coords2D & self );
  94                         Concrete::Length norm( ) const;
  95                         Concrete::LengthSquared normSquared( ) const { return x_ * x_ + y_ * y_; }
  96                         double normScalar( ) const;
  97                         Concrete::UnitFloatPair direction( ) const;
  98                         Concrete::UnitFloatPair direction( Concrete::Length precomputedNorm ) const;
  99                         double angleCosineSignSquared( const Concrete::Coords2D & b ) const; /* Like the ordinary angle cosine, but avoids the use of sqrt. */
 100                         Concrete::UnitFloatPair normalizedOrthogonal( const Coords2D & b ) const
 101                         {
 102                                 // We return the counter-clockwise orthogonal to the normalized vector from ourselves to b
 103                                 return UnitFloatPair( y_.offtype< 1, 0 >( ) - b.y_.offtype< 1, 0 >( ),
 104                                                                                                                         b.x_.offtype< 1, 0 >( ) - x_.offtype< 1, 0 >( ) );
 105                         }
 106                         Concrete::UnitFloatPair normalizedOrthogonalNoFail( const Coords2D & b ) const
 107                         {
 108                                 // We return the counter-clockwise orthogonal to the normalized vector from ourselves to b
 109                                 try
 110                                         {
 111                                                 return normalizedOrthogonal( b );
 112                                         }
 113                                 catch( Exceptions::InternalError & ball )
 114                                         {
 115                                                 return UnitFloatPair( bool( ) );
 116                                         }
 117                         }
 118                         Concrete::UnitFloatPair unNormalizedOrthogonal( const Coords2D & b ) const
 119                         {
 120                                 // Warning!     This is the same as normalizedOrthogonal, but the result is really not normalized!
 121                                 return UnitFloatPair( y_.offtype< 1, 0 >( ) - b.y_.offtype< 1, 0 >( ),
 122                                                                                                                         b.x_.offtype< 1, 0 >( ) - x_.offtype< 1, 0 >( ),
 123                                                                                                                         bool( ) );
 124                         }
 125                         Concrete::Coords2D quarterTurnCounterClockwise( ) const
 126                         {
 127                                 return Concrete::Coords2D( - y_, x_ );
 128                         }
 129                         Concrete::Length distanceTo( const Concrete::Coords2D & other ) const
 130                         {
 131                                 return hypotPhysical( x_ - other.x_, y_ - other.y_ );
 132                         }
 133                 };
 134
 135                 /* The Bezier template classes requires the interface to a vector space with inner product.
 136                  */
 137                 Concrete::Coords2D operator * ( double scalar, const Concrete::Coords2D & coords );
 138                 Concrete::Coords2D operator * ( const Concrete::Coords2D & coords, double scalar );
 139                 Concrete::Coords2D operator + ( const Concrete::Coords2D & coords1, const Concrete::Coords2D & coords2 );
 140                 Concrete::Coords2D operator - ( const Concrete::Coords2D & coords1, const Concrete::Coords2D & coords2 );
 141
 142                 Concrete::Coords2D operator * ( Concrete::Length length, const Concrete::UnitFloatPair & direction );
 143                 Concrete::Coords2D operator * ( const Concrete::UnitFloatPair & direction, Concrete::Length length );
 144
 145                 //      Concrete::UnitFloatPair operator + ( const Concrete::UnitFloatPair & coords1, const Concrete::UnitFloatPair & coords2 );
 146                 //      Concrete::UnitFloatPair operator - ( const Concrete::UnitFloatPair & coords1, const Concrete::UnitFloatPair & coords2 );
 147
 148                 inline Concrete::Coords2D operator / ( const Concrete::Coords2D & coords, double denom )
 149                 {
 150                         double tmp = 1. / denom;
 151                         return Concrete::Coords2D( tmp * coords.x_, tmp * coords.y_ );
 152                 }
 153                 //       inline Concrete::UnitFloatPair operator / ( const Concrete::Coords2D & coords, const Concrete::Length & denom )
 154                 //       {
 155                 //               double tmp = 1. / denom.x_;
 156                 //               return Concrete::UnitFloatPair( tmp * coords.x_.x_, tmp * coords.y_.x_ );
 157                 //       }
 158                 //       inline Concrete::UnitFloatPair operator / ( const Concrete::UnitFloatPair & coords, double denom )
 159                 //       {
 160                 //               double tmp = 1. / denom;
 161                 //               return Concrete::UnitFloatPair( tmp * coords.x_, tmp * coords.y_ );
 162                 //       }
 163
 164                 inline bool operator == ( const Concrete::Coords2D & coords1, const Concrete::Coords2D & coords2 )
 165                 {
 166                         return coords1.x_ == coords2.x_ && coords1.y_ == coords2.y_;
 167                 }
 168                 inline bool operator != ( const Concrete::Coords2D & coords1, const Concrete::Coords2D & coords2 )
 169                 {
 170                         return ! ( coords1 == coords2 );
 171                 }
 172
 173                 inline double inner( const UnitFloatPair & a, const UnitFloatPair & b )
 174                 {
 175                         return a.x_ * b.x_ + a.y_ * b.y_;
 176                 }
 177
 178                 inline Length inner( const UnitFloatPair & a, const Coords2D & b )
 179                 {
 180                         return a.x_ * b.x_ + a.y_ * b.y_;
 181                 }
 182
 183                 inline Length inner( const Coords2D & a, const UnitFloatPair & b )
 184                 {
 185                         return a.x_ * b.x_ + a.y_ * b.y_;
 186                 }
 187
 188                 inline double innerScalar( const UnitFloatPair & a, const Coords2D & b )
 189                 {
 190                         return ( a.x_ * b.x_ + a.y_ * b.y_ ).offtype< 1, 0 >( );
 191                 }
 192
 193                 inline double innerScalar( const Coords2D & a, const UnitFloatPair & b )
 194                 {
 195                         return ( a.x_ * b.x_ + a.y_ * b.y_ ).offtype< 1, 0 >( );
 196                 }
 197
 198                 inline double innerScalar( const Coords2D & a, const Coords2D & b )
 199                 {
 200                         return ( a.x_ * b.x_ + a.y_ * b.y_ ).offtype< 2, 0 >( );
 201                 }
 202
 203
 204                 class UnitFloatTriple
 205                 {
 206                 public:
 207                         double x_;
 208                         double y_;
 209                         double z_;
 210                         UnitFloatTriple( double x, double y, double z );
 211                         UnitFloatTriple( bool ) : x_( 1 ), y_( 0 ), z_( 0 ) { }
 212                         UnitFloatTriple( bool, bool ) : x_( 0 ), y_( 1 ), z_( 0 ) { }
 213                         UnitFloatTriple( bool, bool, bool ) : x_( 0 ), y_( 0 ), z_( 1 ) { }
 214                         explicit UnitFloatTriple( double x, double y, double z, double precomputedNorm )
 215                         {
 216                                 double tmp = 1. / precomputedNorm;
 217                                 x_ = tmp * x;
 218                                 y_ = tmp * y;
 219                                 z_ = tmp * z;
 220                         }
 221                         explicit UnitFloatTriple( double x, double y, double z, bool )
 222                         {
 223                                 x_ = x;
 224                                 y_ = y;
 225                                 z_ = z;
 226                         }
 227                         double normSquaredThatOughtToBeOne( ) const
 228                         {
 229                                 return x_ * x_ + y_ * y_ + z_ * z_;
 230                         }
 231                         Concrete::UnitFloatTriple reflect( const Concrete::UnitFloatTriple dir ) const; /* Does not normalize the result!  Result will have same normSquaredThatOughtToBeOne( ) as *this.*/
 232                         Concrete::UnitFloatTriple unnormalizedScaling( double a ) const; /* Does not normalize the result!  Result will have same normSquaredThatOughtToBeOne( ) as *this.*/
 233                         Concrete::UnitFloatTriple reverse( ) const /* Does not normalize the result!  Result will have same normSquaredThatOughtToBeOne( ) as *this.*/
 234                         {
 235                                 return Concrete::UnitFloatTriple( - x_, - y_, - z_, bool( ) );
 236                         }
 237                         Concrete::UnitFloatTriple rotate( const Concrete::UnitFloatTriple dir, const double angle ) const;
 238                 };
 239
 240                 class Coords3D
 241                 {
 242                 public:
 243                         Length x_;
 244                         Length y_;
 245                         Length z_;
 246                         Coords3D( Concrete::Length x, Concrete::Length y, Concrete::Length z ) : x_( x ), y_( y ), z_( z ) { }
 247                         Coords3D( double x, double y, double z ) : x_( x ), y_( y ), z_( z ) { }
 248                         Coords3D( const Concrete::Coords3D & orig ) : x_( orig.x_ ), y_( orig.y_ ), z_( orig.z_ ) { }
 249                         Coords3D( const Concrete::Coords2D & orig ) : x_( orig.x_ ), y_( orig.y_ ), z_( 0 ) { }
 250                         Coords3D( const Lang::Coords3D & orig );
 251                         Coords3D( const Concrete::Coords3D & base, const Lang::Coords3D & orig );
 252                         operator Lang::Coords3D ( ) const;
 253                         Coords3D * transformedPtr( const Lang::Transform3D & tf ) const;
 254                         Coords3D transformed( const Lang::Transform3D & tf ) const;
 255                         Coords2D * make2D( Concrete::Length eyez ) const;
 256                         Coords2D make2DAutomatic( Concrete::Length eyez ) const;
 257                         friend std::ostream & operator << ( std::ostream & os, const Coords3D & self );
 258                         void negate( ) { x_ = -x_; y_ = -y_; z_ = -z_; }
 259                         double normalizedInner( const Coords3D & b ) const;
 260                         Concrete::Length norm( ) const;
 261                         Concrete::LengthSquared normSquared( ) const { return x_ * x_ + y_ * y_ + z_ * z_; }
 262                         double normScalar( ) const;
 263                         UnitFloatTriple direction( ) const;
 264                         UnitFloatTriple directionNoFail( ) const; /* This function returns an arbitrary direction when direction would fail. */
 265                         Concrete::UnitFloatTriple direction( Concrete::Length precomputedNorm ) const;
 266                 };
 267
 268                 /* The Bezier template classes requires the interface to a vector space.
 269                  */
 270                 Concrete::Coords3D operator * ( double scalar, const Concrete::Coords3D & coords );
 271                 Concrete::Coords3D operator * ( const Concrete::Coords3D & coords, double scalar );
 272                 Concrete::Coords3D operator + ( const Concrete::Coords3D & coords1, const Concrete::Coords3D & coords2 );
 273                 Concrete::Coords3D operator - ( const Concrete::Coords3D & coords1, const Concrete::Coords3D & coords2 );
 274
 275                 Concrete::Coords3D operator * ( Concrete::Length length, const Concrete::UnitFloatTriple & direction );
 276                 Concrete::Coords3D operator * ( const Concrete::UnitFloatTriple & direction, Concrete::Length length );
 277
 278                 //      Concrete::UnitFloatTriple operator + ( const Concrete::UnitFloatTriple & coords1, const Concrete::UnitFloatTriple & coords2 );
 279                 //      Concrete::UnitFloatTriple operator - ( const Concrete::UnitFloatTriple & coords1, const Concrete::UnitFloatTriple & coords2 );
 280
 281                 inline Concrete::Coords3D operator / ( const Concrete::Coords3D & coords, double denom )
 282                 {
 283                         double tmp = 1. / denom;
 284                         return Concrete::Coords3D( tmp * coords.x_, tmp * coords.y_, tmp * coords.z_ );
 285                 }
 286                 //       inline Concrete::UnitFloatTriple operator / ( const Concrete::Coords3D & coords, const Concrete::Length & denom )
 287                 //       {
 288                 //               double tmp = 1. / denom.x_;
 289                 //               return Concrete::UnitFloatTriple( tmp * coords.x_.x_, tmp * coords.y_.x_, tmp * coords.z_.x_ );
 290                 //       }
 291                 //       inline Concrete::UnitFloatTriple operator / ( const Concrete::UnitFloatTriple & coords, double denom )
 292                 //       {
 293                 //               double tmp = 1. / denom;
 294                 //               return Concrete::UnitFloatTriple( tmp * coords.x_, tmp * coords.y_, tmp * coords.z_ );
 295                 //       }
 296
 297                 inline bool operator == ( const Concrete::Coords3D & coords1, const Concrete::Coords3D & coords2 )
 298                 {
 299                         return coords1.x_ == coords2.x_ && coords1.y_ == coords2.y_ && coords1.z_ == coords2.z_;
 300                 }
 301                 inline bool operator != ( const Concrete::Coords3D & coords1, const Concrete::Coords3D & coords2 )
 302                 {
 303                         return ! ( coords1 == coords2 );
 304                 }
 305
 306                 inline double inner( const UnitFloatTriple & a, const UnitFloatTriple & b )
 307                 {
 308                         return a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_;
 309                 }
 310
 311                 inline Length inner( const UnitFloatTriple & a, const Coords3D & b )
 312                 {
 313                         return a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_;
 314                 }
 315
 316                 inline Length inner( const Coords3D & a, const UnitFloatTriple & b )
 317                 {
 318                         return a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_;
 319                 }
 320
 321                 inline LengthSquared inner( const Coords3D & a, const Coords3D & b )
 322                 {
 323                         return a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_;
 324                 }
 325
 326                 inline double innerScalar( const UnitFloatTriple & a, const Coords3D & b )
 327                 {
 328                         return ( a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_ ).offtype< 1, 0 >( );
 329                 }
 330
 331                 inline double innerScalar( const Coords3D & a, const UnitFloatTriple & b )
 332                 {
 333                         return ( a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_ ).offtype< 1, 0 >( );
 334                 }
 335
 336                 inline double innerScalar( const Coords3D & a, const Coords3D & b )
 337                 {
 338                         return ( a.x_ * b.x_ + a.y_ * b.y_ + a.z_ * b.z_ ).offtype< 2, 0 >( );
 339                 }
 340
 341                 /*       inline Coords3D cross( const Coords3D & a, const Coords3D & b ) */
 342                 /*       { */
 343                 /*               return Coords3D( a.y_ * b.z_ - a.z_ * b.y_, */
 344                 /*                                                                                                                              a.z_ * b.x_ - a.x_ * b.z_, */
 345                 /*                                                                                                                              a.x_ * b.y_ - a.y_ * b.x_ ); */
 346                 /*       } */
 347
 348                 inline UnitFloatTriple crossDirection( const UnitFloatTriple & a, const UnitFloatTriple & b )
 349                 {
 350                         try
 351                                 {
 352                                         return UnitFloatTriple( a.y_ * b.z_ - a.z_ * b.y_,
 353                                                                                                                                         a.z_ * b.x_ - a.x_ * b.z_,
 354                                                                                                                                         a.x_ * b.y_ - a.y_ * b.x_ );
 355                                 }
 356                         catch( const Exceptions::InternalError & ball )
 357                                 {
 358                                         throw NonLocalExit::CrossDirectionOfParallel( );
 359                                 }
 360                 }
 361
 362                 inline Coords3D cross( const UnitFloatTriple & a, const Coords3D & b )
 363                 {
 364                         return Coords3D( a.y_ * b.z_ - a.z_ * b.y_,
 365                                                                                          a.z_ * b.x_ - a.x_ * b.z_,
 366                                                                                          a.x_ * b.y_ - a.y_ * b.x_ );
 367                 }
 368
 369                 inline Coords3D cross( const Coords3D & a, const UnitFloatTriple & b )
 370                 {
 371                         return Coords3D( a.y_ * b.z_ - a.z_ * b.y_,
 372                                                                                          a.z_ * b.x_ - a.x_ * b.z_,
 373                                                                                          a.x_ * b.y_ - a.y_ * b.x_ );
 374                 }
 375
 376                 inline UnitFloatTriple crossDirection( const Coords3D & a, const Coords3D & b )
 377                 {
 378                         try
 379                                 {
 380                                         return UnitFloatTriple( ( a.y_ * b.z_ - a.z_ * b.y_ ).offtype< 2, 0 >( ),
 381                                                                                                                                         ( a.z_ * b.x_ - a.x_ * b.z_ ).offtype< 2, 0 >( ),
 382                                                                                                                                         ( a.x_ * b.y_ - a.y_ * b.x_ ).offtype< 2, 0 >( ) );
 383                                 }
 384                         catch( const Exceptions::InternalError & ball )
 385                                 {
 386                                         throw NonLocalExit::CrossDirectionOfParallel( );
 387                                 }
 388                 }
 389
 390                 //       inline Coords3D crossCoords3D( const Coords3D & a, const Coords3D & b )
 391                 //       {
 392                 //               return Coords3D( a.y_.x_ * b.z_.x_ - a.z_.x_ * b.y_.x_,
 393                 //                                                                                                                              a.z_.x_ * b.x_.x_ - a.x_.x_ * b.z_.x_,
 394                 //                                                                                                                              a.x_.x_ * b.y_.x_ - a.y_.x_ * b.x_.x_ );
 395                 //       }
 396
 397                 Area crossMagnitude( const Coords3D & a, const Coords3D & b );
 398
 399         }
 400 }