src/HYDROData/HYDROData_LinearInterpolator.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HYDROData_LinearInterpolator.h"
  24 #include <cstdlib>
  25 #include <cmath>
  26
  27 HYDROData_LinearInterpolator::HYDROData_LinearInterpolator()
  28 {
  29 }
  30
  31 HYDROData_LinearInterpolator::~HYDROData_LinearInterpolator()
  32 {
  33 }
  34
  35 TCollection_AsciiString HYDROData_LinearInterpolator::GetDescription() const
  36 {
  37   TCollection_AsciiString aDescription( "Simple linear interpolator" );
  38
  39   return aDescription;
  40 }
  41
  42 void HYDROData_LinearInterpolator::Calculate()
  43 {
  44   // Reset result data
  45   ClearResults();
  46
  47   // Check input data
  48   std::vector<double> aProfile1 = GetFirstProfileCoords();
  49   std::vector<double> aProfile2 = GetSecondProfileCoords();
  50   int aNbProfilesToCompute = GetNbProfilesToCompute();
  51
  52   int aSize1 = aProfile1.size();
  53   int aSize2 = aProfile2.size();
  54
  55   div_t aDiv1 = std::div( aSize1, 3 );
  56   div_t aDiv2 = std::div( aSize2, 3 );
  57
  58   if ( aNbProfilesToCompute < 1 ||
  59        aSize1 < 2 || aDiv1.rem != 0 ||
  60        aSize2 < 2 || aDiv2.rem != 0 ) {
  61     SetErrorCode( InvalidParametersError );
  62     return;
  63   }
  64
  65   bool isSame = false;
  66   if ( aSize1 < aSize2 ) {
  67     isSame = std::equal( aProfile1.begin(), aProfile1.end(), aProfile2.begin() );
  68   } else {
  69     isSame = std::equal( aProfile2.begin(), aProfile2.end(), aProfile1.begin() );
  70   }
  71
  72   if ( isSame ) {
  73     SetErrorCode( InvalidParametersError );
  74     return;
  75   }
  76
  77   // Linear interpolation
  78   InterpolationError aStatus = OK;
  79
  80   // the first profile should have the equal or less number of points than the second profile
  81   int aNbPoints1 = aDiv1.quot;
  82   int aNbPoints2 = aDiv2.quot;
  83   if ( aNbPoints1 > aNbPoints2 ) {
  84     aProfile1.swap( aProfile2 );
  85     std::swap( aNbPoints1, aNbPoints2 );
  86   }
  87
  88   for ( int k = 0; k <= aNbProfilesToCompute - 1; k++ ) {
  89     std::vector<double> aResultProfile; ///< the k-th computed profile
  90
  91     double aRelParam = (double)( k + 1 ) / ( aNbProfilesToCompute + 1 );
  92     for ( int i = 0; i <= aNbPoints1 - 1; i++ ) {
  93       double aRel = ( aNbPoints2 - 1 ) / ( aNbPoints1 - 1 );
  94       int aPointIndex = (int) std::floor( aRel * i ); ///< point index in the second profile
  95       int anXindex1 = i * 3;
  96       int anXindex2 = aPointIndex * 3;
  97
  98       double anXi = ( aProfile1[anXindex1] * ( 1 - aRelParam ) ) + ( aProfile2[anXindex2] * aRelParam );
  99       double anYi = ( aProfile1[anXindex1 + 1] * ( 1 - aRelParam ) ) + ( aProfile2[anXindex2 + 1] * aRelParam );
 100       double aZi = ( aProfile1[anXindex1 + 2] * ( 1 - aRelParam ) ) + ( aProfile2[anXindex2 + 2] * aRelParam );
 101
 102       aResultProfile.push_back( anXi );
 103       aResultProfile.push_back( anYi );
 104       aResultProfile.push_back( aZi );
 105     }
 106
 107     InsertResultProfile( aResultProfile );
 108   }
 109
 110   SetErrorCode( aStatus );
 111 }