source/facetnormals.cc

   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 Henrik Tidefelt
  17  */
  18
  19 #include "drawabletypes.h"
  20 #include "lighttypes.h"
  21 #include "astclass.h"
  22 #include "constructorrepresentation.h"
  23 #include "globals.h"
  24
  25 using namespace Shapes;
  26
  27
  28 Lang::FacetNormalGray::FacetNormalGray( const Concrete::Coords3D & position,
  29                                                                                                                                                                 const RefCountPtr< const Lang::SpecularReflection > & reflections,
  30                                                                                                                                                                 const Concrete::UnitFloatTriple & reflectionUnitNormal,
  31                                                                                                                                                                 const Concrete::Gray lightMultiply,
  32                                                                                                                                                                 const RefCountPtr< const Lang::SpecularReflection > & autoScattering,
  33                                                                                                                                                                 const Concrete::Gray autoIntensity )
  34         : position_( position ),
  35                 reflections_( reflections ),
  36                 reflectionUnitNormal_( reflectionUnitNormal ),
  37                 lightMultiply_( lightMultiply ),
  38                 autoScattering_( autoScattering ),
  39                 autoIntensity_( autoIntensity )
  40 { }
  41
  42 Lang::FacetNormalGray::~FacetNormalGray( )
  43 { }
  44
  45 RefCountPtr< const Lang::Class > Lang::FacetNormalGray::TypeID( new Lang::SystemFinalClass( strrefdup( "FacetNormalGray" ) ) );
  46 TYPEINFOIMPL( FacetNormalGray );
  47
  48 Concrete::Gray
  49 Lang::FacetNormalGray::compute( const Concrete::Coords3D & point, const Concrete::UnitFloatTriple normal, const Lang::Transform3D & tf, const Concrete::Length eyez, const std::list< RefCountPtr< const Lang::LightSource > > & lights ) const
  50 {
  51         // Here, we must remember that the lights, but only the lights, are already transformed!
  52
  53         Concrete::UnitFloatTriple tf_normal = tf.transformPlaneUnitNormal( reflectionUnitNormal_ );
  54         //      Concrete::Coords3D tf_position( position_.transformed( tf ) );
  55         double grayVal = autoScattering_->illuminate( point, eyez, tf_normal, autoIntensity_ ).gr_;
  56
  57         typedef typeof lights ListType;
  58         for( ListType::const_iterator i = lights.begin( ); i != lights.end( ); ++i )
  59                 {
  60                         double tmp = 0;
  61                         {
  62                                 typedef const Lang::GrayLight LightType;
  63                                 LightType * light = dynamic_cast< LightType * >( i->getPtr( ) );
  64                                 if( light != 0 )
  65                                         {
  66                                                 tmp = reflections_->illuminate( point, eyez, tf_normal, *light ).gr_;
  67                                                 goto foundLightType;
  68                                         }
  69                         }
  70                         {
  71                                 typedef const Lang::RGBLight LightType;
  72                                 LightType * light = dynamic_cast< LightType * >( i->getPtr( ) );
  73                                 if( light != 0 )
  74                                         {
  75                                                 tmp = reflections_->illuminate( point, eyez, tf_normal, *light ).mean( );
  76                                                 goto foundLightType;
  77                                         }
  78                         }
  79                         throw Exceptions::InternalError( "A strange type of light." );
  80                 foundLightType:
  81                         tmp *= lightMultiply_.gr_;
  82                         grayVal += tmp;
  83                         if( grayVal > 1 )
  84                                 {
  85                                         grayVal = 1;
  86                                         break;
  87                                 }
  88                 }
  89
  90         return Concrete::Gray( grayVal );
  91 }
  92
  93 Concrete::Gray
  94 Lang::FacetNormalGray::compute( const Concrete::Coords3D & point, const Lang::Transform3D & tf, const Concrete::Length eyez, const std::list< RefCountPtr< const Lang::LightSource > > & lights ) const
  95 {
  96         return compute( point, reflectionUnitNormal_, tf, eyez, lights );
  97 }
  98
  99 Concrete::Gray
 100 Lang::FacetNormalGray::getDebugColor( ) const
 101 {
 102         return lightMultiply_;
 103 }
 104
 105 RefCountPtr< const Lang::FacetNormalGray >
 106 Lang::FacetNormalGray::transformed( const Lang::Transform3D & tf ) const
 107 {
 108         return RefCountPtr< const Lang::FacetNormalGray >
 109                 ( new Lang::FacetNormalGray( position_.transformed( tf ),
 110                                                                                                                                  reflections_,
 111                                                                                                                                  tf.transformPlaneUnitNormal( reflectionUnitNormal_ ),
 112                                                                                                                                  lightMultiply_,
 113                                                                                                                                  autoScattering_,
 114                                                                                                                                  autoIntensity_ ) );
 115 }
 116
 117
 118 void
 119 Lang::FacetNormalGray::gcMark( Kernel::GCMarkedSet & marked )
 120 {
 121         const_cast< Lang::SpecularReflection * >( reflections_.getPtr( ) )->gcMark( marked );
 122         const_cast< Lang::SpecularReflection * >( autoScattering_.getPtr( ) )->gcMark( marked );
 123 }
 124
 125
 126 Lang::FacetNormalRGB::FacetNormalRGB( const Concrete::Coords3D & position,
 127                                                                                                                                                         const RefCountPtr< const Lang::SpecularReflection > & reflections,
 128                                                                                                                                                         const Concrete::UnitFloatTriple & reflectionUnitNormal,
 129                                                                                                                                                         const Concrete::RGB & lightMultiply,
 130                                                                                                                                                         const RefCountPtr< const Lang::SpecularReflection > & autoScattering,
 131                                                                                                                                                         const Concrete::RGB & autoIntensity )
 132         : position_( position ),
 133                 reflections_( reflections ),
 134                 reflectionUnitNormal_( reflectionUnitNormal ),
 135                 lightMultiply_( lightMultiply ),
 136                 autoScattering_( autoScattering ),
 137                 autoIntensity_( autoIntensity )
 138 { }
 139
 140 Lang::FacetNormalRGB::~FacetNormalRGB( )
 141 { }
 142
 143 RefCountPtr< const Lang::Class > Lang::FacetNormalRGB::TypeID( new Lang::SystemFinalClass( strrefdup( "FacetNormalRGB" ) ) );
 144 TYPEINFOIMPL( FacetNormalRGB );
 145
 146 Concrete::RGB
 147 Lang::FacetNormalRGB::compute( const Concrete::Coords3D & point, const Concrete::UnitFloatTriple normal, const Lang::Transform3D & tf, const Concrete::Length eyez, const std::list< RefCountPtr< const Lang::LightSource > > & lights ) const
 148 {
 149         // Should be analogous to FacetNormalGray::compute.
 150
 151         Concrete::UnitFloatTriple tf_normal = tf.transformPlaneUnitNormal( reflectionUnitNormal_ );
 152         //      Concrete::Coords3D tf_position( position_.transformed( tf ) );
 153         Concrete::RGB val = autoScattering_->illuminate( point, eyez, tf_normal, autoIntensity_ );
 154
 155         typedef typeof lights ListType;
 156         for( ListType::const_iterator i = lights.begin( ); i != lights.end( ); ++i )
 157                 {
 158                         Concrete::RGB tmp( 0, 0, 0);
 159                         {
 160                                 typedef const Lang::GrayLight LightType;
 161                                 LightType * light = dynamic_cast< LightType * >( i->getPtr( ) );
 162                                 if( light != 0 )
 163                                         {
 164                                                 double a = reflections_->illuminate( point, eyez, tf_normal, *light ).gr_;
 165                                                 tmp = Concrete::RGB( a, a, a );
 166                                                 goto foundLightType;
 167                                         }
 168                         }
 169                         {
 170                                 typedef const Lang::RGBLight LightType;
 171                                 LightType * light = dynamic_cast< LightType * >( i->getPtr( ) );
 172                                 if( light != 0 )
 173                                         {
 174                                                 tmp = reflections_->illuminate( point, eyez, tf_normal, *light );
 175                                                 goto foundLightType;
 176                                         }
 177                         }
 178                         throw Exceptions::InternalError( "A strange type of light." );
 179                 foundLightType:
 180                         tmp.r_ *= lightMultiply_.r_;
 181                         tmp.g_ *= lightMultiply_.g_;
 182                         tmp.b_ *= lightMultiply_.b_;
 183                         val = val.addNoCheck( tmp );
 184                         bool doBreak = true;
 185                         if( val.r_ > 1 )
 186                                 {
 187                                         val.r_ = 1;
 188                                 }
 189                         else
 190                                 {
 191                                         doBreak = false;
 192                                 }
 193                         if( val.g_ > 1 )
 194                                 {
 195                                         val.g_ = 1;
 196                                 }
 197                         else
 198                                 {
 199                                         doBreak = false;
 200                                 }
 201                         if( val.b_ > 1 )
 202                                 {
 203                                         val.b_ = 1;
 204                                 }
 205                         else
 206                                 {
 207                                         doBreak = false;
 208                                 }
 209                         if( doBreak )
 210                                 {
 211                                         break;
 212                                 }
 213                 }
 214
 215         return val;
 216 }
 217
 218 Concrete::RGB
 219 Lang::FacetNormalRGB::compute( const Concrete::Coords3D & point, const Lang::Transform3D & tf, const Concrete::Length eyez, const std::list< RefCountPtr< const Lang::LightSource > > & lights ) const
 220 {
 221         return compute( point, reflectionUnitNormal_, tf, eyez, lights );
 222 }
 223
 224 Concrete::RGB
 225 Lang::FacetNormalRGB::getDebugColor( ) const
 226 {
 227         return lightMultiply_;
 228 }
 229
 230 RefCountPtr< const Lang::FacetNormalRGB >
 231 Lang::FacetNormalRGB::transformed( const Lang::Transform3D & tf ) const
 232 {
 233         return RefCountPtr< const Lang::FacetNormalRGB >
 234                 ( new Lang::FacetNormalRGB( position_.transformed( tf ),
 235                                                                                                                                 reflections_,
 236                                                                                                                                 tf.transformPlaneUnitNormal( reflectionUnitNormal_ ),
 237                                                                                                                                 lightMultiply_,
 238                                                                                                                                 autoScattering_,
 239                                                                                                                                 autoIntensity_ ) );
 240 }
 241
 242 void
 243 Lang::FacetNormalRGB::gcMark( Kernel::GCMarkedSet & marked )
 244 {
 245         const_cast< Lang::SpecularReflection * >( reflections_.getPtr( ) )->gcMark( marked );
 246         const_cast< Lang::SpecularReflection * >( autoScattering_.getPtr( ) )->gcMark( marked );
 247 }
 248
 249
 250 Computation::FacetInterpolatorGray::~FacetInterpolatorGray( )
 251 { }
 252
 253
 254 Computation::FacetInterpolatorGray1::~FacetInterpolatorGray1( )
 255 { }
 256
 257 RefCountPtr< const Lang::Gray >
 258 Computation::FacetInterpolatorGray1::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez ) const
 259 {
 260         return RefCountPtr< const Lang::Gray >( new Lang::Gray( n1_->compute( point, tf, eyez, lights ) ) );
 261 }
 262
 263 RefCountPtr< const Lang::Gray >
 264 Computation::FacetInterpolatorGray1::getDebugColor( ) const
 265 {
 266         return RefCountPtr< const Lang::Gray >( new Lang::Gray( n1_->getDebugColor( ) ) );
 267 }
 268
 269 RefCountPtr< const Computation::FacetInterpolatorGray >
 270 Computation::FacetInterpolatorGray1::transformed( const Lang::Transform3D & tf ) const
 271 {
 272         return RefCountPtr< const Computation::FacetInterpolatorGray >
 273                 ( new FacetInterpolatorGray1( n1_->transformed( tf ) ) );
 274 }
 275
 276
 277 void
 278 Computation::FacetInterpolatorGray1::gcMark( Kernel::GCMarkedSet & marked )
 279 {
 280         const_cast< Lang::FacetNormalGray * >( n1_.getPtr( ) )->gcMark( marked );
 281 }
 282
 283
 284 Computation::FacetInterpolatorGray2::FacetInterpolatorGray2( const RefCountPtr< const Lang::FacetNormalGray > & n1,
 285                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalGray > & n2 )
 286         : n1_( n1 ), n2_( n2 ),
 287                 d_( 0 ), t_( 0, 0, 0, bool( ) ), m_( 0 ),
 288                 rotationDirection_( 0, 0, 0, bool( ) )
 289 {
 290         Concrete::Coords3D r = n2_->position( ) - n1_->position( );
 291         d_ = r.norm( );
 292         t_ = r.direction( d_ );
 293         m_ = Concrete::inner( t_, n1_->position( ) );
 294
 295         angle_ = acos( Concrete::inner( n1_->normal( ), n2_->normal( ) ) );
 296         try
 297                 {
 298                         rotationDirection_ = Concrete::crossDirection( n1_->normal( ), n2_->normal( ) );
 299                 }
 300         catch( const NonLocalExit::CrossDirectionOfParallel & ball )
 301                 {
 302                         if( angle_ > 3 )
 303                                 {
 304                                         throw Exceptions::MiscellaneousRequirement( "The facet normals on a facet must not point in opposite directions." );
 305                                 }
 306                         rotationDirection_ = Concrete::UnitFloatTriple( 1, 0, 0, bool( ) );
 307                 }
 308 }
 309
 310 Computation::FacetInterpolatorGray2::~FacetInterpolatorGray2( )
 311 { }
 312
 313 RefCountPtr< const Lang::Gray >
 314 Computation::FacetInterpolatorGray2::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez )
 315  const
 316 {
 317         double w2 = ( Concrete::inner( t_, point ) - m_ ) / d_;
 318         double w1 = 1 - w2;
 319
 320         // The following may not be efficient, but it was easy to code by copying code from Core_rotate3D.
 321         // I've got the feeling that this is much better done using geometric algebra based on Clifford algebra.
 322
 323         Concrete::UnitFloatTriple normal = n1_->normal( ).rotate( rotationDirection_, w2 * angle_ );
 324
 325         Concrete::Gray c1 = n1_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w1 );
 326         Concrete::Gray c2 = n2_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w2 );
 327         Concrete::Gray res = c1.addNoCheck( c2 );
 328         if( res.gr_ < 0 )
 329                 {
 330                         res.gr_ = 0;
 331                 }
 332         else if( res.gr_ > 1 )
 333                 {
 334                         res.gr_ = 1;
 335                 }
 336         return RefCountPtr< const Lang::Gray >( new Lang::Gray( res ) );
 337 }
 338
 339 RefCountPtr< const Lang::Gray >
 340 Computation::FacetInterpolatorGray2::getDebugColor( ) const
 341 {
 342         return RefCountPtr< const Lang::Gray >( new Lang::Gray( n1_->getDebugColor( ) ) );
 343 }
 344
 345 RefCountPtr< const Computation::FacetInterpolatorGray >
 346 Computation::FacetInterpolatorGray2::transformed( const Lang::Transform3D & tf ) const
 347 {
 348         return RefCountPtr< const Computation::FacetInterpolatorGray >
 349                 ( new FacetInterpolatorGray2( n1_->transformed( tf ),
 350                                                                                                                                         n2_->transformed( tf ) ) );
 351 }
 352
 353 void
 354 Computation::FacetInterpolatorGray2::gcMark( Kernel::GCMarkedSet & marked )
 355 {
 356         const_cast< Lang::FacetNormalGray * >( n1_.getPtr( ) )->gcMark( marked );
 357         const_cast< Lang::FacetNormalGray * >( n2_.getPtr( ) )->gcMark( marked );
 358 }
 359
 360
 361 Computation::FacetInterpolatorGray3::FacetInterpolatorGray3( const RefCountPtr< const Lang::FacetNormalGray > & n1,
 362                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalGray > & n2,
 363                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalGray > & n3 )
 364         : n1_( n1 ), n2_( n2 ), n3_( n3 ),
 365                 p1_( n1_->position( ) ), p2_( n2_->position( ) ), p3_( n3_->position( ) ),
 366                 d1_( 0, 0, 0, bool( ) ), d2_( 0, 0, 0, bool( ) ), d3_( 0, 0, 0, bool( ) )
 367 {
 368
 369         Concrete::Coords3D r12 = p2_ - p1_;
 370         Concrete::Coords3D r23 = p3_ - p2_;
 371         Concrete::Coords3D r31 = p1_ - p3_;
 372
 373         d1_ = ( p2_ + r23 * ( Concrete::innerScalar( r23, p1_ - p2_ ) / Concrete::innerScalar( r23, r23 ) ) - p1_ ).direction( );
 374         d2_ = ( p3_ + r31 * ( Concrete::innerScalar( r31, p2_ - p3_ ) / Concrete::innerScalar( r31, r31 ) ) - p2_ ).direction( );
 375         d3_ = ( p1_ + r12 * ( Concrete::innerScalar( r12, p3_ - p1_ ) / Concrete::innerScalar( r12, r12 ) ) - p3_ ).direction( );
 376
 377         // The lengs are obtained by projecting any of the other two points on the direction:
 378         l1_ = Concrete::inner( d1_, p2_ - p1_ );
 379         l2_ = Concrete::inner( d2_, p3_ - p2_ );
 380         l3_ = Concrete::inner( d3_, p1_ - p3_ );
 381 }
 382
 383 Computation::FacetInterpolatorGray3::~FacetInterpolatorGray3( )
 384 { }
 385
 386 RefCountPtr< const Lang::Gray >
 387 Computation::FacetInterpolatorGray3::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez )
 388  const
 389 {
 390         // #define USE_CONVEX_COMBINATION
 391 #ifdef USE_CONVEX_COMBINATION
 392         // This is another way to assign initial weights.       The symmetry is not obvious, but
 393         // follows if we note that the convex combination of the corners
 394         //       (1-(w2+w3)) * p1_ + w2 * p2_ + w3 * p3_
 395         // that result in <point> is unique inside the triangle, and then this is extended to all the plane.
 396         //
 397         Concrete::Coords3D b = point - p1_;
 398         Concrete::Coords3D a1 = p2_ - p1_;
 399         Concrete::Coords3D a2 = p3_ - p1_;
 400         const Physical< 2, 0 > a11 = Concrete::inner( a1, a1 );
 401         const Physical< 2, 0 > a12 = Concrete::inner( a2, a1 );
 402         const Physical< 2, 0 > a22 = Concrete::inner( a2, a2 );
 403         const Physical< 2, 0 > b1 = Concrete::inner( a1, b );
 404         const Physical< 2, 0 > b2 = Concrete::inner( a2, b );
 405         const Physical< 4, 0 > det = ( a11 * a22 - a12 * a12 );
 406         if( det.abs( ) < 1e-8 )
 407                 {
 408                         // The surface has no area, so the color shouldn't matter.
 409                         return Lang::THE_BLACK;
 410                 }
 411         Physical< -4, 0 > invDet = 1. / det;
 412         double w2 = invDet * (   a22 * b1 - a12 * b2 );
 413         double w3 = invDet * ( - a12 * b1 + a11 * b2 );
 414         double w1 = 1 - ( w2 + w3 );
 415         // Inside the triangle, this will do nothing by construction.   However, outside the triangle we want to avoid negative weights, maybe...
 416         // Or maybe negative weights are OK?
 417         w1 = std::max( 0., w1 );
 418         w2 = std::max( 0., w2 );
 419         w3 = std::max( 0., w3 );
 420 #else
 421         // This is one way to assign initial weights.   It is symmetric in construction.
 422         //
 423         double w1 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d1_, point - p1_ ) / l1_ ) );
 424         double w2 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d2_, point - p2_ ) / l2_ ) );
 425         double w3 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d3_, point - p3_ ) / l3_ ) );
 426 #endif
 427
 428
 429         double sumInv = 1 / ( w1 + w2 + w3 );
 430         w1 *= sumInv;
 431         w2 *= sumInv;
 432         w3 *= sumInv;
 433
 434         const Concrete::UnitFloatTriple & n1 = n1_->normal( );
 435         const Concrete::UnitFloatTriple & n2 = n2_->normal( );
 436         const Concrete::UnitFloatTriple & n3 = n3_->normal( );
 437
 438         // Here we do the ugly linear interpolation, and let the UnitFloatTriple constructor normalize the result.
 439         Concrete::UnitFloatTriple normal( w1 * n1.x_ + w2 * n2.x_ + w3 * n3.x_,
 440                                                                                                                                                 w1 * n1.y_ + w2 * n2.y_ + w3 * n3.y_,
 441                                                                                                                                                 w1 * n1.z_ + w2 * n2.z_ + w3 * n3.z_ );
 442
 443         Concrete::Gray c1 = n1_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w1 );
 444         Concrete::Gray c2 = n2_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w2 );
 445         Concrete::Gray c3 = n3_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w3 );
 446
 447         Concrete::Gray res = c1.addNoCheck( c2 ).addNoCheck( c3 );
 448         if( res.gr_ < 0 )
 449                 {
 450                         res.gr_ = 0;
 451                 }
 452         else if( res.gr_ > 1 )
 453                 {
 454                         res.gr_ = 1;
 455                 }
 456         return RefCountPtr< const Lang::Gray >( new Lang::Gray( res ) );
 457 }
 458
 459 RefCountPtr< const Lang::Gray >
 460 Computation::FacetInterpolatorGray3::getDebugColor( ) const
 461 {
 462         return RefCountPtr< const Lang::Gray >( new Lang::Gray( n1_->getDebugColor( ) ) );
 463 }
 464
 465 RefCountPtr< const Computation::FacetInterpolatorGray >
 466 Computation::FacetInterpolatorGray3::transformed( const Lang::Transform3D & tf ) const
 467 {
 468         return RefCountPtr< const Computation::FacetInterpolatorGray >
 469                 ( new FacetInterpolatorGray3( n1_->transformed( tf ),
 470                                                                                                                                         n2_->transformed( tf ),
 471                                                                                                                                         n3_->transformed( tf ) ) );
 472 }
 473
 474 void
 475 Computation::FacetInterpolatorGray3::gcMark( Kernel::GCMarkedSet & marked )
 476 {
 477         const_cast< Lang::FacetNormalGray * >( n1_.getPtr( ) )->gcMark( marked );
 478         const_cast< Lang::FacetNormalGray * >( n2_.getPtr( ) )->gcMark( marked );
 479         const_cast< Lang::FacetNormalGray * >( n3_.getPtr( ) )->gcMark( marked );
 480 }
 481
 482
 483 Computation::FacetInterpolatorRGB::~FacetInterpolatorRGB( )
 484 { }
 485
 486
 487 Computation::FacetInterpolatorRGB1::~FacetInterpolatorRGB1( )
 488 { }
 489
 490 RefCountPtr< const Lang::RGB >
 491 Computation::FacetInterpolatorRGB1::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez ) const
 492 {
 493         return RefCountPtr< const Lang::RGB >( new Lang::RGB( n1_->compute( point, tf, eyez, lights ) ) );
 494 }
 495
 496 RefCountPtr< const Lang::RGB >
 497 Computation::FacetInterpolatorRGB1::getDebugColor( ) const
 498 {
 499         return RefCountPtr< const Lang::RGB >( new Lang::RGB( n1_->getDebugColor( ) ) );
 500 }
 501
 502 RefCountPtr< const Computation::FacetInterpolatorRGB >
 503 Computation::FacetInterpolatorRGB1::transformed( const Lang::Transform3D & tf ) const
 504 {
 505         return RefCountPtr< const Computation::FacetInterpolatorRGB >
 506                 ( new FacetInterpolatorRGB1( n1_->transformed( tf ) ) );
 507 }
 508
 509
 510 void
 511 Computation::FacetInterpolatorRGB1::gcMark( Kernel::GCMarkedSet & marked )
 512 {
 513         const_cast< Lang::FacetNormalRGB * >( n1_.getPtr( ) )->gcMark( marked );
 514 }
 515
 516
 517 Computation::FacetInterpolatorRGB2::FacetInterpolatorRGB2( const RefCountPtr< const Lang::FacetNormalRGB > & n1,
 518                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalRGB > & n2 )
 519         : n1_( n1 ), n2_( n2 ),
 520                 d_( 0 ), t_( 0, 0, 0, bool( ) ), m_( 0 ),
 521                 rotationDirection_( 0, 0, 0, bool( ) )
 522 {
 523         Concrete::Coords3D r = n2_->position( ) - n1_->position( );
 524         d_ = r.norm( );
 525         t_ = r.direction( d_ );
 526         m_ = Concrete::inner( t_, n1_->position( ) );
 527
 528         angle_ = acos( Concrete::inner( n1_->normal( ), n2_->normal( ) ) );
 529         try
 530                 {
 531                         rotationDirection_ = Concrete::crossDirection( n1_->normal( ), n2_->normal( ) );
 532                 }
 533         catch( const NonLocalExit::CrossDirectionOfParallel & ball )
 534                 {
 535                         if( angle_ > 3 )
 536                                 {
 537                                         throw Exceptions::MiscellaneousRequirement( "The facet normals on a facet must not point in opposite directions." );
 538                                 }
 539                         rotationDirection_ = Concrete::UnitFloatTriple( 1, 0, 0, bool( ) );
 540                 }
 541 }
 542
 543 Computation::FacetInterpolatorRGB2::~FacetInterpolatorRGB2( )
 544 { }
 545
 546 RefCountPtr< const Lang::RGB >
 547 Computation::FacetInterpolatorRGB2::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez )
 548  const
 549 {
 550         double w2 = ( Concrete::inner( t_, point ) - m_ ) / d_;
 551         double w1 = 1 - w2;
 552
 553         // The following may not be efficient, but it was easy to code by copying code from Core_rotate3D.
 554         // I've got the feeling that this is much better done using geometric algebra based on Clifford algebra.
 555
 556         Concrete::UnitFloatTriple normal = n1_->normal( ).rotate( rotationDirection_, w2 * angle_ );
 557
 558         Concrete::RGB c1 = n1_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w1 );
 559         Concrete::RGB c2 = n2_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w2 );
 560         Concrete::RGB res = c1.addNoCheck( c2 );
 561         if( res.r_ < 0 )
 562                 {
 563                         res.r_ = 0;
 564                 }
 565         else if( res.r_ > 1 )
 566                 {
 567                         res.r_ = 1;
 568                 }
 569         if( res.g_ < 0 )
 570                 {
 571                         res.g_ = 0;
 572                 }
 573         else if( res.g_ > 1 )
 574                 {
 575                         res.g_ = 1;
 576                 }
 577         if( res.b_ < 0 )
 578                 {
 579                         res.b_ = 0;
 580                 }
 581         else if( res.b_ > 1 )
 582                 {
 583                         res.b_ = 1;
 584                 }
 585         return RefCountPtr< const Lang::RGB >( new Lang::RGB( res ) );
 586 }
 587
 588 RefCountPtr< const Lang::RGB >
 589 Computation::FacetInterpolatorRGB2::getDebugColor( ) const
 590 {
 591         return RefCountPtr< const Lang::RGB >( new Lang::RGB( n1_->getDebugColor( ) ) );
 592 }
 593
 594 RefCountPtr< const Computation::FacetInterpolatorRGB >
 595 Computation::FacetInterpolatorRGB2::transformed( const Lang::Transform3D & tf ) const
 596 {
 597         return RefCountPtr< const Computation::FacetInterpolatorRGB >
 598                 ( new FacetInterpolatorRGB2( n1_->transformed( tf ),
 599                                                                                                                                         n2_->transformed( tf ) ) );
 600 }
 601
 602 void
 603 Computation::FacetInterpolatorRGB2::gcMark( Kernel::GCMarkedSet & marked )
 604 {
 605         const_cast< Lang::FacetNormalRGB * >( n1_.getPtr( ) )->gcMark( marked );
 606         const_cast< Lang::FacetNormalRGB * >( n2_.getPtr( ) )->gcMark( marked );
 607 }
 608
 609
 610 Computation::FacetInterpolatorRGB3::FacetInterpolatorRGB3( const RefCountPtr< const Lang::FacetNormalRGB > & n1,
 611                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalRGB > & n2,
 612                                                                                                                                                                                                                                                  const RefCountPtr< const Lang::FacetNormalRGB > & n3 )
 613         : n1_( n1 ), n2_( n2 ), n3_( n3 ),
 614                 p1_( n1_->position( ) ), p2_( n2_->position( ) ), p3_( n3_->position( ) ),
 615                 d1_( 0, 0, 0, bool( ) ), d2_( 0, 0, 0, bool( ) ), d3_( 0, 0, 0, bool( ) )
 616 {
 617
 618         Concrete::Coords3D r12 = p2_ - p1_;
 619         Concrete::Coords3D r23 = p3_ - p2_;
 620         Concrete::Coords3D r31 = p1_ - p3_;
 621
 622         d1_ = ( p2_ + r23 * ( Concrete::innerScalar( r23, p1_ - p2_ ) / Concrete::innerScalar( r23, r23 ) ) - p1_ ).direction( );
 623         d2_ = ( p3_ + r31 * ( Concrete::innerScalar( r31, p2_ - p3_ ) / Concrete::innerScalar( r31, r31 ) ) - p2_ ).direction( );
 624         d3_ = ( p1_ + r12 * ( Concrete::innerScalar( r12, p3_ - p1_ ) / Concrete::innerScalar( r12, r12 ) ) - p3_ ).direction( );
 625
 626         // The lengs are obtained by projecting any of the other two points on the direction:
 627         l1_ = Concrete::inner( d1_, p2_ - p1_ );
 628         l2_ = Concrete::inner( d2_, p3_ - p2_ );
 629         l3_ = Concrete::inner( d3_, p1_ - p3_ );
 630 }
 631
 632 Computation::FacetInterpolatorRGB3::~FacetInterpolatorRGB3( )
 633 { }
 634
 635 RefCountPtr< const Lang::RGB >
 636 Computation::FacetInterpolatorRGB3::compute( const Lang::Transform3D & tf, const std::list< RefCountPtr< const Lang::LightSource > > & lights, const Concrete::Coords3D & point, const Concrete::Length eyez )
 637  const
 638 {
 639         // #define USE_CONVEX_COMBINATION
 640 #ifdef USE_CONVEX_COMBINATION
 641         // This is another way to assign initial weights.       The symmetry is not obvious, but
 642         // follows if we note that the convex combination of the corners
 643         //       (1-(w2+w3)) * p1_ + w2 * p2_ + w3 * p3_
 644         // that result in <point> is unique inside the triangle, and then this is extended to all the plane.
 645         //
 646         Concrete::Coords3D b = point - p1_;
 647         Concrete::Coords3D a1 = p2_ - p1_;
 648         Concrete::Coords3D a2 = p3_ - p1_;
 649         const Physical< 2, 0 > a11 = Concrete::inner( a1, a1 );
 650         const Physical< 2, 0 > a12 = Concrete::inner( a2, a1 );
 651         const Physical< 2, 0 > a22 = Concrete::inner( a2, a2 );
 652         const Physical< 2, 0 > b1 = Concrete::inner( a1, b );
 653         const Physical< 2, 0 > b2 = Concrete::inner( a2, b );
 654         const Physical< 4, 0 > det = ( a11 * a22 - a12 * a12 );
 655         if( det.abs( ) < 1e-8 )
 656                 {
 657                         // The surface has no area, so the color shouldn't matter.
 658                         return Lang::THE_BLACK;
 659                 }
 660         Physical< -4, 0 > invDet = 1. / det;
 661         double w2 = invDet * (   a22 * b1 - a12 * b2 );
 662         double w3 = invDet * ( - a12 * b1 + a11 * b2 );
 663         double w1 = 1 - ( w2 + w3 );
 664         // Inside the triangle, this will do nothing by construction.   However, outside the triangle we want to avoid negative weights, maybe...
 665         // Or maybe negative weights are OK?
 666         w1 = std::max( 0., w1 );
 667         w2 = std::max( 0., w2 );
 668         w3 = std::max( 0., w3 );
 669 #else
 670         // This is one way to assign initial weights.   It is symmetric in construction.
 671         //
 672         double w1 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d1_, point - p1_ ) / l1_ ) );
 673         double w2 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d2_, point - p2_ ) / l2_ ) );
 674         double w3 = std::max( 0., 1 - static_cast< double >( Concrete::inner( d3_, point - p3_ ) / l3_ ) );
 675 #endif
 676
 677
 678         double sumInv = 1 / ( w1 + w2 + w3 );
 679         w1 *= sumInv;
 680         w2 *= sumInv;
 681         w3 *= sumInv;
 682
 683         const Concrete::UnitFloatTriple & n1 = n1_->normal( );
 684         const Concrete::UnitFloatTriple & n2 = n2_->normal( );
 685         const Concrete::UnitFloatTriple & n3 = n3_->normal( );
 686
 687         // Here we do the ugly linear interpolation, and let the UnitFloatTriple constructor normalize the result.
 688         Concrete::UnitFloatTriple normal( w1 * n1.x_ + w2 * n2.x_ + w3 * n3.x_,
 689                                                                                                                                                 w1 * n1.y_ + w2 * n2.y_ + w3 * n3.y_,
 690                                                                                                                                                 w1 * n1.z_ + w2 * n2.z_ + w3 * n3.z_ );
 691
 692         Concrete::RGB c1 = n1_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w1 );
 693         Concrete::RGB c2 = n2_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w2 );
 694         Concrete::RGB c3 = n3_->compute( point, normal, tf, eyez, lights ).mulNoCheck( w3 );
 695
 696         Concrete::RGB res = c1.addNoCheck( c2 ).addNoCheck( c3 );
 697         if( res.r_ < 0 )
 698                 {
 699                         res.r_ = 0;
 700                 }
 701         else if( res.r_ > 1 )
 702                 {
 703                         res.r_ = 1;
 704                 }
 705         if( res.g_ < 0 )
 706                 {
 707                         res.g_ = 0;
 708                 }
 709         else if( res.g_ > 1 )
 710                 {
 711                         res.g_ = 1;
 712                 }
 713         if( res.b_ < 0 )
 714                 {
 715                         res.b_ = 0;
 716                 }
 717         else if( res.b_ > 1 )
 718                 {
 719                         res.b_ = 1;
 720                 }
 721         return RefCountPtr< const Lang::RGB >( new Lang::RGB( res ) );
 722 }
 723
 724 RefCountPtr< const Lang::RGB >
 725 Computation::FacetInterpolatorRGB3::getDebugColor( ) const
 726 {
 727         return RefCountPtr< const Lang::RGB >( new Lang::RGB( n1_->getDebugColor( ) ) );
 728 }
 729
 730 RefCountPtr< const Computation::FacetInterpolatorRGB >
 731 Computation::FacetInterpolatorRGB3::transformed( const Lang::Transform3D & tf ) const
 732 {
 733         return RefCountPtr< const Computation::FacetInterpolatorRGB >
 734                 ( new FacetInterpolatorRGB3( n1_->transformed( tf ),
 735                                                                                                                                         n2_->transformed( tf ),
 736                                                                                                                                         n3_->transformed( tf ) ) );
 737 }
 738
 739 void
 740 Computation::FacetInterpolatorRGB3::gcMark( Kernel::GCMarkedSet & marked )
 741 {
 742         const_cast< Lang::FacetNormalRGB * >( n1_.getPtr( ) )->gcMark( marked );
 743         const_cast< Lang::FacetNormalRGB * >( n2_.getPtr( ) )->gcMark( marked );
 744         const_cast< Lang::FacetNormalRGB * >( n3_.getPtr( ) )->gcMark( marked );
 745 }