Commit e696a29c authored by reichenbach's avatar reichenbach
Browse files

[CHANGE] anisotropic filter module now filters multi-image datasets

[CHANGE] anisotropic filter module now also aborts and restarts computation when parameters change
parent 76bca782
......@@ -40,7 +40,7 @@ W_LOADABLE_MODULE( WMAnisotropicFiltering )
WMAnisotropicFiltering::WMAnisotropicFiltering():
WModule()
{
m_k = 1.0;
m_k = 1.0;
}
WMAnisotropicFiltering::~WMAnisotropicFiltering()
......@@ -84,15 +84,15 @@ void WMAnisotropicFiltering::properties()
{
m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
m_iterations = m_properties->addProperty( "#Iterations", "Number of iterations.", 10, m_propCondition );
m_iterations->setMax( 1000 );
m_iterations = m_properties->addProperty( "#Iterations", "Number of iterations.", 10, m_propCondition );
m_iterations->setMax( 1000 );
m_Kcoeff = m_properties->addProperty( "K", "The diffusion weighting coefficient. Increase this to better preserve edges.", 9.0, m_propCondition );
m_Kcoeff->setMin( 0.1 );
m_Kcoeff->setMax( 1000.0 );
m_Kcoeff = m_properties->addProperty( "K", "The diffusion weighting coefficient. Increase this to better preserve edges.", 9.0, m_propCondition );
m_Kcoeff->setMin( 0.1 );
m_Kcoeff->setMax( 1000.0 );
m_delta = m_properties->addProperty( "delta", "The time step for integration.", 0.1, m_propCondition );
m_delta->setMax( 10.0 );
m_delta = m_properties->addProperty( "delta", "The time step for integration.", 0.1, m_propCondition );
m_delta->setMax( 10.0 );
WModule::properties();
}
......@@ -131,14 +131,14 @@ void WMAnisotropicFiltering::moduleMain()
WAssert( m_dataSet, "" );
WAssert( m_dataSet->getValueSet(), "" );
infoLog() << "Calculating.";
infoLog() << "Calculating.";
m_k = m_Kcoeff->get( true );
m_d = m_delta->get( true );
m_k = m_Kcoeff->get( true );
m_d = m_delta->get( true );
calcSmoothedImages( m_iterations->get( true ) );
calcSmoothedImages( m_iterations->get( true ) );
infoLog() << "Calculation complete.";
infoLog() << "Calculation complete.";
}
}
}
......@@ -149,143 +149,182 @@ void WMAnisotropicFiltering::moduleMain()
void WMAnisotropicFiltering::calcSmoothedImages( int iterations )
{
if( iterations < 1 )
return;
boost::shared_ptr< std::vector< double > > smoothed( new std::vector< double >( m_dataSet->getGrid()->size() ) );
boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_dataSet->getGrid() );
// fill in data from dataset
copyData( smoothed, grid );
// the 3 derivatives ( actually this is the gradient )
std::vector< double > deriv( 3 * m_dataSet->getGrid()->size() );
// the diffusion coeff
std::vector< double > coeff( m_dataSet->getGrid()->size() );
for( int i = 0; i < iterations; ++i )
{
calcDeriv( deriv, smoothed, grid );
calcCoeff( coeff, deriv, grid );
diffusion( deriv, coeff, smoothed, grid );
}
// create dataset and update connector
boost::shared_ptr< WValueSet< double > > vs( new WValueSet< double >(
m_dataSet->getValueSet()->order(),
m_dataSet->getValueSet()->dimension(),
smoothed,
W_DT_DOUBLE ) );
boost::shared_ptr< WDataSetSingle > ds = m_dataSet->clone( vs );
m_output->updateData( ds );
if( iterations < 1 )
return;
std::size_t numImages = m_dataSet->getValueSet()->rawSize() / m_dataSet->getGrid()->size();
infoLog() << "Images: " << numImages;
boost::shared_ptr< std::vector< double > > smoothed( new std::vector< double >( m_dataSet->getValueSet()->rawSize() ) );
boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_dataSet->getGrid() );
// fill in data from dataset
copyData( smoothed, grid );
// the 3 derivatives ( actually this is the gradient )
std::vector< double > deriv( 3 * m_dataSet->getGrid()->size() );
// the diffusion coeff
std::vector< double > coeff( m_dataSet->getGrid()->size() );
boost::shared_ptr< WProgress > prog( new WProgress( "Smoothing images", numImages ) );
m_progress->addSubProgress( prog );
for( std::size_t k = 0; k < numImages; ++k )
{
for( int i = 0; i < iterations; ++i )
{
calcDeriv( deriv, smoothed, grid, k, numImages );
if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
{
prog->finish();
return;
}
calcCoeff( coeff, deriv, grid );
if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
{
prog->finish();
return;
}
diffusion( deriv, coeff, smoothed, grid, k, numImages );
if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
{
prog->finish();
return;
}
}
++*prog;
}
prog->finish();
// create dataset and update connector
boost::shared_ptr< WValueSet< double > > vs( new WValueSet< double >(
m_dataSet->getValueSet()->order(),
m_dataSet->getValueSet()->dimension(),
smoothed,
W_DT_DOUBLE ) );
boost::shared_ptr< WDataSetSingle > ds = m_dataSet->clone( vs );
m_output->updateData( ds );
}
std::size_t WMAnisotropicFiltering::coordsToIndex( boost::shared_ptr< WGridRegular3D > const& grid, std::size_t x, std::size_t y, std::size_t z )
std::size_t WMAnisotropicFiltering::coordsToIndex( boost::shared_ptr< WGridRegular3D > const& grid,
std::size_t x, std::size_t y, std::size_t z )
{
return x + y * grid->getNbCoordsX() + z * grid->getNbCoordsX() * grid->getNbCoordsY();
return x + y * grid->getNbCoordsX() + z * grid->getNbCoordsX() * grid->getNbCoordsY();
}
void WMAnisotropicFiltering::copyData( boost::shared_ptr< std::vector< double > >& smoothed,
void WMAnisotropicFiltering::copyData( boost::shared_ptr< std::vector< double > >& smoothed, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid )
{
for( std::size_t k = 0; k < grid->size(); ++k )
{
( *smoothed )[ k ] = m_dataSet->getValueAt( k );
}
for( std::size_t k = 0; k < m_dataSet->getValueSet()->rawSize(); ++k )
{
( *smoothed )[ k ] = m_dataSet->getValueSet()->getScalarDouble( k );
}
}
void WMAnisotropicFiltering::calcDeriv( std::vector< double >& deriv, boost::shared_ptr< std::vector< double > > const& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid )
void WMAnisotropicFiltering::calcDeriv( std::vector< double >& deriv, // NOLINT non-const ref
boost::shared_ptr< std::vector< double > > const& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages )
{
std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
std::size_t d[] = { grid->getOffsetX(), grid->getOffsetY(), grid->getOffsetZ() };
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ] = ( ( *smoothed )[ coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) ]
- ( *smoothed )[ coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) ] )
/ ( 2.0 * d[ 0 ] );
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ] = ( ( *smoothed )[ coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) ]
- ( *smoothed )[ coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) ] )
/ ( 2.0 * d[ 1 ] );
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ] = ( ( *smoothed )[ coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) ]
- ( *smoothed )[ coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) ] )
/ ( 2.0 * d[ 2 ] );
}
}
}
std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
std::size_t d[] = { fabs( grid->getOffsetX() ), fabs( grid->getOffsetY() ), fabs( grid->getOffsetZ() ) };
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
double const& x1 = ( *smoothed )[ numImages * coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) + image ];
double const& x2 = ( *smoothed )[ numImages * coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) + image ];
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ] = ( x1 - x2 ) / ( 2.0 * d[ 0 ] );
double const& y1 = ( *smoothed )[ numImages * coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) + image ];
double const& y2 = ( *smoothed )[ numImages * coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) + image ];
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ] = ( y1 - y2 ) / ( 2.0 * d[ 1 ] );
double const& z1 = ( *smoothed )[ numImages * coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) + image ];
double const& z2 = ( *smoothed )[ numImages * coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) + image ];
deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ] = ( z1 - z2 ) / ( 2.0 * d[ 2 ] );
}
}
}
}
void WMAnisotropicFiltering::calcCoeff( std::vector< double >& coeff, std::vector< double > const& deriv,
void WMAnisotropicFiltering::calcCoeff( std::vector< double >& coeff, // NOLINT non-const ref
std::vector< double > const& deriv,
boost::shared_ptr< WGridRegular3D > const& grid )
{
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
// coeff = exp( -sqr( |I|/K ) )
double gradIAbsSquared = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ] * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ]
+ deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ] * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ]
+ deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ] * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
coeff[ coordsToIndex( grid, x, y, z ) ] = 1.0 / exp( gradIAbsSquared / ( m_k * m_k ) );
}
}
}
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
// coeff = exp( -sqr( |I|/K ) )
double gradIAbsSquared = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ]
* deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ]
+ deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ]
* deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ]
+ deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ]
* deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
coeff[ coordsToIndex( grid, x, y, z ) ] = 1.0 / exp( gradIAbsSquared / ( m_k * m_k ) );
}
}
}
}
void WMAnisotropicFiltering::diffusion( std::vector< double > const& deriv, std::vector< double > const& coeff,
boost::shared_ptr< std::vector< double > >& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid )
boost::shared_ptr< std::vector< double > >& smoothed, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages )
{
std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
std::size_t d[] = { grid->getOffsetX(), grid->getOffsetY(), grid->getOffsetZ() };
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
// first deriv of the image intensity
double const& dIx = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ];
double const& dIy = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ];
double const& dIz = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
// first deriv of the diffusion coeff
double dcx = ( coeff[ coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) ]
- coeff[ coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) ] )
/ ( 2.0 * d[ 0 ] );
double dcy = ( coeff[ coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) ]
- coeff[ coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) ] )
/ ( 2.0 * d[ 1 ] );
double dcz = ( coeff[ coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) ]
- coeff[ coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) ] )
/ ( 2.0 * d[ 2 ] );
// diffusion coeff
double const& c = coeff[ coordsToIndex( grid, x, y, z ) ];
// 2nd derivations in x, y, and z of the image intensity
double ddIxx = ( deriv[ 3 * coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) + 0 ]
- deriv[ 3 * coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) + 0 ] )
/ ( 2 * d[ 0 ] );
double ddIyy = ( deriv[ 3 * coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) + 1 ]
- deriv[ 3 * coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) + 1 ] )
/ ( 2 * d[ 1 ] );
double ddIzz = ( deriv[ 3 * coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) + 2 ]
- deriv[ 3 * coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) + 2 ] )
/ ( 2 * d[ 2 ] );
// d * ( grad I * grad c + c * grad grad I )
( *smoothed )[ coordsToIndex( grid, x, y, z ) ] += m_d * ( dIx * dcx + dIy * dcy + dIz * dcz + c * ( ddIxx + ddIyy + ddIzz ) );
}
}
}
std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
std::size_t d[] = { fabs( grid->getOffsetX() ), fabs( grid->getOffsetY() ), fabs( grid->getOffsetZ() ) };
for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
{
for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
{
for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
{
// first deriv of the image intensity
double const& dIx = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ];
double const& dIy = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ];
double const& dIz = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
// first deriv of the diffusion coeff
double dcx = ( coeff[ coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) ]
- coeff[ coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) ] )
/ ( 2.0 * d[ 0 ] );
double dcy = ( coeff[ coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) ]
- coeff[ coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) ] )
/ ( 2.0 * d[ 1 ] );
double dcz = ( coeff[ coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) ]
- coeff[ coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) ] )
/ ( 2.0 * d[ 2 ] );
// diffusion coeff
double const& c = coeff[ coordsToIndex( grid, x, y, z ) ];
// 2nd derivations in x, y, and z of the image intensity
double ddIxx = ( deriv[ 3 * coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) + 0 ]
- deriv[ 3 * coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) + 0 ] )
/ ( 2 * d[ 0 ] );
double ddIyy = ( deriv[ 3 * coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) + 1 ]
- deriv[ 3 * coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) + 1 ] )
/ ( 2 * d[ 1 ] );
double ddIzz = ( deriv[ 3 * coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) + 2 ]
- deriv[ 3 * coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) + 2 ] )
/ ( 2 * d[ 2 ] );
// d * ( grad I * grad c + c * grad grad I )
( *smoothed )[ numImages * coordsToIndex( grid, x, y, z ) + image ] += m_d * ( dIx * dcx + dIy * dcy + dIz * dcz
+ c * ( ddIxx + ddIyy + ddIzz ) );
}
}
}
}
......@@ -26,6 +26,7 @@
#define WMANISOTROPICFILTERING_H
#include <string>
#include <vector>
#include "core/kernel/WModule.h"
#include "core/kernel/WModuleInputData.h"
......@@ -103,67 +104,71 @@ protected:
private:
/**
* Calculates the resulting smoothed image.
*
* \param iterations The number of diffusion time steps.
*/
void calcSmoothedImages( int iterations );
/**
* Calculates grid indices from voxel coords.
*
* \param grid The grid.
* \param x The x-coord.
* \param y The y-coord.
* \param z The z-coord.
*
* \return The voxel index.
*/
std::size_t coordsToIndex( boost::shared_ptr< WGridRegular3D > const& grid,
std::size_t x, std::size_t y, std::size_t z );
/**
* Copy the datasets image data to a temp array.
*
* \param smoothed The temp array to copy to.
* \param grid The grid.
*/
void copyData( boost::shared_ptr< std::vector< double > >& smoothed,
/**
* Calculates the resulting smoothed image.
*
* \param iterations The number of diffusion time steps.
*/
void calcSmoothedImages( int iterations );
/**
* Calculates grid indices from voxel coords.
*
* \param grid The grid.
* \param x The x-coord.
* \param y The y-coord.
* \param z The z-coord.
*
* \return The voxel index.
*/
std::size_t coordsToIndex( boost::shared_ptr< WGridRegular3D > const& grid,
std::size_t x, std::size_t y, std::size_t z );
/**
* Copy the datasets image data to a temp array.
*
* \param smoothed The temp array to copy to.
* \param grid The grid.
*/
void copyData( boost::shared_ptr< std::vector< double > >& smoothed, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid );
/**
* Calculates an array containing the derivations in x, y and z directions of the image
* intensity (i.e. contains the intensity gradient).
*
* \param deriv The memory used for the gradient values.
* \param smoothed The intensity data.
* \param grid The grid.
*/
void calcDeriv( std::vector< double >& deriv, boost::shared_ptr< std::vector< double > > const& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid );
/**
* Calculates an array containing the derivations in x, y and z directions of the image
* intensity (i.e. contains the intensity gradient).
*
* \param deriv The memory used for the gradient values.
* \param smoothed The intensity data.
* \param grid The grid.
* \param image The index of the image to calc the gradient from.
* \param numImages The number of images in this dataset.
*/
void calcDeriv( std::vector< double >& deriv, boost::shared_ptr< std::vector< double > > const& smoothed, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages );
/**
* Calculates the diffusion coeff for every voxel.
*
* \param coeff The memory used for the coeff data.
* \param deriv The gradient data.
* \param grid The grid.
*/
void calcCoeff( std::vector< double >& coeff, std::vector< double > const& deriv,
/**
* Calculates the diffusion coeff for every voxel.
*
* \param coeff The memory used for the coeff data.
* \param deriv The gradient data.
* \param grid The grid.
*/
void calcCoeff( std::vector< double >& coeff, std::vector< double > const& deriv, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid );
/**
* Do the diffusion.
*
* \param deriv The intensity gradient data.
* \param coeff The diffusion coeffs.
* \param smoothed The new smoothed image.
* \param grid The grid.
*/
void diffusion( std::vector< double > const& deriv, std::vector< double > const& coeff,
boost::shared_ptr< std::vector< double > >& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid );
/**
* Do the diffusion.
*
* \param deriv The intensity gradient data.
* \param coeff The diffusion coeffs.
* \param smoothed The new smoothed image.
* \param grid The grid.
* \param image The index of the image to calc the gradient from.
* \param numImages The number of images in this dataset.
*/
void diffusion( std::vector< double > const& deriv, std::vector< double > const& coeff,
boost::shared_ptr< std::vector< double > >& smoothed, // NOLINT non-const ref
boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages );
/**
* An input connector that accepts multi image datasets.
......
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