Commit 4697b1d7 authored by Alexander Wiebel's avatar Alexander Wiebel
Browse files

[ADD] Trilinear interpolation in data given on regular grids without

transformation matrix

Therefore also added the follwing
 * function to make a matrix the identity matrix.
 * cell locator for regular grids without transformation matrix
 * function that gives the vertex ids of a cell
parent 463a7118
......@@ -57,6 +57,29 @@ public:
m_nbCols = newMatrix.m_nbCols;
}
/**
* Makes the matix contain the identity matrix, i.e. 1 on the diagonal.
*/
WMatrix& makeIdentity()
{
size_t nbRows = this->size() / m_nbCols;
for( size_t i = 0; i < nbRows; ++i )
{
for( size_t j = 0; j < m_nbCols; ++j )
{
if( i == j )
{
(*this)( i, j ) = 1;
}
else
{
(*this)( i, j ) = 0;
}
}
}
return *this;
}
/**
* Get number of rows.
*/
......
......@@ -23,11 +23,13 @@
//---------------------------------------------------------------------------
#include <string>
#include <vector>
#include "WDataTexture3D.h"
#include "WValueSet.h"
#include "WGrid.h"
#include "../common/WPrototyped.h"
#include "../common/WException.h"
#include "WDataSetSingle.h"
......@@ -141,3 +143,56 @@ double WDataSetSingle::getValueAt( int x, int y, int z )
return getValueAt( id );
}
double WDataSetSingle::interpolate( wmath::WPosition pos )
{
boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_grid );
// TODO(wiebel): change this to eassert.
if( !grid )
{
throw WException( "This data set has a grid whose type is not yet supported for interpolation." );
}
// TODO(wiebel): change this to eassert.
if( grid->getTransformationMatrix() != wmath::WMatrix<double>( 4, 4 ).makeIdentity() )
{
throw WException( "Only feasible for untranslated grid so far." );
}
// TODO(wiebel): change this to eassert.
if( !( m_valueSet->order() == 0 && m_valueSet->dimension() == 1 ) )
{
throw WException( "Only implemented for scalar values so far." );
}
std::vector< size_t > vertexIds = grid->getCellVertexIds( grid->getCellId( pos ) );
wmath::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 );
double result = 0;
for( size_t i = 0; i < 8; ++i )
{
result += h[i] * getValueAt( vertexIds[i] );
}
return result;
}
......@@ -101,6 +101,15 @@ public:
*/
template< typename T > T getValueAt( size_t id );
/**
* Interpolate the value fo the valueset at the given position
*
* \param pos The position for wich we would like to get a value.
*
* \return Scalar value for that given position
*/
double interpolate( wmath::WPosition pos );
/**
* Get the value stored at a certain grid position of the data set in type double.
*
......
......@@ -288,6 +288,44 @@ wmath::WValue< int > WGridRegular3D::getVoxelCoord( const wmath::WPosition& pos
return result;
}
size_t WGridRegular3D::getCellId( const wmath::WPosition& pos ) const
{
// TODO(wiebel): change this to eassert.
if( m_matrix != wmath::WMatrix<double>( 4, 4 ).makeIdentity() )
{
throw WException( "Only feasible for untranslated grid so far." );
}
wmath::WPosition posRelativeToOrigin = pos - m_origin;
size_t xCellId = floor( posRelativeToOrigin[0] / m_offsetX );
size_t yCellId = floor( posRelativeToOrigin[1] / m_offsetY );
size_t zCellId = floor( posRelativeToOrigin[2] / m_offsetZ );
return xCellId + yCellId * ( m_nbPosX - 1 ) + zCellId * ( m_nbPosX - 1 ) * ( m_nbPosY - 1 );
}
std::vector< size_t > WGridRegular3D::getCellVertexIds( size_t cellId ) const
{
std::vector< size_t > vertices( 8 );
size_t minVertexIdZ = cellId / ( ( m_nbPosX - 1 ) * ( m_nbPosY - 1 ) );
size_t remainderXY = cellId - minVertexIdZ * ( ( m_nbPosX - 1 ) * ( m_nbPosY - 1 ) );
size_t minVertexIdY = remainderXY / ( m_nbPosX - 1 );
size_t minVertexIdX = minVertexIdY % ( m_nbPosX - 1 );
size_t minVertexId = minVertexIdX + minVertexIdY * m_nbPosX + minVertexIdZ * m_nbPosX * m_nbPosY;
vertices[0] = minVertexId;
vertices[1] = minVertexId + 1;
vertices[2] = minVertexId + m_nbPosX;
vertices[3] = minVertexId + m_nbPosX +1;
vertices[4] = minVertexId + m_nbPosX * m_nbPosY;
vertices[5] = minVertexId + m_nbPosX * m_nbPosY + 1;
vertices[6] = minVertexId + m_nbPosX * m_nbPosY + m_nbPosX;
vertices[7] = minVertexId + m_nbPosX * m_nbPosY + m_nbPosX +1;
return vertices;
}
boost::shared_ptr< std::vector< wmath::WPosition > > WGridRegular3D::getVoxelVertices( const wmath::WPosition& point, const double margin ) const
{
typedef boost::shared_ptr< std::vector< wmath::WPosition > > ReturnType;
......
......@@ -312,6 +312,33 @@ public:
*/
wmath::WValue< int > getVoxelCoord( const wmath::WPosition& pos ) const;
/**
* Computes the id of the cell containing the position pos.
*
* \param pos The position selecting the cell.
*/
size_t getCellId( const wmath::WPosition& pos ) const;
/**
* Computes the ids of the vertices of a cell given by its id.
*
* \param cellId The id of the cell we want to know ther vertices of.
* \verbatim
z-axis y-axis
| /
| 6___/_7
|/: /|
4_:___5 |
| :...|.|
|.2 | 3
|_____|/ ____x-axis
0 1
\endverbatim
*
*/
std::vector< size_t > getCellVertexIds( const size_t cellId ) const;
/**
* Computes the vertices for a voxel cuboid around the given point:
*
......
......@@ -86,6 +86,34 @@ public:
TS_ASSERT_EQUALS( dataSetSingle.getGrid(), gridDummy );
TS_ASSERT_DIFFERS( dataSetSingle.getGrid(), other );
}
/**
* Test if the interpolate function works reasonable.
*/
void testInterpolate( void )
{
// create dummies, since they are needed in almost every test
boost::shared_ptr< WGrid > grid = boost::shared_ptr< WGrid >( new WGridRegular3D( 5, 3, 3, 1, 1, 1 ) );
std::vector< double > data( grid->size() );
for( size_t i = 0; i < grid->size(); ++i )
{
data[i] = i;
}
boost::shared_ptr< WValueSet< double > > valueSet = boost::shared_ptr< WValueSet< double > >( new WValueSet< double >( 0, 1, data, W_DT_DOUBLE ) );
WDataSetSingle ds( valueSet, grid );
TS_ASSERT_EQUALS( ds.interpolate( wmath::WPosition() ), data[0] );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 1, 0, 0 ) ), data[1], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 0, 1, 0 ) ), data[5], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 1, 1, 0 ) ), data[6], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 0, 0, 1 ) ), data[15], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 1, 0, 1 ) ), data[16], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 0, 1, 1 ) ), data[20], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 1, 1, 1 ) ), data[21], 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 0.3, 0.4, 0.5 ) ), 9.8, 1e-9 );
TS_ASSERT_DELTA( ds.interpolate( wmath::WPosition( 0.5, 0.5, 0.5 ) ), 10.5, 1e-9 );
}
};
#endif // WDATASETSINGLE_TEST_H
......@@ -486,6 +486,38 @@ public:
"This point: 27 is not part of this grid: nbPosX: 3 nbPosY: 3 nbPosZ: 3" );
}
/**
* Check whether we get the right Ids.
*/
void testGetCellVertexIds( void )
{
WGridRegular3D g( 5, 3, 3, 0, 0, 0, 1, 1, 1 );
std::vector< size_t > vertexIds = g.getCellVertexIds( 5 );
TS_ASSERT_EQUALS( vertexIds.size(), 8 );
TS_ASSERT_EQUALS( vertexIds[0], 6 );
TS_ASSERT_EQUALS( vertexIds[1], 7 );
TS_ASSERT_EQUALS( vertexIds[2], 11 );
TS_ASSERT_EQUALS( vertexIds[3], 12 );
TS_ASSERT_EQUALS( vertexIds[4], 21 );
TS_ASSERT_EQUALS( vertexIds[5], 22 );
TS_ASSERT_EQUALS( vertexIds[6], 26 );
TS_ASSERT_EQUALS( vertexIds[7], 27 );
}
/**
* Check whether we get the right cellId.
*/
void testGetCellId( void )
{
WGridRegular3D g( 5, 3, 3, 0, 0, 0, 1, 1, 1 );
size_t cellId = g.getCellId( wmath::WPosition( 3.3, 1.75, 0.78 ) );
TS_ASSERT_EQUALS( cellId, 7 );
// We should get an exception if the grid is not the origin, as this is not supported so far
WGridRegular3D g1( 5, 3, 3, 1.4, 2.5, 6.4, 1, 1, 1 );
TS_ASSERT_THROWS_ANYTHING( g1.getCellId( wmath::WPosition( 3.3, 1.75, 0.78 ) ) );
}
private:
double m_delta; //!< Maximum amount to values are allowed to differ.
};
......
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