diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp
index a14e4a983f..e9142ce4a0 100644
--- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp
@@ -1,133 +1,145 @@
/*===================================================================
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 <mitkImageCast.h>
#include <mitkImage.h>
#include <mitkIOUtil.h>
#include "mitkCommandLineParser.h"
#include <mitkMultiModalRigidDefaultRegistrationAlgorithm.h>
#include <mitkDiffusionPropertyHelper.h>
#include <mitkAlgorithmHelper.h>
#include <itkExtractDwiChannelFilter.h>
#include <mitkRegistrationHelper.h>
#include <mitkImageMappingHelper.h>
#include <itksys/SystemTools.hxx>
typedef mitk::DiffusionPropertyHelper DPH;
/*!
\brief Copies transformation matrix of one image to another
*/
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("DREG");
parser.setCategory("Preprocessing Tools");
parser.setDescription("TEMPORARY: Rigid registration of two images");
parser.setContributor("MIC");
parser.setArgumentPrefix("--", "-");
parser.addArgument("", "f", mitkCommandLineParser::InputFile, "Fixed:", "fixed image", us::Any(), false);
parser.addArgument("", "m", mitkCommandLineParser::InputFile, "Moving:", "moving image", us::Any(), false);
+ parser.addArgument("", "c", mitkCommandLineParser::StringList, "Coregister:", "apply transform to these images too", us::Any(), false);
parser.addArgument("", "o", mitkCommandLineParser::OutputFile, "Output:", "output image", us::Any(), false);
std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
// mandatory arguments
std::string f = us::any_cast<std::string>(parsedArgs["f"]);
std::string m = us::any_cast<std::string>(parsedArgs["m"]);
std::string o = us::any_cast<std::string>(parsedArgs["o"]);
+ mitkCommandLineParser::StringContainerType coregister;
+ if (parsedArgs.count("c"))
+ coregister = us::any_cast<mitkCommandLineParser::StringContainerType>(parsedArgs["c"]);
+
try
{
typedef itk::Image< float, 3 > ItkFloatImageType;
mitk::Image::Pointer fixed = mitk::IOUtil::Load<mitk::Image>(f);
mitk::Image::Pointer moving = mitk::IOUtil::Load<mitk::Image>(m);
mitk::Image::Pointer fixed_single = fixed;
mitk::Image::Pointer moving_single = moving;
mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::Pointer algo = mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::New();
mitk::MITKAlgorithmHelper helper(algo);
if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(fixed))
{
DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New();
mitk::CastToItkImage(fixed, itkVectorImagePointer);
itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New();
filter->SetInput( itkVectorImagePointer);
filter->SetChannelIndex(0);
filter->Update();
fixed_single = mitk::Image::New();
fixed_single->InitializeByItk( filter->GetOutput() );
fixed_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() );
}
if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(moving))
{
DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New();
mitk::CastToItkImage(moving, itkVectorImagePointer);
itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New();
filter->SetInput( itkVectorImagePointer);
filter->SetChannelIndex(0);
filter->Update();
moving_single = mitk::Image::New();
moving_single->InitializeByItk( filter->GetOutput() );
moving_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() );
}
helper.SetData(moving_single, fixed_single);
mitk::MAPRegistrationWrapper::Pointer reg = helper.GetMITKRegistrationWrapper();
mitk::Image::Pointer registered_image = mitk::ImageMappingHelper::refineGeometry(moving, reg, true);
+ for (auto co : coregister)
+ {
+ std::string out = itksys::SystemTools::GetFilenamePath(o) + "/" + itksys::SystemTools::GetFilenameName(co);
+ mitk::Image::Pointer co_image = mitk::IOUtil::Load<mitk::Image>(co);
+ mitk::IOUtil::Save(mitk::ImageMappingHelper::refineGeometry(co_image, reg, true), out);
+ }
+
if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(registered_image))
{
mitk::DiffusionPropertyHelper propertyHelper( registered_image );
propertyHelper.InitializeImage();
std::string file_extension = itksys::SystemTools::GetFilenameExtension(o);
if (file_extension==".nii" || file_extension==".nii.gz")
mitk::IOUtil::Save(registered_image, "application/vnd.mitk.nii.gz", o);
else
mitk::IOUtil::Save(registered_image, o);
}
else
mitk::IOUtil::Save(registered_image, o);
}
catch (itk::ExceptionObject e)
{
std::cout << e;
return EXIT_FAILURE;
}
catch (std::exception e)
{
std::cout << e.what();
return EXIT_FAILURE;
}
catch (...)
{
std::cout << "ERROR!?!";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
index 376a5d10cd..5e8a7ebb13 100755
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
@@ -1,1727 +1,1729 @@
/*===================================================================
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 <boost/progress.hpp>
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkPoints.h>
#include <vtkPolyLine.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkResampleImageFilter.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <itkBSplineInterpolateImageFunction.h>
#include <itkCastImageFilter.h>
#include <itkImageFileWriter.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkWindowedSincInterpolateImageFunction.h>
#include <itkResampleDwiImageFilter.h>
#include <itkKspaceImageFilter.h>
#include <itkDftImageFilter.h>
#include <itkAddImageFilter.h>
#include <itkConstantPadImageFilter.h>
#include <itkCropImageFilter.h>
#include <mitkAstroStickModel.h>
#include <vtkTransform.h>
#include <iostream>
#include <fstream>
#include <exception>
#include <itkImageDuplicator.h>
#include <itksys/SystemTools.hxx>
#include <mitkIOUtil.h>
#include <itkDiffusionTensor3DReconstructionImageFilter.h>
#include <itkDiffusionTensor3D.h>
#include <itkTractDensityImageFilter.h>
#include <mitkLexicalCast.h>
#include <itkTractsToVectorImageFilter.h>
#include <itkInvertIntensityImageFilter.h>
#include <itkShiftScaleImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkResampleDwiImageFilter.h>
#include <boost/algorithm/string/replace.hpp>
#include <omp.h>
namespace itk
{
template< class PixelType >
TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter()
: m_StatusText("")
, m_UseConstantRandSeed(false)
, m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New())
{
m_RandGen->SetSeed();
m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New();
m_NullDir.Fill(0);
}
template< class PixelType >
TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter()
{
}
template< class PixelType >
TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::
SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images )
{
unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size();
// create slice object
ImageRegion<2> sliceRegion;
sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]);
sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]);
Vector< double, 2 > sliceSpacing;
sliceSpacing[0] = m_WorkingSpacing[0];
sliceSpacing[1] = m_WorkingSpacing[1];
DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0);
auto magnitudeDwiImage = DoubleDwiType::New();
magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() );
magnitudeDwiImage->Allocate();
magnitudeDwiImage->FillBuffer(nullPix);
m_PhaseImage = DoubleDwiType::New();
m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() );
m_PhaseImage->Allocate();
m_PhaseImage->FillBuffer(nullPix);
m_KspaceImage = DoubleDwiType::New();
m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils );
m_KspaceImage->Allocate();
m_KspaceImage->FillBuffer(nullPix);
std::vector< unsigned int > spikeVolume;
for (unsigned int i=0; i<m_Parameters.m_SignalGen.m_Spikes; i++)
spikeVolume.push_back(m_RandGen->GetIntegerVariate()%(compartment_images.at(0)->GetVectorLength()));
std::sort (spikeVolume.begin(), spikeVolume.end());
std::reverse (spikeVolume.begin(), spikeVolume.end());
// calculate coil positions
double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0];
double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1];
double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2];
double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m
m_CoilPointset = mitk::PointSet::New();
std::vector< itk::Vector<double, 3> > coilPositions;
itk::Vector<double, 3> pos; pos.Fill(0.0); pos[1] = -diagonal/2;
itk::Vector<double, 3> center;
center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
for (int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
{
coilPositions.push_back(pos);
m_CoilPointset->InsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center );
double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * itk::Math::pi/180;
vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
rotZ[0][0] = cos(rz);
rotZ[1][1] = rotZ[0][0];
rotZ[0][1] = -sin(rz);
rotZ[1][0] = -rotZ[0][1];
pos.SetVnlVector(rotZ*pos.GetVnlVector());
}
PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
unsigned long lastTick = 0;
boost::progress_display disp(compartment_images.at(0)->GetVectorLength()*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2));
#pragma omp parallel for
for (int g=0; g<(int)compartment_images.at(0)->GetVectorLength(); g++)
{
if (this->GetAbortGenerateData())
continue;
std::vector< unsigned int > spikeSlice;
#pragma omp critical
while (!spikeVolume.empty() && (int)spikeVolume.back()==g)
{
spikeSlice.push_back(m_RandGen->GetIntegerVariate()%compartment_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; z<compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2); z++)
{
std::vector< Float2DImageType::Pointer > compartment_slices;
std::vector< float > t2Vector;
std::vector< float > t1Vector;
for (unsigned int i=0; i<compartment_images.size(); i++)
{
DiffusionSignalModel<double>* signalModel;
if (i<numFiberCompartments)
signalModel = m_Parameters.m_FiberModelList.at(i);
else
signalModel = m_Parameters.m_NonFiberModelList.at(i-numFiberCompartments);
auto slice = Float2DImageType::New();
slice->SetLargestPossibleRegion( 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; y<compartment_images.at(0)->GetLargestPossibleRegion().GetSize(1); y++)
for (unsigned int x=0; x<compartment_images.at(0)->GetLargestPossibleRegion().GetSize(0); x++)
{
Float2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y;
DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]);
}
compartment_slices.push_back(slice);
t2Vector.push_back(signalModel->GetT2());
t1Vector.push_back(signalModel->GetT1());
}
int numSpikes = 0;
while (!spikeSlice.empty() && spikeSlice.back()==z)
{
numSpikes++;
spikeSlice.pop_back();
}
int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils;
if (this->GetAbortGenerateData())
continue;
for (int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
{
// create k-sapce (inverse fourier transform slices)
auto idft = itk::KspaceImageFilter< Float2DImageType::PixelType >::New();
idft->SetCompartmentImages(compartment_slices);
idft->SetT2(t2Vector);
idft->SetT1(t1Vector);
idft->SetUseConstantRandSeed(m_UseConstantRandSeed);
idft->SetParameters(&m_Parameters);
idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)
-compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0);
idft->SetZidx(z);
idft->SetCoilPosition(coilPositions.at(c));
idft->SetFiberBundle(m_FiberBundle);
idft->SetTranslation(m_Translations.at(g));
idft->SetRotation(m_Rotations.at(g));
idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g));
if (c==spikeCoil)
idft->SetSpikesPerSlice(numSpikes);
idft->Update();
#pragma omp critical
if (c==spikeCoil && numSpikes>0)
{
m_SpikeLog += "Volume " + boost::lexical_cast<std::string>(g) + " Coil " + boost::lexical_cast<std::string>(c) + "\n";
m_SpikeLog += idft->GetSpikeLog();
}
Complex2DImageType::Pointer fSlice;
fSlice = idft->GetOutput();
// fourier transform slice
Complex2DImageType::Pointer newSlice;
auto dft = itk::DftImageFilter< Float2DImageType::PixelType >::New();
dft->SetInput(fSlice);
dft->SetParameters(m_Parameters);
dft->Update();
newSlice = dft->GetOutput();
// put slice back into channel g
for (unsigned int y=0; y<fSlice->GetLargestPossibleRegion().GetSize(1); y++)
for (unsigned int x=0; x<fSlice->GetLargestPossibleRegion().GetSize(0); x++)
{
DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y;
Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D);
double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag());
double phase = 0;
if (cPix.real()!=0)
phase = atan( cPix.imag()/cPix.real() );
DoubleDwiType::PixelType real_pix = m_OutputImagesReal.at(c)->GetPixel(index3D);
real_pix[g] = cPix.real();
m_OutputImagesReal.at(c)->SetPixel(index3D, real_pix);
DoubleDwiType::PixelType imag_pix = m_OutputImagesImag.at(c)->GetPixel(index3D);
imag_pix[g] = cPix.imag();
m_OutputImagesImag.at(c)->SetPixel(index3D, imag_pix);
DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D);
DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
{
dwiPix[g] += magn*magn;
phasePix[g] += phase*phase;
}
else
{
dwiPix[g] = magn;
phasePix[g] = phase;
}
//#pragma omp critical
{
magnitudeDwiImage->SetPixel(index3D, dwiPix);
m_PhaseImage->SetPixel(index3D, phasePix);
// k-space image
if (g==0)
{
DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D);
kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D);
m_KspaceImage->SetPixel(index3D, kspacePix);
}
}
}
}
if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
{
for (int y=0; y<static_cast<int>(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++)
for (int x=0; x<static_cast<int>(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++)
{
DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D);
magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
//#pragma omp critical
{
magnitudeDwiImage->SetPixel(index3D, magPix);
m_PhaseImage->SetPixel(index3D, phasePix);
}
}
}
++disp;
unsigned long newTick = 50*disp.count()/disp.expected_count();
for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
PrintToLog("*", false, false, false);
lastTick = newTick;
}
}
PrintToLog("\n", false);
return magnitudeDwiImage;
}
template< class PixelType >
TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >::
NormalizeInsideMask(ItkDoubleImgType::Pointer image)
{
double max = itk::NumericTraits< double >::min();
double min = itk::NumericTraits< double >::max();
itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion());
while(!it.IsAtEnd())
{
if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0)
{
it.Set(0.0);
++it;
continue;
}
if (it.Get()>max)
max = it.Get();
if (it.Get()<min)
min = it.Get();
++it;
}
auto scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
scaler->SetInput(image);
scaler->SetShift(-min);
scaler->SetScale(1.0/(max-min));
scaler->Update();
return scaler->GetOutput();
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages()
{
m_UseRelativeNonFiberVolumeFractions = false;
// check for fiber volume fraction maps
unsigned int fibVolImages = 0;
for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size(); i++)
{
if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage().IsNotNull())
{
PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast<std::string>(i+1));
fibVolImages++;
}
}
// check for non-fiber volume fraction maps
unsigned int nonfibVolImages = 0;
for (std::size_t i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
{
if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage().IsNotNull())
{
PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast<std::string>(i+1));
nonfibVolImages++;
}
}
// not all fiber compartments are using volume fraction maps
// --> non-fiber volume fractions are assumed to be relative to the
// non-fiber volume and not absolute voxel-volume fractions.
// this means if two non-fiber compartments are used but only one of them
// has an associated volume fraction map, the repesctive other volume fraction map
// can be determined as inverse (1-val) of the present volume fraction map-
if ( fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages==1 && m_Parameters.m_NonFiberModelList.size()==2)
{
PrintToLog("Calculating missing non-fiber volume fraction image by inverting the other.\n"
"Assuming non-fiber volume fraction images to contain values relative to the"
" remaining non-fiber volume, not absolute values.");
auto inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
inverter->SetMaximum(1.0);
if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull()
&& m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() )
{
// m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(
// NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) );
inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() );
inverter->Update();
m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput());
}
else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull()
&& m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() )
{
// m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(
// NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) );
inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() );
inverter->Update();
m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput());
}
else
{
itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!"
" Did you use two non fiber compartments but only one volume fraction image?"
" Then it should work and this error is really strange.");
}
m_UseRelativeNonFiberVolumeFractions = true;
nonfibVolImages++;
}
// Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically
if (m_Parameters.m_FiberModelList.size()>2
&& fibVolImages!=m_Parameters.m_FiberModelList.size())
itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!");
// One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically
if (m_Parameters.m_NonFiberModelList.size()>1
&& nonfibVolImages!=m_Parameters.m_NonFiberModelList.size())
itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!");
if (fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages>0)
{
PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n"
"Assuming non-fiber volume fraction images to contain values relative to the"
" remaining non-fiber volume, not absolute values.");
m_UseRelativeNonFiberVolumeFractions = true;
// itk::ImageFileWriter<ItkDoubleImgType>::Pointer wr = itk::ImageFileWriter<ItkDoubleImgType>::New();
// wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage());
// wr->SetFileName("/local/volumefraction.nrrd");
// wr->Update();
}
// initialize the images that store the output volume fraction of each compartment
m_VolumeFractions.clear();
for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size(); i++)
{
auto doubleImg = ItkDoubleImgType::New();
doubleImg->SetSpacing( m_WorkingSpacing );
doubleImg->SetOrigin( m_WorkingOrigin );
doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion );
doubleImg->SetBufferedRegion( m_WorkingImageRegion );
doubleImg->SetRequestedRegion( m_WorkingImageRegion );
doubleImg->Allocate();
doubleImg->FillBuffer(0);
m_VolumeFractions.push_back(doubleImg);
}
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::InitializeData()
{
m_Rotations.clear();
m_Translations.clear();
m_MotionLog = "";
m_SpikeLog = "";
// initialize output dwi image
m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion;
m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)
*m_Parameters.m_SignalGen.m_CroppingFactor);
itk::Point<double,3> shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)
-m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
m_OutputImage = OutputImageType::New();
m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
m_OutputImage->SetOrigin( shiftedOrigin );
m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
m_OutputImage->Allocate();
typename OutputImageType::PixelType temp;
temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
temp.Fill(0.0);
m_OutputImage->FillBuffer(temp);
// images containing real and imaginary part of the dMRI signal for each coil
m_OutputImagesReal.clear();
m_OutputImagesImag.clear();
for (int i=0; i<m_Parameters.m_SignalGen.m_NumberOfCoils; ++i)
{
typename DoubleDwiType::Pointer outputImageReal = DoubleDwiType::New();
outputImageReal->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
outputImageReal->SetOrigin( shiftedOrigin );
outputImageReal->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
outputImageReal->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageReal->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageReal->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageReal->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
outputImageReal->Allocate();
outputImageReal->FillBuffer(temp);
m_OutputImagesReal.push_back(outputImageReal);
typename DoubleDwiType::Pointer outputImageImag = DoubleDwiType::New();
outputImageImag->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
outputImageImag->SetOrigin( shiftedOrigin );
outputImageImag->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
outputImageImag->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageImag->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageImag->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
outputImageImag->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
outputImageImag->Allocate();
outputImageImag->FillBuffer(temp);
m_OutputImagesImag.push_back(outputImageImag);
}
// Apply in-plane upsampling for Gibbs ringing artifact
double upsampling = 1;
if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
upsampling = 2;
m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing;
m_WorkingSpacing[0] /= upsampling;
m_WorkingSpacing[1] /= upsampling;
m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion;
m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling);
m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling);
m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
m_WorkingOrigin[0] += m_WorkingSpacing[0]/2;
m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
m_WorkingOrigin[1] += m_WorkingSpacing[1]/2;
m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
m_WorkingOrigin[2] += m_WorkingSpacing[2]/2;
m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2];
// generate double images to store the individual compartment signals
m_CompartmentImages.clear();
int numFiberCompartments = m_Parameters.m_FiberModelList.size();
int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
for (int i=0; i<numFiberCompartments+numNonFiberCompartments; i++)
{
auto doubleDwi = DoubleDwiType::New();
doubleDwi->SetSpacing( m_WorkingSpacing );
doubleDwi->SetOrigin( m_WorkingOrigin );
doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion );
doubleDwi->SetBufferedRegion( m_WorkingImageRegion );
doubleDwi->SetRequestedRegion( m_WorkingImageRegion );
doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
doubleDwi->Allocate();
DoubleDwiType::PixelType pix;
pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
pix.Fill(0.0);
doubleDwi->FillBuffer(pix);
m_CompartmentImages.push_back(doubleDwi);
}
if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull())
{
m_CompartmentImages.clear();
m_Parameters.m_SignalGen.m_DoAddMotion = false;
m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false;
PrintToLog("Simulating acquisition for input diffusion-weighted image.", false);
auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New();
caster->SetInput(m_InputImage);
caster->Update();
if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
{
PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false);
auto resampler = itk::ResampleDwiImageFilter< double >::New();
resampler->SetInput(caster->GetOutput());
itk::Vector< double, 3 > samplingFactor;
samplingFactor[0] = upsampling;
samplingFactor[1] = upsampling;
samplingFactor[2] = 1;
resampler->SetSamplingFactor(samplingFactor);
resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc);
resampler->Update();
m_CompartmentImages.push_back(resampler->GetOutput());
}
else
m_CompartmentImages.push_back(caster->GetOutput());
for (unsigned int g=0; g<m_Parameters.m_SignalGen.GetNumVolumes(); g++)
{
m_Rotation.Fill(0.0);
m_Translation.Fill(0.0);
m_Rotations.push_back(m_Rotation);
m_Translations.push_back(m_Translation);
}
}
// resample mask image and frequency map to fit upsampled geometry
if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_WorkingImageRegion)
{
PrintToLog("Resampling tissue mask", false);
// rescale mask image (otherwise there are problems with the resampling)
auto rescaler = itk::RescaleIntensityImageFilter<ItkUcharImgType,ItkUcharImgType>::New();
rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage);
rescaler->SetOutputMaximum(100);
rescaler->SetOutputMinimum(0);
rescaler->Update();
// resample mask image
auto resampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
resampler->SetInput(rescaler->GetOutput());
- resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage);
resampler->SetSize(m_WorkingImageRegion.GetSize());
resampler->SetOutputSpacing(m_WorkingSpacing);
resampler->SetOutputOrigin(m_WorkingOrigin);
+ resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection);
+ resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() );
auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
resampler->SetInterpolator(nn_interpolator);
resampler->Update();
m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput();
}
// resample frequency map
if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion()!=m_WorkingImageRegion)
{
PrintToLog("Resampling frequency map", false);
auto resampler = itk::ResampleImageFilter<ItkFloatImgType, ItkFloatImgType>::New();
resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap);
- resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_FrequencyMap);
resampler->SetSize(m_WorkingImageRegion.GetSize());
resampler->SetOutputSpacing(m_WorkingSpacing);
resampler->SetOutputOrigin(m_WorkingOrigin);
+ resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection);
+ resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() );
auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkFloatImgType, double>::New();
resampler->SetInterpolator(nn_interpolator);
resampler->Update();
m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput();
}
m_MaskImageSet = true;
if (m_Parameters.m_SignalGen.m_MaskImage.IsNull())
{
// no input tissue mask is set -> create default
PrintToLog("No tissue mask set", false);
m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New();
m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing );
m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin );
m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion );
m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion );
m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion );
m_Parameters.m_SignalGen.m_MaskImage->Allocate();
m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100);
m_MaskImageSet = false;
}
else
{
PrintToLog("Using tissue mask", false);
}
if (m_Parameters.m_SignalGen.m_DoAddMotion)
{
if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
{
PrintToLog("Random motion artifacts:", false);
PrintToLog("Maximum rotation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
PrintToLog("Maximum translation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
}
else
{
PrintToLog("Linear motion artifacts:", false);
PrintToLog("Maximum rotation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
PrintToLog("Maximum translation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
}
}
if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() )
{
// no motion in first volume
m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
// motion in all other volumes
while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() )
{
m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true);
}
}
// we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false
while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()<m_Parameters.m_SignalGen.GetNumVolumes() )
{
m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
}
m_NumMotionVolumes = 0;
for (unsigned int i=0; i<m_Parameters.m_SignalGen.GetNumVolumes(); ++i)
{
if (m_Parameters.m_SignalGen.m_MotionVolumes[i])
++m_NumMotionVolumes;
}
m_MotionCounter = 0;
// creat image to hold transformed mask (motion artifact)
m_TransformedMaskImage = ItkUcharImgType::New();
auto duplicator = itk::ImageDuplicator<ItkUcharImgType>::New();
duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage);
duplicator->Update();
m_TransformedMaskImage = duplicator->GetOutput();
// second upsampling needed for motion artifacts
ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion;
DoubleVectorType upsampledSpacing = m_WorkingSpacing;
upsampledSpacing[0] /= 4;
upsampledSpacing[1] /= 4;
upsampledSpacing[2] /= 4;
upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4);
upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4);
upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4);
itk::Point<double,3> upsampledOrigin = m_WorkingOrigin;
upsampledOrigin[0] -= m_WorkingSpacing[0]/2;
upsampledOrigin[0] += upsampledSpacing[0]/2;
upsampledOrigin[1] -= m_WorkingSpacing[1]/2;
upsampledOrigin[1] += upsampledSpacing[1]/2;
upsampledOrigin[2] -= m_WorkingSpacing[2]/2;
upsampledOrigin[2] += upsampledSpacing[2]/2;
m_UpsampledMaskImage = ItkUcharImgType::New();
auto upsampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage);
upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage);
upsampler->SetSize(upsampledImageRegion.GetSize());
upsampler->SetOutputSpacing(upsampledSpacing);
upsampler->SetOutputOrigin(upsampledOrigin);
auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
upsampler->SetInterpolator(nn_interpolator);
upsampler->Update();
m_UpsampledMaskImage = upsampler->GetOutput();
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::InitializeFiberData()
{
m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000;
auto caster = itk::CastImageFilter< itk::Image<unsigned char, 3>, itk::Image<float, 3> >::New();
caster->SetInput(m_TransformedMaskImage);
caster->Update();
vtkSmartPointer<vtkFloatArray> weights = m_FiberBundle->GetFiberWeights();
float mean_weight = 0;
for (int i=0; i<weights->GetSize(); i++)
mean_weight += weights->GetValue(i);
mean_weight /= weights->GetSize();
if (mean_weight>0.000001)
for (int i=0; i<weights->GetSize(); i++)
m_FiberBundle->SetFiberWeight(i, weights->GetValue(i)/mean_weight);
else
PrintToLog("\nWarning: streamlines have VERY low weights. Average weight: " + boost::lexical_cast<std::string>(mean_weight) + ". Possible source of calculation errors.", false, true, true);
auto density_calculator = itk::TractDensityImageFilter< itk::Image<float, 3> >::New();
density_calculator->SetFiberBundle(m_FiberBundle);
density_calculator->SetInputImage(caster->GetOutput());
density_calculator->SetBinaryOutput(false);
density_calculator->SetUseImageGeometry(true);
density_calculator->SetOutputAbsoluteValues(true);
density_calculator->Update();
float max_density = density_calculator->GetMaxDensity();
float voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2];
if (m_mmRadius>0)
{
std::stringstream stream;
stream << std::fixed << setprecision(2) << itk::Math::pi*m_mmRadius*m_mmRadius*max_density;
std::string s = stream.str();
PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true);
{
float full_radius = 1000*std::sqrt(voxel_volume/(max_density*itk::Math::pi));
std::stringstream stream;
stream << std::fixed << setprecision(2) << full_radius;
std::string s = stream.str();
PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true);
}
}
else
{
m_mmRadius = std::sqrt(voxel_volume/(max_density*itk::Math::pi));
std::stringstream stream;
stream << std::fixed << setprecision(2) << m_mmRadius*1000;
std::string s = stream.str();
PrintToLog("\nSetting fiber radius to " + s + "µm to obtain full voxel.", false, true, true);
}
// a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy
m_FiberBundleTransformed = m_FiberBundle;
}
template< class PixelType >
bool TractsToDWIImageFilter< PixelType >::PrepareLogFile()
{
assert( ! m_Logfile.is_open() );
std::string filePath;
std::string fileName;
// Get directory name:
if (m_Parameters.m_Misc.m_OutputPath.size() > 0)
{
filePath = m_Parameters.m_Misc.m_OutputPath;
if( *(--(filePath.cend())) != '/')
{
filePath.push_back('/');
}
}
else
{
filePath = mitk::IOUtil::GetTempPath() + '/';
}
// check if directory exists, else use /tmp/:
if( itksys::SystemTools::FileIsDirectory( filePath ) )
{
while( *(--(filePath.cend())) == '/')
{
filePath.pop_back();
}
filePath = filePath + '/';
}
else
{
filePath = mitk::IOUtil::GetTempPath() + '/';
}
// Get file name:
if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() )
{
fileName = m_Parameters.m_Misc.m_ResultNode->GetName();
}
else
{
fileName = "";
}
if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() )
{
fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName;
}
else
{
fileName = "fiberfox";
}
// check if file already exists and DO NOT overwrite existing files:
std::string NameTest = fileName;
int c = 0;
while( itksys::SystemTools::FileExists( filePath + '/' + fileName + ".log" )
&& c <= std::numeric_limits<int>::max() )
{
fileName = NameTest + "_" + boost::lexical_cast<std::string>(c);
++c;
}
try
{
m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() );
}
catch (const std::ios_base::failure &fail)
{
MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what()
<< " while trying to open file" << filePath << '/' << fileName << ".log";
return false;
}
if ( m_Logfile.is_open() )
{
PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false );
return true;
}
else
{
m_StatusText += "Logfile could not be opened!\n";
MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Logfile could not be opened!";
return false;
}
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::GenerateData()
{
// prepare logfile
if ( ! PrepareLogFile() )
{
this->SetAbortGenerateData( true );
return;
}
m_TimeProbe.Start();
// check input data
if (m_FiberBundle.IsNull() && m_InputImage.IsNull())
itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!");
if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull())
itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted"
" image is nullptr! At least one fiber compartment is necessary to simulate diffusion.");
if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull())
itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted"
" image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion.");
if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // no partial volume? remove all but first fiber compartment
while (m_Parameters.m_FiberModelList.size()>1)
m_Parameters.m_FiberModelList.pop_back();
// int baselineIndex = m_Parameters.m_SignalGen.GetFirstBaselineIndex();
// if (baselineIndex<0) { itkExceptionMacro("No baseline index found!"); }
if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed
m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false;
if (m_UseConstantRandSeed) // always generate the same random numbers?
m_RandGen->SetSeed(0);
else
m_RandGen->SetSeed();
InitializeData();
if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation
{
CheckVolumeFractionImages();
InitializeFiberData();
int numFiberCompartments = m_Parameters.m_FiberModelList.size();
int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
double maxVolume = 0;
unsigned long lastTick = 0;
int signalModelSeed = m_RandGen->GetIntegerVariate();
PrintToLog("\n", false, false);
PrintToLog("Generating " + boost::lexical_cast<std::string>(numFiberCompartments+numNonFiberCompartments)
+ "-compartment diffusion-weighted signal.");
std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues();
PrintToLog("b-values: ", false, false, true);
for (auto v : bVals)
PrintToLog(boost::lexical_cast<std::string>(v) + " ", false, false, true);
PrintToLog("\n", false, false, true);
PrintToLog("\n", false, false, true);
int numFibers = m_FiberBundle->GetNumFibers();
boost::progress_display disp(numFibers*m_Parameters.m_SignalGen.GetNumVolumes());
if (m_FiberBundle->GetMeanFiberLength()<5.0)
omp_set_num_threads(2);
PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
for (unsigned int g=0; g<m_Parameters.m_SignalGen.GetNumVolumes(); ++g)
{
// move fibers
SimulateMotion(g);
// Set signal model random generator seeds to get same configuration in each voxel
for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size(); i++)
m_Parameters.m_FiberModelList.at(i)->SetSeed(signalModelSeed);
for (std::size_t i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
m_Parameters.m_NonFiberModelList.at(i)->SetSeed(signalModelSeed);
// storing voxel-wise intra-axonal volume in mm³
auto intraAxonalVolumeImage = ItkDoubleImgType::New();
intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing );
intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin );
intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion );
intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion );
intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion );
intraAxonalVolumeImage->Allocate();
intraAxonalVolumeImage->FillBuffer(0);
maxVolume = 0;
if (this->GetAbortGenerateData())
continue;
vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData();
// generate fiber signal (if there are any fiber models present)
if (!m_Parameters.m_FiberModelList.empty())
{
#pragma omp parallel for
for( int i=0; i<numFibers; i++ )
{
if (this->GetAbortGenerateData())
continue;
float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i);
int numPoints = -1;
std::vector< itk::Vector<double, 3> > points_copy;
#pragma omp critical
{
vtkCell* cell = fiberPolyData->GetCell(i);
numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
for (int j=0; j<numPoints; j++)
points_copy.push_back(GetItkVector(points->GetPoint(j)));
}
if (numPoints<2)
continue;
for( int j=0; j<numPoints - 1; j++)
{
if (this->GetAbortGenerateData())
{
j=numPoints;
continue;
}
itk::Vector<double> v = points_copy.at(j);
itk::Vector<double, 3> dir = points_copy.at(j+1)-v;
if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] )
continue;
dir.Normalize();
itk::Point<float, 3> startVertex = points_copy.at(j);
itk::Index<3> startIndex;
itk::ContinuousIndex<float, 3> startIndexCont;
m_TransformedMaskImage->TransformPhysicalPointToIndex(startVertex, startIndex);
m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont);
itk::Point<float, 3> endVertex = points_copy.at(j+1);
itk::Index<3> endIndex;
itk::ContinuousIndex<float, 3> endIndexCont;
m_TransformedMaskImage->TransformPhysicalPointToIndex(endVertex, endIndex);
m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont);
std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(m_WorkingSpacing, startIndex, endIndex, startIndexCont, endIndexCont);
for (std::pair< itk::Index<3>, double > seg : segments)
{
if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_TransformedMaskImage->GetPixel(seg.first)<=0)
continue;
// generate signal for each fiber compartment
for (int k=0; k<numFiberCompartments; k++)
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(seg.first);
float seg_volume = seg.second*fiberWeight*itk::Math::pi*m_mmRadius*m_mmRadius;
pix[g] += seg_volume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g, dir);
m_CompartmentImages.at(k)->SetPixel(seg.first, pix);
if (k==0)
{
// update fiber volume image
double vol = intraAxonalVolumeImage->GetPixel(seg.first) + seg_volume;
intraAxonalVolumeImage->SetPixel(seg.first, vol);
if (vol>maxVolume) { maxVolume = vol; }
}
}
}
}
#pragma omp critical
{
// progress report
++disp;
unsigned long newTick = 50*disp.count()/disp.expected_count();
for (unsigned int tick = 0; tick<(newTick-lastTick); ++tick)
PrintToLog("*", false, false, false);
lastTick = newTick;
}
}
}
// axon radius not manually defined --> set fullest voxel (maxVolume) to full fiber voxel
double density_correctiony_global = 1.0;
if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001)
density_correctiony_global = m_VoxelVolume/maxVolume;
// generate non-fiber signal
ImageRegionIterator<ItkUcharImgType> it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion());
while(!it3.IsAtEnd())
{
if (it3.Get()>0)
{
DoubleDwiType::IndexType index = it3.GetIndex();
double iAxVolume = intraAxonalVolumeImage->GetPixel(index);
// get non-transformed point (remove headmotion tranformation)
// this point lives in the volume fraction image space
itk::Point<double, 3> volume_fraction_point;
if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_MotionVolumes[g] )
volume_fraction_point = GetMovedPoint(index, false);
else
m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, volume_fraction_point);
if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)
{
if (iAxVolume>0.0001) // scale fiber compartment to voxel
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
pix[g] *= m_VoxelVolume/iAxVolume;
m_CompartmentImages.at(0)->SetPixel(index, pix);
if (g==0)
m_VolumeFractions.at(0)->SetPixel(index, 1);
}
else
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
pix[g] = 0;
m_CompartmentImages.at(0)->SetPixel(index, pix);
SimulateExtraAxonalSignal(index, volume_fraction_point, 0, g);
}
}
else
{
// manually defined axon radius and voxel overflow --> rescale to voxel volume
if ( m_Parameters.m_SignalGen.m_AxonRadius>=0.0001 && iAxVolume>m_VoxelVolume )
{
for (int i=0; i<numFiberCompartments; ++i)
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] *= m_VoxelVolume/iAxVolume;
m_CompartmentImages.at(i)->SetPixel(index, pix);
}
iAxVolume = m_VoxelVolume;
}
// if volume fraction image is set use it, otherwise use global scaling factor
double density_correction_voxel = density_correctiony_global;
if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 )
{
m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage());
double volume_fraction = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator);
if (volume_fraction<0)
mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!";
density_correction_voxel = m_VoxelVolume*volume_fraction/iAxVolume; // remove iAxVolume sclaing and scale to volume_fraction
}
else if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr)
density_correction_voxel = 0.0;
// adjust intra-axonal compartment volume by density correction factor
DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
pix[g] *= density_correction_voxel;
m_CompartmentImages.at(0)->SetPixel(index, pix);
// normalize remaining fiber volume fractions (they are rescaled in SimulateExtraAxonalSignal)
if (iAxVolume>0.0001)
{
for (int i=1; i<numFiberCompartments; i++)
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] /= iAxVolume;
m_CompartmentImages.at(i)->SetPixel(index, pix);
}
}
else
{
for (int i=1; i<numFiberCompartments; i++)
{
DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] = 0;
m_CompartmentImages.at(i)->SetPixel(index, pix);
}
}
iAxVolume = density_correction_voxel*iAxVolume; // new intra-axonal volume = old intra-axonal volume * correction factor
// simulate other compartments
SimulateExtraAxonalSignal(index, volume_fraction_point, iAxVolume, g);
}
}
++it3;
}
}
PrintToLog("\n", false);
}
if (this->GetAbortGenerateData())
{
PrintToLog("\n", false, false);
PrintToLog("Simulation aborted");
return;
}
DoubleDwiType::Pointer doubleOutImage;
double signalScale = m_Parameters.m_SignalGen.m_SignalScale;
if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff
{
PrintToLog("\n", false, false);
PrintToLog("Simulating k-space acquisition using "
+boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_NumberOfCoils)
+" coil(s)");
switch (m_Parameters.m_SignalGen.m_AcquisitionType)
{
case SignalGenerationParameters::SingleShotEpi:
{
PrintToLog("Acquisition type: single shot EPI", false);
break;
}
case SignalGenerationParameters::SpinEcho:
{
PrintToLog("Acquisition type: classic spin echo with cartesian k-space trajectory", false);
break;
}
default:
{
PrintToLog("Acquisition type: single shot EPI", false);
break;
}
}
if (m_Parameters.m_SignalGen.m_NoiseVariance>0 && m_Parameters.m_Misc.m_CheckAddNoiseBox)
PrintToLog("Simulating complex Gaussian noise", false);
if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation)
PrintToLog("Simulating signal relaxation", false);
if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull())
PrintToLog("Simulating distortions", false);
if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
PrintToLog("Simulating ringing artifacts", false);
if (m_Parameters.m_SignalGen.m_EddyStrength>0)
PrintToLog("Simulating eddy currents", false);
if (m_Parameters.m_SignalGen.m_Spikes>0)
PrintToLog("Simulating spikes", false);
if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0)
PrintToLog("Simulating aliasing artifacts", false);
if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0)
PrintToLog("Simulating ghosts", false);
doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages);
signalScale = 1; // already scaled in SimulateKspaceAcquisition()
}
else // don't do k-space stuff, just sum compartments
{
PrintToLog("Summing compartments");
doubleOutImage = m_CompartmentImages.at(0);
for (unsigned int i=1; i<m_CompartmentImages.size(); i++)
{
auto adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New();
adder->SetInput1(doubleOutImage);
adder->SetInput2(m_CompartmentImages.at(i));
adder->Update();
doubleOutImage = adder->GetOutput();
}
}
if (this->GetAbortGenerateData())
{
PrintToLog("\n", false, false);
PrintToLog("Simulation aborted");
return;
}
PrintToLog("Finalizing image");
if (signalScale>1)
PrintToLog(" Scaling signal", false);
if (m_Parameters.m_NoiseModel)
PrintToLog(" Adding noise", false);
unsigned int window = 0;
unsigned int min = itk::NumericTraits<unsigned int>::max();
ImageRegionIterator<OutputImageType> it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion());
DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels());
PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
int lastTick = 0;
while(!it4.IsAtEnd())
{
if (this->GetAbortGenerateData())
{
PrintToLog("\n", false, false);
PrintToLog("Simulation aborted");
return;
}
++disp2;
unsigned long newTick = 50*disp2.count()/disp2.expected_count();
for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
PrintToLog("*", false, false, false);
lastTick = newTick;
typename OutputImageType::IndexType index = it4.GetIndex();
signal = doubleOutImage->GetPixel(index)*signalScale;
if (m_Parameters.m_NoiseModel)
m_Parameters.m_NoiseModel->AddNoise(signal);
for (unsigned int i=0; i<signal.Size(); i++)
{
if (signal[i]>0)
signal[i] = floor(signal[i]+0.5);
else
signal[i] = ceil(signal[i]-0.5);
if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window)
window = signal[i];
if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]<min)
min = signal[i];
}
it4.Set(signal);
++it4;
}
window -= min;
unsigned int level = window/2 + min;
m_LevelWindow.SetLevelWindow(level, window);
this->SetNthOutput(0, m_OutputImage);
PrintToLog("\n", false);
PrintToLog("Finished simulation");
m_TimeProbe.Stop();
if (m_Parameters.m_SignalGen.m_DoAddMotion)
{
PrintToLog("\nHead motion log:", false);
PrintToLog(m_MotionLog, false, false);
}
if (m_Parameters.m_SignalGen.m_Spikes>0)
{
PrintToLog("\nSpike log:", false);
PrintToLog(m_SpikeLog, false, false);
}
if (m_Logfile.is_open()) m_Logfile.close();
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut)
{
// timestamp
if (addTime)
{
m_Logfile << this->GetTime() << " > ";
m_StatusText += this->GetTime() + " > ";
if (stdOut)
std::cout << this->GetTime() << " > ";
}
// message
if (m_Logfile.is_open())
m_Logfile << m;
m_StatusText += m;
if (stdOut)
std::cout << m;
// new line
if (linebreak)
{
if (m_Logfile.is_open())
m_Logfile << "\n";
m_StatusText += "\n";
if (stdOut)
std::cout << "\n";
}
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g)
{
// is motion artifact enabled?
// is the current volume g affected by motion?
if ( m_Parameters.m_SignalGen.m_DoAddMotion
&& m_Parameters.m_SignalGen.m_MotionVolumes[g]
&& g<static_cast<int>(m_Parameters.m_SignalGen.GetNumVolumes()) )
{
if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion )
{
// either undo last transform or work on fresh copy of untransformed fibers
m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy();
m_Rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2)
-m_Parameters.m_SignalGen.m_Rotation[0];
m_Rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2)
-m_Parameters.m_SignalGen.m_Rotation[1];
m_Rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2)
-m_Parameters.m_SignalGen.m_Rotation[2];
m_Translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2)
-m_Parameters.m_SignalGen.m_Translation[0];
m_Translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2)
-m_Parameters.m_SignalGen.m_Translation[1];
m_Translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2)
-m_Parameters.m_SignalGen.m_Translation[2];
}
else
{
m_Rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes;
m_Translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes;
m_MotionCounter++;
}
// move mask image
if (m_MaskImageSet)
{
ImageRegionIterator<ItkUcharImgType> maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion());
m_TransformedMaskImage->FillBuffer(0);
while(!maskIt.IsAtEnd())
{
if (maskIt.Get()<=0)
{
++maskIt;
continue;
}
DoubleDwiType::IndexType index = maskIt.GetIndex();
m_TransformedMaskImage->TransformPhysicalPointToIndex(GetMovedPoint(index, true), index);
if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index))
m_TransformedMaskImage->SetPixel(index,100);
++maskIt;
}
}
if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
{
m_Rotations.push_back(m_Rotation);
m_Translations.push_back(m_Translation);
m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0])
+ "," + boost::lexical_cast<std::string>(m_Rotation[1])
+ "," + boost::lexical_cast<std::string>(m_Rotation[2]) + ";";
m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0])
+ "," + boost::lexical_cast<std::string>(m_Translation[1])
+ "," + boost::lexical_cast<std::string>(m_Translation[2]) + "\n";
}
else
{
m_Rotations.push_back(m_Rotation*m_MotionCounter);
m_Translations.push_back(m_Translation*m_MotionCounter);
m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0]*m_MotionCounter)
+ "," + boost::lexical_cast<std::string>(m_Rotation[1]*m_MotionCounter)
+ "," + boost::lexical_cast<std::string>(m_Rotation[2]*m_MotionCounter) + ";";
m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0]*m_MotionCounter)
+ "," + boost::lexical_cast<std::string>(m_Translation[1]*m_MotionCounter)
+ "," + boost::lexical_cast<std::string>(m_Translation[2]*m_MotionCounter) + "\n";
}
m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]);
}
else
{
m_Rotation.Fill(0.0);
m_Translation.Fill(0.0);
m_Rotations.push_back(m_Rotation);
m_Translations.push_back(m_Translation);
m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0])
+ "," + boost::lexical_cast<std::string>(m_Rotation[1])
+ "," + boost::lexical_cast<std::string>(m_Rotation[2]) + ";";
m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0])
+ "," + boost::lexical_cast<std::string>(m_Translation[1])
+ "," + boost::lexical_cast<std::string>(m_Translation[2]) + "\n";
}
}
template< class PixelType >
itk::Point<double, 3> TractsToDWIImageFilter< PixelType >::GetMovedPoint(itk::Index<3>& index, bool forward)
{
itk::Point<double, 3> transformed_point;
if (forward)
{
m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, transformed_point);
if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
{
transformed_point = m_FiberBundle->TransformPoint(transformed_point.GetVnlVector(),
m_Rotation[0],m_Rotation[1],m_Rotation[2],
m_Translation[0],m_Translation[1],m_Translation[2]);
}
else
{
transformed_point = m_FiberBundle->TransformPoint(transformed_point.GetVnlVector(),
m_Rotation[0]*m_MotionCounter,m_Rotation[1]*m_MotionCounter,m_Rotation[2]*m_MotionCounter,
m_Translation[0]*m_MotionCounter,m_Translation[1]*m_MotionCounter,m_Translation[2]*m_MotionCounter);
}
}
else
{
m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, transformed_point);
if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
{
transformed_point = m_FiberBundle->TransformPoint( transformed_point.GetVnlVector(),
-m_Rotation[0], -m_Rotation[1], -m_Rotation[2],
-m_Translation[0], -m_Translation[1], -m_Translation[2] );
}
else
{
transformed_point = m_FiberBundle->TransformPoint( transformed_point.GetVnlVector(),
-m_Rotation[0]*m_MotionCounter, -m_Rotation[1]*m_MotionCounter, -m_Rotation[2]*m_MotionCounter,
-m_Translation[0]*m_MotionCounter, -m_Translation[1]*m_MotionCounter, -m_Translation[2]*m_MotionCounter );
}
}
return transformed_point;
}
template< class PixelType >
void TractsToDWIImageFilter< PixelType >::
SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point<double, 3>& volume_fraction_point, double intraAxonalVolume, int g)
{
int numFiberCompartments = m_Parameters.m_FiberModelList.size();
int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)
{
// simulate signal for largest non-fiber compartment
int max_compartment_index = 0;
double max_fraction = 0;
if (numNonFiberCompartments>1)
{
for (int i=0; i<numNonFiberCompartments; ++i)
{
m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
double compartment_fraction = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator);
if (compartment_fraction<0)
mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!";
if (compartment_fraction>max_fraction)
{
max_fraction = compartment_fraction;
max_compartment_index = i;
}
}
}
DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(max_compartment_index+numFiberCompartments);
DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index);
pix[g] += m_Parameters.m_NonFiberModelList[max_compartment_index]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume;
doubleDwi->SetPixel(index, pix);
if (g==0)
m_VolumeFractions.at(max_compartment_index+numFiberCompartments)->SetPixel(index, 1);
}
else
{
std::vector< double > fractions;
if (g==0)
m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume);
double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume
if (extraAxonalVolume<0)
{
if (extraAxonalVolume<-0.001)
MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume;
extraAxonalVolume = 0;
}
double interAxonalVolume = 0;
if (numFiberCompartments>1)
interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment
double nonFiberVolume = extraAxonalVolume - interAxonalVolume; // rest of compartment
if (nonFiberVolume<0)
{
if (nonFiberVolume<-0.001)
MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!";
nonFiberVolume = 0;
interAxonalVolume = extraAxonalVolume;
}
double compartmentSum = intraAxonalVolume;
fractions.push_back(intraAxonalVolume/m_VoxelVolume);
// rescale extra-axonal fiber signal
for (int i=1; i<numFiberCompartments; ++i)
{
if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr)
{
m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage());
interAxonalVolume = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume;
if (interAxonalVolume<0)
mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!";
}
DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
pix[g] *= interAxonalVolume;
m_CompartmentImages.at(i)->SetPixel(index, pix);
compartmentSum += interAxonalVolume;
fractions.push_back(interAxonalVolume/m_VoxelVolume);
if (g==0)
m_VolumeFractions.at(i)->SetPixel(index, interAxonalVolume/m_VoxelVolume);
}
for (int i=0; i<numNonFiberCompartments; ++i)
{
double volume = nonFiberVolume;
if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr)
{
m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
volume = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume;
if (volume<0)
mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!";
if (m_UseRelativeNonFiberVolumeFractions)
volume *= nonFiberVolume/m_VoxelVolume;
}
DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index);
pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*volume;
m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix);
compartmentSum += volume;
fractions.push_back(volume/m_VoxelVolume);
if (g==0)
m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, volume/m_VoxelVolume);
}
if (compartmentSum/m_VoxelVolume>1.05)
{
MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")";
for (auto val : fractions)
MITK_ERROR << val;
}
}
}
template< class PixelType >
itk::Point<float, 3> TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3])
{
itk::Point<float, 3> itkPoint;
itkPoint[0] = point[0];
itkPoint[1] = point[1];
itkPoint[2] = point[2];
return itkPoint;
}
template< class PixelType >
itk::Vector<double, 3> TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3])
{
itk::Vector<double, 3> itkVector;
itkVector[0] = point[0];
itkVector[1] = point[1];
itkVector[2] = point[2];
return itkVector;
}
template< class PixelType >
vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3])
{
vnl_vector_fixed<double, 3> vnlVector;
vnlVector[0] = point[0];
vnlVector[1] = point[1];
vnlVector[2] = point[2];
return vnlVector;
}
template< class PixelType >
vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector<float,3>& vector)
{
vnl_vector_fixed<double, 3> 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<std::string>(hours));
out.append(":");
out.append(boost::lexical_cast<std::string>(minutes));
out.append(":");
out.append(boost::lexical_cast<std::string>(seconds));
m_TimeProbe.Start();
return out;
}
}