32 #include <boost/array.hpp>
34 #include "../exceptions/WOutOfBounds.h"
35 #include "../WAssert.h"
36 #include "../WLimits.h"
37 #include "../WStringUtils.h"
39 #include "WPolynomialEquationSolvers.h"
40 #include "linearAlgebra/WLinearAlgebra.h"
56 throw WOutOfBounds( std::string(
"There is no midpoint for an empty line." ) );
58 return line[( line.
size() - 1 ) / 2];
66 double pathLength(
const WLine& line )
70 for(
size_t i = 1; i < line.
size(); ++i )
72 len += length( line[i - 1] - line[i] );
81 if(
size() != numPoints &&
size() > 1 && numPoints > 0 )
83 const double pathL = pathLength( *
this );
84 double newSegmentLength = pathL / ( numPoints - 1 );
85 const double delta = newSegmentLength * 1.0e-10;
86 double remainingLength = 0.0;
88 for(
size_t i = 0; i < (
size() - 1 ); ++i )
90 remainingLength += length(
at( i ) -
at( i + 1 ) );
91 while( ( remainingLength > newSegmentLength ) || std::abs( remainingLength - newSegmentLength ) < delta )
93 remainingLength -= newSegmentLength;
97 WPosition newPoint =
at( i + 1 ) + remainingLength * normalize(
at( i ) -
at( i + 1 ) );
109 else if(
size() == 1 &&
size() < numPoints )
111 for(
size_t i = 0; i < numPoints; ++i )
116 if(
size() != numPoints )
156 for(
size_t i = 1; i <
size(); )
158 if( length( newLine.
back() - ( *this )[i] ) > newSegmentLength )
160 const WPosition& pred = ( *this )[i - 1];
161 if( pred == newLine.
back() )
165 newLine.
push_back( newLine.
back() + normalize( ( *
this )[i] - pred ) * newSegmentLength );
173 double alpha = ( ( *this )[i][0] - pred[0] ) * ( ( *
this )[i][0] - pred[0] ) +
174 ( ( *
this )[i][1] - pred[1] ) * ( ( *
this )[i][1] - pred[1] ) +
175 ( ( *
this )[i][2] - pred[2] ) * ( ( *
this )[i][2] - pred[2] );
177 double beta = 2.0 * ( ( *this )[i][0] - pred[0] ) * ( pred[0] - newLine.
back()[0] ) +
178 2.0 * ( ( *
this )[i][1] - pred[1] ) * ( pred[1] - newLine.
back()[1] ) +
179 2.0 * ( ( *
this )[i][2] - pred[2] ) * ( pred[2] - newLine.
back()[2] );
181 double gamma = ( pred[0] - newLine.
back()[0] ) * ( pred[0] - newLine.
back()[0] ) +
182 ( pred[1] - newLine.
back()[1] ) * ( pred[1] - newLine.
back()[1] ) +
183 ( pred[2] - newLine.
back()[2] ) * ( pred[2] - newLine.
back()[2] ) - newSegmentLength * newSegmentLength;
185 typedef std::pair< std::complex< double >, std::complex< double > > ComplexPair;
186 ComplexPair solution = solveRealQuadraticEquation( alpha, beta, gamma );
188 WAssert( std::imag( solution.first ) == 0.0,
"Invalid quadratic equation while computing resamplingBySegmentLength" );
189 WAssert( std::imag( solution.second ) == 0.0,
"Invalid quadratic equation while computing resamplingBySegmentLength" );
191 if( std::real( solution.first ) > 0.0 )
193 pointOfIntersection = pred + std::real( solution.first ) * ( ( *this )[i] - pred );
197 pointOfIntersection = pred + std::real( solution.second ) * ( ( *this )[i] - pred );
199 newLine.
push_back( pointOfIntersection );
204 if( length( newLine.
back() - ( *this )[
size() - 1] ) > newSegmentLength / 2.0 )
207 newLine.
push_back( newLine.
back() + direction * newSegmentLength );
212 int equalsDelta(
const WLine& line,
const WLine& other,
double delta )
214 size_t pts = ( std::min )( other.
size(), line.
size() );
216 bool sameLines =
true;
217 for( diffPos = 0; ( diffPos < pts ) && sameLines; ++diffPos )
219 for(
int x = 0; x < 3; ++x )
221 sameLines = sameLines && ( std::abs( line[diffPos][x] - other[diffPos][x] ) <= delta );
224 if( sameLines && ( line.
size() == other.
size() ) )
235 double maxSegmentLength(
const WLine& line )
242 for(
size_t i = 0; i < line.
size() - 1; ++i )
244 result = std::max( result, static_cast< double >( length( line[i] - line[i+1] ) ) );
251 const size_t numBasePoints = 4;
252 boost::array< WPosition, numBasePoints > m;
253 boost::array< WPosition, numBasePoints > n;
255 double distance = 0.0;
256 double inverseDistance = 0.0;
257 for(
size_t i = 0; i < numBasePoints; ++i )
259 m[i] = other.
at( ( other.
size() - 1 ) * static_cast< double >( i ) / ( numBasePoints - 1 ) );
260 n[i] =
at( (
size() - 1 ) * static_cast< double >( i ) / ( numBasePoints - 1 ) );
261 distance += length2( m[i] - n[i] );
262 inverseDistance += length2( m[i] -
at( (
size() - 1 ) * static_cast< double >( numBasePoints - 1 - i ) / ( numBasePoints - 1 ) ) );
264 distance /=
static_cast< double >( numBasePoints );
265 inverseDistance /=
static_cast< double >( numBasePoints );
267 if( inverseDistance < distance )