Commit 76bca782 authored by Andre's avatar Andre
Browse files

[ADD] a module implementing a simple anisotropic diffusion filter, needs to be...

[ADD] a module implementing a simple anisotropic diffusion filter, needs to be adapted to multi-image data
[REMOVE] berman tracking removed from toolbox
parent cd562f78
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include <algorithm>
#include <string>
#include <sstream>
#include <vector>
#include "core/kernel/WKernel.h"
#include "core/dataHandler/WDataHandler.h"
#include "core/common/WPropertyHelper.h"
#include "WMAnisotropicFiltering.xpm"
#include "WMAnisotropicFiltering.h"
// This line is needed by the module loader to actually find your module.
W_LOADABLE_MODULE( WMAnisotropicFiltering )
WMAnisotropicFiltering::WMAnisotropicFiltering():
WModule()
{
m_k = 1.0;
}
WMAnisotropicFiltering::~WMAnisotropicFiltering()
{
}
boost::shared_ptr< WModule > WMAnisotropicFiltering::factory() const
{
return boost::shared_ptr< WModule >( new WMAnisotropicFiltering() );
}
const std::string WMAnisotropicFiltering::getName() const
{
return "Anisotropic Filter";
}
const std::string WMAnisotropicFiltering::getDescription() const
{
return "Filters MRI data, preserving edges.";
}
const char** WMAnisotropicFiltering::getXPMIcon() const
{
return anisotropicFiltering_xpm;
}
void WMAnisotropicFiltering::connectors()
{
m_input = boost::shared_ptr< WModuleInputData < WDataSetSingle > >(
new WModuleInputData< WDataSetSingle >( shared_from_this(), "in", "The input dataset." ) );
addConnector( m_input );
m_output = boost::shared_ptr< WModuleOutputData < WDataSetSingle > >(
new WModuleOutputData< WDataSetSingle >( shared_from_this(), "out", "The extracted image." ) );
addConnector( m_output );
WModule::connectors();
}
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_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 );
WModule::properties();
}
void WMAnisotropicFiltering::activate()
{
WModule::activate();
}
void WMAnisotropicFiltering::moduleMain()
{
m_moduleState.setResetable( true, true );
m_moduleState.add( m_input->getDataChangedCondition() );
m_moduleState.add( m_propCondition );
ready();
while( !m_shutdownFlag() )
{
m_moduleState.wait();
if( m_shutdownFlag() )
{
break;
}
boost::shared_ptr< WDataSetSingle > newDataSet = m_input->getData();
bool dataChanged = ( m_dataSet != newDataSet );
bool dataValid = ( newDataSet );
if( dataValid )
{
if( dataChanged || m_iterations->changed() || m_Kcoeff->changed() )
{
m_dataSet = newDataSet;
WAssert( m_dataSet, "" );
WAssert( m_dataSet->getValueSet(), "" );
infoLog() << "Calculating.";
m_k = m_Kcoeff->get( true );
m_d = m_delta->get( true );
calcSmoothedImages( m_iterations->get( true ) );
infoLog() << "Calculation complete.";
}
}
}
debugLog() << "Shutting down...";
debugLog() << "Finished! Good Bye!";
}
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 );
}
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();
}
void WMAnisotropicFiltering::copyData( boost::shared_ptr< std::vector< double > >& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid )
{
for( std::size_t k = 0; k < grid->size(); ++k )
{
( *smoothed )[ k ] = m_dataSet->getValueAt( k );
}
}
void WMAnisotropicFiltering::calcDeriv( std::vector< double >& deriv, boost::shared_ptr< std::vector< double > > const& smoothed,
boost::shared_ptr< WGridRegular3D > const& grid )
{
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 ] );
}
}
}
}
void WMAnisotropicFiltering::calcCoeff( std::vector< double >& coeff, 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 ) );
}
}
}
}
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 )
{
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 ) );
}
}
}
}
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WMANISOTROPICFILTERING_H
#define WMANISOTROPICFILTERING_H
#include <string>
#include "core/kernel/WModule.h"
#include "core/kernel/WModuleInputData.h"
#include "core/kernel/WModuleOutputData.h"
#include "core/dataHandler/WDataSetSingle.h"
#include "core/dataHandler/WValueSet.h"
/**
* This module smoothes images of a dataset while preserving edges.
*
* \ingroup modules
*/
class WMAnisotropicFiltering : public WModule
{
public:
/**
* Standard constructor.
*/
WMAnisotropicFiltering();
/**
* Destructor.
*/
virtual ~WMAnisotropicFiltering();
/**
* Gives back the name of this module.
* \return the module's name.
*/
virtual const std::string getName() const;
/**
* Gives back a description of this module.
* \return description to module.
*/
virtual const std::string getDescription() const;
/**
* Due to the prototype design pattern used to build modules, this method returns a new instance of this method. NOTE: it
* should never be initialized or modified in some other way. A simple new instance is required.
*
* \return the prototype used to create every module in OpenWalnut.
*/
virtual boost::shared_ptr< WModule > factory() const;
/**
* Return an icon for this module.
*
* \return the icon
*/
virtual const char** getXPMIcon() const;
protected:
/**
* Entry point after loading the module. Runs in separate thread.
*/
virtual void moduleMain();
/**
* Initialize the connectors this module is using.
*/
virtual void connectors();
/**
* Initialize the properties for this module.
*/
virtual void properties();
/**
* Callback for m_active.
*/
virtual void activate();
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,
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 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,
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 );
/**
* An input connector that accepts multi image datasets.
*/
boost::shared_ptr< WModuleInputData< WDataSetSingle > > m_input;
/**
* An output connector for the output dataset.
*/
boost::shared_ptr< WModuleOutputData< WDataSetSingle > > m_output;
/**
* This is a pointer to the dataset the module is currently working on.
*/
boost::shared_ptr< WDataSetSingle > m_dataSet;
/**
* A condition used to notify about changes in several properties.
*/
boost::shared_ptr< WCondition > m_propCondition;
//! The edge preservation parameter used for diffusion coeff calculation.
double m_k;
//! The size of the time steps.
double m_d;
//! A property for the number of iterations.
WPropInt m_iterations;
//! A property for the edge preservation parameter.
WPropDouble m_Kcoeff;
//! A property for the time step size.
WPropDouble m_delta;
};
#endif // WMANISOTROPICFILTERING_H
/* XPM */
static const char * anisotropicFiltering_xpm[] = {
"32 32 5 1",
" c #FFFFFF",
". c #D2D2D2",
"+ c #9C9C9C",
"@ c #6C6C6C",
"# c #3C3C3C",
" ",
" .............. ",
" .............. ",
" .............. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" ++++++++++++++. ",
" +++++++++@@@@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ##############@@@@@@@@@@@ ",
" ############## "};
ADD_MODULE( HARDIToSphericalHarmonics )
ADD_MODULE( anisotropicFiltering )
ADD_MODULE( applyMask )
ADD_MODULE( atlasCreator )
ADD_MODULE( atlasSurfaces )
......@@ -34,7 +35,6 @@ ADD_MODULE( webglSupport )
IF( Teem_FOUND )
ADD_MODULE( teemGlyphs ${Teem_LIBRARIES} )
# ADD_MODULE( bermanTracking ${Teem_LIBRARIES} )
ENDIF( Teem_FOUND )
# demonstration modules
......
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