Commit 57140131 authored by Mathias Goldau's avatar Mathias Goldau
Browse files

[ADD] New interpolation for eigen vector fields

parent d158258a
......@@ -27,6 +27,8 @@
#include <string>
#include <vector>
#include <boost/array.hpp>
#include "../common/WAssert.h"
#include "WDataSetSingle.h"
#include "WDataSetVector.h"
......@@ -78,56 +80,95 @@ boost::shared_ptr< WPrototyped > WDataSetVector::getPrototype()
return m_prototype;
}
WVector3D WDataSetVector::interpolate( const WPosition& pos, bool *success ) const
namespace
{
boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_grid );
boost::array< double, 8 > computePrefactors( const WPosition& pos, boost::shared_ptr< const WGrid > i_grid,
boost::shared_ptr< const WValueSetBase > i_valueSet, bool *success, boost::shared_ptr< std::vector< size_t > > vertexIds )
{
boost::shared_ptr< const WGridRegular3D > grid = boost::shared_dynamic_cast< const WGridRegular3D >( i_grid );
WAssert( grid, "This data set has a grid whose type is not yet supported for interpolation." );
WAssert( grid->isNotRotated(), "Only feasible for grids that are only translated or scaled so far." );
WAssert( ( m_valueSet->order() == 1 && m_valueSet->dimension() == 3 ),
"Only implemented for 3D Vectors so far." );
WAssert( grid, "This data set has a grid whose type is not yet supported for interpolation." );
WAssert( grid->isNotRotated(), "Only feasible for grids that are only translated or scaled so far." );
WAssert( ( i_valueSet->order() == 1 && i_valueSet->dimension() == 3 ),
"Only implemented for 3D Vectors so far." );
boost::array< double, 8 > h;
bool isInside = true;
size_t cellId = grid->getCellId( pos, &isInside );
bool isInside = true;
size_t cellId = grid->getCellId( pos, &isInside );
if( !isInside )
{
*success = false;
return h;
}
*success = true; // set it here, before the real work is done, because this cannot fail anymore.
*vertexIds = grid->getCellVertexIds( cellId );
WPosition localPos = pos - grid->getPosition( ( *vertexIds )[0] );
double lambdaX = localPos[0] / grid->getOffsetX();
double lambdaY = localPos[1] / grid->getOffsetY();
double lambdaZ = localPos[2] / grid->getOffsetZ();
// lZ lY
// | /
// | 6___/_7
// |/: /|
// 4_:___5 |
// | :...|.|
// |.2 | 3
// |_____|/ ____lX
// 0 1
h[0] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[1] = ( lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[2] = ( 1 - lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[3] = ( lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[4] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[5] = ( lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[6] = ( 1 - lambdaX ) * ( lambdaY ) * ( lambdaZ );
h[7] = ( lambdaX ) * ( lambdaY ) * ( lambdaZ );
return h;
}
}
WVector3D WDataSetVector::interpolate( const WPosition& pos, bool *success ) const
{
boost::shared_ptr< std::vector< size_t > > vertexIds( new std::vector< size_t > );
boost::array< double, 8 > h = computePrefactors( pos, m_grid, m_valueSet, success, vertexIds );
WVector3D result( 0.0, 0.0, 0.0 );
if( !isInside )
if( *success ) // only if pos was iniside the grid, we proivde a result different to 0.0, 0.0, 0.0
{
*success = false;
return WVector3D();
for( size_t i = 0; i < 8; ++i )
{
result += h[i] * getVectorAt( ( *vertexIds )[i] );
}
}
std::vector< size_t > vertexIds = grid->getCellVertexIds( cellId );
WPosition localPos = pos - grid->getPosition( vertexIds[0] );
double lambdaX = localPos[0] / grid->getOffsetX();
double lambdaY = localPos[1] / grid->getOffsetY();
double lambdaZ = localPos[2] / grid->getOffsetZ();
std::vector< double > h( 8 );
// lZ lY
// | /
// | 6___/_7
// |/: /|
// 4_:___5 |
// | :...|.|
// |.2 | 3
// |_____|/ ____lX
// 0 1
h[0] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[1] = ( lambdaX ) * ( 1 - lambdaY ) * ( 1 - lambdaZ );
h[2] = ( 1 - lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[3] = ( lambdaX ) * ( lambdaY ) * ( 1 - lambdaZ );
h[4] = ( 1 - lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[5] = ( lambdaX ) * ( 1 - lambdaY ) * ( lambdaZ );
h[6] = ( 1 - lambdaX ) * ( lambdaY ) * ( lambdaZ );
h[7] = ( lambdaX ) * ( lambdaY ) * ( lambdaZ );
WVector3D result( 0, 0, 0 );
for( size_t i = 0; i < 8; ++i )
return result;
}
WVector3D WDataSetVector::eigenVectorInterpolate( const WPosition& pos, bool *success ) const
{
boost::shared_ptr< std::vector< size_t > > vertexIds( new std::vector< size_t > );
boost::array< double, 8 > h = computePrefactors( pos, m_grid, m_valueSet, success, vertexIds );
WVector3D result( 0.0, 0.0, 0.0 );
if( *success ) // only if pos was iniside the grid, we proivde a result different to 0.0, 0.0, 0.0
{
result += h[i] * getVectorAt( vertexIds[i] );
for( size_t i = 0; i < 8; ++i )
{
double sign = 1.0;
if( getVectorAt( ( *vertexIds )[0] ).dotProduct( getVectorAt( ( *vertexIds )[i] ) ) < 0.0 )
{
sign = -1.0;
}
result += h[i] * sign * getVectorAt( ( *vertexIds )[i] );
}
}
*success = true;
return result;
}
......
......@@ -100,6 +100,17 @@ public:
*/
WVector3D interpolate( const WPosition &pos, bool *success ) const;
/**
* Interpolates the very same way as \ref interpolate but it assures that all vecs are aligned to point into the same
* half-space. This is useful for eigenvector fields, where -v, and v both are eigenvectors.
*
* \param pos Position to interpolate a vector for
* \param success return parameter which is true if pos was inside of the grid, false otherwise.
*
* \return The resulting interpolated vector.
*/
WVector3D eigenVectorInterpolate( const WPosition &pos, bool *success ) const;
/**
* Get the vector on the given position in value set.
* \note currently only implmented for WVector3D
......
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