diff --git a/Applications/Diffusion/CMakeLists.txt b/Applications/Diffusion/CMakeLists.txt
index a07c306ded..4aad065532 100644
--- a/Applications/Diffusion/CMakeLists.txt
+++ b/Applications/Diffusion/CMakeLists.txt
@@ -1,80 +1,81 @@
project(MitkDiffusion)
set(DIFFUSIONAPP_NAME MitkDiffusion)
set(_app_options)
if(MITK_SHOW_CONSOLE_WINDOW)
list(APPEND _app_options SHOW_CONSOLE)
endif()
# Create a cache entry for the provisioning file which is used to export
# the file name in the MITKConfig.cmake file. This will keep external projects
# which rely on this file happy.
set(DIFFUSIONIMAGINGAPP_PROVISIONING_FILE "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${DIFFUSIONAPP_NAME}.provisioning" CACHE INTERNAL "${DIFFUSIONAPP_NAME} provisioning file" FORCE)
set(_plugins
org.commontk.configadmin
org.commontk.eventadmin
org.blueberry.osgi
org.blueberry.compat
org.blueberry.core.runtime
org.blueberry.core.expressions
org.blueberry.solstice.common
org.blueberry.core.commands
org.blueberry.ui
org.blueberry.ui.qt
org.blueberry.ui.qt.log
org.blueberry.ui.qt.help
org.mitk.core.services
org.mitk.gui.common
org.mitk.planarfigure
org.mitk.core.ext
org.mitk.diffusionimaging
org.mitk.gui.qt.application
org.mitk.gui.qt.ext
org.mitk.gui.qt.diffusionimagingapp
org.mitk.gui.qt.common
org.mitk.gui.qt.stdmultiwidgeteditor
org.mitk.gui.qt.common.legacy
org.mitk.gui.qt.datamanager
org.mitk.gui.qt.measurementtoolbox
org.mitk.gui.qt.segmentation
org.mitk.gui.qt.volumevisualization
org.mitk.gui.qt.diffusionimaging
org.mitk.gui.qt.imagenavigator
org.mitk.gui.qt.moviemaker
org.mitk.gui.qt.basicimageprocessing
org.mitk.gui.qt.registration
org.mitk.gui.qt.properties
+ org.mitk.gui.qt.candystore
)
# Plug-ins listed below will not be
# - added as a build-time dependency to the executable
# - listed in the provisioning file for the executable
# - installed if they are external plug-ins
set(_exclude_plugins
org.blueberry.test
org.blueberry.uitest
org.mitk.gui.qt.coreapplication
org.mitk.gui.qt.extapplication
)
FunctionCreateBlueBerryApplication(
NAME ${DIFFUSIONAPP_NAME}
DESCRIPTION "MITK Diffusion"
PLUGINS ${_plugins}
EXCLUDE_PLUGINS ${_exclude_plugins}
${_app_options}
)
mitk_use_modules(TARGET ${DIFFUSIONAPP_NAME} MODULES qtsingleapplication)
# Add meta dependencies (e.g. on auto-load modules from depending modules)
if(ALL_META_DEPENDENCIES)
add_dependencies(${DIFFUSIONAPP_NAME} ${ALL_META_DEPENDENCIES})
endif()
# Add a build time dependency to legacy BlueBerry bundles.
if(MITK_MODULES_ENABLED_PLUGINS)
add_dependencies(${DIFFUSIONAPP_NAME} ${MITK_MODULES_ENABLED_PLUGINS})
endif()
diff --git a/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user b/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user
index d698c9843a..5280950531 100644
--- a/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user
+++ b/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user
@@ -1,63 +1,63 @@
-
+
ProjectExplorer.Project.ActiveTarget
-1
ProjectExplorer.Project.EditorSettings
true
false
true
Cpp
CppGlobal
QmlJS
QmlJSGlobal
2
UTF-8
false
4
false
true
1
true
0
true
0
8
true
1
true
true
true
true
ProjectExplorer.Project.PluginSettings
ProjectExplorer.Project.TargetCount
0
ProjectExplorer.Project.Updater.EnvironmentId
{3dea2ab9-fb51-4719-a902-a8416c2a1947}
ProjectExplorer.Project.Updater.FileVersion
15
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
index 11afc7b593..cdb43eec2e 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
@@ -1,99 +1,100 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
/*===================================================================
This file is based heavily on a corresponding ITK filter.
===================================================================*/
#ifndef __itkDwiNormilzationFilter_h_
#define __itkDwiNormilzationFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
#include
#include
#include
#include
namespace itk{
/** \class DwiNormilzationFilter
- * \brief Max-Normalizes the data vectors either using the global baseline maximum or the voxelwise baseline value.
+ * \brief Normalizes the data vectors either using the specified reference value or the voxelwise baseline value.
*/
template< class TInPixelType >
class DwiNormilzationFilter :
public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > >
{
public:
typedef DwiNormilzationFilter Self;
typedef SmartPointer Pointer;
typedef SmartPointer ConstPointer;
typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(DwiNormilzationFilter, ImageToImageFilter)
typedef itk::Image< double, 3 > DoubleImageType;
typedef itk::Image< unsigned char, 3 > UcharImageType;
typedef itk::Image< unsigned short, 3 > BinImageType;
typedef itk::Image< TInPixelType, 3 > TInPixelImageType;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
typedef mitk::DiffusionImage< short >::GradientDirectionType GradientDirectionType;
typedef mitk::DiffusionImage< short >::GradientDirectionContainerType::Pointer GradientContainerType;
typedef itk::LabelStatisticsImageFilter< TInPixelImageType,BinImageType > StatisticsFilterType;
typedef itk::Statistics::Histogram< typename TInPixelImageType::PixelType > HistogramType;
typedef typename HistogramType::MeasurementType MeasurementType;
typedef itk::ShiftScaleImageFilter ShiftScaleImageFilterType;
itkSetMacro( GradientDirections, GradientContainerType )
- itkSetMacro( NewMax, TInPixelType )
- itkSetMacro( UseGlobalMax, bool )
+ itkSetMacro( ScalingFactor, TInPixelType )
+ itkSetMacro( UseGlobalReference, bool )
itkSetMacro( MaskImage, UcharImageType::Pointer )
+ itkSetMacro( Reference, double )
protected:
DwiNormilzationFilter();
~DwiNormilzationFilter() {}
void PrintSelf(std::ostream& os, Indent indent) const;
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType);
UcharImageType::Pointer m_MaskImage;
GradientContainerType m_GradientDirections;
int m_B0Index;
- TInPixelType m_NewMax;
- bool m_UseGlobalMax;
- double m_GlobalMax;
+ TInPixelType m_ScalingFactor;
+ bool m_UseGlobalReference;
+ double m_Reference;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkDwiNormilzationFilter.txx"
#endif
#endif //__itkDwiNormilzationFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.txx
index ff11de92ab..59faf8479b 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.txx
@@ -1,278 +1,127 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef __itkDwiNormilzationFilter_txx
#define __itkDwiNormilzationFilter_txx
#include
#include
#include
#define _USE_MATH_DEFINES
#include
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkImageRegionIterator.h"
#include
namespace itk {
template< class TInPixelType >
DwiNormilzationFilter< TInPixelType>::DwiNormilzationFilter()
:m_B0Index(-1)
- ,m_NewMax(1000)
- ,m_UseGlobalMax(false)
+ ,m_ScalingFactor(1000)
+ ,m_UseGlobalReference(false)
{
this->SetNumberOfRequiredInputs( 1 );
}
template< class TInPixelType >
void DwiNormilzationFilter< TInPixelType>::BeforeThreadedGenerateData()
{
typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
m_B0Index = -1;
for (unsigned int i=0; iGetVectorLength(); i++)
{
GradientDirectionType g = m_GradientDirections->GetElement(i);
if (g.magnitude()<0.001)
m_B0Index = i;
}
if (m_B0Index==-1)
itkExceptionMacro(<< "DwiNormilzationFilter: No b-Zero indecies found!");
-// // cleanup
-// itk::ImageRegionIterator inIt(inputImagePointer, inputImagePointer->GetLargestPossibleRegion());
-// while ( !inIt.IsAtEnd() )
-// {
-// typename InputImageType::PixelType pix = inIt.Get();
-// for (unsigned int i=0; iGetVectorLength(); i++)
-// {
-// if (pix[i]<0)
-// pix[i] = 0;
-// else if(pix[i]>pix[m_B0Index])
-// pix[i] = pix[m_B0Index];
-// }
-// inIt.Set(pix);
-// ++inIt;
-// }
-
-// if (m_MaskImage.IsNull())
-// {
-// m_MaskImage = UcharImageType::New();
-// m_MaskImage->SetRegions(inputImagePointer->GetLargestPossibleRegion());
-// m_MaskImage->SetOrigin(inputImagePointer->GetOrigin());
-// m_MaskImage->SetSpacing(inputImagePointer->GetSpacing());
-// m_MaskImage->SetDirection(inputImagePointer->GetDirection());
-// m_MaskImage->Allocate();
-// m_MaskImage->FillBuffer(1);
-// }
-
-// BinImageType::Pointer outputBinMaskImage = BinImageType::New();
-// outputBinMaskImage->SetRegions(m_MaskImage->GetLargestPossibleRegion());
-// outputBinMaskImage->SetOrigin(m_MaskImage->GetOrigin());
-// outputBinMaskImage->SetSpacing(m_MaskImage->GetSpacing());
-// outputBinMaskImage->SetDirection(m_MaskImage->GetDirection());
-// outputBinMaskImage->Allocate();
-
-// typename TInPixelImageType::Pointer b0Image = TInPixelImageType::New();
-// b0Image->SetRegions(m_MaskImage->GetLargestPossibleRegion());
-// b0Image->SetOrigin(m_MaskImage->GetOrigin());
-// b0Image->SetSpacing(m_MaskImage->GetSpacing());
-// b0Image->SetDirection(m_MaskImage->GetDirection());
-// b0Image->Allocate();
-// b0Image->FillBuffer(0);
-
-// itk::ImageRegionIterator it(inputImagePointer, inputImagePointer->GetLargestPossibleRegion());
-// while( !it.IsAtEnd() )
-// {
-// if (it.Get()[m_B0Index]>0)
-// b0Image->SetPixel(it.GetIndex(), it.Get()[m_B0Index]);
-// if( m_MaskImage->GetPixel(it.GetIndex()) > 0)
-// outputBinMaskImage->SetPixel(it.GetIndex(), 1);
-// else
-// outputBinMaskImage->SetPixel(it.GetIndex(), 0);
-
-// ++it;
-// }
-
-// typename StatisticsFilterType::Pointer statistics = StatisticsFilterType::New();
-// statistics->SetLabelInput( outputBinMaskImage );
-// statistics->SetInput( b0Image );
-// statistics->Update();
-
-// typename HistogramType::Pointer histogram = HistogramType::New();
-// typename HistogramType::SizeType sz;
-// typename HistogramType::MeasurementVectorType lowerBound;
-// typename HistogramType::MeasurementVectorType upperBound;
-// sz.SetSize(1);
-// lowerBound.SetSize(1);
-// upperBound.SetSize(1);
-// histogram->SetMeasurementVectorSize(1);
-
-// // Numver of Bins
-// sz[0] = 256;
-
-// lowerBound.Fill(statistics->GetMinimum(1));
-// upperBound.Fill(statistics->GetMaximum(1));
-
-// histogram->Initialize(sz, lowerBound, upperBound);
-// histogram->SetToZero();
-
-// typename HistogramType::MeasurementVectorType measurement;
-// measurement.SetSize(1); measurement[0] = 0;
-
-// itk::ImageRegionIterator b0It(b0Image, b0Image->GetLargestPossibleRegion());
-// while ( !b0It.IsAtEnd() )
-// {
-// if (m_MaskImage->GetPixel(b0It.GetIndex()) > 0)
-// {
-// TInPixelType value = b0It.Value();
-// measurement[0] = value;
-// histogram->IncreaseFrequencyOfMeasurement(measurement, 1);
-// }
-// ++b0It;
-// }
-
-// // Find bin with max frequency
-// double maxFreq =0 ;
-// unsigned int index = 0;
-// for (unsigned int ii =0 ; ii < sz[0];++ii)
-// {
-// if( maxFreq < histogram->GetFrequency(ii))
-// {
-// maxFreq = histogram->GetFrequency(ii);
-// index = ii;
-// }
-// }
-
-// typename StatisticsFilterType::Pointer statisticsImageFilter = StatisticsFilterType::New();
-// statisticsImageFilter->SetInput(b0Image);
-// statisticsImageFilter->SetLabelInput(outputBinMaskImage);
-// statisticsImageFilter->Update();
-// double variance = statisticsImageFilter->GetVariance(1);
-
-// typename ShiftScaleImageFilterType::Pointer shiftScaleImageFilter = ShiftScaleImageFilterType::New();
-// MITK_INFO << "******************* " << histogram->GetBinMax(0,index);
-// m_GlobalMax = histogram->GetBinMax(0,index);
-// m_UseGlobalMax = true;
-// MITK_INFO << "******************* " << maxFreq;
-// shiftScaleImageFilter->SetShift( 0 );//-histogram->GetBinMax(0,index));
-// shiftScaleImageFilter->SetScale( ( (double) 1 ) / std::sqrt(variance));
-
-// for (unsigned int i=0; iGetVectorLength(); i++)
-// {
-// typename TInPixelImageType::Pointer channel = TInPixelImageType::New();
-// channel->SetSpacing( this->GetInput()->GetSpacing() );
-// channel->SetOrigin( this->GetInput()->GetOrigin() );
-// channel->SetDirection( this->GetInput()->GetDirection() );
-// channel->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion() );
-// channel->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
-// channel->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
-// channel->Allocate();
-
-// ImageRegionIterator it(channel, channel->GetLargestPossibleRegion());
-// while(!it.IsAtEnd())
-// {
-// typename InputImageType::PixelType pix = inputImagePointer->GetPixel(it.GetIndex());
-// it.Set(pix[i]);
-// ++it;
-// }
-
-// shiftScaleImageFilter->SetInput( channel );
-// shiftScaleImageFilter->Update();
-// DoubleImageType::Pointer normalized = shiftScaleImageFilter->GetOutput();
-
-// ImageRegionIterator it2(normalized, normalized->GetLargestPossibleRegion());
-// while(!it2.IsAtEnd())
-// {
-// typename InputImageType::PixelType pix = inputImagePointer->GetPixel(it2.GetIndex());
-// pix[i] = m_NewMax*it2.Get();
-// inputImagePointer->SetPixel(it2.GetIndex(), pix);
-// ++it2;
-// }
-// }
-
+ if (m_UseGlobalReference)
+ MITK_INFO << "Using global reference value: " << m_Reference;
}
template< class TInPixelType >
void DwiNormilzationFilter< TInPixelType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType )
{
- MITK_INFO << "8";
typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
oit.GoToBegin();
typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
typename OutputImageType::PixelType nullPix;
nullPix.SetSize(inputImagePointer->GetVectorLength());
nullPix.Fill(0);
InputIteratorType git( inputImagePointer, outputRegionForThread );
git.GoToBegin();
while( !git.IsAtEnd() )
{
typename InputImageType::PixelType pix = git.Get();
typename OutputImageType::PixelType outPix;
outPix.SetSize(inputImagePointer->GetVectorLength());
double S0 = pix[m_B0Index];
- if (m_UseGlobalMax)
- S0 = m_GlobalMax;
+ if (m_UseGlobalReference)
+ S0 = m_Reference;
- if (S0>0.1)
+ if (S0>0.1 && pix[m_B0Index]>0)
{
for (unsigned int i=0; iGetVectorLength(); i++)
{
double val = (double)pix[i];
if (val!=val || val<0)
val = 0;
else
{
val /= S0;
- val *= (double)m_NewMax;
+ val *= (double)m_ScalingFactor;
}
outPix[i] = (TInPixelType)val;
}
oit.Set(outPix);
}
else
oit.Set(nullPix);
++oit;
++git;
}
std::cout << "One Thread finished calculation" << std::endl;
}
template< class TInPixelType >
void
DwiNormilzationFilter< TInPixelType>
::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os,indent);
}
}
#endif // __itkDwiNormilzationFilter_txx
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h
new file mode 100644
index 0000000000..5db210b559
--- /dev/null
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h
@@ -0,0 +1,78 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+/*===================================================================
+
+This file is based heavily on a corresponding ITK filter.
+
+===================================================================*/
+#ifndef __itkExtractDwiChannelFilter_h_
+#define __itkExtractDwiChannelFilter_h_
+
+#include "itkImageToImageFilter.h"
+#include "itkVectorImage.h"
+#include
+#include
+#include
+
+namespace itk{
+/** \class ExtractDwiChannelFilter
+ * \brief Remove spcified channels from diffusion-weighted image.
+ */
+
+template< class TInPixelType >
+class ExtractDwiChannelFilter :
+ public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > >
+{
+
+public:
+
+ typedef ExtractDwiChannelFilter Self;
+ typedef SmartPointer Pointer;
+ typedef SmartPointer ConstPointer;
+ typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > > Superclass;
+
+ /** Method for creation through the object factory. */
+ itkFactorylessNewMacro(Self)
+ itkCloneMacro(Self)
+
+ /** Runtime information support. */
+ itkTypeMacro(ExtractDwiChannelFilter, ImageToImageFilter)
+
+ typedef typename Superclass::InputImageType InputImageType;
+ typedef typename Superclass::OutputImageType OutputImageType;
+ typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
+
+ itkSetMacro( ChannelIndex, unsigned int )
+
+ protected:
+ ExtractDwiChannelFilter();
+ ~ExtractDwiChannelFilter() {}
+
+ void BeforeThreadedGenerateData();
+ void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType );
+
+ unsigned int m_ChannelIndex;
+};
+
+}
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkExtractDwiChannelFilter.txx"
+#endif
+
+#endif //__itkExtractDwiChannelFilter_h_
+
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.txx
new file mode 100644
index 0000000000..b6a38a6076
--- /dev/null
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.txx
@@ -0,0 +1,73 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef __itkExtractDwiChannelFilter_txx
+#define __itkExtractDwiChannelFilter_txx
+
+#include
+#include
+#include
+
+#define _USE_MATH_DEFINES
+#include
+
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+
+namespace itk {
+
+
+template< class TInPixelType >
+ExtractDwiChannelFilter< TInPixelType>::ExtractDwiChannelFilter()
+{
+ this->SetNumberOfRequiredInputs( 1 );
+}
+
+template< class TInPixelType >
+void ExtractDwiChannelFilter< TInPixelType>::BeforeThreadedGenerateData()
+{
+ typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
+ if ( inputImagePointer->GetVectorLength()<=m_ChannelIndex )
+ itkExceptionMacro("Index out of bounds!");
+}
+
+template< class TInPixelType >
+void ExtractDwiChannelFilter< TInPixelType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType )
+{
+ typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
+
+ ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
+ oit.GoToBegin();
+
+ typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
+ typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
+
+ InputIteratorType git( inputImagePointer, outputRegionForThread );
+ git.GoToBegin();
+ while( !git.IsAtEnd() )
+ {
+ oit.Set( git.Get()[m_ChannelIndex] );
+ ++oit;
+ ++git;
+ }
+
+ std::cout << "One Thread finished calculation" << std::endl;
+}
+
+}
+
+#endif // __itkExtractDwiChannelFilter_txx
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h
new file mode 100644
index 0000000000..7e18435d36
--- /dev/null
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h
@@ -0,0 +1,85 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+/*===================================================================
+
+This file is based heavily on a corresponding ITK filter.
+
+===================================================================*/
+#ifndef __itkRemoveDwiChannelFilter_h_
+#define __itkRemoveDwiChannelFilter_h_
+
+#include "itkImageToImageFilter.h"
+#include "itkVectorImage.h"
+#include
+#include
+#include
+
+namespace itk{
+/** \class RemoveDwiChannelFilter
+ * \brief Remove spcified channels from diffusion-weighted image.
+ */
+
+template< class TInPixelType >
+class RemoveDwiChannelFilter :
+ public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > >
+{
+
+public:
+
+ typedef RemoveDwiChannelFilter Self;
+ typedef SmartPointer Pointer;
+ typedef SmartPointer ConstPointer;
+ typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass;
+
+ /** Method for creation through the object factory. */
+ itkFactorylessNewMacro(Self)
+ itkCloneMacro(Self)
+
+ /** Runtime information support. */
+ itkTypeMacro(RemoveDwiChannelFilter, ImageToImageFilter)
+
+ typedef typename Superclass::InputImageType InputImageType;
+ typedef typename Superclass::OutputImageType OutputImageType;
+ typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
+ typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionType DirectionType;
+ typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionContainerType DirectionContainerType;
+
+ void SetChannelIndices( std::vector< unsigned int > indices ){ m_ChannelIndices = indices; }
+ void SetDirections( typename DirectionContainerType::Pointer directions ){ m_Directions = directions; }
+ typename DirectionContainerType::Pointer GetNewDirections(){ return m_NewDirections; }
+
+ protected:
+ RemoveDwiChannelFilter();
+ ~RemoveDwiChannelFilter() {}
+ void PrintSelf(std::ostream& os, Indent indent) const;
+
+ void BeforeThreadedGenerateData();
+ void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType id );
+
+ std::vector< unsigned int > m_ChannelIndices;
+ typename DirectionContainerType::Pointer m_Directions;
+ typename DirectionContainerType::Pointer m_NewDirections;
+};
+
+}
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkRemoveDwiChannelFilter.txx"
+#endif
+
+#endif //__itkRemoveDwiChannelFilter_h_
+
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.txx
new file mode 100644
index 0000000000..cc87ff38dc
--- /dev/null
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.txx
@@ -0,0 +1,133 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center,
+Division of Medical and Biological Informatics.
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#ifndef __itkRemoveDwiChannelFilter_txx
+#define __itkRemoveDwiChannelFilter_txx
+
+#include
+#include
+#include
+
+#define _USE_MATH_DEFINES
+#include
+
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+
+namespace itk {
+
+
+template< class TInPixelType >
+RemoveDwiChannelFilter< TInPixelType>::RemoveDwiChannelFilter()
+{
+ this->SetNumberOfRequiredInputs( 1 );
+}
+
+template< class TInPixelType >
+void RemoveDwiChannelFilter< TInPixelType>::BeforeThreadedGenerateData()
+{
+ typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
+ if ( inputImagePointer->GetVectorLength()-m_ChannelIndices.size()<=0 )
+ itkExceptionMacro("No channels remaining!");
+
+ typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
+ outputImage->SetSpacing( inputImagePointer->GetSpacing() );
+ outputImage->SetOrigin( inputImagePointer->GetOrigin() );
+ outputImage->SetDirection( inputImagePointer->GetDirection() );
+ outputImage->SetLargestPossibleRegion( inputImagePointer->GetLargestPossibleRegion() );
+ outputImage->SetBufferedRegion( inputImagePointer->GetLargestPossibleRegion() );
+ outputImage->SetRequestedRegion( inputImagePointer->GetLargestPossibleRegion() );
+ outputImage->Allocate();
+ outputImage->SetVectorLength( inputImagePointer->GetVectorLength()-m_ChannelIndices.size() );
+ typename OutputImageType::PixelType nullPix;
+ nullPix.SetSize(outputImage->GetVectorLength());
+ nullPix.Fill(0);
+ outputImage->FillBuffer(nullPix);
+ this->SetNthOutput(0, outputImage);
+
+ m_NewDirections = DirectionContainerType::New();
+
+ int chIdx = 0;
+ for (unsigned int i=0; iGetVectorLength(); i++)
+ {
+ bool use = true;
+ for (unsigned int j=0; jInsertElement(chIdx, m_Directions->GetElement(i));
+ ++chIdx;
+ MITK_INFO << "Using channel " << i;
+ }
+ }
+}
+
+template< class TInPixelType >
+void RemoveDwiChannelFilter< TInPixelType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType id )
+{
+ typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
+
+ ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
+ oit.GoToBegin();
+
+ typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
+ typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
+
+ InputIteratorType git( inputImagePointer, outputRegionForThread );
+ git.GoToBegin();
+ while( !git.IsAtEnd() )
+ {
+ int chIdx = 0;
+ typename OutputImageType::PixelType pix = oit.Get();
+ for (unsigned int i=0; iGetVectorLength(); i++)
+ {
+ bool use = true;
+ for (unsigned int j=0; j
+void
+RemoveDwiChannelFilter< TInPixelType>
+::PrintSelf(std::ostream& os, Indent indent) const
+{
+ Superclass::PrintSelf(os,indent);
+}
+
+}
+
+#endif // __itkRemoveDwiChannelFilter_txx
diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake
index c7b3f04aed..1e2a79653d 100644
--- a/Modules/DiffusionImaging/DiffusionCore/files.cmake
+++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake
@@ -1,129 +1,131 @@
set(CPP_FILES
# DicomImport
DicomImport/mitkDicomDiffusionImageReader.cpp
# DicomImport/mitkGroupDiffusionHeadersFilter.cpp
DicomImport/mitkDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp
# DataStructures -> DWI
IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp
IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp
IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp
IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp
IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp
# DataStructures -> QBall
IODataStructures/QBallImages/mitkQBallImageSource.cpp
IODataStructures/QBallImages/mitkQBallImage.cpp
# DataStructures -> Tensor
IODataStructures/TensorImages/mitkTensorImage.cpp
#IODataStructures/mitkRegistrationObject.cpp
# Rendering
Rendering/vtkMaskedProgrammableGlyphFilter.cpp
Rendering/mitkVectorImageVtkGlyphMapper3D.cpp
Rendering/vtkOdfSource.cxx
Rendering/vtkThickPlane.cxx
Rendering/mitkOdfNormalizationMethodProperty.cpp
Rendering/mitkOdfScaleByProperty.cpp
# Algorithms
Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp
Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp
Algorithms/itkDwiGradientLengthCorrectionFilter.cpp
# Registration Algorithms & Co.
Algorithms/Registration/mitkRegistrationWrapper.cpp
Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp
# MultishellProcessing
Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.cpp
# Function Collection
mitkDiffusionFunctionCollection.cpp
)
set(H_FILES
# function Collection
mitkDiffusionFunctionCollection.h
# Rendering
Rendering/mitkDiffusionImageMapper.h
Rendering/mitkOdfVtkMapper2D.h
# Reconstruction
Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h
Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
Algorithms/Reconstruction/itkPointShell.h
Algorithms/Reconstruction/itkOrientationDistributionFunction.h
Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h
# MultishellProcessing
Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h
Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h
# IO Datastructures
IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h
# Algorithms
Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h
Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h
Algorithms/itkTensorDerivedMeasurementsFilter.h
Algorithms/itkBrainMaskExtractionImageFilter.h
Algorithms/itkB0ImageExtractionImageFilter.h
Algorithms/itkB0ImageExtractionToSeparateImageFilter.h
Algorithms/itkTensorImageToDiffusionImageFilter.h
Algorithms/itkTensorToL2NormImageFilter.h
Algorithms/itkGaussianInterpolateImageFunction.h
Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h
Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h
Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h
Algorithms/itkCartesianToPolarVectorImageFilter.h
Algorithms/itkPolarToCartesianVectorImageFilter.h
Algorithms/itkDistanceMapFilter.h
Algorithms/itkProjectionFilter.h
Algorithms/itkResidualImageFilter.h
Algorithms/itkExtractChannelFromRgbaImageFilter.h
Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h
Algorithms/itkMergeDiffusionImagesFilter.h
Algorithms/itkDwiPhantomGenerationFilter.h
Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h
Algorithms/itkMrtrixPeakImageConverter.h
Algorithms/itkFslPeakImageConverter.h
Algorithms/itkShCoefficientImageImporter.h
Algorithms/itkShCoefficientImageExporter.h
Algorithms/itkOdfMaximaExtractionFilter.h
Algorithms/itkResampleDwiImageFilter.h
Algorithms/itkDwiGradientLengthCorrectionFilter.h
Algorithms/itkAdcImageFilter.h
Algorithms/itkDwiNormilzationFilter.h
Algorithms/itkSplitDWImageFilter.h
+ Algorithms/itkRemoveDwiChannelFilter.h
+ Algorithms/itkExtractDwiChannelFilter.h
Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h
Algorithms/mitkDiffusionImageToDiffusionImageFilter.h
Algorithms/itkNonLocalMeansDenoisingFilter.h
Algorithms/itkVectorImageToImageFilter.h
)
set( TOOL_FILES
)
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp
index 2da65f6632..152251cc54 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp
@@ -1,227 +1,237 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Coindex[1]right (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "itkTractDensityImageFilter.h"
// VTK
#include
#include
#include
// misc
#include
#include
namespace itk{
template< class OutputImageType >
TractDensityImageFilter< OutputImageType >::TractDensityImageFilter()
: m_InvertImage(false)
, m_FiberBundle(NULL)
, m_UpsamplingFactor(1)
, m_InputImage(NULL)
, m_BinaryOutput(false)
, m_UseImageGeometry(false)
, m_OutputAbsoluteValues(false)
+ , m_UseTrilinearInterpolation(false)
{
}
template< class OutputImageType >
TractDensityImageFilter< OutputImageType >::~TractDensityImageFilter()
{
}
template< class OutputImageType >
itk::Point TractDensityImageFilter< OutputImageType >::GetItkPoint(double point[3])
{
itk::Point itkPoint;
itkPoint[0] = point[0];
itkPoint[1] = point[1];
itkPoint[2] = point[2];
return itkPoint;
}
template< class OutputImageType >
void TractDensityImageFilter< OutputImageType >::GenerateData()
{
// generate upsampled image
mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
typename OutputImageType::Pointer outImage = this->GetOutput();
// calculate new image parameters
itk::Vector newSpacing;
mitk::Point3D newOrigin;
itk::Matrix newDirection;
ImageRegion<3> upsampledRegion;
if (m_UseImageGeometry && !m_InputImage.IsNull())
{
MITK_INFO << "TractDensityImageFilter: using image geometry";
newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor;
upsampledRegion = m_InputImage->GetLargestPossibleRegion();
newOrigin = m_InputImage->GetOrigin();
typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize();
size[0] *= m_UpsamplingFactor;
size[1] *= m_UpsamplingFactor;
size[2] *= m_UpsamplingFactor;
upsampledRegion.SetSize(size);
newDirection = m_InputImage->GetDirection();
}
else
{
MITK_INFO << "TractDensityImageFilter: using fiber bundle geometry";
newSpacing = geometry->GetSpacing()/m_UpsamplingFactor;
newOrigin = geometry->GetOrigin();
mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
newOrigin[0] += bounds.GetElement(0);
newOrigin[1] += bounds.GetElement(2);
newOrigin[2] += bounds.GetElement(4);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
newDirection[j][i] = geometry->GetMatrixColumn(i)[j];
upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor);
upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor);
upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor);
}
typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize();
// apply new image parameters
outImage->SetSpacing( newSpacing );
outImage->SetOrigin( newOrigin );
outImage->SetDirection( newDirection );
outImage->SetRegions( upsampledRegion );
outImage->Allocate();
outImage->FillBuffer(0.0);
int w = upsampledSize[0];
int h = upsampledSize[1];
int d = upsampledSize[2];
// set/initialize output
OutPixelType* outImageBufferPointer = (OutPixelType*)outImage->GetBufferPointer();
// resample fiber bundle
float minSpacing = 1;
if(newSpacing[0]GetDeepCopy();
- m_FiberBundle->ResampleFibers(minSpacing);
+ m_FiberBundle->ResampleFibers(minSpacing/10);
MITK_INFO << "TractDensityImageFilter: starting image generation";
vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData();
vtkSmartPointer vLines = fiberPolyData->GetLines();
vLines->InitTraversal();
int numFibers = m_FiberBundle->GetNumFibers();
boost::progress_display disp(numFibers);
for( int i=0; iGetNextCell ( numPoints, points );
// fill output image
for( int j=0; j vertex = GetItkPoint(fiberPolyData->GetPoint(points[j]));
itk::Index<3> index;
itk::ContinuousIndex contIndex;
outImage->TransformPhysicalPointToIndex(vertex, index);
outImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
+ if (!m_UseTrilinearInterpolation)
+ {
+ if (m_BinaryOutput)
+ outImage->SetPixel(index, 1);
+ else
+ outImage->SetPixel(index, outImage->GetPixel(index)+0.01);
+ continue;
+ }
+
float frac_x = contIndex[0] - index[0];
float frac_y = contIndex[1] - index[1];
float frac_z = contIndex[2] - index[2];
if (frac_x<0)
{
index[0] -= 1;
frac_x += 1;
}
if (frac_y<0)
{
index[1] -= 1;
frac_y += 1;
}
if (frac_z<0)
{
index[2] -= 1;
frac_z += 1;
}
frac_x = 1-frac_x;
frac_y = 1-frac_y;
frac_z = 1-frac_z;
// int coordinates inside image?
if (index[0] < 0 || index[0] >= w-1)
continue;
if (index[1] < 0 || index[1] >= h-1)
continue;
if (index[2] < 0 || index[2] >= d-1)
continue;
if (m_BinaryOutput)
{
outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] = 1;
outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] = 1;
outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] = 1;
outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] = 1;
outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] = 1;
outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] = 1;
outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] = 1;
outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] = 1;
}
else
{
outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] += ( frac_x)*( frac_y)*( frac_z);
outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] += ( frac_x)*(1-frac_y)*( frac_z);
outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] += ( frac_x)*( frac_y)*(1-frac_z);
outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] += ( frac_x)*(1-frac_y)*(1-frac_z);
outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] += (1-frac_x)*( frac_y)*( frac_z);
outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] += (1-frac_x)*( frac_y)*(1-frac_z);
outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] += (1-frac_x)*(1-frac_y)*( frac_z);
outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] += (1-frac_x)*(1-frac_y)*(1-frac_z);
}
}
}
if (!m_OutputAbsoluteValues && !m_BinaryOutput)
{
MITK_INFO << "TractDensityImageFilter: max-normalizing output image";
OutPixelType max = 0;
for (int i=0; i0)
for (int i=0; i
#include
#include
#include
#include
namespace itk{
/**
* \brief Generates tract density images from input fiberbundles (Calamante 2010). */
template< class OutputImageType >
class TractDensityImageFilter : public ImageSource< OutputImageType >
{
public:
typedef TractDensityImageFilter Self;
typedef ProcessObject Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
typedef typename OutputImageType::PixelType OutPixelType;
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
itkTypeMacro( TractDensityImageFilter, ImageSource )
itkSetMacro( UpsamplingFactor, float) ///< use higher resolution for ouput image
itkGetMacro( UpsamplingFactor, float) ///< use higher resolution for ouput image
itkSetMacro( InvertImage, bool) ///< voxelvalue = 1-voxelvalue
itkGetMacro( InvertImage, bool) ///< voxelvalue = 1-voxelvalue
itkSetMacro( BinaryOutput, bool) ///< generate binary fiber envelope
itkGetMacro( BinaryOutput, bool) ///< generate binary fiber envelope
itkSetMacro( OutputAbsoluteValues, bool) ///< output absolute values of the number of fibers per voxel
itkGetMacro( OutputAbsoluteValues, bool) ///< output absolute values of the number of fibers per voxel
itkSetMacro( UseImageGeometry, bool) ///< use input image geometry to initialize output image
itkGetMacro( UseImageGeometry, bool) ///< use input image geometry to initialize output image
itkSetMacro( FiberBundle, mitk::FiberBundleX::Pointer) ///< input fiber bundle
itkSetMacro( InputImage, typename OutputImageType::Pointer) ///< use input image geometry to initialize output image
+ itkSetMacro( UseTrilinearInterpolation, bool )
void GenerateData();
protected:
itk::Point GetItkPoint(double point[3]);
TractDensityImageFilter();
virtual ~TractDensityImageFilter();
typename OutputImageType::Pointer m_InputImage; ///< use input image geometry to initialize output image
mitk::FiberBundleX::Pointer m_FiberBundle; ///< input fiber bundle
float m_UpsamplingFactor; ///< use higher resolution for ouput image
bool m_InvertImage; ///< voxelvalue = 1-voxelvalue
bool m_BinaryOutput; ///< generate binary fiber envelope
bool m_UseImageGeometry; ///< use input image geometry to initialize output image
bool m_OutputAbsoluteValues; ///< do not normalize image values to 0-1
+ bool m_UseTrilinearInterpolation;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTractDensityImageFilter.cpp"
#endif
#endif // __itkTractDensityImageFilter_h__
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
index af5223aafd..b8ec972e16 100755
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp
@@ -1,1258 +1,1204 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "itkTractsToDWIImageFilter.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
+#include
#include
+#include
namespace itk
{
template< class PixelType >
TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter()
: m_FiberBundle(NULL)
, m_StatusText("")
, m_UseConstantRandSeed(false)
, m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New())
- , m_NoAcquisitionSimulation(false)
{
m_RandGen->SetSeed();
}
template< class PixelType >
TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter()
{
}
template< class PixelType >
TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images )
{
+ int numFiberCompartments = m_Parameters.m_FiberModelList.size();
// create slice object
ImageRegion<2> sliceRegion;
sliceRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]);
sliceRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]);
Vector< double, 2 > sliceSpacing;
sliceSpacing[0] = m_UpsampledSpacing[0];
sliceSpacing[1] = m_UpsampledSpacing[1];
// frequency map slice
SliceType::Pointer fMapSlice = NULL;
if (m_Parameters.m_FrequencyMap.IsNotNull())
{
fMapSlice = SliceType::New();
ImageRegion<2> region;
region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]);
region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]);
fMapSlice->SetLargestPossibleRegion( region );
fMapSlice->SetBufferedRegion( region );
fMapSlice->SetRequestedRegion( region );
fMapSlice->Allocate();
fMapSlice->FillBuffer(0.0);
}
DoubleDwiType::Pointer newImage = DoubleDwiType::New();
newImage->SetSpacing( m_Parameters.m_ImageSpacing );
newImage->SetOrigin( m_Parameters.m_ImageOrigin );
newImage->SetDirection( m_Parameters.m_ImageDirection );
newImage->SetLargestPossibleRegion( m_Parameters.m_ImageRegion );
newImage->SetBufferedRegion( m_Parameters.m_ImageRegion );
newImage->SetRequestedRegion( m_Parameters.m_ImageRegion );
newImage->SetVectorLength( images.at(0)->GetVectorLength() );
newImage->Allocate();
std::vector< unsigned int > spikeVolume;
for (unsigned int i=0; iGetIntegerVariate()%images.at(0)->GetVectorLength());
std::sort (spikeVolume.begin(), spikeVolume.end());
std::reverse (spikeVolume.begin(), spikeVolume.end());
m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n";
m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
unsigned long lastTick = 0;
boost::progress_display disp(2*images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2));
for (unsigned int g=0; gGetVectorLength(); g++)
{
std::vector< unsigned int > spikeSlice;
while (!spikeVolume.empty() && spikeVolume.back()==g)
{
spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2));
spikeVolume.pop_back();
}
std::sort (spikeSlice.begin(), spikeSlice.end());
std::reverse (spikeSlice.begin(), spikeSlice.end());
for (unsigned int z=0; zGetLargestPossibleRegion().GetSize(2); z++)
{
std::vector< SliceType::Pointer > compartmentSlices;
std::vector< double > t2Vector;
for (unsigned int i=0; i* signalModel;
- if (iSetLargestPossibleRegion( sliceRegion );
slice->SetBufferedRegion( sliceRegion );
slice->SetRequestedRegion( sliceRegion );
slice->SetSpacing(sliceSpacing);
slice->Allocate();
slice->FillBuffer(0.0);
// extract slice from channel g
for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++)
for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++)
{
SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]);
if (fMapSlice.IsNotNull() && i==0)
fMapSlice->SetPixel(index2D, m_Parameters.m_FrequencyMap->GetPixel(index3D));
}
compartmentSlices.push_back(slice);
t2Vector.push_back(signalModel->GetT2());
}
if (this->GetAbortGenerateData())
return NULL;
// create k-sapce (inverse fourier transform slices)
itk::Size<2> outSize; outSize.SetElement(0, m_Parameters.m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_Parameters.m_ImageRegion.GetSize(1));
itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New();
idft->SetCompartmentImages(compartmentSlices);
idft->SetT2(t2Vector);
idft->SetUseConstantRandSeed(m_UseConstantRandSeed);
idft->SetParameters(m_Parameters);
idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0);
idft->SetDiffusionGradientDirection(m_Parameters.GetGradientDirection(g));
idft->SetFrequencyMapSlice(fMapSlice);
idft->SetOutSize(outSize);
int numSpikes = 0;
while (!spikeSlice.empty() && spikeSlice.back()==z)
{
numSpikes++;
spikeSlice.pop_back();
}
idft->SetSpikesPerSlice(numSpikes);
idft->Update();
ComplexSliceType::Pointer fSlice;
fSlice = idft->GetOutput();
++disp;
unsigned long newTick = 50*disp.count()/disp.expected_count();
for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
m_StatusText += "*";
lastTick = newTick;
// fourier transform slice
SliceType::Pointer newSlice;
itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New();
dft->SetInput(fSlice);
dft->Update();
newSlice = dft->GetOutput();
// put slice back into channel g
for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++)
for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++)
{
DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D);
pix3D[g] = newSlice->GetPixel(index2D);
newImage->SetPixel(index3D, pix3D);
}
++disp;
newTick = 50*disp.count()/disp.expected_count();
for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
m_StatusText += "*";
lastTick = newTick;
}
}
m_StatusText += "\n\n";
return newImage;
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::GenerateData()
{
m_TimeProbe.Start();
m_StatusText = "Starting simulation\n";
// check input data
if (m_FiberBundle.IsNull())
itkExceptionMacro("Input fiber bundle is NULL!");
+ if (m_Parameters.m_FiberModelList.empty())
+ itkExceptionMacro("No diffusion model for fiber compartments defined!");
+
if (m_Parameters.m_NonFiberModelList.empty())
itkExceptionMacro("No diffusion model for non-fiber compartments defined!");
int baselineIndex = m_Parameters.GetFirstBaselineIndex();
if (baselineIndex<0)
itkExceptionMacro("No baseline index found!");
+ if (!m_Parameters.m_SimulateKspaceAcquisition)
+ m_Parameters.m_DoAddGibbsRinging = false;
+
if (m_UseConstantRandSeed) // always generate the same random numbers?
m_RandGen->SetSeed(0);
else
m_RandGen->SetSeed();
// initialize output dwi image
ImageRegion<3> croppedRegion = m_Parameters.m_ImageRegion; croppedRegion.SetSize(1, croppedRegion.GetSize(1)*m_Parameters.m_CroppingFactor);
itk::Point shiftedOrigin = m_Parameters.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_ImageRegion.GetSize(1)-croppedRegion.GetSize(1))*m_Parameters.m_ImageSpacing[1]/2;
typename OutputImageType::Pointer outImage = OutputImageType::New();
outImage->SetSpacing( m_Parameters.m_ImageSpacing );
outImage->SetOrigin( shiftedOrigin );
outImage->SetDirection( m_Parameters.m_ImageDirection );
outImage->SetLargestPossibleRegion( croppedRegion );
outImage->SetBufferedRegion( croppedRegion );
outImage->SetRequestedRegion( croppedRegion );
outImage->SetVectorLength( m_Parameters.GetNumVolumes() );
outImage->Allocate();
typename OutputImageType::PixelType temp;
temp.SetSize(m_Parameters.GetNumVolumes());
temp.Fill(0.0);
outImage->FillBuffer(temp);
// ADJUST GEOMETRY FOR FURTHER PROCESSING
// is input slize size a power of two?
unsigned int x=m_Parameters.m_ImageRegion.GetSize(0); unsigned int y=m_Parameters.m_ImageRegion.GetSize(1);
ItkDoubleImgType::SizeType pad; pad[0]=x%2; pad[1]=y%2; pad[2]=0;
m_Parameters.m_ImageRegion.SetSize(0, x+pad[0]);
m_Parameters.m_ImageRegion.SetSize(1, y+pad[1]);
if (m_Parameters.m_FrequencyMap.IsNotNull() && (pad[0]>0 || pad[1]>0))
{
itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New();
zeroPadder->SetInput(m_Parameters.m_FrequencyMap);
zeroPadder->SetConstant(0);
zeroPadder->SetPadUpperBound(pad);
zeroPadder->Update();
m_Parameters.m_FrequencyMap = zeroPadder->GetOutput();
}
if (m_Parameters.m_MaskImage.IsNotNull() && (pad[0]>0 || pad[1]>0))
{
itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New();
zeroPadder->SetInput(m_Parameters.m_MaskImage);
zeroPadder->SetConstant(0);
zeroPadder->SetPadUpperBound(pad);
zeroPadder->Update();
m_Parameters.m_MaskImage = zeroPadder->GetOutput();
}
// Apply in-plane upsampling for Gibbs ringing artifact
double upsampling = 1;
if (m_Parameters.m_DoAddGibbsRinging)
upsampling = 2;
m_UpsampledSpacing = m_Parameters.m_ImageSpacing;
m_UpsampledSpacing[0] /= upsampling;
m_UpsampledSpacing[1] /= upsampling;
m_UpsampledImageRegion = m_Parameters.m_ImageRegion;
m_UpsampledImageRegion.SetSize(0, m_Parameters.m_ImageRegion.GetSize()[0]*upsampling);
m_UpsampledImageRegion.SetSize(1, m_Parameters.m_ImageRegion.GetSize()[1]*upsampling);
m_UpsampledOrigin = m_Parameters.m_ImageOrigin;
m_UpsampledOrigin[0] -= m_Parameters.m_ImageSpacing[0]/2; m_UpsampledOrigin[0] += m_UpsampledSpacing[0]/2;
m_UpsampledOrigin[1] -= m_Parameters.m_ImageSpacing[1]/2; m_UpsampledOrigin[1] += m_UpsampledSpacing[1]/2;
m_UpsampledOrigin[2] -= m_Parameters.m_ImageSpacing[2]/2; m_UpsampledOrigin[2] += m_UpsampledSpacing[2]/2;
// generate double images to store the individual compartment signals
- std::vector< DoubleDwiType::Pointer > compartments;
- for (unsigned int i=0; iSetSpacing( m_UpsampledSpacing );
doubleDwi->SetOrigin( m_UpsampledOrigin );
doubleDwi->SetDirection( m_Parameters.m_ImageDirection );
doubleDwi->SetLargestPossibleRegion( m_UpsampledImageRegion );
doubleDwi->SetBufferedRegion( m_UpsampledImageRegion );
doubleDwi->SetRequestedRegion( m_UpsampledImageRegion );
doubleDwi->SetVectorLength( m_Parameters.GetNumVolumes() );
doubleDwi->Allocate();
DoubleDwiType::PixelType pix;
pix.SetSize(m_Parameters.GetNumVolumes());
pix.Fill(0.0);
doubleDwi->FillBuffer(pix);
- compartments.push_back(doubleDwi);
+ m_CompartmentImages.push_back(doubleDwi);
}
// initialize output volume fraction images
m_VolumeFractions.clear();
- for (unsigned int i=0; iSetSpacing( m_UpsampledSpacing );
doubleImg->SetOrigin( m_UpsampledOrigin );
doubleImg->SetDirection( m_Parameters.m_ImageDirection );
doubleImg->SetLargestPossibleRegion( m_UpsampledImageRegion );
doubleImg->SetBufferedRegion( m_UpsampledImageRegion );
doubleImg->SetRequestedRegion( m_UpsampledImageRegion );
doubleImg->Allocate();
doubleImg->FillBuffer(0);
m_VolumeFractions.push_back(doubleImg);
}
// get volume fraction images
ItkDoubleImgType::Pointer sumImage = ItkDoubleImgType::New();
bool foundVolumeFractionImage = false;
- for (unsigned int i=0; iGetVolumeFractionImage().IsNotNull())
{
foundVolumeFractionImage = true;
itk::ConstantPadImageFilter::Pointer zeroPadder = itk::ConstantPadImageFilter::New();
zeroPadder->SetInput(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
zeroPadder->SetConstant(0);
zeroPadder->SetPadUpperBound(pad);
zeroPadder->Update();
m_Parameters.m_NonFiberModelList[i]->SetVolumeFractionImage(zeroPadder->GetOutput());
sumImage->SetSpacing( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetSpacing() );
sumImage->SetOrigin( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetOrigin() );
sumImage->SetDirection( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetDirection() );
sumImage->SetLargestPossibleRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() );
sumImage->SetBufferedRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() );
sumImage->SetRequestedRegion( m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion() );
sumImage->Allocate();
sumImage->FillBuffer(0);
break;
}
}
if (!foundVolumeFractionImage)
{
sumImage->SetSpacing( m_UpsampledSpacing );
sumImage->SetOrigin( m_UpsampledOrigin );
sumImage->SetDirection( m_Parameters.m_ImageDirection );
sumImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
sumImage->SetBufferedRegion( m_UpsampledImageRegion );
sumImage->SetRequestedRegion( m_UpsampledImageRegion );
sumImage->Allocate();
sumImage->FillBuffer(0.0);
}
- for (unsigned int i=0; iGetVolumeFractionImage().IsNull())
{
ItkDoubleImgType::Pointer doubleImg = ItkDoubleImgType::New();
doubleImg->SetSpacing( sumImage->GetSpacing() );
doubleImg->SetOrigin( sumImage->GetOrigin() );
doubleImg->SetDirection( sumImage->GetDirection() );
doubleImg->SetLargestPossibleRegion( sumImage->GetLargestPossibleRegion() );
doubleImg->SetBufferedRegion( sumImage->GetLargestPossibleRegion() );
doubleImg->SetRequestedRegion( sumImage->GetLargestPossibleRegion() );
doubleImg->Allocate();
- doubleImg->FillBuffer(1.0/m_Parameters.m_NonFiberModelList.size());
+ doubleImg->FillBuffer(1.0/numNonFiberCompartments);
m_Parameters.m_NonFiberModelList[i]->SetVolumeFractionImage(doubleImg);
}
ImageRegionIterator it(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage(), m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion());
while(!it.IsAtEnd())
{
sumImage->SetPixel(it.GetIndex(), sumImage->GetPixel(it.GetIndex())+it.Get());
++it;
}
}
- for (unsigned int i=0; i it(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage(), m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion());
while(!it.IsAtEnd())
{
if (sumImage->GetPixel(it.GetIndex())>0)
it.Set(it.Get()/sumImage->GetPixel(it.GetIndex()));
++it;
}
}
// resample mask image and frequency map to fit upsampled geometry
if (m_Parameters.m_DoAddGibbsRinging)
{
if (m_Parameters.m_MaskImage.IsNotNull())
{
// rescale mask image (otherwise there are problems with the resampling)
itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New();
rescaler->SetInput(0,m_Parameters.m_MaskImage);
rescaler->SetOutputMaximum(100);
rescaler->SetOutputMinimum(0);
rescaler->Update();
// resample mask image
itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New();
resampler->SetInput(rescaler->GetOutput());
resampler->SetOutputParametersFromImage(m_Parameters.m_MaskImage);
resampler->SetSize(m_UpsampledImageRegion.GetSize());
resampler->SetOutputSpacing(m_UpsampledSpacing);
resampler->SetOutputOrigin(m_UpsampledOrigin);
itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator
= itk::NearestNeighborInterpolateImageFunction::New();
resampler->SetInterpolator(nn_interpolator);
resampler->Update();
m_Parameters.m_MaskImage = resampler->GetOutput();
itk::ImageFileWriter::Pointer w = itk::ImageFileWriter::New();
w->SetFileName("/local/mask_ups.nrrd");
w->SetInput(m_Parameters.m_MaskImage);
w->Update();
}
// resample frequency map
if (m_Parameters.m_FrequencyMap.IsNotNull())
{
itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New();
resampler->SetInput(m_Parameters.m_FrequencyMap);
resampler->SetOutputParametersFromImage(m_Parameters.m_FrequencyMap);
resampler->SetSize(m_UpsampledImageRegion.GetSize());
resampler->SetOutputSpacing(m_UpsampledSpacing);
resampler->SetOutputOrigin(m_UpsampledOrigin);
itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator
= itk::NearestNeighborInterpolateImageFunction::New();
resampler->SetInterpolator(nn_interpolator);
resampler->Update();
m_Parameters.m_FrequencyMap = resampler->GetOutput();
}
}
// no input tissue mask is set -> create default
bool maskImageSet = true;
if (m_Parameters.m_MaskImage.IsNull())
{
m_StatusText += "No tissue mask set\n";
MITK_INFO << "No tissue mask set";
m_Parameters.m_MaskImage = ItkUcharImgType::New();
m_Parameters.m_MaskImage->SetSpacing( m_UpsampledSpacing );
m_Parameters.m_MaskImage->SetOrigin( m_UpsampledOrigin );
m_Parameters.m_MaskImage->SetDirection( m_Parameters.m_ImageDirection );
m_Parameters.m_MaskImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
m_Parameters.m_MaskImage->SetBufferedRegion( m_UpsampledImageRegion );
m_Parameters.m_MaskImage->SetRequestedRegion( m_UpsampledImageRegion );
m_Parameters.m_MaskImage->Allocate();
m_Parameters.m_MaskImage->FillBuffer(1);
maskImageSet = false;
}
else
{
m_StatusText += "Using tissue mask\n";
MITK_INFO << "Using tissue mask";
}
m_Parameters.m_ImageRegion = croppedRegion;
x=m_Parameters.m_ImageRegion.GetSize(0); y=m_Parameters.m_ImageRegion.GetSize(1);
if ( x%2 == 1 )
m_Parameters.m_ImageRegion.SetSize(0, x+1);
if ( y%2 == 1 )
m_Parameters.m_ImageRegion.SetSize(1, y+1);
// resample fiber bundle for sufficient voxel coverage
m_StatusText += "\n"+this->GetTime()+" > Resampling fibers ...\n";
double segmentVolume = 0.0001;
float minSpacing = 1;
if(m_UpsampledSpacing[0]GetDeepCopy();
double volumeAccuracy = 10;
fiberBundle->ResampleFibers(minSpacing/volumeAccuracy);
double mmRadius = m_Parameters.m_AxonRadius/1000;
if (mmRadius>0)
segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/volumeAccuracy;
double maxVolume = 0;
double voxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2];
ofstream logFile;
if (m_Parameters.m_DoAddMotion)
{
std::string fileName = "fiberfox_motion_0.log";
std::string filePath = mitk::IOUtil::GetTempPath();
if (m_Parameters.m_OutputPath.size()>0)
filePath = m_Parameters.m_OutputPath;
int c = 1;
while (itksys::SystemTools::FileExists((filePath+fileName).c_str()))
{
fileName = "fiberfox_motion_";
fileName += boost::lexical_cast(c);
fileName += ".log";
c++;
}
logFile.open((filePath+fileName).c_str());
logFile << "0 rotation: 0,0,0; translation: 0,0,0\n";
if (m_Parameters.m_DoRandomizeMotion)
{
m_StatusText += "Adding random motion artifacts:\n";
m_StatusText += "Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_Rotation) + "°\n";
m_StatusText += "Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_Translation) + "mm\n";
}
else
{
m_StatusText += "Adding linear motion artifacts:\n";
m_StatusText += "Maximum rotation: " + boost::lexical_cast(m_Parameters.m_Rotation) + "°\n";
m_StatusText += "Maximum translation: " + boost::lexical_cast(m_Parameters.m_Translation) + "mm\n";
}
m_StatusText += "Motion logfile: " + (filePath+fileName) + "\n";
MITK_INFO << "Adding motion artifacts";
MITK_INFO << "Maximum rotation: " << m_Parameters.m_Rotation;
MITK_INFO << "Maxmimum translation: " << m_Parameters.m_Translation;
}
maxVolume = 0;
- m_StatusText += "\n"+this->GetTime()+" > Generating " + boost::lexical_cast(m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size()) + "-compartment diffusion-weighted signal.\n";
+ m_StatusText += "\n"+this->GetTime()+" > Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal.\n";
int numFibers = m_FiberBundle->GetNumFibers();
boost::progress_display disp(numFibers*m_Parameters.GetNumVolumes());
// get transform for motion artifacts
FiberBundleType fiberBundleTransformed = fiberBundle;
- VectorType rotation = m_Parameters.m_Rotation/m_Parameters.GetNumVolumes();
- VectorType translation = m_Parameters.m_Translation/m_Parameters.GetNumVolumes();
+ DoubleVectorType rotation = m_Parameters.m_Rotation/m_Parameters.GetNumVolumes();
+ DoubleVectorType translation = m_Parameters.m_Translation/m_Parameters.GetNumVolumes();
// creat image to hold transformed mask (motion artifact)
ItkUcharImgType::Pointer tempTissueMask = ItkUcharImgType::New();
itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New();
duplicator->SetInputImage(m_Parameters.m_MaskImage);
duplicator->Update();
tempTissueMask = duplicator->GetOutput();
// second upsampling needed for motion artifacts
ImageRegion<3> upsampledImageRegion = m_UpsampledImageRegion;
- itk::Vector upsampledSpacing = m_UpsampledSpacing;
+ DoubleVectorType upsampledSpacing = m_UpsampledSpacing;
upsampledSpacing[0] /= 4;
upsampledSpacing[1] /= 4;
upsampledSpacing[2] /= 4;
upsampledImageRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]*4);
upsampledImageRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]*4);
upsampledImageRegion.SetSize(2, m_UpsampledImageRegion.GetSize()[2]*4);
itk::Point upsampledOrigin = m_UpsampledOrigin;
upsampledOrigin[0] -= m_UpsampledSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2;
upsampledOrigin[1] -= m_UpsampledSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2;
upsampledOrigin[2] -= m_UpsampledSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2;
ItkUcharImgType::Pointer upsampledTissueMask = ItkUcharImgType::New();
itk::ResampleImageFilter::Pointer upsampler = itk::ResampleImageFilter::New();
upsampler->SetInput(m_Parameters.m_MaskImage);
upsampler->SetOutputParametersFromImage(m_Parameters.m_MaskImage);
upsampler->SetSize(upsampledImageRegion.GetSize());
upsampler->SetOutputSpacing(upsampledSpacing);
upsampler->SetOutputOrigin(upsampledOrigin);
itk::NearestNeighborInterpolateImageFunction::Pointer nn_interpolator
= itk::NearestNeighborInterpolateImageFunction::New();
upsampler->SetInterpolator(nn_interpolator);
upsampler->Update();
upsampledTissueMask = upsampler->GetOutput();
unsigned long lastTick = 0;
- if (true)
+ switch (m_Parameters.m_DiffusionDirectionMode)
+ {
+ case(FiberfoxParameters<>::FIBER_TANGENT_DIRECTIONS):
{
-
m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n";
m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
for (unsigned int g=0; gSetSpacing( m_UpsampledSpacing );
intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin );
intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection );
intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion );
intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion );
intraAxonalVolumeImage->Allocate();
intraAxonalVolumeImage->FillBuffer(0);
vtkPolyData* fiberPolyData = fiberBundleTransformed->GetFiberPolyData();
// generate fiber signal (if there are any fiber models present)
if (!m_Parameters.m_FiberModelList.empty())
for( int i=0; iGetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (numPoints<2)
continue;
for( int j=0; jGetAbortGenerateData())
{
m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
return;
}
double* temp = points->GetPoint(j);
itk::Point vertex = GetItkPoint(temp);
itk::Vector v = GetItkVector(temp);
itk::Vector dir(3);
if (jGetPoint(j+1))-v;
else
dir = v-GetItkVector(points->GetPoint(j-1));
if (dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2])
continue;
itk::Index<3> idx;
itk::ContinuousIndex contIndex;
tempTissueMask->TransformPhysicalPointToIndex(vertex, idx);
tempTissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
if (!tempTissueMask->GetLargestPossibleRegion().IsInside(idx) || tempTissueMask->GetPixel(idx)<=0)
continue;
// generate signal for each fiber compartment
- for (unsigned int k=0; kSetFiberDirection(dir);
- DoubleDwiType::PixelType pix = compartments.at(k)->GetPixel(idx);
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(idx);
pix[g] += segmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g);
- compartments.at(k)->SetPixel(idx, pix);
+
+ m_CompartmentImages.at(k)->SetPixel(idx, pix);
}
// update fiber volume image
double vol = intraAxonalVolumeImage->GetPixel(idx) + segmentVolume;
intraAxonalVolumeImage->SetPixel(idx, vol);
if (g==0 && vol>maxVolume)
maxVolume = vol;
}
// progress report
++disp;
unsigned long newTick = 50*disp.count()/disp.expected_count();
for (unsigned int tick = 0; tick<(newTick-lastTick); tick++)
m_StatusText += "*";
lastTick = newTick;
}
// generate non-fiber signal
ImageRegionIterator it3(tempTissueMask, tempTissueMask->GetLargestPossibleRegion());
double fact = 1;
if (m_Parameters.m_AxonRadius<0.0001 || maxVolume>voxelVolume)
fact = voxelVolume/maxVolume;
while(!it3.IsAtEnd())
{
if (it3.Get()>0)
{
DoubleDwiType::IndexType index = it3.GetIndex();
// get fiber volume fraction
double intraAxonalVolume = intraAxonalVolumeImage->GetPixel(index)*fact;
- for (unsigned int i=0; iGetPixel(index);
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] *= fact;
- compartments.at(i)->SetPixel(index, pix);
+ m_CompartmentImages.at(i)->SetPixel(index, pix);
}
if (intraAxonalVolume>0.0001 && m_Parameters.m_DoDisablePartialVolume) // only fiber in voxel
{
- DoubleDwiType::PixelType pix = compartments.at(0)->GetPixel(index);
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
pix[g] *= voxelVolume/intraAxonalVolume;
- compartments.at(0)->SetPixel(index, pix);
+ m_CompartmentImages.at(0)->SetPixel(index, pix);
m_VolumeFractions.at(0)->SetPixel(index, 1);
- for (unsigned int i=1; iGetPixel(index);
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] = 0;
- compartments.at(i)->SetPixel(index, pix);
+ m_CompartmentImages.at(i)->SetPixel(index, pix);
}
}
else
{
m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/voxelVolume);
double extraAxonalVolume = voxelVolume-intraAxonalVolume; // non-fiber volume
double interAxonalVolume = 0;
- if (m_Parameters.m_FiberModelList.size()>1)
+ if (numFiberCompartments>1)
interAxonalVolume = extraAxonalVolume * intraAxonalVolume/voxelVolume; // inter-axonal fraction of non fiber compartment scales linearly with f
double other = extraAxonalVolume - interAxonalVolume; // rest of compartment
- double singleinter = interAxonalVolume/(m_Parameters.m_FiberModelList.size()-1);
+ double singleinter = interAxonalVolume/(numFiberCompartments-1);
// adjust non-fiber and intra-axonal signal
- for (unsigned int i=1; iGetPixel(index);
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
if (intraAxonalVolume>0) // remove scaling by intra-axonal volume from inter-axonal compartment
pix[g] /= intraAxonalVolume;
pix[g] *= singleinter;
- compartments.at(i)->SetPixel(index, pix);
+ m_CompartmentImages.at(i)->SetPixel(index, pix);
m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume);
}
- for (unsigned int i=0; i point;
- tempTissueMask->TransformIndexToPhysicalPoint(index, point);
-
- if (m_Parameters.m_DoAddMotion)
- {
- if (m_Parameters.m_DoRandomizeMotion && g>0)
- point = fiberBundle->TransformPoint(point.GetVnlVector(), -rotation[0],-rotation[1],-rotation[2],-translation[0],-translation[1],-translation[2]);
- else
- point = fiberBundle->TransformPoint(point.GetVnlVector(), -rotation[0]*g,-rotation[1]*g,-rotation[2]*g,-translation[0]*g,-translation[1]*g,-translation[2]*g);
- }
+ itk::Point point;
+ tempTissueMask->TransformIndexToPhysicalPoint(index, point);
+ if (m_Parameters.m_DoAddMotion)
+ {
+ if (m_Parameters.m_DoRandomizeMotion && g>0)
+ point = fiberBundle->TransformPoint(point.GetVnlVector(), -rotation[0],-rotation[1],-rotation[2],-translation[0],-translation[1],-translation[2]);
+ else
+ point = fiberBundle->TransformPoint(point.GetVnlVector(), -rotation[0]*g,-rotation[1]*g,-rotation[2]*g,-translation[0]*g,-translation[1]*g,-translation[2]*g);
+ }
+ for (int i=0; i1)
+ if (numNonFiberCompartments>1)
{
DoubleDwiType::IndexType newIndex;
m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->TransformPhysicalPointToIndex(point, newIndex);
if (!m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetLargestPossibleRegion().IsInside(newIndex))
- {
- MITK_WARN << "Volume fraction image is too small for the chosen motion artifacts! Due to motion a volume fraction outside of the specified image volume is requested.";
continue;
- }
weight = m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()->GetPixel(newIndex);
}
- DoubleDwiType::Pointer doubleDwi = compartments.at(i+m_Parameters.m_FiberModelList.size());
+ DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(i+numFiberCompartments);
DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index);
pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*other*weight;
doubleDwi->SetPixel(index, pix);
- m_VolumeFractions.at(i+m_Parameters.m_FiberModelList.size())->SetPixel(index, other/voxelVolume*weight);
+ m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, other/voxelVolume*weight);
}
}
}
++it3;
}
// move fibers
if (m_Parameters.m_DoAddMotion && gGetDeepCopy();
rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[0]*2)-m_Parameters.m_Rotation[0];
rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[1]*2)-m_Parameters.m_Rotation[1];
rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Rotation[2]*2)-m_Parameters.m_Rotation[2];
translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[0]*2)-m_Parameters.m_Translation[0];
translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[1]*2)-m_Parameters.m_Translation[1];
translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_Translation[2]*2)-m_Parameters.m_Translation[2];
}
// rotate mask image
if (maskImageSet)
{
ImageRegionIterator maskIt(upsampledTissueMask, upsampledTissueMask->GetLargestPossibleRegion());
tempTissueMask->FillBuffer(0);
while(!maskIt.IsAtEnd())
{
if (maskIt.Get()<=0)
{
++maskIt;
continue;
}
DoubleDwiType::IndexType index = maskIt.GetIndex();
itk::Point point;
upsampledTissueMask->TransformIndexToPhysicalPoint(index, point);
if (m_Parameters.m_DoRandomizeMotion)
point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]);
else
point = fiberBundle->TransformPoint(point.GetVnlVector(), rotation[0]*(g+1),rotation[1]*(g+1),rotation[2]*(g+1),translation[0]*(g+1),translation[1]*(g+1),translation[2]*(g+1));
tempTissueMask->TransformPhysicalPointToIndex(point, index);
if (tempTissueMask->GetLargestPossibleRegion().IsInside(index))
tempTissueMask->SetPixel(index,100);
++maskIt;
}
}
// rotate fibers
if (logFile.is_open())
{
logFile << g+1 << " rotation: " << rotation[0] << "," << rotation[1] << "," << rotation[2] << ";";
logFile << " translation: " << translation[0] << "," << translation[1] << "," << translation[2] << "\n";
}
fiberBundleTransformed->TransformFibers(rotation[0],rotation[1],rotation[2],translation[0],translation[1],translation[2]);
}
}
+ break;
}
- else
+ case (FiberfoxParameters<>::MAIN_FIBER_DIRECTIONS):
{
- m_NoAcquisitionSimulation = true;
- m_Parameters.m_SignalScale = 1;
- std::vector< RawShModel* > wm_models;
- std::vector< RawShModel* > gm_models;
+ typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType;
+ typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType;
+
+ itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New();
+ fOdfFilter->SetFiberBundle(fiberBundle);
+ fOdfFilter->SetMaskImage(tempTissueMask);
+ fOdfFilter->SetAngularThreshold(cos(45*M_PI/180));
+ fOdfFilter->SetNormalizeVectors(false);
+ fOdfFilter->SetUseWorkingCopy(false);
+ fOdfFilter->SetSizeThreshold(0);
+ fOdfFilter->SetMaxNumDirections(3);
+ fOdfFilter->Update();
+ ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer();
- if (m_InputDwi.IsNotNull())
{
- typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, double > TensorReconstructionImageFilterType;
- TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New();
- filter->SetGradientImage( m_InputDwi->GetDirections(), m_InputDwi->GetVectorImage() );
- filter->SetBValue(m_InputDwi->GetReferenceBValue());
- filter->Update();
- itk::Image< itk::DiffusionTensor3D< double >, 3 >::Pointer tensorImage = filter->GetOutput();
-
- const int shOrder = 4;
- const int NumCoeffs = (shOrder*shOrder + shOrder + 2)/2 + shOrder;
- typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType;
- typename QballFilterType::Pointer qballfilter = QballFilterType::New();
- qballfilter->SetGradientImage( m_InputDwi->GetDirections(), m_InputDwi->GetVectorImage() );
- qballfilter->SetBValue(m_InputDwi->GetReferenceBValue());
- qballfilter->SetLambda(0.006);
- qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL);
- qballfilter->Update();
- QballFilterType::CoefficientImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage();
-
- int b0Index;
- for (unsigned int i=0; iGetDirectionsWithoutMeasurementFrame()->Size(); i++)
- if ( m_InputDwi->GetDirectionsWithoutMeasurementFrame()->GetElement(i).magnitude()<0.001 )
- {
- b0Index = i;
- break;
- }
+ ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage();
+ typedef itk::ImageFileWriter< ItkUcharImgType > WriterType;
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName("/local/NumDirections.nrrd");
+ writer->SetInput(numDirImage);
+ writer->Update();
+ }
- m_StatusText += "Sampling signal kernels.\n";
- ImageRegionIterator< itk::Image< itk::DiffusionTensor3D< double >, 3 > > it(tensorImage, tensorImage->GetLargestPossibleRegion());
- while(!it.IsAtEnd())
+ ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New();
+ intraAxonalVolumeImage->SetSpacing( m_UpsampledSpacing );
+ intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin );
+ intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection );
+ intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
+ intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion );
+ intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion );
+ intraAxonalVolumeImage->Allocate();
+ intraAxonalVolumeImage->FillBuffer(0);
+
+ itk::TractDensityImageFilter< ItkDoubleImgType >::Pointer generator = itk::TractDensityImageFilter< ItkDoubleImgType >::New();
+ generator->SetFiberBundle(fiberBundle);
+ generator->SetBinaryOutput(false);
+ generator->SetOutputAbsoluteValues(false);
+ generator->SetInputImage(intraAxonalVolumeImage);
+ generator->SetUseImageGeometry(true);
+ generator->Update();
+ intraAxonalVolumeImage = generator->GetOutput();
+
+ m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n";
+ m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
+ boost::progress_display disp(tempTissueMask->GetLargestPossibleRegion().GetNumberOfPixels());
+ ImageRegionIterator< ItkUcharImgType > it(tempTissueMask, tempTissueMask->GetLargestPossibleRegion());
+ while(!it.IsAtEnd())
+ {
+ ++disp;
+ unsigned long newTick = 50*disp.count()/disp.expected_count();
+ for (unsigned int tick = 0; tick<(newTick-lastTick); tick++)
+ m_StatusText += "*";
+ lastTick = newTick;
+
+ if (this->GetAbortGenerateData())
{
- bool valid = true;
- for (int i=0; iGetNumberOfChannels(); i++)
- {
- if (m_InputDwi->GetVectorImage()->GetPixel(it.GetIndex())[i]<=0 || m_InputDwi->GetVectorImage()->GetPixel(it.GetIndex())[i]>m_InputDwi->GetVectorImage()->GetPixel(it.GetIndex())[b0Index])
- valid = false;
- }
- if (valid && tempTissueMask->GetPixel(it.GetIndex())>0)
+ m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
+ return;
+ }
+
+ if (it.Get()>0)
+ {
+ int count = 0;
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(it.GetIndex());
+ for (unsigned int i=0; iSize(); i++)
{
- itk::DiffusionTensor3D< double > tensor = it.Get();
- double FA = tensor.GetFractionalAnisotropy();
- if (FA>0.7 && FA<0.9)
- {
- RawShModel* model = new RawShModel();
- model->SetGradientList( m_Parameters.GetGradientDirections() );
- itk::Vector< float, NumCoeffs > itkv = itkFeatureImage->GetPixel(it.GetIndex());
- vnl_vector_fixed< double, NumCoeffs > coeffs;
- for (unsigned int c=0; cSetB0Signal( m_InputDwi->GetVectorImage()->GetPixel(it.GetIndex())[b0Index] );
- if (!model->SetShCoefficients( coeffs ))
- {
- ++it;
- continue;
- }
- wm_models.push_back(model);
- MITK_INFO << "WM KERNEL: " << it.GetIndex();
- }
- else if (FA>0.0 && FA<0.15)
+ itk::Vector< float, 3> dir = directionImageContainer->GetElement(i)->GetPixel(it.GetIndex());
+ double norm = dir.GetNorm();
+ if (norm>0.0001)
{
- RawShModel* model = new RawShModel();
- model->SetGradientList( m_Parameters.GetGradientDirections() );
- itk::Vector< float, NumCoeffs > itkv = itkFeatureImage->GetPixel(it.GetIndex());
- vnl_vector_fixed< double, NumCoeffs > coeffs;
- for (unsigned int c=0; cSetB0Signal( m_InputDwi->GetVectorImage()->GetPixel(it.GetIndex())[b0Index] );
- if (!model->SetShCoefficients( coeffs ))
- {
- ++it;
- continue;
- }
- gm_models.push_back(model);
- MITK_INFO << "GM/CSF KERNEL: " << it.GetIndex();
+ int modelIndex = m_RandGen->GetIntegerVariate(m_Parameters.m_FiberModelList.size()-1);
+ m_Parameters.m_FiberModelList.at(modelIndex)->SetFiberDirection(dir);
+ pix += m_Parameters.m_FiberModelList.at(modelIndex)->SimulateMeasurement()*norm;
+ count++;
}
+ }
+ if (count>0)
+ pix /= count;
- if (wm_models.size()>=100 && gm_models.size()>=100)
- break;
+ pix *= intraAxonalVolumeImage->GetPixel(it.GetIndex());
+ // GM/CSF
+ {
+ int modelIndex = m_RandGen->GetIntegerVariate(m_Parameters.m_NonFiberModelList.size()-1);
+ pix += (1-intraAxonalVolumeImage->GetPixel(it.GetIndex()))*m_Parameters.m_NonFiberModelList.at(modelIndex)->SimulateMeasurement();
}
- ++it;
+
+ m_CompartmentImages.at(0)->SetPixel(it.GetIndex(), pix);
}
- MITK_INFO << "Using pool of " << wm_models.size() << " WM and " << gm_models.size() << " GM/CSF kernels";
+ ++it;
}
-
+ break;
+ }
+ case (FiberfoxParameters<>::RANDOM_DIRECTIONS):
+ {
ItkUcharImgType::Pointer numDirectionsImage = ItkUcharImgType::New();
numDirectionsImage->SetSpacing( m_UpsampledSpacing );
numDirectionsImage->SetOrigin( m_UpsampledOrigin );
numDirectionsImage->SetDirection( m_Parameters.m_ImageDirection );
numDirectionsImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
numDirectionsImage->SetBufferedRegion( m_UpsampledImageRegion );
numDirectionsImage->SetRequestedRegion( m_UpsampledImageRegion );
numDirectionsImage->Allocate();
numDirectionsImage->FillBuffer(0);
- ItkDoubleImgType::Pointer intraAxonalVolumeImage = ItkDoubleImgType::New();
- intraAxonalVolumeImage->SetSpacing( m_UpsampledSpacing );
- intraAxonalVolumeImage->SetOrigin( m_UpsampledOrigin );
- intraAxonalVolumeImage->SetDirection( m_Parameters.m_ImageDirection );
- intraAxonalVolumeImage->SetLargestPossibleRegion( m_UpsampledImageRegion );
- intraAxonalVolumeImage->SetBufferedRegion( m_UpsampledImageRegion );
- intraAxonalVolumeImage->SetRequestedRegion( m_UpsampledImageRegion );
- intraAxonalVolumeImage->Allocate();
- intraAxonalVolumeImage->FillBuffer(0);
-
m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n";
m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
boost::progress_display disp(tempTissueMask->GetLargestPossibleRegion().GetNumberOfPixels());
ImageRegionIterator it(tempTissueMask, tempTissueMask->GetLargestPossibleRegion());
while(!it.IsAtEnd())
{
++disp;
-
unsigned long newTick = 50*disp.count()/disp.expected_count();
for (unsigned int tick = 0; tick<(newTick-lastTick); tick++)
m_StatusText += "*";
lastTick = newTick;
if (this->GetAbortGenerateData())
{
m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
return;
}
if (it.Get()>0)
{
int numFibs = m_RandGen->GetIntegerVariate(2)+1;
- DoubleDwiType::PixelType pix = compartments.at(0)->GetPixel(it.GetIndex());
+ DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(it.GetIndex());
+ double volume = m_RandGen->GetVariateWithClosedRange(0.3);
+
+ double sum = 0;
+ std::vector< double > fractions;
+ for (int i=0; iGetVariateWithClosedRange(0.5));
+ sum += fractions.at(i);
+ }
+ for (int i=0; i > directions;
for (int i=0; iGetIntegerVariate(wm_models.size()-1);
- itk::Vector fib;
+ DoubleVectorType fib;
fib[0] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
fib[1] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
fib[2] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
fib.Normalize();
double min = 0;
for (unsigned int d=0; dmin)
min = angle;
}
- if (min<0.7)
+ if (min<0.5)
{
- wm_models.at(modelIndex)->SetFiberDirection(fib);
- pix += wm_models.at(modelIndex)->SimulateMeasurement()/numFibs;
+ m_Parameters.m_FiberModelList.at(0)->SetFiberDirection(fib);
+ pix += m_Parameters.m_FiberModelList.at(0)->SimulateMeasurement()*fractions[i];
directions.push_back(fib);
}
else
i--;
}
- compartments.at(0)->SetPixel(it.GetIndex(), pix);
+ pix *= (1-volume);
+
+ // CSF/GM
+ {
+// int modelIndex = m_RandGen->GetIntegerVariate(m_Parameters.m_NonFiberModelList.size()-1);
+ pix += volume*m_Parameters.m_NonFiberModelList.at(0)->SimulateMeasurement();
+ }
+
+ m_CompartmentImages.at(0)->SetPixel(it.GetIndex(), pix);
numDirectionsImage->SetPixel(it.GetIndex(), numFibs);
- // int numDirs = 0;
- // double volume = m_RandGen->GetVariateWithClosedRange(voxelVolume);
- // int numFibs = m_RandGen->GetIntegerVariate(2)+1;
- // std::vector< double > fractions;
- // double sum = 0;
- // for (int i=0; iGetVariateWithClosedRange(1));
- // sum += fractions.at(i);
- // }
-
- // std::vector< itk::Vector > directions;
- // for (int i=0; i0)
- // fractions[i] /= sum;
-
- // itk::Vector dir;
- // dir[0] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
- // dir[1] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
- // dir[2] = m_RandGen->GetVariateWithClosedRange(2)-1.0;
- // dir.Normalize();
-
- // if (fractions.at(i)*numFibs > 0.3 && volume/voxelVolume>0.01)
- // {
- // numDirs += 1;
- // for (unsigned int d=0; dSetFiberDirection(dir);
- // DoubleDwiType::PixelType pix = compartments.at(k)->GetPixel(it2.GetIndex());
- // pix += volume*fractions.at(i)*m_Parameters.m_FiberModelList[k]->SimulateMeasurement();
- // compartments.at(k)->SetPixel(it2.GetIndex(), pix);
- // }
- // }
-
- // intraAxonalVolumeImage->SetPixel(it2.GetIndex(), volume);
- // if (volume>maxVolume)
- // maxVolume = volume;
-
- // // for (unsigned int i=0; i::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New();
- wr->SetInput(numDirectionsImage);
- wr->SetFileName(mitk::IOUtil::GetTempPath()+"/NumDirections.nrrd");
- wr->Update();
}
+
+ itk::ImageFileWriter< ItkUcharImgType >::Pointer wr = itk::ImageFileWriter< ItkUcharImgType >::New();
+ wr->SetInput(numDirectionsImage);
+ wr->SetFileName("/local/NumDirections.nrrd");
+ wr->Update();
+ }
}
if (logFile.is_open())
{
logFile << "DONE";
logFile.close();
}
m_StatusText += "\n\n";
if (this->GetAbortGenerateData())
{
m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
return;
}
// do k-space stuff
DoubleDwiType::Pointer doubleOutImage;
- if ( !m_NoAcquisitionSimulation && (m_Parameters.m_Spikes>0 || m_Parameters.m_FrequencyMap.IsNotNull() || m_Parameters.m_KspaceLineOffset>0 || m_Parameters.m_DoSimulateRelaxation || m_Parameters.m_EddyStrength>0 || m_Parameters.m_DoAddGibbsRinging || m_Parameters.m_CroppingFactor<1.0) )
+ if ( m_Parameters.m_SimulateKspaceAcquisition )
{
m_StatusText += this->GetTime()+" > Adjusting complex signal\n";
MITK_INFO << "Adjusting complex signal:";
if (m_Parameters.m_DoSimulateRelaxation)
m_StatusText += "Simulating signal relaxation\n";
if (m_Parameters.m_FrequencyMap.IsNotNull())
m_StatusText += "Simulating distortions\n";
if (m_Parameters.m_DoAddGibbsRinging)
m_StatusText += "Simulating ringing artifacts\n";
if (m_Parameters.m_EddyStrength>0)
m_StatusText += "Simulating eddy currents\n";
if (m_Parameters.m_Spikes>0)
m_StatusText += "Simulating spikes\n";
if (m_Parameters.m_CroppingFactor<1.0)
m_StatusText += "Simulating aliasing artifacts\n";
if (m_Parameters.m_KspaceLineOffset>0)
m_StatusText += "Simulating ghosts\n";
- doubleOutImage = DoKspaceStuff(compartments);
- m_Parameters.m_SignalScale = 1;
+ doubleOutImage = DoKspaceStuff(m_CompartmentImages);
+ m_Parameters.m_SignalScale = 1; // already scaled in DoKspaceStuff
}
else
{
m_StatusText += this->GetTime()+" > Summing compartments\n";
MITK_INFO << "Summing compartments";
- doubleOutImage = compartments.at(0);
+ doubleOutImage = m_CompartmentImages.at(0);
- for (unsigned int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New();
adder->SetInput1(doubleOutImage);
- adder->SetInput2(compartments.at(i));
+ adder->SetInput2(m_CompartmentImages.at(i));
adder->Update();
doubleOutImage = adder->GetOutput();
}
}
if (this->GetAbortGenerateData())
{
m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
return;
}
m_StatusText += this->GetTime()+" > Finalizing image\n";
MITK_INFO << "Finalizing image";
if (m_Parameters.m_SignalScale>1)
m_StatusText += " Scaling signal\n";
if (m_Parameters.m_NoiseModel!=NULL)
m_StatusText += " Adding noise\n";
unsigned int window = 0;
unsigned int min = itk::NumericTraits::max();
ImageRegionIterator it4 (outImage, outImage->GetLargestPossibleRegion());
DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.GetNumVolumes());
boost::progress_display disp2(outImage->GetLargestPossibleRegion().GetNumberOfPixels());
m_StatusText += "0% 10 20 30 40 50 60 70 80 90 100%\n";
m_StatusText += "|----|----|----|----|----|----|----|----|----|----|\n*";
lastTick = 0;
while(!it4.IsAtEnd())
{
if (this->GetAbortGenerateData())
{
m_StatusText += "\n"+this->GetTime()+" > Simulation aborted\n";
return;
}
++disp2;
unsigned long newTick = 50*disp2.count()/disp2.expected_count();
for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
m_StatusText += "*";
lastTick = newTick;
typename OutputImageType::IndexType index = it4.GetIndex();
signal = doubleOutImage->GetPixel(index)*m_Parameters.m_SignalScale;
if (m_Parameters.m_NoiseModel!=NULL)
- {
- DoubleDwiType::PixelType accu = signal; accu.Fill(0.0);
- for (unsigned int i=0; iAddNoise(temp);
- accu += temp;
- }
- signal = accu/m_Parameters.m_Repetitions;
- }
+ m_Parameters.m_NoiseModel->AddNoise(signal);
for (unsigned int i=0; i0)
signal[i] = floor(signal[i]+0.5);
else
signal[i] = ceil(signal[i]-0.5);
- if (!m_Parameters.IsBaselineIndex(i) && signal[i]>window)
+ if ( (!m_Parameters.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window)
window = signal[i];
- if (!m_Parameters.IsBaselineIndex(i) && signal[i]SetNthOutput(0, outImage);
m_StatusText += "\n\n";
m_StatusText += "Finished simulation\n";
m_StatusText += "Simulation time: "+GetTime();
m_TimeProbe.Stop();
}
template< class PixelType >
itk::Point TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3])
{
itk::Point itkPoint;
itkPoint[0] = point[0];
itkPoint[1] = point[1];
itkPoint[2] = point[2];
return itkPoint;
}
template< class PixelType >
itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3])
{
itk::Vector itkVector;
itkVector[0] = point[0];
itkVector[1] = point[1];
itkVector[2] = point[2];
return itkVector;
}
template< class PixelType >
vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3])
{
vnl_vector_fixed vnlVector;
vnlVector[0] = point[0];
vnlVector[1] = point[1];
vnlVector[2] = point[2];
return vnlVector;
}
template< class PixelType >
vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector)
{
vnl_vector_fixed vnlVector;
vnlVector[0] = vector[0];
vnlVector[1] = vector[1];
vnlVector[2] = vector[2];
return vnlVector;
}
template< class PixelType >
double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) {
return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
}
template< class PixelType >
std::string TractsToDWIImageFilter< PixelType >::GetTime()
{
m_TimeProbe.Stop();
unsigned long total = RoundToNearest(m_TimeProbe.GetTotal());
unsigned long hours = total/3600;
unsigned long minutes = (total%3600)/60;
unsigned long seconds = total%60;
std::string out = "";
out.append(boost::lexical_cast(hours));
out.append(":");
out.append(boost::lexical_cast(minutes));
out.append(":");
out.append(boost::lexical_cast(seconds));
m_TimeProbe.Start();
return out;
}
}
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h
index cf7303c01c..dfb640ad23 100755
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.h
@@ -1,112 +1,111 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef __itkTractsToDWIImageFilter_h__
#define __itkTractsToDWIImageFilter_h__
#include
#include
#include
#include
#include
#include
#include
#include
#include
namespace itk
{
/**
* \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model.
* See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details.
*/
template< class PixelType >
class TractsToDWIImageFilter : public ImageSource< itk::VectorImage< PixelType, 3 > >
{
public:
typedef TractsToDWIImageFilter Self;
typedef ImageSource< itk::VectorImage< PixelType, 3 > > Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
typedef typename Superclass::OutputImageType OutputImageType;
typedef itk::Image ItkDoubleImgType;
typedef itk::Image ItkUcharImgType;
typedef mitk::FiberBundleX::Pointer FiberBundleType;
typedef itk::VectorImage< double, 3 > DoubleDwiType;
typedef itk::Matrix MatrixType;
typedef itk::Image< double, 2 > SliceType;
typedef itk::VnlForwardFFTImageFilter::OutputImageType ComplexSliceType;
- typedef itk::Vector< double,3> VectorType;
+ typedef itk::Vector< double,3> DoubleVectorType;
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
itkTypeMacro( TractsToDWIImageFilter, ImageSource )
/** Input */
itkSetMacro( FiberBundle, FiberBundleType ) ///< Input fiber bundle
itkSetMacro( UseConstantRandSeed, bool ) ///< Seed for random generator.
- itkSetMacro( InputDwi, mitk::DiffusionImage::Pointer )
void SetParameters( FiberfoxParameters param ) ///< Simulation parameters.
{ m_Parameters = param; }
/** Output */
FiberfoxParameters GetParameters(){ return m_Parameters; }
std::vector< ItkDoubleImgType::Pointer > GetVolumeFractions() ///< one double image for each compartment containing the corresponding volume fraction per voxel
{ return m_VolumeFractions; }
mitk::LevelWindow GetLevelWindow(){ return m_LevelWindow; }
itkGetMacro( StatusText, std::string )
void GenerateData();
protected:
TractsToDWIImageFilter();
virtual ~TractsToDWIImageFilter();
itk::Point GetItkPoint(double point[3]);
itk::Vector GetItkVector(double point[3]);
vnl_vector_fixed GetVnlVector(double point[3]);
vnl_vector_fixed GetVnlVector(Vector< float, 3 >& vector);
double RoundToNearest(double num);
std::string GetTime();
/** Transform generated image compartment by compartment, channel by channel and slice by slice using DFT and add k-space artifacts. */
DoubleDwiType::Pointer DoKspaceStuff(std::vector< DoubleDwiType::Pointer >& images);
mitk::FiberfoxParameters m_Parameters;
itk::Vector m_UpsampledSpacing;
itk::Point m_UpsampledOrigin;
ImageRegion<3> m_UpsampledImageRegion;
FiberBundleType m_FiberBundle;
mitk::LevelWindow m_LevelWindow;
std::vector< ItkDoubleImgType::Pointer > m_VolumeFractions;
std::string m_StatusText;
itk::TimeProbe m_TimeProbe;
bool m_UseConstantRandSeed;
- mitk::DiffusionImage::Pointer m_InputDwi;
- bool m_NoAcquisitionSimulation;
itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen;
+
+ std::vector< DoubleDwiType::Pointer > m_CompartmentImages;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTractsToDWIImageFilter.cpp"
#endif
#endif
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
index ef8b897311..6c437a400b 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp
@@ -1,835 +1,556 @@
#include "itkTractsToVectorImageFilter.h"
// VTK
#include
#include
#include
// ITK
#include
#include
// misc
#define _USE_MATH_DEFINES
#include
#include
namespace itk{
static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2)
{
return (v1.magnitude()>v2.magnitude());
}
template< class PixelType >
TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter():
m_AngularThreshold(0.7),
m_Epsilon(0.999),
m_MaskImage(NULL),
m_NormalizeVectors(false),
m_UseWorkingCopy(true),
- m_UseTrilinearInterpolation(false),
m_MaxNumDirections(3),
- m_Thres(0.5),
+ m_SizeThreshold(0.2),
m_NumDirectionsImage(NULL)
{
this->SetNumberOfRequiredOutputs(1);
}
template< class PixelType >
TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter()
{
}
template< class PixelType >
vnl_vector_fixed TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[])
{
vnl_vector_fixed vnlVector;
vnlVector[0] = point[0];
vnlVector[1] = point[1];
vnlVector[2] = point[2];
return vnlVector;
}
template< class PixelType >
itk::Point TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[])
{
itk::Point itkPoint;
itkPoint[0] = point[0];
itkPoint[1] = point[1];
itkPoint[2] = point[2];
return itkPoint;
}
template< class PixelType >
void TractsToVectorImageFilter< PixelType >::GenerateData()
{
- mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
+ mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
// calculate new image parameters
itk::Vector spacing;
itk::Point origin;
itk::Matrix direction;
ImageRegion<3> imageRegion;
if (!m_MaskImage.IsNull())
{
spacing = m_MaskImage->GetSpacing();
imageRegion = m_MaskImage->GetLargestPossibleRegion();
origin = m_MaskImage->GetOrigin();
direction = m_MaskImage->GetDirection();
}
else
{
spacing = geometry->GetSpacing();
origin = geometry->GetOrigin();
mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds();
origin[0] += bounds.GetElement(0);
origin[1] += bounds.GetElement(2);
origin[2] += bounds.GetElement(4);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
direction[j][i] = geometry->GetMatrixColumn(i)[j];
imageRegion.SetSize(0, geometry->GetExtent(0));
imageRegion.SetSize(1, geometry->GetExtent(1));
imageRegion.SetSize(2, geometry->GetExtent(2));
m_MaskImage = ItkUcharImgType::New();
m_MaskImage->SetSpacing( spacing );
m_MaskImage->SetOrigin( origin );
m_MaskImage->SetDirection( direction );
m_MaskImage->SetRegions( imageRegion );
m_MaskImage->Allocate();
m_MaskImage->FillBuffer(1);
}
OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize();
m_OutImageSpacing = m_MaskImage->GetSpacing();
m_ClusteredDirectionsContainer = ContainerType::New();
- // initialize crossings image
- m_CrossingsImage = ItkUcharImgType::New();
- m_CrossingsImage->SetSpacing( spacing );
- m_CrossingsImage->SetOrigin( origin );
- m_CrossingsImage->SetDirection( direction );
- m_CrossingsImage->SetRegions( imageRegion );
- m_CrossingsImage->Allocate();
- m_CrossingsImage->FillBuffer(0);
-
// initialize num directions image
m_NumDirectionsImage = ItkUcharImgType::New();
m_NumDirectionsImage->SetSpacing( spacing );
m_NumDirectionsImage->SetOrigin( origin );
m_NumDirectionsImage->SetDirection( direction );
m_NumDirectionsImage->SetRegions( imageRegion );
m_NumDirectionsImage->Allocate();
m_NumDirectionsImage->FillBuffer(0);
+ // initialize direction images
+ m_DirectionImageContainer = DirectionImageContainerType::New();
+
// resample fiber bundle
double minSpacing = 1;
if(m_OutImageSpacing[0]GetDeepCopy();
// resample fiber bundle for sufficient voxel coverage
- m_FiberBundle->ResampleFibers(minSpacing/3);
+ m_FiberBundle->ResampleFibers(minSpacing/10);
// iterate over all fibers
vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData();
- vtkSmartPointer vLines = fiberPolyData->GetLines();
- vLines->InitTraversal();
int numFibers = m_FiberBundle->GetNumFibers();
- itk::TimeProbe clock;
m_DirectionsContainer = ContainerType::New();
- if (m_UseTrilinearInterpolation)
- MITK_INFO << "Generating directions from tractogram (trilinear interpolation)";
- else
- MITK_INFO << "Generating directions from tractogram";
-
+ MITK_INFO << "Generating directions from tractogram";
boost::progress_display disp(numFibers);
for( int i=0; iGetNextCell ( numPoints, points );
+ vtkCell* cell = fiberPolyData->GetCell(i);
+ int numPoints = cell->GetNumberOfPoints();
+ vtkPoints* points = cell->GetPoints();
if (numPoints<2)
continue;
- itk::Index<3> index; index.Fill(0);
- itk::ContinuousIndex contIndex;
- vnl_vector_fixed dir, wDir;
+ vnl_vector_fixed