Commit 64c943e3 authored by Cyril Poupon's avatar Cyril Poupon
Browse files

adding addTortuosity() to LightCurve3d

parent aaf6a693
#include <gkg-processing-container/LightCurve3d.h>
#include <gkg-processing-numericalanalysis/CubicSpline3d.h>
#include <gkg-processing-numericalanalysis/NumericalAnalysisSelector.h>
#include <gkg-processing-numericalanalysis/RandomGenerator.h>
#include <gkg-processing-transform/Rotation3dFunctions.h>
#include <gkg-processing-algobase/Math.h>
#include <gkg-core-io/TypeOf.h>
#include <gkg-core-exception/Exception.h>
#include <algorithm>
......@@ -1409,6 +1413,251 @@ void gkg::LightCurve3d< T >::reverse()
}
template < class T >
T gkg::LightCurve3d< T >::addTortuosity(
T magnitude,
T angularDispersionInDegrees,
T waveLength,
const gkg::Vector3d< T >& referenceCurveOrientation,
const gkg::BoundingBox< T >& fieldOfView,
int32_t excludedExtremityPointCount,
int32_t maximumIterationCount )
{
try
{
////////////////////////////////////////////////////////////////////////////
// sanity check
////////////////////////////////////////////////////////////////////////////
int32_t pointCount = ( int32_t )_points.size();
if ( pointCount < 3 )
{
throw std::runtime_error( "point count should be at least 3" );
}
if ( excludedExtremityPointCount < 0 )
{
throw std::runtime_error( "excluded point count must be positive" );
}
if ( 2 * excludedExtremityPointCount >= pointCount )
{
throw std::runtime_error( "excluded point count is two large with respect"
" to the curve pint count" );
}
////////////////////////////////////////////////////////////////////////////
// converting target angular dispersion in radians
////////////////////////////////////////////////////////////////////////////
T angularDispersionInRadians = angularDispersionInDegrees * M_PI / 180.0;
////////////////////////////////////////////////////////////////////////////
// pointing to the num. analysis factory and creating a random generator
////////////////////////////////////////////////////////////////////////////
gkg::NumericalAnalysisImplementationFactory*
factory = gkg::NumericalAnalysisSelector::getInstance().
getImplementationFactory();
gkg::RandomGenerator randomGenerator( gkg::RandomGenerator::Taus );
////////////////////////////////////////////////////////////////////////////
// computing the neighborhood size
////////////////////////////////////////////////////////////////////////////
T curveResolution = ( _points[ 1 ] - _points[ 0 ] ).getNorm();
int32_t neighborhoodSize = ( int32_t )( ( 2.0f * waveLength ) /
curveResolution );
int32_t iterationCount = 0;
T inducedAngularDispersionInRadians = 0.0f;
int32_t selectedPointIndex = 0;
int32_t lowerPointIndex = 0;
int32_t upperPointIndex = 0;
int32_t lowerStartingPointIndex = 0;
int32_t upperStartingPointIndex = 0;
T oldAngularDispersionInRadians = 0.0;
T newAngularDispersionInRadians = 0.0;
int32_t p = 0;
gkg::Vector3d< T > firstAxis;
gkg::Vector3d< T > secondAxis;
gkg::Vector3d< T > thirdAxis;
T displacementAlongSecondAxis = 0.0f;
T displacementAlongThirdAxis = 0.0f;
double offset = 0.0f;
bool exitingFieldOfView = false;
gkg::Vector3d< T > displacement;
int32_t neighborPointIndex = 0;
int32_t displacementIndex = 0;
while ( ( inducedAngularDispersionInRadians <
angularDispersionInRadians ) &&
( iterationCount < maximumIterationCount ) )
{
// randomly selecting a point
//// a modifier pour mettre les limites du tirage dans le constructeur
selectedPointIndex = ( int32_t )factory->getUniformRandomUInt32(
randomGenerator,
( uint32_t )( pointCount / 2 ) ) +
pointCount / 4;
// computing lower and upper bounds of point indices to be moved
lowerPointIndex = std::max( selectedPointIndex - neighborhoodSize,
excludedExtremityPointCount );
upperPointIndex = std::min( selectedPointIndex + neighborhoodSize,
pointCount - 1 -
excludedExtremityPointCount );
// computing lower and upper indices of all impacted point(s)
lowerStartingPointIndex = std::max( lowerPointIndex - 1,
excludedExtremityPointCount );
upperStartingPointIndex = upperPointIndex;
if ( upperPointIndex == pointCount - 1 - excludedExtremityPointCount )
{
-- upperStartingPointIndex;
}
// computing the old angular deviation of the selected curve portion
oldAngularDispersionInRadians = 0.0f;
for ( p = lowerStartingPointIndex; p < upperStartingPointIndex; p++ )
{
oldAngularDispersionInRadians += gkg::getLineAngles(
_points[ p + 1 ] - _points[ p ],
referenceCurveOrientation );
}
// computing the orthonormal trieder at the selected point
if ( selectedPointIndex < pointCount - 1 )
{
firstAxis = _points[ selectedPointIndex + 1 ] -
_points[ selectedPointIndex ];
firstAxis.normalize();
}
else
{
firstAxis = _points[ selectedPointIndex ] -
_points[ selectedPointIndex - 1 ];
firstAxis.normalize();
}
gkg::getRandomOrthonormalTriederFromFirstAxis( firstAxis,
secondAxis,
thirdAxis );
// computing displacement along the second and third axis
displacementAlongSecondAxis =
( T )factory->getGaussianRandomNumber( randomGenerator,
0.0,
magnitude );
displacementAlongThirdAxis =
( T )factory->getGaussianRandomNumber( randomGenerator,
0.0,
magnitude );
// computing local deformation
exitingFieldOfView = false;
std::vector< gkg::Vector3d< T > >
displacements( upperPointIndex - lowerPointIndex + 1 );
for ( neighborPointIndex = lowerPointIndex;
neighborPointIndex <= upperPointIndex;
neighborPointIndex++ )
{
const gkg::Vector3d< T >&
neighborPoint = _points[ neighborPointIndex ];
offset = ( double )std::abs( neighborPointIndex -
selectedPointIndex );
displacement = ( secondAxis * displacementAlongSecondAxis +
thirdAxis * displacementAlongThirdAxis ) *
( T )std::exp( -( offset * offset ) /
( double )neighborhoodSize );
if ( !fieldOfView.contains( neighborPoint + displacement ) )
{
exitingFieldOfView = true;
break;
}
displacements[ displacementIndex ] = displacement;
++ displacementIndex;
}
if ( !exitingFieldOfView )
{
displacementIndex = 0;
for ( neighborPointIndex = lowerPointIndex;
neighborPointIndex <= upperPointIndex;
neighborPointIndex ++ )
{
// moving each sphere
_points[ neighborPointIndex ] += displacements[ displacementIndex ];
++ displacementIndex;
}
// computing the new angular deviation of the selected curve portion
newAngularDispersionInRadians = 0.0f;
for ( p = lowerStartingPointIndex; p < upperStartingPointIndex; p++ )
{
newAngularDispersionInRadians += gkg::getLineAngles(
_points[ p + 1 ] - _points[ p ],
referenceCurveOrientation );
}
inducedAngularDispersionInRadians =
std::abs( newAngularDispersionInRadians -
oldAngularDispersionInRadians ) /
pointCount;
}
++ iterationCount;
}
return inducedAngularDispersionInRadians * 180.0 / M_PI;
}
GKG_CATCH( "template < class T > "
"T gkg::LightCurve3d< T >::addTortuosity( "
"T magnitude, "
"T angularDispersionInDegrees, "
"T waveLength, "
"const gkg::Vector3d< T >& referenceCurveOrientation, "
"const gkg::BoundingBox< T >& fieldOfView, "
"int32_t excludedExtremityPointCount, "
"int32_t maximumIterationCount )" );
}
std::ostream& std::operator<<( std::ostream& os,
const gkg::LightCurve3d< float >& thing )
{
......
......@@ -3,6 +3,7 @@
#include <gkg-processing-coordinates/Vector3d.h>
#include <gkg-processing-container/BoundingBox.h>
#include <gkg-processing-numericalanalysis/Matrix.h>
#include <vector>
#include <list>
......@@ -107,6 +108,18 @@ class LightCurve3d
void reverse();
// return the induced angular dispersion in degrees
// excludedExtremityPointCount: do not move points outside of the range of
// indices [ excludedExtremityPointCount;
// pointCount - 1 - excludedExtremityPointCount ]
T addTortuosity( T magnitude,
T angularDispersionInDegrees,
T waveLength,
const Vector3d< T >& referenceCurveOrientation,
const BoundingBox< T >& fieldOfView,
int32_t excludedExtremityPointCount = 0,
int32_t maximumIterationCount = 100000 );
protected:
std::vector< Vector3d< T > > _points;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment